! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE BINREAD ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE binread(nx,ny,nz,nzsoil,nstyps,grdbas,inch,time,x,y,z,zp, & 2,5 zpsoil,uprt, vprt, wprt, ptprt, pprt, & qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, & ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, & soiltyp,stypfrct,vegtyp,lai,roufns,veg, & tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & ireturn) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Read in binary data set created by ARPS using history dump format ! No. 1. ! All data read in are located at the original staggered grid points ! !----------------------------------------------------------------------- ! ! AUTHOR: Ming Xue ! 2/27/92. ! ! MODIFICATION HISTORY: ! ! 6/08/92 ! Added full documentation (K. Brewster) ! ! 7/14/92 (K. Brewster) ! Added runname, comment and version number reading ! ! 8/20/92 (M. Xue) ! Added data reading of computational z coordinate array z. ! ! 4/23/93 (M. Xue) ! New data format. ! ! 02/06/95 (Y. Liu) ! Added map projection parameters into the binary dumping ! ! 03/26/96 (G. Bassett) ! Backwards compatibility added for ARPS 3.2 and ARPS 4.0 binary ! history dump formats. ! ! 12/09/1998 (Donghai Wang) ! Added the snow cover. ! ! 05/15/2002 (J. Brotzge) ! Added additional variables to allow for multiple soil schemes ! ! 1 June 2002 Eric Kemp ! Bug fixes for new soil variables. ! !----------------------------------------------------------------------- ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! nzsoil Number of grid points in the soil ! ! grdbas Data read flag. ! =1, only grid and base state arrays will be read ! =0, all arrays will be read based on data ! parameter setting. ! inch Channel number for binary reading. ! This channel must be opened for unformatted reading ! by the calling routine. ! ! OUTPUT: ! ! time Time in seconds of data read from "filename" ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp z coordinate of grid points in physical space (m) ! zpsoil z coordinate of grid points in the soil (m) ! ! uprt x component of perturbation velocity (m/s) ! vprt y component of perturbation velocity (m/s) ! wprt Vertical component of perturbation velocity in ! Cartesian coordinates (m/s). ! ptprt Perturbation potential temperature (K) ! pprt Perturbation pressure (Pascal) ! ! qvprt Perturbation water vapor mixing ratio (kg/kg) ! qc Cloud water mixing ratio (kg/kg) ! qr Rainwater mixing ratio (kg/kg) ! qi Cloud ice mixing ratio (kg/kg) ! qs Snow mixing ratio (kg/kg) ! qh Hail mixing ratio (kg/kg) ! tke Turbulent Kinetic Energy ((m/s)**2) ! ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! ! ubar Base state x velocity component (m/s) ! vbar Base state y velocity component (m/s) ! wbar Base state z velocity component (m/s) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! rhobar Base state air density (kg/m**3) ! qvbar Base state water vapor mixing ratio (kg/kg) ! ! soiltyp Soil type ! stypfrct Soil type fraction ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil Soil temperature (K) ! qsoil Soil moisture (m**3/m**3) ! wetcanp Canopy water amount ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! prcrate Precipitation rates ! ! radfrc Radiation forcing (K/s) ! radsw Solar radiation reaching the surface ! rnflx Net radiation flux absorbed by surface ! radswnet Net shortwave radiation ! radlwin Incoming longwave radiation ! ! usflx Surface flux of u-momentum (kg/(m*s**2)) ! vsflx Surface flux of v-momentum (kg/(m*s**2)) ! ptsflx Surface heat flux (K*kg/(m**2 * s )) ! qvsflx Surface moisture flux of (kg/(m**2 * s)) ! ! ireturn Return status indicator ! =0, successful read of all data ! =1, error reading data ! =2, end-of-file reached during read attempt ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Number of grid points in 3 directions INTEGER :: nzsoil ! Number of grid points in the soil REAL :: x (nx) ! x-coord. of the physical and compu ! -tational grid. Defined at u-point(m). REAL :: y (ny) ! y-coord. of the physical and compu ! -tational grid. Defined at v-point(m). REAL :: z (nz) ! z-coord. of the computational grid. ! Defined at w-point on the staggered ! grid(m). REAL :: zp (nx,ny,nz) ! Physical height coordinate defined at ! w-point of the staggered grid(m). REAL :: zpsoil(nx,ny,nzsoil) ! Physical height coordinate defined at ! w-point of the soil (m) INTEGER :: grdbas ! Data read flag. INTEGER :: inch ! Channel number for binary reading REAL :: time ! Time in seconds of data read ! from "filename" REAL :: uprt (nx,ny,nz) ! Perturbation u-velocity (m/s) REAL :: vprt (nx,ny,nz) ! Perturbation v-velocity (m/s) REAL :: wprt (nx,ny,nz) ! Perturbation w-velocity (m/s) REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K) REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal) REAL :: qvprt (nx,ny,nz) ! Perturbation water vapor mixing ! ratio (kg/kg) REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg) REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s) REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s) REAL :: wbar (nx,ny,nz) ! Base state w-velocity (m/s) REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K) REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal) REAL :: rhobar(nx,ny,nz) ! Base state air density (kg/m**3) REAL :: qvbar (nx,ny,nz) ! Base state water vapor mixing ratio INTEGER :: nstyps ! Number of soil type INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type REAL :: stypfrct(nx,ny,nstyps) ! Soil type fraction INTEGER :: vegtyp(nx,ny) ! Vegetation type REAL :: lai (nx,ny) ! Leaf Area Index REAL :: roufns (nx,ny) ! Surface roughness REAL :: veg (nx,ny) ! Vegetation fraction REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K) REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount REAL :: snowdpth(nx,ny) ! Snow depth (m) REAL :: raing(nx,ny) ! Grid supersaturation rain REAL :: rainc(nx,ny) ! Cumulus convective rain REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s)) ! prcrate(1,1,1) = total precip. rate ! prcrate(1,1,2) = grid scale precip. rate ! prcrate(1,1,3) = cumulus precip. rate ! prcrate(1,1,4) = microphysics precip. rate REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s) REAL :: radsw (nx,ny) ! Solar radiation reaching the surface REAL :: rnflx (nx,ny) ! Net radiation flux absorbed by surface REAL :: radswnet(nx,ny) ! Net shortwave radiation REAL :: radlwin(nx,ny) ! Incoming longwave radiation REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2)) REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2)) REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m**2*s)) REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s)) INTEGER :: ireturn ! Return status indicator ! !----------------------------------------------------------------------- ! ! Parameters describing routine that wrote the gridded data ! !----------------------------------------------------------------------- ! ! 06/28/2002 Zuwen He ! ! fmtver??: to label each data a version. ! intver??: an integer to allow faster comparison than fmtver??, ! which are strings. ! ! Verion 5.00: significant change in soil variables since version 4.10. ! CHARACTER (LEN=40) :: fmtver320,fmtver400,fmtver410,fmtver500 INTEGER :: intver,intver320,intver400,intver410,intver500 PARAMETER (fmtver320='003.20 Binary Data',intver320=320) PARAMETER (fmtver400='004.00 Binary Data',intver400=400) PARAMETER (fmtver410='004.10 Binary Data',intver410=410) PARAMETER (fmtver500='005.00 Binary Data',intver500=500) CHARACTER (LEN=40) :: fmtverin ! CHARACTER (LEN=10) :: tmunit ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! INTEGER :: lchanl PARAMETER (lchanl=6) ! Channel number for formatted printing. INTEGER :: i,j,k,is INTEGER :: nstyp1 CHARACTER (LEN=12) :: label ! INTEGER :: nxin,nyin,nzin INTEGER :: nxin,nyin,nzin,nzsoilin INTEGER :: bgrdin,bbasin,bvarin,bicein,btkein,btrbin INTEGER :: idummy REAL :: rdummy ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'indtflg.inc' INCLUDE 'globcst.inc' INCLUDE 'grid.inc' ! Grid & map parameters. INCLUDE 'mp.inc' ! mpi parameters. ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !----------------------------------------------------------------------- ! ! Read header info ! !----------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) fmtverin IF (fmtverin == fmtver320) THEN intver=intver320 ELSE IF (fmtverin == fmtver400) THEN intver=intver400 ELSE IF (fmtverin == fmtver410) THEN intver=intver410 ELSE IF (fmtverin == fmtver500) THEN intver=intver500 ELSE IF (myproc == 0) WRITE(6,'(/1x,a,a,a/)') & 'Incoming data format, fmtverin=',fmtverin, & ', not found. The Job stopped.' CALL arpsstop('arpstop called from binread. ',1) END IF IF (myproc == 0) WRITE(6,'(/1x,a,a/)') & 'Incoming data format, fmtverin=',fmtverin READ(inch,ERR=110,END=120) runname READ(inch,ERR=110,END=120) nocmnt IF( nocmnt > 0 ) THEN DO i=1,nocmnt READ(inch,ERR=110,END=120) cmnt(i) END DO END IF IF (myproc == 0) & WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') runname IF( nocmnt > 0 ) THEN IF(myproc == 0)THEN DO i=1,nocmnt WRITE(6,'(1x,a)') cmnt(i) END DO END IF END IF READ(inch,ERR=110,END=120) time,tmunit ! !----------------------------------------------------------------------- ! ! Get dimensions of data in binary file and check against ! the dimensions passed to BINREAD ! !----------------------------------------------------------------------- ! IF (intver <= intver410) THEN READ(inch,ERR=110,END=120) nxin, nyin, nzin nzsoilin = 2 ! for version prior to 410, it is a two-layer soil model ELSE IF (intver >= intver500) THEN READ(inch,ERR=110,END=120) nxin, nyin, nzin,nzsoilin END IF ! ! Data validation: dimensions ! IF( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz) THEN IF(myproc == 0)THEN WRITE(6,'(1x,a)') & ' Dimensions in BINREAD inconsistent with data.' WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin WRITE(6,'(1x,a,3I15)') ' Expected: ', nx, ny, nz WRITE(6,'(1x,a)') ' Program aborted in BINREAD.' END IF CALL arpsstop('arpstop called from binread nx-ny-nz read ',1) END IF IF(nzsoilin /= nzsoil) THEN IF (intver <= intver410) THEN IF(myproc == 0) WRITE(6,'(1x,a,a/,2(1x,a/))') & ' The incoming data version is ', fmtverin, & ' In the input file, nzsoil must be set to 2. ', & ' Program aborted in BINREAD.' ELSE IF (intver >= intver500) THEN IF(myproc == 0)THEN WRITE(6,'(1x,a)') & ' Dimensions in BINREAD inconsistent with data.' WRITE(6,'(1x,a,I15)') ' Read were: ', nzsoilin WRITE(6,'(1x,a,I15)') ' Expected: ', nzsoil WRITE(6,'(1x,a)') ' Program aborted in BINREAD.' END IF END IF CALL arpsstop('arpstop called from binread nzsoil read ',1) END IF ! !----------------------------------------------------------------------- ! ! Read in flags for different data groups ! !----------------------------------------------------------------------- ! IF( grdbas == 1 ) THEN ! Read grid and base state arrays IF (myproc == 0) & WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)') & 'To read grid and base state data at time ', time, & ' secs = ',(time/60.),' mins.' READ(inch,ERR=110,END=120) & bgrdin,bbasin,bvarin,mstin,bicein, & btrbin,idummy,idummy,landin,totin, & btkein,idummy,idummy,mapproj,month, & day, year, hour, minute, second ELSE ! Normal data reading IF (myproc == 0) & WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:', & time,' secs = ',(time/60.),' mins.' READ(inch,ERR=110,END=120) & grdin,basin,varin,mstin,icein, & trbin, sfcin,rainin,landin,totin, & tkein,idummy,idummy,mapproj, month, & day, year, hour, minute, second END IF READ(inch,ERR=110,END=120) & umove,vmove,xgrdorg,ygrdorg,trulat1, & trulat2,trulon,sclfct,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy, & tstop,thisdmp,latitud,ctrlat,ctrlon IF ( totin /= 0 ) THEN ! !----------------------------------------------------------------------- ! ! Read in additional parameters for ARPS history dump 4.0 or later ! version. ! !----------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) & nstyp1, prcin, radin, flxin,snowcin, & snowin,idummy,idummy,idummy,idummy, & idummy,idummy,idummy,idummy,idummy, & idummy,idummy,idummy,idummy,idummy IF ( nstyp1 < 1 ) THEN nstyp1 = 1 END IF READ(inch,ERR=110,END=120) & rdummy,rdummy,rdummy,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy END IF ! !----------------------------------------------------------------------- ! ! Read in x, y, z and zp arrays. ! !---------------------------------------------------------------------- ! IF( grdin == 1 .OR. grdbas == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) x IF (myproc == 0) WRITE(lchanl,910) label,' x.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) y IF (myproc == 0) WRITE(lchanl,910) label,' y.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) z IF (myproc == 0) WRITE(lchanl,910) label,' z.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) zp IF (myproc == 0) WRITE(lchanl,910) label,' zp.' IF (intver <= intver410) THEN ! ! nzsoil must equal to 2, 06/28/2002, Zuwen ! assume zpsoil(,,2) is one meter below the surface. ! DO j=1,ny DO i=1,nx zpsoil(i,j,1)=zp(i,j,2) zpsoil(i,j,2)=zpsoil(i,j,1)-1. END DO END DO IF (myproc == 0) THEN WRITE(lchanl,910) ' Assign zpsoil. ' WRITE(lchanl,910) ' Assume zpsoil(,,1) is zp(,,2). ' WRITE(lchanl,910) ' Assume zpsoil(,,2) is zp(,,2)-1. ' END IF ELSE IF (intver >= intver500) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) zpsoil IF (myproc == 0) WRITE(lchanl,910) label,' zpsoil.' END IF ! intver END IF ! grdin ! !----------------------------------------------------------------------- ! ! Read in base state fields ! !---------------------------------------------------------------------- ! IF( basin == 1 .OR. grdbas == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ubar IF (myproc == 0) & WRITE(lchanl,910) label,' ubar.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) vbar IF (myproc == 0) & WRITE(lchanl,910) label,' vbar.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) wbar IF (myproc == 0) & WRITE(lchanl,910) label,' wbar.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ptbar IF (myproc == 0) & WRITE(lchanl,910) label,' ptbar.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) pbar IF (myproc == 0) & WRITE(lchanl,910) label,' pbar.' IF( mstin == 1) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) qvbar IF (myproc == 0) & WRITE(lchanl,910) label,' qvbar.' END IF IF (landin == 1) THEN IF (nstyp1 <= 1) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((soiltyp(i,j,1),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' soiltyp.' ELSE DO is=1,nstyp1 IF (is <= nstyps) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((soiltyp(i,j,is),i=1,nx),j=1,ny) IF (myproc == 0) & WRITE(lchanl,910) label,' soiltyp.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((stypfrct(i,j,is),i=1,nx),j=1,ny) IF (myproc == 0) & WRITE(lchanl,910) label,' stypfrct.' ELSE READ(inch,ERR=110,END=120) label IF (myproc == 0) & WRITE(lchanl,910) label,'skipping soiltyp' READ(inch,ERR=110,END=120) READ(inch,ERR=110,END=120) label IF (myproc == 0) & WRITE(lchanl,910) label,'skipping stypfrct.' READ(inch,ERR=110,END=120) ENDIF END DO END IF CALL fix_stypfrct_nstyp(nx,ny,nstyp1,nstyp,stypfrct) READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) vegtyp IF (myproc == 0) WRITE(lchanl,910) label,' vegtyp.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) lai IF (myproc == 0) WRITE(lchanl,910) label,' lai.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) roufns IF (myproc == 0) WRITE(lchanl,910) label,' roufns.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) veg IF (myproc == 0) WRITE(lchanl,910) label,' veg.' END IF END IF IF( grdbas == 1 ) GO TO 930 IF( varin == 1 ) THEN IF ( totin == 0 ) THEN ! !----------------------------------------------------------------------- ! ! Read in perturbations from history dump ! !----------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) uprt IF (myproc == 0) WRITE(lchanl,910) label,' uprt.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) vprt IF (myproc == 0) WRITE(lchanl,910) label,' vprt.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) wprt IF (myproc == 0) WRITE(lchanl,910) label,' wprt.' ! !----------------------------------------------------------------------- ! ! Read in scalars ! !---------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ptprt IF (myproc == 0) WRITE(lchanl,910) label,' ptprt.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) pprt IF (myproc == 0) WRITE(lchanl,910) label,' pprt.' ELSE ! !----------------------------------------------------------------------- ! ! Read in total values of variables from history dump ! !---------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) uprt IF (myproc == 0) WRITE(lchanl,910) label,' u.' DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx uprt(i,j,k) = uprt(i,j,k) - ubar(i,j,k) END DO END DO END DO READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) vprt IF (myproc == 0) WRITE(lchanl,910) label,' v.' DO k=1,nz-1 DO j=1,ny DO i=1,nx-1 vprt(i,j,k) = vprt(i,j,k) - vbar(i,j,k) END DO END DO END DO READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) wprt IF (myproc == 0) WRITE(lchanl,910) label,' w.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ptprt IF (myproc == 0) WRITE(lchanl,910) label,' pt.' DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 ptprt(i,j,k) = ptprt(i,j,k) - ptbar(i,j,k) END DO END DO END DO READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) pprt IF (myproc == 0) WRITE(lchanl,910) label,' p.' DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 pprt(i,j,k) = pprt(i,j,k) - pbar(i,j,k) END DO END DO END DO END IF END IF ! !----------------------------------------------------------------------- ! ! Read in moisture variables ! !----------------------------------------------------------------------- ! IF( mstin == 1 ) THEN IF ( totin == 0 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) qvprt IF (myproc == 0) WRITE(lchanl,910) label,' qvprt.' ELSE READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) qvprt IF (myproc == 0) WRITE(lchanl,910) label,' qv.' DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 qvprt(i,j,k) = qvprt(i,j,k) - qvbar(i,j,k) END DO END DO END DO END IF READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) qc IF (myproc == 0) WRITE(lchanl,910) label,' qc.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) qr IF (myproc == 0) WRITE(lchanl,910) label,' qr.' IF( rainin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) raing IF (myproc == 0) WRITE(lchanl,910) label,' raing.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) rainc IF (myproc == 0) WRITE(lchanl,910) label,' rainc.' END IF IF( prcin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((prcrate(i,j,1),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' prcrate1.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((prcrate(i,j,2),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' prcrate2.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((prcrate(i,j,3),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' prcrate3.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((prcrate(i,j,4),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' prcrate4.' END IF IF( icein == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) qi IF (myproc == 0) WRITE(lchanl,910) label,' qi.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) qs IF (myproc == 0) WRITE(lchanl,910) label,' qs.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) qh IF (myproc == 0) WRITE(lchanl,910) label,' qh.' END IF END IF IF( tkein == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) tke IF (myproc == 0) WRITE(lchanl,910) label,' tke.' END IF IF( trbin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) kmh IF (myproc == 0) WRITE(lchanl,910) label,' kmh.' IF ( intver >= intver410 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) kmv IF (myproc == 0) WRITE(lchanl,910) label,' kmv.' END IF ! intver END IF IF( sfcin == 1 ) THEN IF (nstyp1 <= 1) THEN IF (intver <= intver410) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((tsoil(i,j,1,0),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' tsfc.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((tsoil(i,j,2,0),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((qsoil(i,j,1,0),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' wetsfc.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((qsoil(i,j,2,0),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' wetdp.' ELSE IF (intver >= intver500) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((tsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) IF (myproc == 0) WRITE(lchanl,910) label,' qsoil.' END IF ! intver READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((wetcanp(i,j,0),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' wetcanp.' ELSE DO is=0,nstyp1 IF (is <= nstyps) THEN IF (intver <= intver410) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((tsoil(i,j,1,is),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' tsfc.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((tsoil(i,j,2,is),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((qsoil(i,j,1,is),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' wetsfc.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((qsoil(i,j,2,is),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' wetdp.' ELSE IF (intver >= intver500) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((tsoil(i,j,k,is),i=1,nx),j=1,ny),k=1,nzsoil) IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qsoil(i,j,k,is),i=1,nx),j=1,ny),k=1,nzsoil) IF (myproc == 0) WRITE(lchanl,910) label,' qsoil.' END IF ! intver READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((wetcanp(i,j,is),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' wetcanp.' ELSE READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) IF (myproc == 0) & WRITE(lchanl,910) label,'skipping tsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) IF (myproc == 0) & WRITE(lchanl,910) label,'skipping qsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) IF (myproc == 0) & WRITE(lchanl,910) label,'skipping wetcanp.' ENDIF END DO END IF CALL fix_soil_nstyp(nx,ny,nzsoil,nstyp1,nstyp,tsoil,qsoil,wetcanp) IF (snowcin == 1) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) IF (myproc == 0) WRITE(lchanl,910) label,' snowcvr -- discarding.' END IF IF (snowin == 1) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((snowdpth(i,j),i=1,nx),j=1,ny) IF (myproc == 0) WRITE(lchanl,910) label,' snowdpth.' END IF END IF IF( radin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) radfrc IF (myproc == 0) WRITE(lchanl,910) label,' radfrc.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) radsw IF (myproc == 0) WRITE(lchanl,910) label,' radsw.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) rnflx IF (myproc == 0) WRITE(lchanl,910) label,' rnflx.' IF (intver <= intver410) THEN ! ! 06/28/2002 Zuwen He ! ! Do not know how to initialized radswnet, radlwin, yet, ! at this moment. ! radswnet = 0. radlwin = 0. ELSE IF (intver >= intver500) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) radswnet WRITE(lchanl,910) label,' radswnet.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) radlwin WRITE(lchanl,910) label,' radlwin.' END IF ! intver END IF IF( flxin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) usflx IF (myproc == 0) & WRITE(lchanl,910) label,' usflx.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) vsflx IF (myproc == 0) & WRITE(lchanl,910) label,' vsflx.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ptsflx IF (myproc == 0) & WRITE(lchanl,910) label,' ptsflx.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) qvsflx IF (myproc == 0) & WRITE(lchanl,910) label,' qvsflx.' END IF 910 FORMAT(1X,'Field ',a12,' was read into array',a) ! !----------------------------------------------------------------------- ! ! Friendly exit message ! !---------------------------------------------------------------------- ! 930 CONTINUE IF (myproc == 0) & WRITE(6,'(/a,F8.1,a/)') & ' Data at time=', time/60,' (min) were successfully read.' ireturn = 0 RETURN ! !----------------------------------------------------------------------- ! ! Error during read ! !---------------------------------------------------------------------- ! 110 CONTINUE WRITE(6,'(/a/)') ' Error reading data in BINREAD' ireturn=1 RETURN ! !----------------------------------------------------------------------- ! ! End-of-file during read ! !---------------------------------------------------------------------- ! 120 CONTINUE WRITE(6,'(/a/)') ' End of file reached in BINREAD' ireturn=2 RETURN END SUBROUTINE binread ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE BINREADSPLIT ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE binreadsplit(nx,ny,nz,nzsoil,nstyps,grdbas,inch,time,x,y,z, & 2,120 zp,zpsoil,uprt, vprt, wprt, ptprt, pprt, & qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, & ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, & soiltyp,stypfrct,vegtyp,lai,roufns,veg, & tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & ireturn) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Read in binary data set created by ARPS using history dump format ! No. 1. ! All data read in are located at the original staggered grid points ! !----------------------------------------------------------------------- ! ! AUTHOR: Yunheng Wang ! 9/03/02. ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! nzsoil Number of grid points in the soil ! ! grdbas Data read flag. ! =1, only grid and base state arrays will be read ! =0, all arrays will be read based on data ! parameter setting. ! inch Channel number for binary reading. ! This channel must be opened for unformatted reading ! by the calling routine. ! ! OUTPUT: ! ! time Time in seconds of data read from "filename" ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp z coordinate of grid points in physical space (m) ! zpsoil z coordinate of grid points in the soil (m) ! ! uprt x component of perturbation velocity (m/s) ! vprt y component of perturbation velocity (m/s) ! wprt Vertical component of perturbation velocity in ! Cartesian coordinates (m/s). ! ptprt Perturbation potential temperature (K) ! pprt Perturbation pressure (Pascal) ! ! qvprt Perturbation water vapor mixing ratio (kg/kg) ! qc Cloud water mixing ratio (kg/kg) ! qr Rainwater mixing ratio (kg/kg) ! qi Cloud ice mixing ratio (kg/kg) ! qs Snow mixing ratio (kg/kg) ! qh Hail mixing ratio (kg/kg) ! tke Turbulent Kinetic Energy ((m/s)**2) ! ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! ! ubar Base state x velocity component (m/s) ! vbar Base state y velocity component (m/s) ! wbar Base state z velocity component (m/s) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! rhobar Base state air density (kg/m**3) ! qvbar Base state water vapor mixing ratio (kg/kg) ! ! soiltyp Soil type ! stypfrct Soil type fraction ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil Soil temperature (K) ! qsoil Soil moisture (m**3/m**3) ! wetcanp Canopy water amount ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! prcrate Precipitation rates ! ! radfrc Radiation forcing (K/s) ! radsw Solar radiation reaching the surface ! rnflx Net radiation flux absorbed by surface ! radswnet Net shortwave radiation ! radlwin Incoming longwave radiation ! ! usflx Surface flux of u-momentum (kg/(m*s**2)) ! vsflx Surface flux of v-momentum (kg/(m*s**2)) ! ptsflx Surface heat flux (K*kg/(m**2 * s )) ! qvsflx Surface moisture flux of (kg/(m**2 * s)) ! ! ireturn Return status indicator ! =0, successful read of all data ! =1, error reading data ! =2, end-of-file reached during read attempt ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Number of grid points in 3 directions INTEGER :: nzsoil ! Number of grid points in the soil REAL :: x (nx) ! x-coord. of the physical and compu ! -tational grid. Defined at u-point(m). REAL :: y (ny) ! y-coord. of the physical and compu ! -tational grid. Defined at v-point(m). REAL :: z (nz) ! z-coord. of the computational grid. ! Defined at w-point on the staggered ! grid(m). REAL :: zp (nx,ny,nz) ! Physical height coordinate defined at ! w-point of the staggered grid(m). REAL :: zpsoil(nx,ny,nzsoil) ! Physical height coordinate defined at ! w-point of the soil (m) INTEGER :: grdbas ! Data read flag. INTEGER :: inch ! Channel number for binary reading REAL :: time ! Time in seconds of data read ! from "filename" REAL :: uprt (nx,ny,nz) ! Perturbation u-velocity (m/s) REAL :: vprt (nx,ny,nz) ! Perturbation v-velocity (m/s) REAL :: wprt (nx,ny,nz) ! Perturbation w-velocity (m/s) REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K) REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal) REAL :: qvprt (nx,ny,nz) ! Perturbation water vapor mixing ! ratio (kg/kg) REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg) REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s) REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s) REAL :: wbar (nx,ny,nz) ! Base state w-velocity (m/s) REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K) REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal) REAL :: rhobar(nx,ny,nz) ! Base state air density (kg/m**3) REAL :: qvbar (nx,ny,nz) ! Base state water vapor mixing ratio INTEGER :: nstyps ! Number of soil type INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type REAL :: stypfrct(nx,ny,nstyps) ! Soil type fraction INTEGER :: vegtyp(nx,ny) ! Vegetation type REAL :: lai (nx,ny) ! Leaf Area Index REAL :: roufns (nx,ny) ! Surface roughness REAL :: veg (nx,ny) ! Vegetation fraction REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K) REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount REAL :: snowdpth(nx,ny) ! Snow depth (m) REAL :: raing(nx,ny) ! Grid supersaturation rain REAL :: rainc(nx,ny) ! Cumulus convective rain REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s)) ! prcrate(1,1,1) = total precip. rate ! prcrate(1,1,2) = grid scale precip. rate ! prcrate(1,1,3) = cumulus precip. rate ! prcrate(1,1,4) = microphysics precip. rate REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s) REAL :: radsw (nx,ny) ! Solar radiation reaching the surface REAL :: rnflx (nx,ny) ! Net radiation flux absorbed by surface REAL :: radswnet(nx,ny) ! Net shortwave radiation REAL :: radlwin(nx,ny) ! Incoming longwave radiation REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2)) REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2)) REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m**2*s)) REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s)) INTEGER :: ireturn ! Return status indicator ! !----------------------------------------------------------------------- ! ! Parameters describing routine that wrote the gridded data ! !----------------------------------------------------------------------- ! ! fmtver??: to label each data a version. ! intver??: an integer to allow faster comparison than fmtver??, ! which are strings. ! ! Verion 5.00: significant change in soil variables since version 4.10. ! CHARACTER (LEN=40) :: fmtver320,fmtver400,fmtver410,fmtver500 INTEGER :: intver,intver320,intver400,intver410,intver500 PARAMETER (fmtver320='003.20 Binary Data',intver320=320) PARAMETER (fmtver400='004.00 Binary Data',intver400=400) PARAMETER (fmtver410='004.10 Binary Data',intver410=410) PARAMETER (fmtver500='005.00 Binary Data',intver500=500) CHARACTER (LEN=40) :: fmtverin ! CHARACTER (LEN=10) :: tmunit ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! INTEGER :: lchanl PARAMETER (lchanl=6) ! Channel number for formatted printing. INTEGER :: i,j,k,is INTEGER :: nstyp1 CHARACTER (LEN=12) :: label ! INTEGER :: nxin,nyin,nzin INTEGER :: nxin,nyin,nzin,nzsoilin INTEGER :: bgrdin,bbasin,bvarin,bicein,btkein,btrbin INTEGER :: idummy REAL :: rdummy INTEGER :: nxlg, nylg, nzlg, nzsoillg, n3d INTEGER, ALLOCATABLE :: var3di(:,:,:) REAL, ALLOCATABLE :: var2d(:,:), var3d(:,:,:) ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'indtflg.inc' INCLUDE 'globcst.inc' INCLUDE 'grid.inc' ! Grid & map parameters. INCLUDE 'mp.inc' ! mpi parameters. ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! nxlg = (nx-3)*nproc_x+3 nylg = (ny-3)*nproc_y+3 nzlg = nz nzsoillg = nzsoil n3d = MAX(nz, nzsoil, nstyps+1, 4) ALLOCATE(var2d(nxlg,nylg)) ALLOCATE(var3d(nxlg,nylg,n3d)) ALLOCATE(var3di(nxlg,nylg,n3d)) IF (myproc == 0) THEN !----------------------------------------------------------------------- ! ! Read header info ! !----------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) fmtverin IF (fmtverin == fmtver320) THEN intver=intver320 ELSE IF (fmtverin == fmtver400) THEN intver=intver400 ELSE IF (fmtverin == fmtver410) THEN intver=intver410 ELSE IF (fmtverin == fmtver500) THEN intver=intver500 ELSE WRITE(6,'(/1x,a,a,a/)') & 'Incoming data format, fmtverin=',fmtverin, & ', not found. The Job stopped.' CALL arpsstop('arpstop called from BINREADSPLIT. ',1) END IF WRITE(6,'(/1x,a,a/)') & 'Incoming data format, fmtverin=',fmtverin READ(inch,ERR=110,END=120) runname READ(inch,ERR=110,END=120) nocmnt IF( nocmnt > 0 ) THEN DO i=1,nocmnt READ(inch,ERR=110,END=120) cmnt(i) END DO END IF WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') runname IF( nocmnt > 0 ) THEN DO i=1,nocmnt WRITE(6,'(1x,a)') cmnt(i) END DO END IF READ(inch,ERR=110,END=120) time,tmunit ! !----------------------------------------------------------------------- ! ! Get dimensions of data in binary file and check against ! the dimensions passed to BINREADSPLIT ! !----------------------------------------------------------------------- ! IF (intver <= intver410) THEN READ(inch,ERR=110,END=120) nxin, nyin, nzin nzsoilin = 2 ! for version prior to 410, it is a two-layer soil model ELSE IF (intver >= intver500) THEN READ(inch,ERR=110,END=120) nxin, nyin, nzin,nzsoilin END IF ! ! Data validation: dimensions ! IF( nxin /= nxlg .OR. nyin /= nylg .OR. nzin /= nzlg) THEN WRITE(6,'(1x,a)') & ' Dimensions in BINREADSPLIT inconsistent with data.' WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin WRITE(6,'(1x,a,3I15)') ' Expected: ', nx, ny, nz WRITE(6,'(1x,a)') ' Program aborted in BINREADSPLIT.' CALL arpsstop('arpstop called from BINREADSPLIT nx-ny-nz read ',1) END IF IF(nzsoilin /= nzsoillg) THEN IF (intver <= intver410) THEN WRITE(6,'(1x,a,a/,2(1x,a/))') & ' The incoming data version is ', fmtverin, & ' In the input file, nzsoil must be set to 2. ', & ' Program aborted in BINREADSPLIT.' ELSE IF (intver >= intver500) THEN WRITE(6,'(1x,a)') & ' Dimensions in BINREADSPLIT inconsistent with data.' WRITE(6,'(1x,a,I15)') ' Read were: ', nzsoilin WRITE(6,'(1x,a,I15)') ' Expected: ', nzsoillg WRITE(6,'(1x,a)') ' Program aborted in BINREADSPLIT.' END IF CALL arpsstop('arpstop called from BINREADSPLIT nzsoil read ',1) END IF nxin = (nxin-3)/nproc_x+3 nyin = (nyin-3)/nproc_y+3 END IF CALL mpupdatei(intver, 1) CALL mpupdatec(runname, 40) CALL mpupdater(time, 1) CALL mpupdatec(tmunit, 10) CALL mpupdatei(nxin, 1) CALL mpupdatei(nyin, 1) CALL mpupdatei(nzin, 1) CALL mpupdatei(nzsoilin, 1) ! !----------------------------------------------------------------------- ! ! Read in flags for different data groups ! !----------------------------------------------------------------------- ! IF (myproc == 0) THEN IF( grdbas == 1 ) THEN ! Read grid and base state arrays WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)') & 'To read grid and base state data at time ', time, & ' secs = ',(time/60.),' mins.' READ(inch,ERR=110,END=120) & bgrdin,bbasin,bvarin,mstin,bicein, & btrbin,idummy,idummy,landin,totin, & btkein,idummy,idummy,mapproj,month, & day, year, hour, minute, second ELSE ! Normal data reading WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:', & time,' secs = ',(time/60.),' mins.' READ(inch,ERR=110,END=120) & grdin,basin,varin,mstin,icein, & trbin, sfcin,rainin,landin,totin, & tkein,idummy,idummy,mapproj, month, & day, year, hour, minute, second END IF READ(inch,ERR=110,END=120) & umove,vmove,xgrdorg,ygrdorg,trulat1, & trulat2,trulon,sclfct,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy, & tstop,thisdmp,latitud,ctrlat,ctrlon IF ( totin /= 0 ) THEN ! !----------------------------------------------------------------------- ! ! Read in additional parameters for ARPS history dump 4.0 or later ! version. ! !----------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) & nstyp1, prcin, radin, flxin,snowcin, & snowin,idummy,idummy,idummy,idummy, & idummy,idummy,idummy,idummy,idummy, & idummy,idummy,idummy,idummy,idummy IF ( nstyp1 < 1 ) THEN nstyp1 = 1 END IF READ(inch,ERR=110,END=120) & rdummy,rdummy,rdummy,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy END IF END IF ! myproc == 0 CALL mpupdatei(mstin,1) CALL mpupdatei(landin,1) CALL mpupdatei(totin,1) CALL mpupdatei(mapproj,1) CALL mpupdatei(month,1) CALL mpupdatei(day,1) CALL mpupdatei(year,1) CALL mpupdatei(hour,1) CALL mpupdatei(minute,1) CALL mpupdatei(second,1) IF(grdbas == 1) THEN CALL mpupdatei(bgrdin,1) CALL mpupdatei(bbasin,1) CALL mpupdatei(bvarin,1) CALL mpupdatei(btrbin,1) CALL mpupdatei(btkein,1) ELSE CALL mpupdatei(grdin,1) CALL mpupdatei(basin,1) CALL mpupdatei(varin,1) CALL mpupdatei(trbin,1) CALL mpupdatei(tkein,1) CALL mpupdatei(icein,1) CALL mpupdatei(sfcin,1) CALL mpupdatei(rainin,1) END IF CALL mpupdater(umove,1) CALL mpupdater(vmove,1) CALL mpupdater(xgrdorg,1) CALL mpupdater(ygrdorg,1) CALL mpupdater(trulat1,1) CALL mpupdater(trulat2,1) CALL mpupdater(trulon,1) CALL mpupdater(sclfct,1) CALL mpupdater(tstop,1) CALL mpupdater(thisdmp,1) CALL mpupdater(latitud,1) CALL mpupdater(ctrlat,1) CALL mpupdater(ctrlon,1) IF(totin /= 0) THEN CALL mpupdatei(nstyp1,1) CALL mpupdatei(prcin,1) CALL mpupdatei(radin,1) CALL mpupdatei(flxin,1) CALL mpupdatei(snowcin,1) CALL mpupdatei(snowin,1) END IF ! !----------------------------------------------------------------------- ! ! Read in x, y, z and zp arrays. ! !---------------------------------------------------------------------- ! IF( grdin == 1 .OR. grdbas == 1 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (var2d(i,1),i=1,nxlg) WRITE(lchanl,910) label,' x.' END IF CALL mpisplit1dx(var2d(:,1),nx,x) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (var2d(1,j),j=1,nylg) WRITE(lchanl,910) label,' y.' END IF CALL mpisplit1dy(var2d(1,:),ny,y) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) z WRITE(lchanl,910) label,' z.' END IF CALL mpupdater(z,nzlg) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' zp.' END IF CALL mpisplit3d(var3d,nx,ny,nz,zp) IF (intver <= intver410) THEN ! ! nzsoil must equal to 2, 06/28/2002, Zuwen ! assume zpsoil(,,2) is one meter below the surface. ! DO j=1,ny DO i=1,nx zpsoil(i,j,1)=zp(i,j,2) zpsoil(i,j,2)=zpsoil(i,j,1)-1. END DO END DO IF (myproc == 0) THEN WRITE(lchanl,'(1x,a)') ' Assign zpsoil. ' WRITE(lchanl,'(1x,a)') ' Assume zpsoil(,,1) is zp(:,:,2). ' WRITE(lchanl,'(1x,a)') ' Assume zpsoil(,,2) is zp(:,:,2)-1. ' END IF ELSE IF (intver >= intver500) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg) WRITE(lchanl,910) label,' zpsoil.' END IF CALL mpisplit3d(var3d,nx,ny,nzsoil,zpsoil) END IF ! intver END IF ! grdin ! !----------------------------------------------------------------------- ! ! Read in base state fields ! !---------------------------------------------------------------------- ! IF( basin == 1 .OR. grdbas == 1 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' ubar.' END IF CALL mpisplit3d(var3d,nx,ny,nz,ubar) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' vbar.' END IF CALL mpisplit3d(var3d,nx,ny,nz,vbar) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' wbar.' END IF CALL mpisplit3d(var3d,nx,ny,nz,wbar) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' ptbar.' END IF CALL mpisplit3d(var3d,nx,ny,nz,ptbar) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' pbar.' END IF CALL mpisplit3d(var3d,nx,ny,nz,pbar) IF( mstin == 1) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' qvbar.' END IF CALL mpisplit3d(var3d,nx,ny,nz,qvbar) END IF IF (landin == 1) THEN IF (nstyp1 <= 1) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((var3di(i,j,1),i=1,nxlg),j=1,nylg) WRITE(lchanl,910) label,' soiltyp.' END IF CALL mpisplit3di(var3di,nx,ny,1,soiltyp(:,:,1)) ELSE IF(myproc == 0) THEN var3di(:,:,:) = 0 var3d(:,:,:) = 0.0 DO is=1,nstyp1 IF (is <= nstyps) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((var3di(i,j,is),i=1,nxlg),j=1,nylg) WRITE(lchanl,910) label,' soiltyp.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((var3d(i,j,is),i=1,nxlg),j=1,nylg) WRITE(lchanl,910) label,' stypfrct.' ELSE READ(inch,ERR=110,END=120) label WRITE(lchanl,910) label,'skipping soiltyp' READ(inch,ERR=110,END=120) READ(inch,ERR=110,END=120) label WRITE(lchanl,910) label,'skipping stypfrct.' READ(inch,ERR=110,END=120) END IF END DO ENDIF CALL mpisplit3di(var3di,nx,ny,nstyps,soiltyp) CALL mpisplit3d(var3d,nx,ny,nstyps,stypfrct) END IF CALL fix_stypfrct_nstyp(nx,ny,nstyp1,nstyp,stypfrct) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((var3di(i,j,1),i=1,nxlg),j=1,nylg) WRITE(lchanl,910) label,' vegtyp.' END IF CALL mpisplit3di(var3di,nx,ny,1,vegtyp) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' lai.' END IF CALL mpisplit3d(var2d,nx,ny,1,lai) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' roufns.' END IF CALL mpisplit3d(var2d,nx,ny,1,roufns) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' veg.' END IF CALL mpisplit3d(var2d,nx,ny,1,veg) END IF END IF IF( grdbas == 1 ) GO TO 930 IF( varin == 1 ) THEN IF ( totin == 0 ) THEN ! !----------------------------------------------------------------------- ! ! Read in perturbations from history dump ! !----------------------------------------------------------------------- ! IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' uprt.' END IF CALL mpisplit3d(var3d,nx,ny,nz,uprt) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' vprt.' END IF CALL mpisplit3d(var3d,nx,ny,nz,vprt) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' wprt.' END IF CALL mpisplit3d(var3d,nx,ny,nz,wprt) ! !----------------------------------------------------------------------- ! ! Read in scalars ! !---------------------------------------------------------------------- ! IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' ptprt.' END IF CALL mpisplit3d(var3d,nx,ny,nz,ptprt) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' pprt.' END IF CALL mpisplit3d(var3d,nx,ny,nz,pprt) ELSE ! !----------------------------------------------------------------------- ! ! Read in total values of variables from history dump ! !---------------------------------------------------------------------- ! IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' u.' END IF CALL mpisplit3d(var3d,nx,ny,nz,uprt) DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx uprt(i,j,k) = uprt(i,j,k) - ubar(i,j,k) END DO END DO END DO IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' v.' END IF CALL mpisplit3d(var3d,nx,ny,nz,vprt) DO k=1,nz-1 DO j=1,ny DO i=1,nx-1 vprt(i,j,k) = vprt(i,j,k) - vbar(i,j,k) END DO END DO END DO IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' w.' END IF CALL mpisplit3d(var3d,nx,ny,nz,wprt) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' pt.' END IF CALL mpisplit3d(var3d,nx,ny,nz,ptprt) DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 ptprt(i,j,k) = ptprt(i,j,k) - ptbar(i,j,k) END DO END DO END DO IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' p.' END IF CALL mpisplit3d(var3d,nx,ny,nz,pprt) DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 pprt(i,j,k) = pprt(i,j,k) - pbar(i,j,k) END DO END DO END DO END IF END IF ! !----------------------------------------------------------------------- ! ! Read in moisture variables ! !----------------------------------------------------------------------- ! IF( mstin == 1 ) THEN IF ( totin == 0 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' qvprt.' END IF CALL mpisplit3d(var3d,nx,ny,nz,qvprt) ELSE IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' qv.' END IF CALL mpisplit3d(var3d,nx,ny,nz,qvprt) DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 qvprt(i,j,k) = qvprt(i,j,k) - qvbar(i,j,k) END DO END DO END DO END IF IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' qc.' END IF CALL mpisplit3d(var3d,nx,ny,nz,qc) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' qr.' END IF CALL mpisplit3d(var3d,nx,ny,nz,qr) IF( rainin == 1 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' raing.' END IF CALL mpisplit3d(var2d,nx,ny,1,raing) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' rainc.' END IF CALL mpisplit3d(var2d,nx,ny,1,rainc) END IF IF( prcin == 1 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((var3d(i,j,1),i=1,nxlg),j=1,nylg) WRITE(lchanl,910) label,' prcrate1.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((var3d(i,j,2),i=1,nxlg),j=1,nylg) WRITE(lchanl,910) label,' prcrate2.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((var3d(i,j,3),i=1,nxlg),j=1,nylg) WRITE(lchanl,910) label,' prcrate3.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((var3d(i,j,4),i=1,nxlg),j=1,nylg) WRITE(lchanl,910) label,' prcrate4.' END IF CALL mpisplit3d(var3d,nx,ny,4,prcrate) END IF IF( icein == 1 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' qi.' END IF CALL mpisplit3d(var3d,nx,ny,nz,qi) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' qs.' END IF CALL mpisplit3d(var3d,nx,ny,nz,qs) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' qh.' END IF CALL mpisplit3d(var3d,nx,ny,nz,qh) END IF END IF IF( tkein == 1 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' tke.' END IF CALL mpisplit3d(var3d,nx,ny,nz,tke) END IF IF( trbin == 1 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' kmh.' END IF CALL mpisplit3d(var3d,nx,ny,nz,kmh) IF ( intver >= intver410 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' kmv.' END IF CALL mpisplit3d(var3d,nx,ny,nz,kmv) END IF ! intver END IF IF( sfcin == 1 ) THEN IF (nstyp1 <= 1) THEN IF (intver <= intver410) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' tsfc.' END IF CALL mpisplit3d(var2d,nx,ny,1,tsoil(:,:,1,0)) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' tsoil.' END IF CALL mpisplit3d(var2d,nx,ny,1,tsoil(:,:,2,0)) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' wetsfc.' END IF CALL mpisplit3d(var2d,nx,ny,1,qsoil(:,:,1,0)) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' wetdp.' END IF CALL mpisplit3d(var2d,nx,ny,1,qsoil(:,:,2,0)) ELSE IF (intver >= intver500) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg) WRITE(lchanl,910) label,' tsoil.' END IF CALL mpisplit3d(var3d,nx,ny,nzsoil,tsoil(:,:,:,0)) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg) WRITE(lchanl,910) label,' qsoil.' END IF CALL mpisplit3d(var3d,nx,ny,nzsoil,qsoil(:,:,:,0)) END IF ! intver IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' wetcanp.' END IF CALL mpisplit3d(var2d,nx,ny,1,wetcanp(:,:,0)) ELSE DO is=0,nstyp1 IF (is <= nstyps) THEN IF (intver <= intver410) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' tsfc.' END IF CALL mpisplit3d(var2d,nx,ny,1,tsoil(:,:,1,is)) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' tsoil.' END IF CALL mpisplit3d(var2d,nx,ny,1,tsoil(:,:,2,is)) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' wetsfc.' END IF CALL mpisplit3d(var2d,nx,ny,1,qsoil(:,:,1,is)) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' wetdp.' END IF CALL mpisplit3d(var2d,nx,ny,1,qsoil(:,:,2,is)) ELSE IF (intver >= intver500) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg) WRITE(lchanl,910) label,' tsoil.' END IF CALL mpisplit3d(var3d,nx,ny,nzsoil,tsoil(:,:,:,is)) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg) WRITE(lchanl,910) label,' qsoil.' END IF CALL mpisplit3d(var3d,nx,ny,nzsoil,qsoil(:,:,:,is)) END IF ! intver IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' wetcanp.' END IF CALL mpisplit3d(var2d,nx,ny,1,wetcanp(:,:,is)) ELSE IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,'skipping tsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,'skipping qsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,'skipping wetcanp.' END IF ENDIF END DO END IF CALL fix_soil_nstyp(nx,ny,nzsoil,nstyp1,nstyp,tsoil,qsoil,wetcanp) IF (snowcin == 1) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,' snowcvr -- discarding.' END IF END IF IF (snowin == 1) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' snowdpth.' END IF CALL mpisplit3d(var2d,nx,ny,1,snowdpth) END IF END IF IF( radin == 1 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg) WRITE(lchanl,910) label,' radfrc.' END IF CALL mpisplit3d(var3d,nx,ny,nz,radfrc) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' radsw.' END IF CALL mpisplit3d(var2d,nx,ny,1,radsw) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' rnflx.' END IF CALL mpisplit3d(var2d,nx,ny,1,rnflx) IF (intver <= intver410) THEN ! ! Do not know how to initialized radswnet, radlwin, yet, ! at this moment. ! radswnet = 0. radlwin = 0. ELSE IF (intver >= intver500) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' radswnet.' END IF CALL mpisplit3d(var2d,nx,ny,1,radswnet) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' radlwin.' END IF CALL mpisplit3d(var2d,nx,ny,1,radlwin) END IF ! intver END IF IF( flxin == 1 ) THEN IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' usflx.' END IF CALL mpisplit3d(var2d,nx,ny,1,usflx) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' vsflx.' END IF CALL mpisplit3d(var2d,nx,ny,1,vsflx) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' ptsflx.' END IF CALL mpisplit3d(var2d,nx,ny,1,ptsflx) IF(myproc == 0) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) var2d WRITE(lchanl,910) label,' qvsflx.' END IF CALL mpisplit3d(var2d,nx,ny,1,qvsflx) END IF DEALLOCATE(var2d, var3d, var3di) 910 FORMAT(1X,'Field ',a12,' was read into array',a) ! !----------------------------------------------------------------------- ! ! Friendly exit message ! !---------------------------------------------------------------------- ! 930 CONTINUE IF (myproc == 0) & WRITE(6,'(/a,F8.1,a/)') & ' Data at time=', time/60,' (min) were successfully read.' ireturn = 0 RETURN ! !----------------------------------------------------------------------- ! ! Error during read ! !---------------------------------------------------------------------- ! 110 CONTINUE WRITE(6,'(/a/)') ' Error reading data in BINREADSPLIT' ireturn=1 RETURN ! !----------------------------------------------------------------------- ! ! End-of-file during read ! !---------------------------------------------------------------------- ! 120 CONTINUE WRITE(6,'(/a/)') ' End of file reached in BINREADSPLIT' ireturn=2 RETURN END SUBROUTINE binreadsplit ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE BN2READ ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE bn2read(nx,ny,nz,nzsoil,nstyps,grdbas,inch,time,x,y,z,zp, & 2,4 zpsoil,uprt, vprt, wprt, ptprt, pprt, & qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, & ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, & soiltyp,stypfrct,vegtyp,lai,roufns,veg, & tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & ireturn) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Read in binary data set created by ARPS using history dump format ! No.2. ! All data read in are located at the original staggered grid points ! !----------------------------------------------------------------------- ! ! AUTHOR: Ming Xue ! 2/27/92. ! ! MODIFICATION HISTORY: ! ! 6/08/92 Added full documentation (K. Brewster) ! ! 7/14/92 (K. Brewster) ! Added runname, comment and version number reading ! ! 8/20/92 (M. Xue) ! Added data reading of computational z coordinate array z. ! ! 4/23/93 (M. Xue) ! New data format. ! ! 02/06/95 (Y. Liu) ! Added map projection parameters into the second binary dumping ! ! 12/09/1998 (Donghai Wang) ! Added the snow cover. ! ! 05/15/2002 (J. Brotzge) ! Added variables to allow for multiple soil schemes. ! !----------------------------------------------------------------------- ! ! INPUT : ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! nzsoil Number of grid points in the soil ! ! grdbas Data read flag. ! =1, only grid and base state arrays will be read ! =0, all arrays will be read based on data ! parameter setting. ! inch Channel number for binary reading. ! This channel must be opened for unformatted reading ! by the calling routine. ! ! OUTPUT: ! ! time Time in seconds of data read from "filename" ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp z coordinate of grid points in physical space (m) ! zpsoil z coordinate of grid points in the soil (m) ! ! uprt x component of perturbation velocity (m/s) ! vprt y component of perturbation velocity (m/s) ! wprt Vertical component of perturbation velocity in ! Cartesian coordinates (m/s). ! ptprt Perturbation potential temperature (K) ! pprt Perturbation pressure (Pascal) ! ! qvprt Perturbation water vapor mixing ratio (kg/kg) ! qc Cloud water mixing ratio (kg/kg) ! qr Rainwater mixing ratio (kg/kg) ! qi Cloud ice mixing ratio (kg/kg) ! qs Snow mixing ratio (kg/kg) ! qh Hail mixing ratio (kg/kg) ! tke Turbulent Kinetic Energy ((m/s)**2) ! ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! ! ubar Base state x velocity component (m/s) ! vbar Base state y velocity component (m/s) ! wbar Base state z velocity component (m/s) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! rhobar Base state air density (kg/m**3) ! qvbar Base state water vapor mixing ratio (kg/kg) ! ! soiltyp Soil type ! stypfrct Soil type fraction ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil Soil temperature (K) ! qsoil Soil moisture (m**3/m**3) ! wetcanp Canopy water amount ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! prcrate Precipitation rates ! ! radfrc Radiation forcing (K/s) ! radsw Solar radiation reaching the surface ! rnflx Net radiation flux absorbed by surface ! radswnet Net shortwave radiation ! radlwin Incoming longwave radiation ! ! usflx Surface flux of u-momentum (kg/(m*s**2)) ! vsflx Surface flux of v-momentum (kg/(m*s**2)) ! ptsflx Surface heat flux (K*kg/(m**2 * s )) ! qvsflx Surface moisture flux of (kg/(m**2 * s)) ! ! ireturn Return status indicator ! =0 successful read of all data ! =1 error reading data ! =2 end-of-file reached during read attempt ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Number of grid points in 3 directions INTEGER :: nzsoil ! Number of grid points in 3 directions REAL :: x (nx) ! x-coord. of the physical and compu ! -tational grid. Defined at u-point(m). REAL :: y (ny) ! y-coord. of the physical and compu ! -tational grid. Defined at v-point(m). REAL :: z (nz) ! z-coord. of the computational grid. ! Defined at w-point on the staggered ! grid(m). REAL :: zp (nx,ny,nz) ! Physical height coordinate defined at ! w-point of the staggered grid(m). REAL :: zpsoil(nx,ny,nzsoil) ! Physical height coordinate defined at ! w-point of the soil (m) INTEGER :: grdbas ! Data read flag. INTEGER :: inch ! Channel number for binary reading REAL :: time ! Time in seconds of data read ! from "filename" REAL :: uprt (nx,ny,nz) ! Perturbation u-velocity (m/s) REAL :: vprt (nx,ny,nz) ! Perturbation v-velocity (m/s) REAL :: wprt (nx,ny,nz) ! Perturbation w-velocity (m/s) REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K) REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal) REAL :: qvprt (nx,ny,nz) ! Perturbation water vapor mixing ! ratio (kg/kg) REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg) REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s) REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s) REAL :: wbar (nx,ny,nz) ! Base state w-velocity (m/s) REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K) REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal) REAL :: rhobar(nx,ny,nz) ! Base state air density (kg/m**3) REAL :: qvbar (nx,ny,nz) ! Base state water vapor mixing ratio INTEGER :: nstyps ! Number of soil type INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type REAL :: stypfrct(nx,ny,nstyps) ! Soil type INTEGER :: vegtyp(nx,ny) ! Vegetation type REAL :: lai (nx,ny) ! Leaf Area Index REAL :: roufns (nx,ny) ! Surface roughness REAL :: veg (nx,ny) ! Vegetation fraction REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K) REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount REAL :: snowdpth(nx,ny) ! Snow depth (m) REAL :: raing(nx,ny) ! Grid supersaturation rain REAL :: rainc(nx,ny) ! Cumulus convective rain REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s)) ! prcrate(1,1,1) = total precip. rate ! prcrate(1,1,2) = grid scale precip. rate ! prcrate(1,1,3) = cumulus precip. rate ! prcrate(1,1,4) = microphysics precip. rate REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s) REAL :: radsw (nx,ny) ! Solar radiation reaching the surface REAL :: rnflx (nx,ny) ! Net radiation flux absorbed by surface REAL :: radswnet(nx,ny) ! Net shortwave radiation REAL :: radlwin(nx,ny) ! Incoming longwave radiation REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2)) REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2)) REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m**2*s)) REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s)) INTEGER :: ireturn ! Return status indicator ! !----------------------------------------------------------------------- ! ! Parameters describing routine that wrote the gridded data ! !----------------------------------------------------------------------- ! CHARACTER (LEN=40) :: fmtver0,fmtver1,fmtverin PARAMETER (fmtver0='004.10 2nd Binary Data') PARAMETER (fmtver1='004.10 2nd Binary Data') CHARACTER (LEN=10) :: tmunit ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! INTEGER :: lchanl PARAMETER (lchanl=6) ! Channel number for formatted printing. INTEGER :: i,j,k,is INTEGER :: nstyp1 CHARACTER (LEN=12) :: label ! INTEGER :: nxin,nyin,nzin INTEGER :: nxin,nyin,nzin,nzsoilin INTEGER :: bgrdin,bbasin,bvarin,bicein,btkein,btrbin,idummy REAL :: rdummy ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'indtflg.inc' INCLUDE 'globcst.inc' INCLUDE 'grid.inc' ! Grid & map parameters. ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !----------------------------------------------------------------------- ! ! Read header info ! !----------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) fmtverin IF( fmtverin /= fmtver0 .AND. fmtverin /= fmtver1 ) THEN WRITE(6,'(/1x,a/1x,2a/1x,2a/1x,2a/1x,a)') & 'Data format incompatible with the data reader.', & 'Format of data is ',fmtverin,' Format of reader is ',fmtver1, & 'compitable to ',fmtver0, '. Job stopped.' CALL arpsstop('arpstop called from bn2read header read ',1) END IF READ(inch,ERR=110,END=120) runname READ(inch,ERR=110,END=120) nocmnt IF( nocmnt > 0 ) THEN DO i=1,nocmnt READ(inch,ERR=110,END=120) cmnt(i) END DO END IF WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') runname IF( nocmnt > 0 ) THEN DO i=1,nocmnt WRITE(6,'(1x,a)') cmnt(i) END DO END IF READ(inch,ERR=110,END=120) time,tmunit ! !----------------------------------------------------------------------- ! ! Get dimensions of data in binary file and check against ! the dimensions passed to BN2READ ! !----------------------------------------------------------------------- ! ! READ(inch,ERR=110,END=120) nxin, nyin, nzin READ(inch,ERR=110,END=120) nxin, nyin, nzin, nzsoilin ! IF( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN IF( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz .OR. nzsoilin /= nzsoil) THEN WRITE(6,'(1x,a)') & ' Dimensions in BN2READ inconsistent with data.' ! WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin WRITE(6,'(1x,a,4I15)') ' Read were: ', nxin, nyin, nzin, nzsoilin WRITE(6,'(1x,a)') & ' Program aborted in BN2READ.' ! CALL arpsstop('arpstop called from bn2read nx-ny-nz read',1) CALL arpsstop('arpstop called from bn2read nx-ny-nz-nzsoil read',1) END IF ! !----------------------------------------------------------------------- ! ! Read in flags for different data groups ! !----------------------------------------------------------------------- ! IF( grdbas == 1 ) THEN ! Read grid and base state arrays WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)') & 'To read grid and base state data at time ', time, & ' secs = ',(time/60.),' mins.' READ(inch,ERR=110,END=120) & bgrdin,bbasin,bvarin,mstin,bicein, & btrbin,idummy,idummy,landin,totin, & btkein,idummy,idummy,mapproj,month, & day,year,hour,minute,second ELSE ! Normal data reading WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:', & time,' secs = ',(time/60.),' mins.' READ(inch,ERR=110,END=120) & grdin,basin,varin,mstin,icein, & trbin,sfcin,rainin,landin,totin, & tkein,idummy,idummy,mapproj,month, & day,year,hour,minute,second END IF READ(inch,ERR=110,END=120) & umove,vmove,xgrdorg,ygrdorg,trulat1, & trulat2,trulon,sclfct,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy, & tstop,thisdmp,latitud,ctrlat,ctrlon IF ( totin /= 0 ) THEN ! !----------------------------------------------------------------------- ! ! Read in additional parameters for ARPS history dump 4.1 or later ! version. ! !----------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) & nstyp1, prcin, radin, flxin,snowcin, & snowin,idummy,idummy,idummy,idummy, & idummy,idummy,idummy,idummy,idummy, & idummy,idummy,idummy,idummy,idummy IF ( nstyp1 < 1 ) THEN nstyp1 = 1 END IF READ(inch,ERR=110,END=120) & rdummy,rdummy,rdummy,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy, & rdummy,rdummy,rdummy,rdummy,rdummy END IF ! !----------------------------------------------------------------------- ! ! Read in x,y and z at grid cell centers (scalar points). ! !---------------------------------------------------------------------- ! IF( grdin == 1 .OR. grdbas == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (x(i),i=1,nx) WRITE(lchanl,910) label,' x.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (y(j),j=1,ny) WRITE(lchanl,910) label,' y.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (z(k),k=1,nz) WRITE(lchanl,910) label,' z.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((zp(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' zp.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((zpsoil(i,j,k),i=1,nx),j=1,ny),k=1,nzsoil) WRITE(lchanl,910) label,' zpsoil.' END IF ! grdin ! !----------------------------------------------------------------------- ! ! Read in base state fields ! !---------------------------------------------------------------------- ! IF( basin == 1 .OR. grdbas == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((ubar(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' ubar.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((vbar(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' vbar.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((wbar(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' wbar.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((ptbar(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' ptbar.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((pbar(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' pbar.' IF( mstin == 1) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qvbar(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' qvbar.' END IF IF(landin == 1) THEN IF( nstyp1 <= 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((soiltyp(i,j,1),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' soiltyp.' ELSE DO is=1,nstyp1 IF (is <= nstyps) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((soiltyp(i,j,is),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' soiltyp.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((stypfrct(i,j,is),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' stypfrct.' ELSE READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,'skipping soiltyp.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,'skipping stypfrct.' ENDIF END DO END IF CALL fix_stypfrct_nstyp(nx,ny,nstyp1,nstyp,stypfrct) READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((vegtyp (i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' vegtyp.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((lai (i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' lai.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((roufns (i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' roufns.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((veg (i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' veg.' END IF END IF IF( grdbas == 1 ) GO TO 930 IF( varin == 1 ) THEN IF ( totin == 0 ) THEN ! !----------------------------------------------------------------------- ! ! Read in uprt, vprt, and wprt ! !---------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((uprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' uprt.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((vprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' vprt.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((wprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' wprt.' ! !----------------------------------------------------------------------- ! ! Read in scalars ! !---------------------------------------------------------------------- ! READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((ptprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' ptprt.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((pprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' pprt.' ELSE READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((uprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' u.' DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx uprt(i,j,k) = uprt(i,j,k) - ubar(i,j,k) END DO END DO END DO READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((vprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' v.' DO k=1,nz-1 DO j=1,ny DO i=1,nx-1 vprt(i,j,k) = vprt(i,j,k) - vbar(i,j,k) END DO END DO END DO READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((wprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' w.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((ptprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' pt.' DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 ptprt(i,j,k) = ptprt(i,j,k) - ptbar(i,j,k) END DO END DO END DO READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((pprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' p.' DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 pprt(i,j,k) = pprt(i,j,k) - pbar(i,j,k) END DO END DO END DO END IF END IF ! !----------------------------------------------------------------------- ! ! Read in moisture variables ! !---------------------------------------------------------------------- ! IF( mstin == 1 ) THEN IF ( totin == 0 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qvprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' qvprt.' ELSE READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qvprt(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' qv.' DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 qvprt(i,j,k) = qvprt(i,j,k) - qvbar(i,j,k) END DO END DO END DO END IF READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qc(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' qc.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qr(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' qr.' IF( rainin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((raing(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' raing.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((rainc(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' rainc.' END IF IF( prcin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((prcrate(i,j,1),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' prcrate1.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((prcrate(i,j,2),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' prcrate2.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((prcrate(i,j,3),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' prcrate3.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & ((prcrate(i,j,4),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' prcrate4.' END IF IF( icein == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qi(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' qi.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qs(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' qs.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((qh(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' qh.' END IF END IF IF( tkein == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((tke(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' tke.' END IF IF( trbin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((kmh(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' kmh.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((kmv(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' kmv.' END IF IF( sfcin == 1 ) THEN IF (nstyp1 <= 1) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((tsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) WRITE(lchanl,910) label,' tsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) (((qsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) WRITE(lchanl,910) label,' qsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((wetcanp(i,j,0),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' wetcanp.' ELSE DO is=0,nstyp1 IF (is <= nstyps) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120)(((tsoil(i,j,k,is),i=1,nx),j=1,ny), & k=1,nzsoil) WRITE(lchanl,910) label,' tsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120)(((qsoil(i,j,k,is),i=1,nx),j=1,ny), & k=1,nzsoil) WRITE(lchanl,910) label,' qsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120)((wetcanp(i,j,is),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' wetcanp.' ELSE READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,'skipping tsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,'skipping qsoil.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,'skipping wetcanp.' ENDIF END DO END IF CALL fix_soil_nstyp(nx,ny,nzsoil,nstyp1,nstyp,tsoil,qsoil,wetcanp) IF(snowcin == 1) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) WRITE(lchanl,910) label,' snowcvr -- discarding.' END IF IF(snowin == 1) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120)((snowdpth(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' snowdpth.' END IF END IF IF( radin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) & (((radfrc(i,j,k),i=1,nx),j=1,ny),k=1,nz) WRITE(lchanl,910) label,' radfrc.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((radsw(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' radsw.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((rnflx(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' rnflx.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((radswnet(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' radswnet.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((radlwin(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' radlwin.' END IF IF( flxin == 1 ) THEN READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((usflx(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' usflx.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((vsflx(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' vsflx.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((ptsflx(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' ptsflx.' READ(inch,ERR=110,END=120) label READ(inch,ERR=110,END=120) ((qvsflx(i,j),i=1,nx),j=1,ny) WRITE(lchanl,910) label,' qvsflx.' END IF 910 FORMAT(1X,'Field ',a12,' was read into array',a) ! !----------------------------------------------------------------------- ! ! Friendly exit message ! !---------------------------------------------------------------------- ! 930 CONTINUE WRITE(6,'(/a,F8.1,a/)') & ' Data at time=', time/60,' (min) were successfully read.' ireturn = 0 RETURN ! !----------------------------------------------------------------------- ! ! Error during read ! !---------------------------------------------------------------------- ! 110 CONTINUE WRITE(6,'(/a/)') ' Error reading data in BN2READ' ireturn=1 RETURN ! !----------------------------------------------------------------------- ! ! End-of-file during read ! !---------------------------------------------------------------------- ! 120 CONTINUE WRITE(6,'(/a/)') ' End of file reached in BN2READ' ireturn=2 RETURN END SUBROUTINE bn2read ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE BINDUMP ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE bindump(nx,ny,nz,nzsoil,nstyps, nchanl, grdbas, & 1,54 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv, & ubar,vbar,ptbar,pbar,rhobar,qvbar, & x,y,z,zp,zpsoil, & soiltyp,stypfrct,vegtyp,lai,roufns,veg, & tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & tem1) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Write history data into channel nchanl as binary data. ! ! All data read in are located at the original staggered grid points. ! ! Note: coordinate fields are dumped as 3 dimensional fields which ! have been converted from meters to kilometers. This is for the ! convenience of the plotting applications. ! ! The last 4 characters of the 12 character label written out ! with each 1-,2-, or 3-d array is used by the splitdump and ! joinfiles subroutines used by the message passing version of the ! ARPS (and also by some auxiliary ARPS I/O routines) ! to determine the data type of the array. ! Key to the labels: ! ! 'nnnnnnn tdds' ! ! n - characters containing the name of the variable. ! t - type of variable: "r" for real and "i" for integer. ! dd - number of dimensions: "1d" "2d" or "3d". ! s - staggered dimension: "0" for centered, ! "1" for staggered in x, ! "2" for staggered in y, ! "3" for staggered in z. ! !----------------------------------------------------------------------- ! ! ! !----------------------------------------------------------------------- ! ! AUTHOR: Ming Xue ! 3/10/92. ! ! MODIFICATION HISTORY: ! ! 6/06/92 (M. Xue) ! Added full documentation. ! ! 7/13/92 (K. Brewster) ! Added runname, comment and version number writing ! ! 8/23/92 (M. Xue) ! Modify to perform the dumping of both base and t-dependent arrays ! and added control on grid staggering. ! ! 4/4/93 (M. Xue) ! Modified, so that data on the original staggered grid are written ! out. Averaging to the volume center is no longer done. ! ! 9/1/94 (Y. Lu) ! Cleaned up documentation. ! ! 02/06/95 (Y. Liu) ! Added map projection parameters into the binary dumping ! ! 03/26/96 (G. Bassett) ! Labels were modified to include information about array type. ! This information is used by splitdump and joinfiles subroutines. ! ! 12/09/1998 (Donghai Wang) ! Added the snow cover. ! ! 05/15/2002 (J. Brotzge) ! Added to allow for multiple soil schemes. ! !----------------------------------------------------------------------- ! ! INPUT: ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! nzsoil Number of grid points in the soil ! ! nchanl FORTRAN I/O channel number for history data output. ! grdbas Flag indicating if this is a call for the data dump ! of grid and base state arrays only. If so, grdbas=1. ! ! u x component of velocity at a given time level (m/s) ! v y component of velocity at a given time level (m/s) ! w Vertical component of Cartesian velocity at a given ! time level (m/s) ! ptprt Perturbation potential temperature at a given time ! level (K) ! pprt Perturbation pressure at a given time level (Pascal) ! qv Water vapor specific humidity at a given time level (kg/kg) ! qc Cloud water mixing ratio at a given time level (kg/kg) ! qr Rainwater mixing ratio at a given time level (kg/kg) ! qi Cloud ice mixing ratio at a given time level (kg/kg) ! qs Snow mixing ratio at a given time level (kg/kg) ! qh Hail mixing ratio at a given time level (kg/kg) ! tke Turbulent Kinetic Energy ((m/s)**2) ! ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! ! ubar Base state zonal velocity component (m/s) ! vbar Base state meridional velocity component (m/s) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! rhobar Base state density (kg/m**3) ! qvbar Base state water vapor specific humidity (kg/kg) ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp Vertical coordinate of grid points in physical space (m) ! zpsoil Vertical coordinate of grid points in the soil (m) ! ! soiltyp Soil type ! stypfrct Soil type fraction ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil Soil temperature (K) ! qsoil Soil moisture (m**3/m**3) ! wetcanp Canopy water amount ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! prcrate Precipitation rates ! ! radfrc Radiation forcing (K/s) ! radsw Solar radiation reaching the surface ! rnflx Net radiation flux absorbed by surface ! radswnet Net shortwave radiation ! radlwin Incoming longwave radiation ! ! usflx Surface flux of u-momentum (kg/(m*s**2)) ! vsflx Surface flux of v-momentum (kg/(m*s**2)) ! ptsflx Surface heat flux (K*kg/(m**2*s)) ! qvsflx Surface moisture flux (kg/(m**2*s)) ! ! OUTPUT: ! ! None. ! ! WORK ARRAY: ! ! tem1 Temporary work array. ! ! !----------------------------------------------------------------------- ! ! The following parameters are passed into this subroutine through ! a common block in globcst.inc, and they determine which ! variables are output. ! ! grdout =0 or 1. If grdout=0, grid variables are not dumped. ! basout =0 or 1. If basout=0, base state variables are not dumped. ! varout =0 or 1. If varout=0, model perturbation variables are not dumped. ! mstout =0 or 1. If mstout=0, water variables are not dumped. ! rainout=0 or 1. If rainout=0, rain variables are not dumped. ! prcout =0 or 1. If prcout=0, precipitation rates are not dumped. ! iceout =0 or 1. If iceout=0, qi, qs and qh are not dumped. ! tkeout =0 or 1. If tkeout=0, tke is not dumped. ! trbout =0 or 1. If trbout=0, turbulence parameter km is not dumped. ! sfcout =0 or 1. If sfcout=0, surface variables are not dumped. ! landout=0 or 1. If landout=0, surface propertty arrays are not dumped. ! radout =0 or 1. If radout =0, radiation arrays are not dumped. ! flxout =0 or 1. If flxout =0, surface flux arrays are not dumped. ! ! These following parameters are also passed in through common ! blocks in globcst.inc. ! ! runname,curtim,umove,vmove,xgrdorg,ygrdorg ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Number of grid points in 3 directions INTEGER :: nzsoil ! Number of grid points in the soil INTEGER :: nchanl ! FORTRAN I/O channel number for output INTEGER :: grdbas ! If this is a grid/base state array dump REAL :: u (nx,ny,nz) ! Total u-velocity (m/s) REAL :: v (nx,ny,nz) ! Total v-velocity (m/s) REAL :: w (nx,ny,nz) ! Total w-velocity (m/s) REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K) REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal) REAL :: qv (nx,ny,nz) ! Water vapor specific humidity (kg/kg) REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg) REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s) REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s) REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K) REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal) REAL :: rhobar(nx,ny,nz) ! Base state air density (kg/m**3) REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific humidity ! (kg/kg) REAL :: x (nx) ! The x-coord. of the physical and ! computational grid. Defined at u-point. REAL :: y (ny) ! The y-coord. of the physical and ! computational grid. Defined at v-point. REAL :: z (nz) ! The z-coord. of the computational grid. ! Defined at w-point on the staggered grid. REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at ! w-point of the staggered grid. REAL :: zpsoil(nx,ny,nzsoil) ! The physical height coordinate defined at ! w-point of the soil INTEGER :: nstyps ! Number of soil types INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type REAL :: stypfrct(nx,ny,nstyps) ! Soil type fractions INTEGER :: vegtyp (nx,ny) ! Vegetation type REAL :: lai (nx,ny) ! Leaf Area Index REAL :: roufns (nx,ny) ! Surface roughness REAL :: veg (nx,ny) ! Vegetation fraction REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K) REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount REAL :: snowdpth(nx,ny) ! Snow depth (m) REAL :: raing(nx,ny) ! Grid supersaturation rain REAL :: rainc(nx,ny) ! Cumulus convective rain REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s)) ! prcrate(1,1,1) = total precip. rate ! prcrate(1,1,2) = grid scale precip. rate ! prcrate(1,1,3) = cumulus precip. rate ! prcrate(1,1,4) = microphysics precip. rate REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s) REAL :: radsw (nx,ny) ! Solar radiation reaching the surface REAL :: rnflx (nx,ny) ! Net radiation flux absorbed by surface REAL :: radswnet(nx,ny) ! Net shortwave radiation REAL :: radlwin(nx,ny) ! Incoming longwave radiation REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2)) REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2)) REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m**2*s)) REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s)) REAL :: tem1 (nx,ny,nz) ! Temporary work array ! !----------------------------------------------------------------------- ! ! Parameters describing this routine ! !----------------------------------------------------------------------- ! ! 06/28/2002 Zuwen He ! ! fmtver??: to label each data a version. ! intver??: an integer to allow faster comparison than fmtver??, ! which are strings. ! ! Verion 5.00: significant change in soil variables since version 4.10. ! CHARACTER (LEN=40) :: fmtver,fmtver410,fmtver500 INTEGER :: intver,intver410,intver500 PARAMETER (fmtver410='004.10 Binary Data',intver410=410) PARAMETER (fmtver500='005.00 Binary Data',intver500=500) CHARACTER (LEN=10) :: tmunit PARAMETER (tmunit='seconds ') ! !----------------------------------------------------------------------- ! ! Misc. local variables: ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k,l,is INTEGER :: idummy REAL :: rdummy ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'grid.inc' ! Grid & map parameters. INCLUDE 'mp.inc' ! mpi parameters. ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! intver = intver500 ! for the time being, in the future, we will ! allow to dump data in the different version ! intver will be assigned from input file IF (intver == intver410) THEN fmtver=fmtver410 ELSE IF (intver == intver500) THEN fmtver=fmtver500 ELSE IF (myproc == 0) WRITE(6,'(/1x,a,i10,a/)') & 'Data format, intver=',intver,', not found. The Job stopped.' CALL arpsstop('arpstop called from bindump.',1) END IF IF (myproc == 0) WRITE(6,'(/1x,a,a,a/)') & 'Data dump format, fmtver=',fmtver,'. ' IF (myproc == 0) & WRITE(6,'(1x,a,f13.3/)') 'Writing history data at time=', curtim ! !----------------------------------------------------------------------- ! ! Write header info ! !----------------------------------------------------------------------- ! WRITE(nchanl) fmtver WRITE(nchanl) runname WRITE(nchanl) nocmnt IF( nocmnt > 0 ) THEN DO l=1,nocmnt WRITE(nchanl) cmnt(l) END DO END IF WRITE(nchanl) curtim,tmunit IF (intver == intver410) THEN WRITE(nchanl) nx,ny,nz ELSE IF (intver == intver500) THEN WRITE(nchanl) nx,ny,nz,nzsoil END IF ! !----------------------------------------------------------------------- ! ! Write the flags for different data groups. ! !----------------------------------------------------------------------- ! idummy = 0 IF( grdbas == 1 ) THEN WRITE(nchanl) 1, 1, 0, mstout, 0, & 0, 0, 0, landout,totout, & 0, idummy, idummy, mapproj, month, & day, year, hour, minute, second ELSE WRITE(nchanl) grdout, basout, varout, mstout, iceout, & trbout, sfcout, rainout,landout,totout, & tkeout, idummy, idummy, mapproj, month, & day, year, hour, minute, second END IF rdummy = 0.0 WRITE(nchanl) umove, vmove, xgrdorg, ygrdorg, trulat1, & trulat2, trulon, sclfct, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy, & tstop, thisdmp, latitud, ctrlat, ctrlon ! !----------------------------------------------------------------------- ! ! If totout=1, write additional parameters to history dump files. ! This is for ARPS version 4.1.2 or later. ! !----------------------------------------------------------------------- ! IF ( totout == 1 ) THEN WRITE(nchanl) nstyp, prcout, radout, flxout, 0, & ! 0 for snowcvr snowout,idummy, idummy, idummy, idummy, & idummy, idummy, idummy, idummy, idummy, & idummy, idummy, idummy, idummy, idummy WRITE(nchanl) rdummy, rdummy, rdummy, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy END IF ! !----------------------------------------------------------------------- ! ! If grdout=1 or grdbas=1, write out grid variables ! !----------------------------------------------------------------------- ! IF(grdout == 1 .OR. grdbas == 1 ) THEN WRITE(nchanl) 'x coord r1d1' WRITE(nchanl) x WRITE(nchanl) 'y coord r1d2' WRITE(nchanl) y WRITE(nchanl) 'z coord r1d3' WRITE(nchanl) z CALL edgfill(zp,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz) WRITE(nchanl) 'zp coor r3d0' WRITE(nchanl) zp IF (intver == intver500) THEN CALL edgfill(zpsoil,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nzsoil,1,nzsoil) WRITE(nchanl) 'zpsoil rs3d0' WRITE(nchanl) zpsoil END IF END IF ! grdout ! !----------------------------------------------------------------------- ! ! If basout=1, write out base state variables. ! !----------------------------------------------------------------------- ! IF(basout == 1 .OR. grdbas == 1 ) THEN CALL edgfill(ubar,1,nx,1,nx, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'ubar r3d1' WRITE(nchanl) ubar CALL edgfill(vbar,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1) WRITE(nchanl) 'vbar r3d2' WRITE(nchanl) vbar DO k=1,nz DO j=1,ny DO i=1,nx tem1(i,j,k) = 0.0 END DO END DO END DO WRITE(nchanl) 'wbar r3d3' WRITE(nchanl) tem1 CALL edgfill(ptbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'ptbar r3d0' WRITE(nchanl) ptbar CALL edgfill(pbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'pbar r3d0' WRITE(nchanl) pbar IF(mstout == 1) THEN CALL edgfill(qvbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'qvbar r3d0' WRITE(nchanl) qvbar END IF IF(landout == 1) THEN IF( nstyp <= 1 ) THEN CALL iedgfill(soiltyp(1,1,1),1,nx,1,nx-1, 1,ny,1,ny-1, & 1,1,1,1) WRITE(nchanl) 'soiltyp i2d0' WRITE(nchanl) ((soiltyp(i,j,1),i=1,nx),j=1,ny) ELSE DO is=1,nstyp CALL iedgfill(soiltyp(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1, & 1,1,1,1) WRITE(nchanl) 'soiltyp i2d0' WRITE(nchanl) ((soiltyp(i,j,is),i=1,nx),j=1,ny) CALL edgfill(stypfrct(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1, & 1,1,1,1) WRITE(nchanl) 'stypfrc r2d0' WRITE(nchanl) ((stypfrct(i,j,is),i=1,nx),j=1,ny) END DO END IF CALL iedgfill(vegtyp ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'vegtyp i2d0' WRITE(nchanl) vegtyp CALL edgfill(lai ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'lai r2d0' WRITE(nchanl) lai CALL edgfill(roufns ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'roufns r2d0' WRITE(nchanl) roufns CALL edgfill(veg ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'veg r2d0' WRITE(nchanl) veg END IF END IF IF ( grdbas == 1 ) RETURN ! !----------------------------------------------------------------------- ! ! If varout = 1, Write out uprt, vprt, wprt, ptprt, pprt. ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Write out u,v and w. ! !----------------------------------------------------------------------- ! IF(varout == 1) THEN IF ( totout == 0 ) THEN ! !----------------------------------------------------------------------- ! ! Write out perturbations to history dump ! !----------------------------------------------------------------------- ! DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx tem1(i,j,k)=u(i,j,k)-ubar(i,j,k) END DO END DO END DO CALL edgfill(tem1,1,nx,1,nx, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'uprt r3d1' WRITE(nchanl) tem1 DO k=1,nz-1 DO i=1,nx-1 DO j=1,ny tem1(i,j,k)=v(i,j,k)-vbar(i,j,k) END DO END DO END DO CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1) WRITE(nchanl) 'vprt r3d2' WRITE(nchanl) tem1 CALL edgfill(w,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz) WRITE(nchanl) 'wprt r3d3' WRITE(nchanl) w ! !----------------------------------------------------------------------- ! ! Write out scalars ! !----------------------------------------------------------------------- ! CALL edgfill(ptprt,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'ptprt r3d0' WRITE(nchanl) ptprt CALL edgfill(pprt,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'pprt r3d0' WRITE(nchanl) pprt ELSE ! !----------------------------------------------------------------------- ! ! Write out total values to history dump ! !----------------------------------------------------------------------- ! CALL edgfill(u,1,nx,1,nx, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'u r3d1' WRITE(nchanl) u CALL edgfill(v,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1) WRITE(nchanl) 'v r3d2' WRITE(nchanl) v CALL edgfill(w,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz) WRITE(nchanl) 'w r3d3' WRITE(nchanl) w ! !----------------------------------------------------------------------- ! ! Write out scalars ! !----------------------------------------------------------------------- ! DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 tem1(i,j,k) = ptbar(i,j,k) + ptprt(i,j,k) END DO END DO END DO CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'pt r3d0' WRITE(nchanl) tem1 DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 tem1(i,j,k) = pbar(i,j,k) + pprt(i,j,k) END DO END DO END DO CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'p r3d0' WRITE(nchanl) tem1 END IF END IF ! varout ! !----------------------------------------------------------------------- ! ! If mstout = 1, write out moisture scalars. ! !----------------------------------------------------------------------- ! IF(mstout == 1) THEN IF( totout == 0 ) THEN ! !----------------------------------------------------------------------- ! ! Write out perturbation to history dump ! !----------------------------------------------------------------------- ! DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 tem1(i,j,k)=qv(i,j,k)-qvbar(i,j,k) END DO END DO END DO CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'qvprt r3d0' WRITE(nchanl) tem1 ELSE ! !----------------------------------------------------------------------- ! ! Write out total values to history dump ! !----------------------------------------------------------------------- ! CALL edgfill(qv,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'qv r3d0' WRITE(nchanl) qv END IF CALL edgfill(qc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'qc r3d0' WRITE(nchanl) qc CALL edgfill(qr,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'qr r3d0' WRITE(nchanl) qr IF(rainout == 1) THEN CALL edgfill(raing, 1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'raing r2d0' WRITE(nchanl) raing CALL edgfill(rainc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'rainc r2d0' WRITE(nchanl) rainc END IF !rainout IF ( prcout == 1 ) THEN CALL edgfill(prcrate,1,nx,1,nx-1, 1,ny,1,ny-1, 1,4,1,4) WRITE(nchanl) 'prcrat1 r2d0' WRITE(nchanl) ((prcrate(i,j,1),i=1,nx),j=1,ny) WRITE(nchanl) 'prcrat2 r2d0' WRITE(nchanl) ((prcrate(i,j,2),i=1,nx),j=1,ny) WRITE(nchanl) 'prcrat3 r2d0' WRITE(nchanl) ((prcrate(i,j,3),i=1,nx),j=1,ny) WRITE(nchanl) 'prcrat4 r2d0' WRITE(nchanl) ((prcrate(i,j,4),i=1,nx),j=1,ny) END IF IF(iceout == 1) THEN CALL edgfill(qi,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'qi r3d0' WRITE(nchanl) qi CALL edgfill(qs,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'qs r3d0' WRITE(nchanl) qs CALL edgfill(qh,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'qh r3d0' WRITE(nchanl) qh END IF !iceout END IF !mstout ! !----------------------------------------------------------------------- ! ! If tkeout = 1, write out tke. ! !----------------------------------------------------------------------- ! IF( tkeout == 1 ) THEN CALL edgfill(tke,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'tke r3d0' WRITE(nchanl) tke END IF ! !----------------------------------------------------------------------- ! ! If trbout = 1, write out the turbulence parameter, km. ! !----------------------------------------------------------------------- ! IF( trbout == 1 ) THEN CALL edgfill(kmh,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'kmh r3d0' WRITE(nchanl) kmh CALL edgfill(kmv,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'kmv r3d0' WRITE(nchanl) kmv END IF ! trbout ! !----------------------------------------------------------------------- ! ! If sfcout = 1, write out the surface variables, ! tsoil, qsoil, and wetcanp. ! !----------------------------------------------------------------------- ! IF( sfcout == 1) THEN IF( nstyp <= 1 ) THEN CALL edgfill(tsoil(1,1,1,0), 1,nx,1,nx-1, 1,ny,1,ny-1, & 1,nzsoil,1,nzsoil) CALL edgfill(qsoil(1,1,1,0), 1,nx,1,nx-1, 1,ny,1,ny-1, & 1,nzsoil,1,nzsoil) IF (intver == intver410) THEN ! ! 06/28/2002 Zuwen He ! ! In version 410, tsoil is the average of tsoil from 2 to nzsoil in later ! version, and wetdp is similar. ! tem1(:,:,1)=0. DO k=2,nzsoil DO j=1,ny DO i=1,nx tem1(i,j,1)=tem1(i,j,1)+tsoil(i,j,k,0) END DO END DO END DO DO j=1,ny DO i=1,nx tem1(i,j,1)=tem1(i,j,1)/float(nzsoil-1) END DO END DO WRITE(nchanl) 'tsfc r2d0' WRITE(nchanl) ((tsoil(i,j,1,0),i=1,nx),j=1,ny) WRITE(nchanl) 'tsoil r2d0' WRITE(nchanl) ((tem1(i,j,1),i=1,nx),j=1,ny) tem1(:,:,1)=0. DO k=2,nzsoil DO j=1,ny DO i=1,nx tem1(i,j,1)=tem1(i,j,1)+qsoil(i,j,k,0) END DO END DO END DO DO j=1,ny DO i=1,nx tem1(i,j,1)=tem1(i,j,1)/float(nzsoil-1) END DO END DO WRITE(nchanl) 'wetsfc r2d0' WRITE(nchanl) ((qsoil(i,j,1,0),i=1,nx),j=1,ny) WRITE(nchanl) 'wetdp r2d0' WRITE(nchanl) ((tem1(i,j,1),i=1,nx),j=1,ny) ELSE IF (intver == intver500) THEN WRITE(nchanl) 'tsoil rs3d0' WRITE(nchanl) (((tsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) WRITE(nchanl) 'qsoil rs3d0' WRITE(nchanl) (((qsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) END IF ! intver CALL edgfill(wetcanp(1,1,0),1,nx,1,nx-1, 1,ny,1,ny-1, & 1,1,1,1) WRITE(nchanl) 'wetcanp r2d0' WRITE(nchanl) ((wetcanp(i,j,0),i=1,nx),j=1,ny) ELSE DO is=0,nstyp CALL edgfill(tsoil(1,1,1,is), 1,nx,1,nx-1, 1,ny,1,ny-1, & 1,nzsoil,1,nzsoil) CALL edgfill(qsoil(1,1,1,is), 1,nx,1,nx-1, 1,ny,1,ny-1, & 1,nzsoil,1,nzsoil) IF (intver == intver410) THEN ! ! 06/28/2002 Zuwen He ! ! In version 410, tsoil is the average of tsoil from 2 to nzsoil in later ! version, and wetdp is similar. ! tem1(:,:,1)=0. DO k=2,nzsoil DO j=1,ny DO i=1,nx tem1(i,j,1)=tem1(i,j,1)+tsoil(i,j,k,is) END DO END DO END DO DO j=1,ny DO i=1,nx tem1(i,j,1)=tem1(i,j,1)/float(nzsoil-1) END DO END DO WRITE(nchanl) 'tsfc r2d0' WRITE(nchanl) ((tsoil(i,j,1,is),i=1,nx),j=1,ny) WRITE(nchanl) 'tsoil r2d0' WRITE(nchanl) ((tem1(i,j,1),i=1,nx),j=1,ny) tem1(:,:,1)=0. DO k=2,nzsoil DO j=1,ny DO i=1,nx tem1(i,j,1)=tem1(i,j,1)+qsoil(i,j,k,is) END DO END DO END DO DO j=1,ny DO i=1,nx tem1(i,j,1)=tem1(i,j,1)/float(nzsoil-1) END DO END DO WRITE(nchanl) 'wetsfc r2d0' WRITE(nchanl) ((qsoil(i,j,1,is),i=1,nx),j=1,ny) WRITE(nchanl) 'wetdp r2d0' WRITE(nchanl) ((tem1(i,j,1),i=1,nx),j=1,ny) ELSE IF (intver == intver500) THEN WRITE(nchanl) 'tsoil rs3d0' WRITE(nchanl) (((tsoil(i,j,k,is),i=1,nx),j=1,ny),k=1,nzsoil) WRITE(nchanl) 'qsoil rs3d0' WRITE(nchanl) (((qsoil(i,j,k,is),i=1,nx),j=1,ny),k=1,nzsoil) END IF ! intver CALL edgfill(wetcanp(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1, & 1,1,1,1) WRITE(nchanl) 'wetcanp r2d0' WRITE(nchanl) ((wetcanp(i,j,is),i=1,nx),j=1,ny) END DO END IF IF (snowout == 1) THEN CALL edgfill(snowdpth,1,nx,1,nx-1, 1,ny,1,ny-1, & 1,1,1,1) WRITE(nchanl) 'snowdpthr2d0' WRITE(nchanl) ((snowdpth(i,j),i=1,nx),j=1,ny) END IF END IF ! !----------------------------------------------------------------------- ! ! If radout = 1, write out the radiation arrays ! !----------------------------------------------------------------------- ! IF( radout == 1 ) THEN CALL edgfill(radfrc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1) WRITE(nchanl) 'radfrc r3d0' WRITE(nchanl) radfrc CALL edgfill(radsw,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'radsw r2d0' WRITE(nchanl) radsw CALL edgfill(rnflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'rnflx r2d0' WRITE(nchanl) rnflx IF (intver >= intver500) THEN CALL edgfill(radswnet,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'radswnetr2d0' WRITE(nchanl) radswnet CALL edgfill(radlwin,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'radlwin r2d0' WRITE(nchanl) radlwin END IF ! intver END IF ! radout ! !----------------------------------------------------------------------- ! ! If flxout = 1, write out the surface fluxes ! !----------------------------------------------------------------------- ! IF( flxout == 1 ) THEN CALL edgfill(usflx,1,nx,1,nx, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'usflx r2d0' WRITE(nchanl) usflx CALL edgfill(vsflx,1,nx,1,nx-1, 1,ny,1,ny, 1,1,1,1) WRITE(nchanl) 'vsflx r2d0' WRITE(nchanl) vsflx CALL edgfill(ptsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'ptsflx r2d0' WRITE(nchanl) ptsflx CALL edgfill(qvsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1) WRITE(nchanl) 'qvsflx r2d0' WRITE(nchanl) qvsflx END IF ! flxout RETURN END SUBROUTINE bindump ! ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE BINJOINDUMP ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE binjoindump(nx,ny,nz,nzsoil,nstyps, nchanl, grdbas, & 1,114 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv, & ubar,vbar,ptbar,pbar,rhobar,qvbar, & x,y,z,zp,zpsoil, & soiltyp,stypfrct,vegtyp,lai,roufns,veg, & tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & tem1) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Write history data into channel nchanl as binary data. ! ! All data read in are located at the original staggered grid points. ! ! Note: coordinate fields are dumped as 3 dimensional fields which ! have been converted from meters to kilometers. This is for the ! convenience of the plotting applications. ! ! The last 4 characters of the 12 character label written out ! with each 1-,2-, or 3-d array is used by the splitdump and ! joinfiles subroutines used by the message passing version of the ! ARPS (and also by some auxiliary ARPS I/O routines) ! to determine the data type of the array. ! Key to the labels: ! ! 'nnnnnnn tdds' ! ! n - characters containing the name of the variable. ! t - type of variable: "r" for real and "i" for integer. ! dd - number of dimensions: "1d" "2d" or "3d". ! s - staggered dimension: "0" for centered, ! "1" for staggered in x, ! "2" for staggered in y, ! "3" for staggered in z. ! !----------------------------------------------------------------------- ! ! ! !----------------------------------------------------------------------- ! ! AUTHOR: Yunheng Wang ! 8/16/02. ! Based on subroutine bindump. ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT: ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! nzsoil Number of grid points in the soil ! ! nchanl FORTRAN I/O channel number for history data output. ! grdbas Flag indicating if this is a call for the data dump ! of grid and base state arrays only. If so, grdbas=1. ! ! u x component of velocity at a given time level (m/s) ! v y component of velocity at a given time level (m/s) ! w Vertical component of Cartesian velocity at a given ! time level (m/s) ! ptprt Perturbation potential temperature at a given time ! level (K) ! pprt Perturbation pressure at a given time level (Pascal) ! qv Water vapor specific humidity at a given time level (kg/kg) ! qc Cloud water mixing ratio at a given time level (kg/kg) ! qr Rainwater mixing ratio at a given time level (kg/kg) ! qi Cloud ice mixing ratio at a given time level (kg/kg) ! qs Snow mixing ratio at a given time level (kg/kg) ! qh Hail mixing ratio at a given time level (kg/kg) ! tke Turbulent Kinetic Energy ((m/s)**2) ! ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! ! ubar Base state zonal velocity component (m/s) ! vbar Base state meridional velocity component (m/s) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! rhobar Base state density (kg/m**3) ! qvbar Base state water vapor specific humidity (kg/kg) ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp Vertical coordinate of grid points in physical space (m) ! zpsoil Vertical coordinate of grid points in the soil (m) ! ! soiltyp Soil type ! stypfrct Soil type fraction ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil Soil temperature (K) ! qsoil Soil moisture (m**3/m**3) ! wetcanp Canopy water amount ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! prcrate Precipitation rates ! ! radfrc Radiation forcing (K/s) ! radsw Solar radiation reaching the surface ! rnflx Net radiation flux absorbed by surface ! radswnet Net shortwave radiation ! radlwin Incoming longwave radiation ! ! usflx Surface flux of u-momentum (kg/(m*s**2)) ! vsflx Surface flux of v-momentum (kg/(m*s**2)) ! ptsflx Surface heat flux (K*kg/(m**2*s)) ! qvsflx Surface moisture flux (kg/(m**2*s)) ! ! OUTPUT: ! ! None. ! ! WORK ARRAY: ! ! tem1 Temporary work array. ! ! !----------------------------------------------------------------------- ! ! The following parameters are passed into this subroutine through ! a common block in globcst.inc, and they determine which ! variables are output. ! ! grdout =0 or 1. If grdout=0, grid variables are not dumped. ! basout =0 or 1. If basout=0, base state variables are not dumped. ! varout =0 or 1. If varout=0, model perturbation variables are not dumped. ! mstout =0 or 1. If mstout=0, water variables are not dumped. ! rainout=0 or 1. If rainout=0, rain variables are not dumped. ! prcout =0 or 1. If prcout=0, precipitation rates are not dumped. ! iceout =0 or 1. If iceout=0, qi, qs and qh are not dumped. ! tkeout =0 or 1. If tkeout=0, tke is not dumped. ! trbout =0 or 1. If trbout=0, turbulence parameter km is not dumped. ! sfcout =0 or 1. If sfcout=0, surface variables are not dumped. ! landout=0 or 1. If landout=0, surface propertty arrays are not dumped. ! radout =0 or 1. If radout =0, radiation arrays are not dumped. ! flxout =0 or 1. If flxout =0, surface flux arrays are not dumped. ! ! These following parameters are also passed in through common ! blocks in globcst.inc. ! ! runname,curtim,umove,vmove,xgrdorg,ygrdorg ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Number of grid points in 3 directions INTEGER :: nzsoil ! Number of grid points in the soil INTEGER :: nchanl ! FORTRAN I/O channel number for output INTEGER :: grdbas ! If this is a grid/base state array dump REAL :: u (nx,ny,nz) ! Total u-velocity (m/s) REAL :: v (nx,ny,nz) ! Total v-velocity (m/s) REAL :: w (nx,ny,nz) ! Total w-velocity (m/s) REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K) REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal) REAL :: qv (nx,ny,nz) ! Water vapor specific humidity (kg/kg) REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg) REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s) REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s) REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K) REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal) REAL :: rhobar(nx,ny,nz) ! Base state air density (kg/m**3) REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific humidity ! (kg/kg) REAL :: x (nx) ! The x-coord. of the physical and ! computational grid. Defined at u-point. REAL :: y (ny) ! The y-coord. of the physical and ! computational grid. Defined at v-point. REAL :: z (nz) ! The z-coord. of the computational grid. ! Defined at w-point on the staggered grid. REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at ! w-point of the staggered grid. REAL :: zpsoil(nx,ny,nzsoil) ! The physical height coordinate defined at ! w-point of the soil INTEGER :: nstyps ! Number of soil types INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type REAL :: stypfrct(nx,ny,nstyps) ! Soil type fractions INTEGER :: vegtyp (nx,ny) ! Vegetation type REAL :: lai (nx,ny) ! Leaf Area Index REAL :: roufns (nx,ny) ! Surface roughness REAL :: veg (nx,ny) ! Vegetation fraction REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K) REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount REAL :: snowdpth(nx,ny) ! Snow depth (m) REAL :: raing(nx,ny) ! Grid supersaturation rain REAL :: rainc(nx,ny) ! Cumulus convective rain REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s)) ! prcrate(1,1,1) = total precip. rate ! prcrate(1,1,2) = grid scale precip. rate ! prcrate(1,1,3) = cumulus precip. rate ! prcrate(1,1,4) = microphysics precip. rate REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s) REAL :: radsw (nx,ny) ! Solar radiation reaching the surface REAL :: rnflx (nx,ny) ! Net radiation flux absorbed by surface REAL :: radswnet(nx,ny) ! Net shortwave radiation REAL :: radlwin(nx,ny) ! Incoming longwave radiation REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2)) REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2)) REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m**2*s)) REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s)) REAL :: tem1 (nx,ny,nz) ! Temporary work array ! !----------------------------------------------------------------------- ! ! Parameters describing this routine ! !----------------------------------------------------------------------- ! ! 06/28/2002 Zuwen He ! ! fmtver??: to label each data a version. ! intver??: an integer to allow faster comparison than fmtver??, ! which are strings. ! ! Verion 5.00: significant change in soil variables since version 4.10. ! CHARACTER (LEN=40) :: fmtver,fmtver410,fmtver500 INTEGER :: intver,intver410,intver500 PARAMETER (fmtver410='004.10 Binary Data',intver410=410) PARAMETER (fmtver500='005.00 Binary Data',intver500=500) CHARACTER (LEN=10) :: tmunit PARAMETER (tmunit='seconds ') ! !----------------------------------------------------------------------- ! ! Misc. local variables: ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k,l,is INTEGER :: idummy REAL :: rdummy INTEGER :: nxlg, nylg, nzlg, nzsoillg INTEGER :: n3d, istat REAL, ALLOCATABLE :: out1d(:) REAL, ALLOCATABLE :: out3d(:,:,:) INTEGER, ALLOCATABLE :: out3di(:,:,:) REAL, ALLOCATABLE :: outtsoil(:,:,:,:), outqsoil(:,:,:,:) ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'grid.inc' ! Grid & map parameters. INCLUDE 'mp.inc' ! mpi parameters. ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! nxlg = nproc_x*(nx-3)+3 nylg = nproc_y*(ny-3)+3 nzlg = nz nzsoillg = nzsoil n3d = MAX(nzlg, nzsoillg, nstyps+1, 4) ! 3rd dimenson for out3d intver = intver500 ! for the time being, in the future, we will ! allow to dump data in the different version ! intver will be assigned from input file IF (intver == intver410) THEN fmtver=fmtver410 ELSE IF (intver == intver500) THEN fmtver=fmtver500 ELSE IF (myproc == 0) WRITE(6,'(/1x,a,i10,a/)') & 'Data format, intver=',intver,', not found. The Job stopped.' CALL arpsstop('arpstop called from binjoindumpp.',1) END IF IF (myproc == 0) THEN WRITE(6,'(/1x,a,a,a/)') & 'Data dump format, fmtver=',fmtver,'. ' WRITE(6,'(1x,a,f13.3/)') 'Writing history data at time=', curtim ! !----------------------------------------------------------------------- ! ! Write header info ! !----------------------------------------------------------------------- ! WRITE(nchanl) fmtver WRITE(nchanl) runname WRITE(nchanl) nocmnt IF( nocmnt > 0 ) THEN DO l=1,nocmnt WRITE(nchanl) cmnt(l) END DO END IF WRITE(nchanl) curtim,tmunit IF (intver == intver410) THEN WRITE(nchanl) nxlg,nylg,nzlg ELSE IF (intver == intver500) THEN WRITE(nchanl) nxlg,nylg,nzlg,nzsoillg END IF ! !----------------------------------------------------------------------- ! ! Write the flags for different data groups. ! !----------------------------------------------------------------------- ! idummy = 0 IF( grdbas == 1 ) THEN WRITE(nchanl) 1, 1, 0, mstout, 0, & 0, 0, 0, landout,totout, & 0, idummy, idummy, mapproj, month, & day, year, hour, minute, second ELSE WRITE(nchanl) grdout, basout, varout, mstout, iceout, & trbout, sfcout, rainout,landout,totout, & tkeout, idummy, idummy, mapproj, month, & day, year, hour, minute, second END IF rdummy = 0.0 WRITE(nchanl) umove, vmove, xgrdorg, ygrdorg, trulat1, & trulat2, trulon, sclfct, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy, & tstop, thisdmp, latitud, ctrlat, ctrlon ! !----------------------------------------------------------------------- ! ! If totout=1, write additional parameters to history dump files. ! This is for ARPS version 4.1.2 or later. ! !----------------------------------------------------------------------- ! IF ( totout == 1 ) THEN WRITE(nchanl) nstyp, prcout, radout, flxout, 0, & ! 0 for snowcvr snowout,idummy, idummy, idummy, idummy, & idummy, idummy, idummy, idummy, idummy, & idummy, idummy, idummy, idummy, idummy WRITE(nchanl) rdummy, rdummy, rdummy, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy END IF END IF ! myproc == 0 ALLOCATE (out1d( MAX(nxlg,nylg,nzlg) ),stat=istat) IF (istat /= 0) & CALL arpsstop("BINJOINDUMP: ERROR allocating out1d, returning",1) ALLOCATE (out3d( nxlg,nylg, n3d ),stat=istat) IF (istat /= 0) & CALL arpsstop("BINJOINDUMP: ERROR allocating out3d, returning",1) ALLOCATE (out3di( nxlg,nylg, nstyps ),stat=istat) IF (istat /= 0) & CALL arpsstop("BINJOINDUMP: ERROR allocating out3di, returning",1) ALLOCATE (outtsoil( nxlg,nylg, nzsoillg, 0:nstyps ),stat=istat) IF (istat /= 0) & CALL arpsstop("BINJOINDUMP: ERROR allocating outtsoil, returning",1) ALLOCATE (outqsoil( nxlg,nylg, nzsoillg, 0:nstyps ),stat=istat) IF (istat /= 0) & CALL arpsstop("BINJOINDUMP: ERROR allocating outqsoil, returning",1) !----------------------------------------------------------------------- ! ! If grdout=1 or grdbas=1, write out grid variables ! !----------------------------------------------------------------------- ! IF(grdout == 1 .OR. grdbas == 1 ) THEN CALL mpimerge1dx(x,nx,out1d) IF (myproc == 0) THEN WRITE(nchanl) 'x coord r1d1' WRITE(nchanl) out1d(1:nxlg) END IF CALL mpimerge1dy(y,ny,out1d) IF (myproc == 0) THEN WRITE(nchanl) 'y coord r1d2' WRITE(nchanl) out1d(1:nylg) END IF IF (myproc == 0) THEN WRITE(nchanl) 'z coord r1d3' WRITE(nchanl) z END IF CALL mpimerge3d(zp,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg) WRITE(nchanl) 'zp coor r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF IF (intver == intver500) THEN CALL mpimerge3d(zpsoil,nx,ny,nzsoil,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzsoillg,1,nzsoillg) WRITE(nchanl) 'zpsoil rs3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzsoillg) END IF END IF END IF ! grdout ! !----------------------------------------------------------------------- ! ! If basout=1, write out base state variables. ! !----------------------------------------------------------------------- ! IF(basout == 1 .OR. grdbas == 1 ) THEN CALL mpimerge3d(ubar,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'ubar r3d1' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(vbar,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'vbar r3d2' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF IF (myproc == 0) THEN out3d(:,:,:) = 0.0 WRITE(nchanl) 'wbar r3d3' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(ptbar,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'ptbar r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(pbar,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'pbar r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF IF(mstout == 1) THEN CALL mpimerge3d(qvbar,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'qvbar r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF END IF IF(landout == 1) THEN IF( nstyp <= 1 ) THEN CALL mpimerge3di(soiltyp(:,:,1),nx,ny,1,out3di) IF (myproc == 0) THEN CALL iedgfill(out3di(1,1,1),1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,1,1,1) WRITE(nchanl) 'soiltyp i2d0' WRITE(nchanl) out3di(1:nxlg,1:nylg,1) END IF ELSE CALL mpimerge3di(soiltyp,nx,ny,nstyp,out3di) CALL mpimerge3d(stypfrct,nx,ny,nstyp,out3d) IF (myproc == 0) THEN DO is=1,nstyp CALL iedgfill(out3di(1,1,is),1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,1,1,1) WRITE(nchanl) 'soiltyp i2d0' WRITE(nchanl) ((out3di(i,j,is),i=1,nxlg),j=1,nylg) CALL edgfill(out3d(1,1,is),1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,1,1,1) WRITE(nchanl) 'stypfrc r2d0' WRITE(nchanl) ((out3d(i,j,is),i=1,nxlg),j=1,nylg) END DO END IF END IF CALL mpimerge3di(vegtyp,nx,ny,1,out3di) IF (myproc == 0) THEN CALL iedgfill(out3di ,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'vegtyp i2d0' WRITE(nchanl) out3di(1:nxlg,1:nylg,1) END IF CALL mpimerge3d(lai,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'lai r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF CALL mpimerge3d(roufns,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d ,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'roufns r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF CALL mpimerge3d(veg,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d ,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'veg r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF END IF END IF IF ( grdbas == 1 ) RETURN ! !----------------------------------------------------------------------- ! ! If varout = 1, Write out uprt, vprt, wprt, ptprt, pprt. ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Write out u,v and w. ! !----------------------------------------------------------------------- ! IF(varout == 1) THEN IF ( totout == 0 ) THEN ! !----------------------------------------------------------------------- ! ! Write out perturbations to history dump ! !----------------------------------------------------------------------- ! DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx tem1(i,j,k)=u(i,j,k)-ubar(i,j,k) END DO END DO END DO CALL mpimerge3d(tem1,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'uprt r3d1' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF DO k=1,nz-1 DO i=1,nx-1 DO j=1,ny tem1(i,j,k)=v(i,j,k)-vbar(i,j,k) END DO END DO END DO CALL mpimerge3d(tem1,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'vprt r3d2' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(w,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg) WRITE(nchanl) 'wprt r3d3' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF ! !----------------------------------------------------------------------- ! ! Write out scalars ! !----------------------------------------------------------------------- ! CALL mpimerge3d(ptprt,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'ptprt r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(pprt,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'pprt r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF ELSE ! !----------------------------------------------------------------------- ! ! Write out total values to history dump ! !----------------------------------------------------------------------- ! CALL mpimerge3d(u,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'u r3d1' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(v,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'v r3d2' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(w,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg) WRITE(nchanl) 'w r3d3' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF ! !----------------------------------------------------------------------- ! ! Write out scalars ! !----------------------------------------------------------------------- ! DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 tem1(i,j,k) = ptbar(i,j,k) + ptprt(i,j,k) END DO END DO END DO CALL mpimerge3d(tem1,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'pt r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 tem1(i,j,k) = pbar(i,j,k) + pprt(i,j,k) END DO END DO END DO CALL mpimerge3d(tem1,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'p r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF END IF END IF ! varout ! !----------------------------------------------------------------------- ! ! If mstout = 1, write out moisture scalars. ! !----------------------------------------------------------------------- ! IF(mstout == 1) THEN IF( totout == 0 ) THEN ! !----------------------------------------------------------------------- ! ! Write out perturbation to history dump ! !----------------------------------------------------------------------- ! DO k=1,nz-1 DO j=1,ny-1 DO i=1,nx-1 tem1(i,j,k)=qv(i,j,k)-qvbar(i,j,k) END DO END DO END DO CALL mpimerge3d(tem1,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'qvprt r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF ELSE ! !----------------------------------------------------------------------- ! ! Write out total values to history dump ! !----------------------------------------------------------------------- ! CALL mpimerge3d(qv,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'qv r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF END IF CALL mpimerge3d(qc,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'qc r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(qr,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'qr r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF IF(rainout == 1) THEN CALL mpimerge3d(raing,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d, 1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'raing r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF CALL mpimerge3d(rainc,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'rainc r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF END IF !rainout IF ( prcout == 1 ) THEN CALL mpimerge3d(prcrate,nx,ny,4,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,4,1,4) WRITE(nchanl) 'prcrat1 r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) WRITE(nchanl) 'prcrat2 r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,2) WRITE(nchanl) 'prcrat3 r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,3) WRITE(nchanl) 'prcrat4 r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,4) END IF END IF IF(iceout == 1) THEN CALL mpimerge3d(qi,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'qi r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(qs,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'qs r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(qh,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'qh r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF END IF !iceout END IF !mstout ! !----------------------------------------------------------------------- ! ! If tkeout = 1, write out tke. ! !----------------------------------------------------------------------- ! IF( tkeout == 1 ) THEN CALL mpimerge3d(tke,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'tke r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF END IF ! !----------------------------------------------------------------------- ! ! If trbout = 1, write out the turbulence parameter, km. ! !----------------------------------------------------------------------- ! IF( trbout == 1 ) THEN CALL mpimerge3d(kmh,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'kmh r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(kmv,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'kmv r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF END IF ! trbout ! !----------------------------------------------------------------------- ! ! If sfcout = 1, write out the surface variables, ! tsoil, qsoil, and wetcanp. ! !----------------------------------------------------------------------- ! IF( sfcout == 1) THEN IF( nstyp <= 1 ) THEN CALL mpimerge3d(tsoil(:,:,:,0),nx,ny,nzsoil,out3d) IF (myproc == 0) THEN CALL edgfill(out3d, 1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,nzsoillg,1,nzsoillg) WRITE(nchanl) 'tsoil rs3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzsoillg) END IF CALL mpimerge3d(qsoil(:,:,:,0),nx,ny,nzsoil,out3d) IF (myproc == 0) THEN CALL edgfill(out3d, 1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,nzsoillg,1,nzsoillg) WRITE(nchanl) 'qsoil rs3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzsoillg) END IF CALL mpimerge3d(wetcanp(:,:,0),nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,1,1,1) WRITE(nchanl) 'wetcanp r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF ELSE CALL mpimerge4d(tsoil,nx,ny,nzsoil,nstyp+1,outtsoil) CALL mpimerge4d(qsoil,nx,ny,nzsoil,nstyp+1,outqsoil) CALL mpimerge3d(wetcanp,nx,ny,nstyp+1,out3d) IF (myproc == 0) THEN DO is=0,nstyp CALL edgfill(outtsoil(1,1,1,is), 1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,nzsoillg,1,nzsoillg) CALL edgfill(outqsoil(1,1,1,is), 1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,nzsoillg,1,nzsoillg) WRITE(nchanl) 'tsoil rs3d0' WRITE(nchanl) outtsoil(1:nxlg,1:nylg,1:nzsoillg,is) WRITE(nchanl) 'qsoil rs3d0' WRITE(nchanl) outqsoil(1:nxlg,1:nylg,1:nzsoillg,is) CALL edgfill(out3d(1,1,is+1),1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,1,1,1) WRITE(nchanl) 'wetcanp r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,is+1) END DO END IF END IF IF (snowout == 1) THEN CALL mpimerge3d(snowdpth,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, & 1,1,1,1) WRITE(nchanl) 'snowdpthr2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF END IF END IF ! !----------------------------------------------------------------------- ! ! If radout = 1, write out the radiation arrays ! !----------------------------------------------------------------------- ! IF( radout == 1 ) THEN CALL mpimerge3d(radfrc,nx,ny,nz,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1) WRITE(nchanl) 'radfrc r3d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg) END IF CALL mpimerge3d(radsw,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'radsw r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF CALL mpimerge3d(rnflx,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'rnflx r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF IF (intver >= intver500) THEN CALL mpimerge3d(radswnet,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'radswnetr2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF CALL mpimerge3d(radlwin,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'radlwin r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF END IF ! intver END IF ! radout ! !----------------------------------------------------------------------- ! ! If flxout = 1, write out the surface fluxes ! !----------------------------------------------------------------------- ! IF( flxout == 1 ) THEN CALL mpimerge3d(usflx,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'usflx r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF CALL mpimerge3d(vsflx,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg, 1,1,1,1) WRITE(nchanl) 'vsflx r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF CALL mpimerge3d(ptsflx,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'ptsflx r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF CALL mpimerge3d(qvsflx,nx,ny,1,out3d) IF (myproc == 0) THEN CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1) WRITE(nchanl) 'qvsflx r2d0' WRITE(nchanl) out3d(1:nxlg,1:nylg,1) END IF END IF ! flxout DEALLOCATE(out1d) DEALLOCATE(out3d, out3di) DEALLOCATE(outtsoil, outqsoil) RETURN END SUBROUTINE binjoindump ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE BN2DUMP ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE bn2dump(nx,ny,nz,nzsoil,nstyps, nchanl, grdbas, & 1 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv, & ubar,vbar,ptbar,pbar,rhobar,qvbar, & x,y,z,zp,zpsoil, & soiltyp,stypfrct,vegtyp,lai,roufns,veg, & tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & tem1) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Write history data into channel nchanl as binary data. ! ! This routine can dump the data arrays in a model subdomain and ! at selected data points. ! ! All data read in are located at the original staggered grid points. ! ! Note: coordinate fields are dumped as 3 dimensional fields which ! have been converted from meters to kilometers. This is for the ! convenience of the plotting applications. ! !----------------------------------------------------------------------- ! ! AUTHOR: Ming Xue ! 3/10/92. ! ! MODIFICATION HISTORY: ! ! 4/4/93 (M. Xue) ! Modified, so that data on the original staggered grid are written ! out. Averaging to the volume center is no longer done. ! ! 9/1/94 (Y. Lu) ! Cleaned up documentation. ! ! 02/06/95 (Y. Liu) ! Added map projection parameters into the second binary dumping ! ! 12/09/1998 (Donghai Wang) ! Added the snow cover. ! ! 05/15/2002 (J. Brotzge) ! Added to allow for multiple soil schemes. ! !----------------------------------------------------------------------- ! ! INPUT: ! ! nx Number of grid points in the x-direction (east/west) ! ny Number of grid points in the y-direction (north/south) ! nz Number of grid points in the vertical ! nzsoil Number of grid points in the soil ! ! nchanl FORTRAN I/O channel number for history data output. ! grdbas Flag indicating if this is a call for the data dump ! of grid and base state arrays only. If so, grdbas=1. ! ! u x component of velocity at a given time level (m/s) ! v y component of velocity at a given time level (m/s) ! w Vertical component of Cartesian velocity at a given ! time level (m/s) ! ptprt Perturbation potential temperature at a given time ! level (K) ! pprt Perturbation pressure at a given time level (Pascal) ! qv Water vapor specific humidity at a given time level (kg/kg) ! qc Cloud water mixing ratio at a given time level (kg/kg) ! qr Rainwater mixing ratio at a given time level (kg/kg) ! qi Cloud ice mixing ratio at a given time level (kg/kg) ! qs Snow mixing ratio at a given time level (kg/kg) ! qh Hail mixing ratio at a given time level (kg/kg) ! tke Turbulent Kinetic Energy ((m/s)**2) ! ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! ! ubar Base state zonal velocity component (m/s) ! vbar Base state meridional velocity component (m/s) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! rhobar Base state density (kg/m**3) ! qvbar Base state water vapor specific humidity (kg/kg) ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp Vertical coordinate of grid points in physical space (m) ! zpsoil Vertical coordinate of grid points in soil (m) ! ! soiltyp Soil type ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil Soil temperature (K) ! qsoil Soil moisture (m**3/m**3) ! wetcanp Canopy water amount ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! prcrate Precipitation rates ! ! radfrc Radiation forcing (K/s) ! radsw Solar radiation reaching the surface ! rnflx Net radiation flux absorbed by surface ! radswnet Net shortwave radiation ! radlwin Incoming longwave radiation ! ! usflx Surface flux of u-momentum (kg/(m*s**2)) ! vsflx Surface flux of v-momentum (kg/(m*s**2)) ! ptsflx Surface heat flux (K*kg/(m**2 * s )) ! qvsflx Surface moisture flux of (kg/(m**2 * s)) ! ! OUTPUT: ! ! None. ! ! WORK ARRAY: ! ! tem1 Temporary work array. ! ! !----------------------------------------------------------------------- ! ! The following parameters are passed into this subroutine through ! a common block in globcst.inc. These parameters determine which ! variables are output. ! ! grdout =0 or 1. If grdout=0, grid variables are not dumped. ! basout =0 or 1. If basout=0, base state variables are not dumped. ! varout =0 or 1. If varout=0, model perturbation variables are not dumped. ! mstout =0 or 1. If mstout=0, water variables are not dumped. ! rainout=0 or 1. If rainout=0, rain variables are not dumped. ! prcout =0 or 1. If prcout=0, precipitation rates are not dumped. ! iceout =0 or 1. If iceout=0, qi, qs and qh are not dumped. ! trbout =0 or 1. If trbout=0, turbulence parameter km is not dumped. ! tkeout =0 or 1. If tkeout=0, tke is not dumped. ! radout =0 or 1. If radout=0, radiation arrays are not dumped. ! flxout =0 or 1. If flxout=0, surface fluxes are not dumped. ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Number of grid points in 3 directions INTEGER :: nzsoil ! Number of grid points in the soil INTEGER :: nchanl ! FORTRAN I/O channel number for output INTEGER :: grdbas ! If this is a grid/base state array dump REAL :: u (nx,ny,nz) ! Total u-velocity (m/s) REAL :: v (nx,ny,nz) ! Total v-velocity (m/s) REAL :: w (nx,ny,nz) ! Total w-velocity (m/s) REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K) REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal) REAL :: qv (nx,ny,nz) ! Water vapor specific humidity (kg/kg) REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg) REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s) REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s) REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K) REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal) REAL :: rhobar(nx,ny,nz) ! Base state air density (kg/m**3) REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific humidity ! (kg/kg) REAL :: x (nx) ! The x-coord. of the physical and ! computational grid. Defined at u-point. REAL :: y (ny) ! The y-coord. of the physical and ! computational grid. Defined at v-point. REAL :: z (nz) ! The z-coord. of the computational grid. ! Defined at w-point on the staggered grid. REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at ! w-point of the staggered grid. REAL :: zpsoil(nx,ny,nzsoil) ! The physical height coordinate defined at ! w-point of the soil INTEGER :: nstyps ! Number of soil types INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type REAL :: stypfrct(nx,ny,nstyps) ! Soil type INTEGER :: vegtyp(nx,ny) ! Vegetation type REAL :: lai (nx,ny) ! Leaf Area Index REAL :: roufns (nx,ny) ! Surface roughness REAL :: veg (nx,ny) ! Vegetation fraction REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K) REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount REAL :: snowdpth(nx,ny) ! Snow depth (m) REAL :: raing(nx,ny) ! Grid supersaturation rain REAL :: rainc(nx,ny) ! Cumulus convective rain REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s)) ! prcrate(1,1,1) = total precip. rate ! prcrate(1,1,2) = grid scale precip. rate ! prcrate(1,1,3) = cumulus precip. rate ! prcrate(1,1,4) = microphysics precip. rate REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s) REAL :: radsw (nx,ny) ! Solar radiation reaching the surface REAL :: rnflx (nx,ny) ! Net radiation flux absorbed by surface REAL :: radswnet(nx,ny) ! Net shortwave radiation REAL :: radlwin(nx,ny) ! Incoming longwave radiation REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2)) REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2)) REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m**2*s)) REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s)) REAL :: tem1 (nx,ny,nz) ! Temporary work array ! !----------------------------------------------------------------------- ! ! Parameters describing this routine ! !----------------------------------------------------------------------- ! CHARACTER (LEN=40) :: fmtver PARAMETER (fmtver='004.10 2nd Binary Data') CHARACTER (LEN=10) :: tmunit PARAMETER (tmunit='seconds ') ! !----------------------------------------------------------------------- ! ! Misc. local variables: ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k,l,is, idummy REAL :: rdummy INTEGER :: nxout,nyout,nzout ! The size of array to be written out. INTEGER :: nzsoilout ! The size of array to be written out. INTEGER :: ist ,ind ,isk ,jst ,jnd ,jsk ,kst ,knd ,ksk INTEGER :: ist1,ind1,isk1,jst1,jnd1,jsk1,kst1,knd1,ksk1 INTEGER :: lst, lnd, lsk INTEGER :: setdomn,setskip SAVE setdomn, setskip SAVE ist,ind,isk,jst,jnd,jsk,kst,knd,ksk,lst,lnd,lsk DATA setdomn, setskip /0,0/ ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'grid.inc' ! Grid & map parameters. ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ IF( setdomn == 0) THEN ! If these parameters are nevers set ... ist = 1 jst = 1 kst = 1 lst = 1 ind = nx jnd = ny knd = nz lnd = nzsoil END IF IF( setskip == 0) THEN ! If these parameters are nevers set ... isk = 1 jsk = 1 ksk = 1 lsk = 1 END IF WRITE(6,'(1x,a,f13.3/)') 'Writing history data at time=', curtim ! !----------------------------------------------------------------------- ! ! Write header info ! !----------------------------------------------------------------------- ! WRITE(nchanl) fmtver WRITE(nchanl) runname WRITE(nchanl) nocmnt IF( nocmnt > 0 ) THEN DO l=1,nocmnt WRITE(nchanl) cmnt(l) END DO END IF WRITE(nchanl) curtim,tmunit nxout = ist+(ind-ist)/isk nyout = jst+(jnd-jst)/jsk nzout = kst+(knd-kst)/ksk nzsoilout = lst+(lnd-lst)/lsk WRITE(nchanl) nxout, nyout, nzout, nzsoilout PRINT*,'nxout= ',nxout,'nyout= ',nyout,'nzout= ',nzout, & 'nzsoilout= ',nzsoilout ! !----------------------------------------------------------------------- ! ! Write the flags for different data groups. ! !----------------------------------------------------------------------- ! idummy = 0 IF( grdbas == 1 ) THEN WRITE(nchanl) 1, 1, 0, mstout, 0, & 0, idummy, idummy, landout, totout, & idummy, idummy, idummy, mapproj, month, & day, year, hour, minute, second ELSE WRITE(nchanl) grdout, basout, varout, mstout, iceout, & trbout, sfcout, rainout,landout,totout, & tkeout, idummy, idummy, mapproj, month, & day, year, hour, minute, second END IF rdummy = 0.0 WRITE(nchanl) umove, vmove, xgrdorg, ygrdorg, trulat1, & trulat2, trulon, sclfct, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy, & tstop, thisdmp, latitud, ctrlat, ctrlon IF ( totout /= 0 ) THEN ! !----------------------------------------------------------------------- ! ! Add new parameters for new version of history data dump ! !----------------------------------------------------------------------- ! WRITE(nchanl) nstyp, prcout, radout, flxout, 0, & ! 0 for snowcvr snowout,idummy, idummy, idummy, idummy, & idummy, idummy, idummy, idummy, idummy, & idummy, idummy, idummy, idummy, idummy WRITE(nchanl) rdummy, rdummy, rdummy, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy, & rdummy, rdummy, rdummy, rdummy, rdummy END IF ! !----------------------------------------------------------------------- ! ! If grdout=1 or grdbas=1, write out grid variables ! !----------------------------------------------------------------------- ! IF(grdout == 1 .OR. grdbas == 1 ) THEN WRITE(nchanl) 'x coordinate' WRITE(nchanl) ( x(i), i=ist,ind,isk) WRITE(nchanl) 'y coordinate' WRITE(nchanl) ( y(j), j=jst,jnd,jsk) WRITE(nchanl) 'z coordinate' WRITE(nchanl) ( z(k), k=kst,knd,ksk) WRITE(nchanl) 'zp coord ' WRITE(nchanl) ((( zp(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'zpsoil coord ' WRITE(nchanl) (((zpsoil(i,j,l),i=ist,ind,isk),j=jst,jnd,jsk), & l=lst,lnd,lsk) END IF ! grdout ! !----------------------------------------------------------------------- ! ! If basout=1, write out base state variables. ! !----------------------------------------------------------------------- ! IF(basout == 1 .OR. grdbas == 1 ) THEN WRITE(nchanl) 'ubar ' WRITE(nchanl) ((( ubar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'vbar ' WRITE(nchanl) ((( vbar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) DO k=kst,knd,ksk DO j=jst,jnd,jsk DO i=ist,ind,isk tem1(i,j,k) = 0.0 END DO END DO END DO WRITE(nchanl) 'wbar ' WRITE(nchanl) (((tem1(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'ptbar ' WRITE(nchanl) (((ptbar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'pbar ' WRITE(nchanl) (((pbar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) IF(mstout == 1) THEN WRITE(nchanl) 'qvbar ' WRITE(nchanl) (((qvbar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) END IF IF(landout == 1) THEN IF( nstyp <= 1 ) THEN WRITE(nchanl) 'soiltyp ' WRITE(nchanl) ((soiltyp(i,j,1),i=ist,ind,isk), & j=jst,jnd,jsk) ELSE DO is=1,nstyp WRITE(nchanl) 'soiltyp ' WRITE(nchanl) ((soiltyp(i,j,is),i=ist,ind,isk), & j=jst,jnd,jsk) WRITE(nchanl) 'stypfrct ' WRITE(nchanl) ((stypfrct(i,j,is),i=ist,ind,isk), & j=jst,jnd,jsk) END DO END IF WRITE(nchanl) 'vegtyp ' WRITE(nchanl) ((vegtyp (i,j),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'lai ' WRITE(nchanl) ((lai (i,j),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'roufns ' WRITE(nchanl) ((roufns (i,j),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'veg ' WRITE(nchanl) ((veg (i,j),i=ist,ind,isk),j=jst,jnd,jsk) END IF END IF IF ( grdbas == 1 ) RETURN ! !----------------------------------------------------------------------- ! ! If varout = 1, Write out uprt, vprt, wprt, ptprt, pprt. ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Write out u, v and w ! !----------------------------------------------------------------------- ! IF(varout == 1) THEN IF ( totout == 0 ) THEN ! !----------------------------------------------------------------------- ! ! Write out perturbatios to history dump ! !----------------------------------------------------------------------- ! WRITE(nchanl) 'uprt ' WRITE(nchanl) ((( u(i,j,k)-ubar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'vprt ' WRITE(nchanl) ((( v(i,j,k)-vbar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'wprt ' WRITE(nchanl) ((( w(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) ! !----------------------------------------------------------------------- ! ! Write out scalars ! !----------------------------------------------------------------------- ! WRITE(nchanl) 'ptprt ' WRITE(nchanl) ((( ptprt(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'pprt ' WRITE(nchanl) ((( pprt(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) ELSE ! !----------------------------------------------------------------------- ! ! Write out total values to history dump ! !----------------------------------------------------------------------- ! WRITE(nchanl) 'u ' WRITE(nchanl) ((( u(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'v ' WRITE(nchanl) ((( v(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'w ' WRITE(nchanl) ((( w(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) ! !----------------------------------------------------------------------- ! ! Write out scalars ! !----------------------------------------------------------------------- ! WRITE(nchanl) 'pt ' WRITE(nchanl) ((( ptprt(i,j,k)+ptbar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'p ' WRITE(nchanl) ((( pprt(i,j,k)+pbar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) END IF END IF ! varout ! !----------------------------------------------------------------------- ! ! If mstout = 1, write out moisture scalars. ! !----------------------------------------------------------------------- ! IF(mstout == 1) THEN IF ( totout == 0 ) THEN ! !----------------------------------------------------------------------- ! ! Write out perturbations to history dump ! !----------------------------------------------------------------------- ! WRITE(nchanl) 'qvprt ' WRITE(nchanl) ((( qv(i,j,k)-qvbar(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) ELSE ! !----------------------------------------------------------------------- ! ! Write out total values to history dump ! !----------------------------------------------------------------------- ! WRITE(nchanl) 'qv ' WRITE(nchanl) ((( qv(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) END IF WRITE(nchanl) 'qc ' WRITE(nchanl) ((( qc(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'qr ' WRITE(nchanl) ((( qr(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) IF(rainout == 1) THEN WRITE(nchanl) 'raing ' WRITE(nchanl) ((raing(i,j),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'rainc ' WRITE(nchanl) ((rainc(i,j),i=ist,ind,isk),j=jst,jnd,jsk) END IF !rainout IF ( prcout == 1 ) THEN WRITE(nchanl) 'prcrate1 ' WRITE(nchanl) ((prcrate(i,j,1),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'prcrate2 ' WRITE(nchanl) ((prcrate(i,j,2),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'prcrate3 ' WRITE(nchanl) ((prcrate(i,j,3),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'prcrate4 ' WRITE(nchanl) ((prcrate(i,j,4),i=ist,ind,isk),j=jst,jnd,jsk) END IF ! prcout IF(iceout == 1) THEN WRITE(nchanl) 'qi ' WRITE(nchanl) ((( qi(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'qs ' WRITE(nchanl) ((( qs(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'qh ' WRITE(nchanl) ((( qh(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) END IF !iceout END IF !mstout ! !----------------------------------------------------------------------- ! ! If tkeout = 1, write out turbulence parameter, km. ! !----------------------------------------------------------------------- ! IF( tkeout == 1 ) THEN WRITE(nchanl) 'tke ' WRITE(nchanl) ((( tke(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) END IF ! tkeout ! !----------------------------------------------------------------------- ! ! If trbout = 1, write out turbulence parameter, km. ! !----------------------------------------------------------------------- ! IF( trbout == 1 ) THEN WRITE(nchanl) 'kmh ' WRITE(nchanl) ((( kmh(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'kmv ' WRITE(nchanl) ((( kmv(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) END IF ! trbout ! !----------------------------------------------------------------------- ! ! If sfcout = 1, write out the surface variables, tsoil, ! qsoil, and wetcanp. ! !----------------------------------------------------------------------- ! IF( sfcout == 1) THEN IF( nstyp <= 1 ) THEN WRITE(nchanl) 'tsoil ' WRITE(nchanl) (((tsoil(i,j,l,0),i=ist,ind,isk),j=jst,jnd,jsk), & l=lst,lnd,lsk) WRITE(nchanl) 'qsoil ' WRITE(nchanl) (((qsoil(i,j,l,0),i=ist,ind,isk),j=jst,jnd,jsk), & l=lst,lnd,lsk) WRITE(nchanl) 'wetcanp ' WRITE(nchanl) ((wetcanp(i,j,0),i=ist,ind,isk),j=jst,jnd,jsk) ELSE DO is=0,nstyp WRITE(nchanl) 'tsoil ' WRITE(nchanl) (((tsoil(i,j,l,is),i=ist,ind,isk), & j=jst,jnd,jsk),l=lst,lnd,lsk) WRITE(nchanl) 'qsoil ' WRITE(nchanl) (((qsoil(i,j,l,is),i=ist,ind,isk), & j=jst,jnd,jsk),l=lst,lnd,lsk) WRITE(nchanl) 'wetcanp ' WRITE(nchanl) ((wetcanp(i,j,is),i=ist,ind,isk), & j=jst,jnd,jsk) END DO END IF IF(snowout == 1) THEN WRITE(nchanl) 'snowdpth ' WRITE(nchanl) ((snowdpth(i,j),i=ist,ind,isk), & j=jst,jnd,jsk) END IF END IF ! sfcout done ! !----------------------------------------------------------------------- ! ! If radout = 1, write out radiation arrays, radfrc, radsw and ! rnflx. ! !----------------------------------------------------------------------- ! IF( radout == 1 ) THEN WRITE(nchanl) 'radfrc ' WRITE(nchanl) ((( radfrc(i,j,k), & i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk) WRITE(nchanl) 'radsw ' WRITE(nchanl) (( radsw(i,j), & i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'rnflx' WRITE(nchanl) (( rnflx(i,j), & i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'radswnet ' WRITE(nchanl) (( radswnet(i,j), & i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'radlwin ' WRITE(nchanl) (( radlwin(i,j), & i=ist,ind,isk),j=jst,jnd,jsk) END IF ! radout ! !----------------------------------------------------------------------- ! ! If flxout = 1, write out surface fluxes ! !----------------------------------------------------------------------- ! IF( flxout == 1 ) THEN WRITE(nchanl) 'usflx ' WRITE(nchanl) (( usflx(i,j),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'vsflx ' WRITE(nchanl) (( vsflx(i,j),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'ptsflx ' WRITE(nchanl) (( ptsflx(i,j),i=ist,ind,isk),j=jst,jnd,jsk) WRITE(nchanl) 'qvsflx ' WRITE(nchanl) (( qvsflx(i,j),i=ist,ind,isk),j=jst,jnd,jsk) END IF ! flxout RETURN ENTRY bdmpdomn(ist1,ind1,jst1,jnd1,kst1,knd1) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! To set the start and end indicies of the model subdomain ! in which the data is dumped out. ! !----------------------------------------------------------------------- ! ist = ist1 jst = jst1 kst = kst1 ind = ind1 jnd = jnd1 knd = knd1 setdomn = 1 RETURN ENTRY bdmpskip(isk1, jsk1, ksk1) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! To set data skip parameters for data dump. ! !----------------------------------------------------------------------- ! isk = isk1 jsk = jsk1 ksk = ksk1 setskip = 1 RETURN END SUBROUTINE bn2dump