!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE DTAREAD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE dtaread(nx,ny,nz,nstyps, & 23,130
hinfmt, nchanl,grdbasfn,lengbf,datafn,lendtf, time, &
x,y,z,zp, 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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1, tem2, tem3)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the reading of history data of various formats.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 2/19/1992.
!
! MODIFICATION HISTORY:
!
! 10/06/1992 (K. Brewster)
! Added ENTRY SETGBRD to allow the reading of
! more than one grid-base file in a program.
!
! 4/23/93 (M. Xue)
! New data format.
!
! 3/30/94 (M.Xue)
! GrADS data format and surface fields.
!
! 12/01/95 (Yuhe Liu)
! Fixed a bug in netread call for dumping grid and base data file.
! The file name was previously given to datafn, instead of grdbasfn.
!
! 12/01/95 (Yuhe Liu)
! Changed the order of reading grid and base file and history files.
! Previously the code read history file first and then checked if
! need to read the grid and base file. In order to cooperate with
! GRIB data format which no longer contains perturbation variables,
! the base fields have to ready so that the perturbation fields can
! be derived from total fields substracting the base fields.
!
! 12/07/95 (Yuhe Liu)
! Added a new data dumping format, GRIB. The parameter definitions
! is mostly in consistance with NMC standard, despite those
! variables which are not in the WMO and NMC's tables (Table 2)
!
! 01/22/1996 (Yuhe Liu)
! Changed the ARPS history dumpings to dump the total values of
! ARPS variables, instead of the perturbated values.
!
! 12/09/1998 (Donghai Wang)
! Added the snow cover.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx,ny,nz The dimension of data arrays
!
! hinfmt The format of the history data dump
! =1, machine dependent unformatted binary dump,
! =2, formatted ascii dump,
! =3, NCSA HDF file format.
! =4, machine dependent packed unformatted binary dump
!
! nchanl FORTRAN I/O channel number for history data output.
!
! grdbasfn Name of the grid/base state array file
! lengbf Length of the grid/base state data file name string
! datafn Name of the other time dependent data file
! lendtf Length of the data file name string
!
! DATA ARRAYS READ IN:
!
! 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)
!
! 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
!
! tsfc Temperature at surface (K)
! tsoil Deep soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! 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
!
! 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:
!
! time The time of the input data (s)
! 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)
!
! 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
!
! tsfc Temperature at surface (K)
! tsoil Deep soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! 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
!
! 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
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: time ! The time of the input data (s)
INTEGER :: hinfmt ! The format of the history data dump
INTEGER :: lengbf ! Length of the grid/base state data
! file name string
INTEGER :: lendtf ! Length of the data file name string
CHARACTER(LEN=*) :: grdbasfn ! Name of the grid/base state array file
CHARACTER(LEN=*) :: datafn ! Name of the other time dependent
! data file
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 :: uprt (nx,ny,nz) ! x component of perturbation velocity
! (m/s)
REAL :: vprt (nx,ny,nz) ! y component of perturbation velocity
! (m/s)
REAL :: wprt (nx,ny,nz) ! vertical component of perturbation
! velocity in Cartesian coordinates
! (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) ! Rainwater 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 x velocity component (m/s)
REAL :: vbar (nx,ny,nz) ! Base state y velocity component (m/s)
REAL :: wbar (nx,ny,nz) ! Base state z velocity component (m/s)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: rhobar(nx,ny,nz) ! Base state air density (kg/m**3)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal)
REAL :: qvbar (nx,ny,nz) ! Base state water vapor mixing ratio
! (kg/kg)
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 :: tsfc (nx,ny,0:nstyps) ! Temperature at surface (K)
REAL :: tsoil (nx,ny,0:nstyps) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture
REAL :: wetdp (nx,ny,0:nstyps) ! Deep soil moisture
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 :: 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*s**2))
REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s))
REAL :: tem1(nx,ny,nz) ! Temporary work array
REAL :: tem2(nx,ny,nz) ! Temporary work array
REAL :: tem3(nx,ny,nz) ! Temporary work array
INTEGER :: grdread,iread
SAVE grdread
INTEGER :: nchanl ! FORTRAN I/O channel number for output
INTEGER :: istat
INTEGER :: ireturn ! Return status indicator
INTEGER :: grdbas ! Wether this is a grid/base state
! array dump
REAL :: btime ! The time of the base state data
REAL :: p0inv, amin, amax
INTEGER :: i,j,k
INTEGER :: ierr
LOGICAL :: fexist
INTEGER :: packed
CHARACTER (LEN=80) :: rname
INTEGER :: is
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'indtflg.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'mp.inc' ! mpi parameters.
INCLUDE 'phycst.inc'
DATA grdread /0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ireturn = 0
p0inv=1./p0
IF (lengbf < 1) lengbf = len_trim(grdbasfn)
IF (lendtf < 1) lendtf = len_trim(datafn)
!
!-----------------------------------------------------------------------
!
! Open and read grid and base state data file depending on the
! values of parameters grdin and basin, which are read in from the
! time dependent data set. If grdin or basin is zero, the grid and
! base state arrays have to be read in from a separate file.
!
!-----------------------------------------------------------------------
!
IF ( hinfmt == 9 ) GO TO 500
IF( grdread == 0 ) THEN
grdbas = 1
rname = runname
INQUIRE(FILE=trim(grdbasfn(1:lengbf)), EXIST = fexist )
IF( fexist ) GO TO 200
INQUIRE(FILE=trim(grdbasfn(1:lengbf))//'.Z', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs
( trim(grdbasfn(1:lengbf))//'.Z' )
GO TO 200
END IF
INQUIRE(FILE=trim(grdbasfn(1:lengbf))//'.gz', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( trim(grdbasfn(1:lengbf))//'.gz' )
GO TO 200
END IF
WRITE(6,'(/1x,a,/1x,a/)') &
'File '//trim(grdbasfn(1:lengbf))// &
' or its compressed version not found.', &
'Program stopped in DTAREAD.'
CALL arpsstop('arpsstop called from dtaread during base state read',1)
200 CONTINUE
!
!-----------------------------------------------------------------------
!
! Read grid and base state fields.
!
!-----------------------------------------------------------------------
!
IF( hinfmt == 1 ) THEN
CALL getunit( nchanl )
!
!-----------------------------------------------------------------------
!
! Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(trim(grdbasfn(1:lengbf)), '-F f77 -N ieee', ierr)
OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 999
CALL binread(nx,ny,nz,nstyps, grdbas, nchanl, btime, x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF( hinfmt == 2 ) THEN
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))), &
STATUS='old',FORM='formatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 999
CALL ascread(nx,ny,nz,nstyps,grdbas,nchanl,btime,x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF( hinfmt == 3 ) THEN
CALL hdfread(nx,ny,nz,nstyps,grdbas, trim(grdbasfn(1:lengbf)), &
btime,x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1)
ELSE IF( hinfmt == 4 ) THEN
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 999
! CALL pakread(nx,ny,nz,nstyps,grdbas,nchanl,btime,x,y,z,zp, &
! 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, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
! raing,rainc,prcrate, &
! radfrc,radsw,rnflx, &
! usflx,vsflx,ptsflx,qvsflx, &
! ireturn, tem1,tem2)
CALL pakread
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF( hinfmt == 6 ) THEN
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 999
CALL bn2read(nx,ny,nz,nstyps,grdbas,nchanl,btime,x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF (hinfmt == 7) THEN ! NetCDF format
packed = 0
! CALL netread(trim(grdbasfn(1:lengbf)),nx,ny,nz,nstyps,packed, &
! grdbas, x, y, z, zp, &
! 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, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
! raing,rainc,prcrate, &
! radfrc,radsw,rnflx, &
! usflx,vsflx,ptsflx,qvsflx, &
! tem1, tem2, tem3)
CALL netread
ELSE IF (hinfmt == 8) THEN ! NetCDF packed format
packed = 1
! CALL netread(trim(grdbasfn(1:lengbf)),nx,ny,nz,nstyps,packed, &
! grdbas, x, y, z, zp, &
! 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, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
! raing,rainc,prcrate, &
! radfrc,radsw,rnflx, &
! usflx,vsflx,ptsflx,qvsflx, &
! tem1, tem2, tem3)
CALL netread
ELSE IF( hinfmt == 10 ) THEN
!
!-----------------------------------------------------------------------
!
! Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(trim(grdbasfn(1:lengbf)), '-F f77 -N ieee', ierr)
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))),STATUS='old', &
FORM='unformatted',IOSTAT= istat )
CALL gribread(nx,ny,nz,nstyps,grdbas,nchanl,btime,x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE
WRITE(6,'(a,i3,a)') &
' Data format flag had an invalid value ', &
hinfmt ,' program stopped.'
CALL arpsstop('arpsstop called from dtaread during base state &
&read',1)
END IF
grdread = 1
runname = rname
END IF
500 CONTINUE
!
!-----------------------------------------------------------------------
!
! Read data fields.
!
!-----------------------------------------------------------------------
!
grdbas = 0
INQUIRE(FILE=trim(datafn(1:lendtf)), EXIST = fexist )
IF( fexist ) GO TO 100
INQUIRE(FILE=trim(datafn(1:lendtf))//'.Z', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( trim(datafn(1:lendtf))//'.Z' )
GO TO 100
END IF
INQUIRE(FILE=trim(datafn(1:lendtf))//'.gz', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( trim(datafn(1:lendtf))//'.gz' )
GO TO 100
END IF
WRITE(6,'(/1x,a,/1x,a/)') &
'File '//trim(datafn(1:lendtf)) &
//' or its compressed version not found.', &
'Program stopped in DTAREAD.'
CALL arpsstop('arpsstop called from dtaread during base read-2',1)
100 CONTINUE
IF( hinfmt == 1 ) THEN
!
!-----------------------------------------------------------------------
!
! Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(trim(datafn(1:lendtf)), '-F f77 -N ieee', ierr)
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 998
CALL binread(nx,ny,nz,nstyps,grdbas, nchanl, time, x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF( hinfmt == 2 ) THEN
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))), &
STATUS='old',FORM='formatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 998
CALL ascread(nx,ny,nz,nstyps,grdbas,nchanl,time,x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF( hinfmt == 3 ) THEN
CALL hdfread(nx,ny,nz,nstyps,grdbas,trim(datafn(1:lendtf)), &
time,x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1)
ELSE IF( hinfmt == 4 ) THEN
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 998
! CALL pakread(nx,ny,nz,nstyps,grdbas,nchanl,time,x,y,z,zp, &
! 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, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
! raing,rainc,prcrate, &
! radfrc,radsw,rnflx, &
! usflx,vsflx,ptsflx,qvsflx, &
! ireturn, tem1,tem2)
CALL pakread
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF( hinfmt == 6 ) THEN
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 998
CALL bn2read(nx,ny,nz,nstyps,grdbas, nchanl, time, x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF (hinfmt == 7) THEN ! NetCDF format
packed = 0
! CALL netread(trim(datafn(1:lendtf)),nx,ny,nz,nstyps,packed, &
! grdbas, x, y, z, zp, &
! 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, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
! raing,rainc,prcrate, &
! radfrc,radsw,rnflx, &
! usflx,vsflx,ptsflx,qvsflx, &
! tem1, tem2, tem3)
CALL netread
ELSE IF (hinfmt == 8) THEN ! NetCDF packed format
packed = 1
! CALL netread(trim(datafn(1:lendtf)),nx,ny,nz,nstyps,packed, &
! grdbas, x, y, z, zp, &
! 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, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
! raing,rainc,prcrate, &
! radfrc,radsw,rnflx, &
! usflx,vsflx,ptsflx,qvsflx, &
! tem1, tem2, tem3)
CALL netread
ELSE IF( hinfmt == 9 ) THEN
!
!-----------------------------------------------------------------------
!
! hinfmt = 9, read GrADS format data file. Since the DTAREAD
! subroutine doesn't define u, v, w, pt, p, qv, we have to use
! those base variables such ubar, vbar, wbar, ptbar, pbar, and
! qvbar as temporary arrays to store u, v, w, pt, p, qv.
!
!-----------------------------------------------------------------------
!
CALL gradsread(nx,ny,nz,nstyps, &
nchanl, trim(datafn(1:lendtf)), time, &
x,y,z,zp, uprt,vprt,wprt,ptprt,pprt, &
qvprt,qc,qr,qi,qs,qh, tke,kmh,kmv, &
ubar,vbar,ptbar,pbar,rhobar,qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1)
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
ubar (i,j,k) = ubar (i,j,k) - uprt (i,j,k)
vbar (i,j,k) = vbar (i,j,k) - vprt (i,j,k)
wbar (i,j,k) = 0.0
pbar (i,j,k) = pbar (i,j,k) - pprt (i,j,k)
ptbar(i,j,k) = ptbar(i,j,k) - ptprt(i,j,k)
qvbar(i,j,k) = qvbar(i,j,k) - qvprt(i,j,k)
END DO
END DO
END DO
grdin = 1
ELSE IF( hinfmt == 10 ) THEN
!
!-----------------------------------------------------------------------
!
! hinfmt = 10, read GRIB format data file.
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(trim(datafn(1:lendtf)), '-F f77 -N ieee', ierr)
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))),STATUS='old', &
FORM='unformatted',IOSTAT= istat )
CALL gribread(nx,ny,nz,nstyps, grdbas, nchanl, time, x,y,z,zp, &
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, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE
WRITE(6,'(a,i3,a)') &
' Data format flag had an invalid value ', &
hinfmt ,' program stopped.'
CALL arpsstop('arpsstop called from dtaread during read',1)
END IF
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
rhobar(i,j,k)= pbar(i,j,k)/ &
( rd * ptbar(i,j,k)*(pbar(i,j,k)*p0inv)**rddcp )
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set dx, dy and dz from read-in data
!
!-----------------------------------------------------------------------
!
dx = x(2) - x(1)
dy = y(2) - y(1)
dz = z(2) - z(1)
!
IF (myproc == 0) &
WRITE(6,'(/1x,a/)') 'Min. and max. of the data arrays read in:'
CALL edgfill(x,1,nx,1,nx,1,1,1,1,1,1,1,1)
CALL a3dmax0lcl(x,1,nx,1,nx,1,1,1,1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(/1x,2(a,e13.6))') 'xmin = ', amin,', xmax =',amax
CALL edgfill(y,1,ny,1,ny,1,1,1,1,1,1,1,1)
CALL a3dmax0lcl(y,1,ny,1,ny,1,1,1,1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'ymin = ', amin,', ymax =',amax
CALL edgfill(z,1,nz,1,nz,1,1,1,1,1,1,1,1)
CALL a3dmax0lcl(z,1,nz,1,nz,1,1,1,1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'zmin = ', amin,', zmax =',amax
CALL edgfill(zp,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz)
CALL a3dmax0lcl(zp,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'zpmin = ', amin,', zpmax =',amax
CALL edgfill(ubar,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(ubar,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'ubarmin = ', amin,', ubarmax =',amax
CALL edgfill(vbar,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1)
CALL a3dmax0lcl(vbar,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'vbarmin = ', amin,', vbarmax =',amax
CALL edgfill(wbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz)
CALL a3dmax0lcl(wbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'wbarmin = ', amin,', wbarmax =',amax
CALL edgfill(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'ptbarmin= ', amin,', ptbarmax=',amax
CALL edgfill(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'pbarmin = ', amin,', pbarmax =',amax
IF( mstin == 1) THEN
CALL edgfill(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'qvbarmin= ', amin, ', qvbarmax=',amax
END IF
CALL edgfill(uprt,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(uprt,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'uprtmin = ', amin,', uprtmax =',amax
CALL edgfill(vprt,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1)
CALL a3dmax0lcl(vprt,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'vprtmin = ', amin,', vprtmax =',amax
CALL edgfill(wprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz)
CALL a3dmax0lcl(wprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'wprtmin = ', amin,', wprtmax =',amax
CALL edgfill(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'ptprtmin= ', amin,', ptprtmax=',amax
CALL edgfill(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'pprtmin = ', amin,', pprtmax =',amax
IF( mstin == 1) THEN
CALL edgfill(qvprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(qvprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'qvprtmin= ', amin, &
', qvprtmax=',amax
CALL edgfill(qc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(qc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'qcmin = ', amin, &
', qcmax =',amax
CALL edgfill(qr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(qr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'qrmin = ', amin, &
', qrmax =',amax
IF( icein == 1) THEN
CALL edgfill(qi,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(qi,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'qimin = ', amin, &
', qimax =',amax
CALL edgfill(qs,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(qs,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'qsmin = ', amin, &
', qsmax =',amax
CALL edgfill(qh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(qh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') 'qhmin = ', amin, &
', qhmax =',amax
END IF
IF(rainin == 1) THEN
CALL edgfill(raing,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(raing,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'raingmin= ', amin,', raingmax=',amax
CALL edgfill(rainc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(rainc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'raincmin= ', amin,', raincmax=',amax
END IF
IF( prcin == 1 ) THEN
CALL edgfill(prcrate(1,1,1),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(prcrate(1,1,1),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'prcr1min= ', amin,', prcr1max=',amax
CALL edgfill(prcrate(1,1,2),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(prcrate(1,1,2),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'prcr2min= ', amin,', prcr2max=',amax
CALL edgfill(prcrate(1,1,3),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(prcrate(1,1,3),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'prcr3min= ', amin,', prcr3max=',amax
CALL edgfill(prcrate(1,1,4),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(prcrate(1,1,4),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'prcr4min= ', amin,', prcr4max=',amax
END IF
END IF
IF (tkein == 1) THEN
CALL edgfill(tke,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(tke,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'tkemin = ', amin,', tkemax =',amax
END IF
IF (trbin == 1) THEN
CALL edgfill(kmh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(kmh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'kmhmin = ', amin,', kmhmax =',amax
CALL edgfill(kmv,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(kmv,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'kmvmin = ', amin,', kmvmax =',amax
END IF
IF( sfcin == 1 ) THEN
DO is = 0, nstyp
IF(myproc == 0) &
PRINT*,'Max/min of soil model variables for type index ', is
CALL edgfill(tsfc(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(tsfc(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF(myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'tsfcmin = ', amin,', tsfcmax =',amax
CALL edgfill(tsoil(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(tsoil(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF(myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'tsoilmin= ', amin,', tsoilmax=',amax
CALL edgfill(wetsfc(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(wetsfc(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF(myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'wetsmin = ', amin,', wetsmax =',amax
CALL edgfill(wetdp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(wetdp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF(myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'wetdmin = ', amin,', wetdmax =',amax
CALL edgfill(wetcanp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(wetcanp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF(myproc == 0) &
WRITE(6,'(1x,2(a,e13.6))') &
'wetcmin = ', amin,', wetcmax =',amax
END DO
END IF
IF ( radin == 1 ) THEN
CALL edgfill(radfrc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0lcl(radfrc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'radfmin = ', amin,', radfmax =',amax
CALL edgfill(radsw,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(radsw,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'radswmin= ', amin,', radswmax=',amax
CALL edgfill(rnflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(rnflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'rnflxmin= ', amin,', rnflxmax=',amax
END IF
IF ( flxin == 1 ) THEN
CALL edgfill(usflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(usflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'usflxmin= ', amin,', usflxmax=',amax
CALL edgfill(vsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(vsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'vsflxmin= ', amin,', vsflxmax=',amax
CALL edgfill(ptsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(ptsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'ptflxmin= ', amin,', ptflxmax=',amax
CALL edgfill(qvsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL a3dmax0lcl(qvsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
IF (myproc == 0) &
WRITE(6,'(1x,2(a,e13.4))') &
'qvflxmin= ', amin,', qvflxmax=',amax
END IF
RETURN
!##################################################################
!##################################################################
!###### ######
!###### ENTRY SETGBRD ######
!###### ######
!##################################################################
!##################################################################
ENTRY setgbrd (iread)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Set the gridread parameter so that the grid/base state file
! will (iread=0) or will not (iread=1) be read. This is useful
! if more than one history file is to be accessed during a run
! (primarily in data assimilation experiments).
!
! As a default, the grid/base file is read on the first call to
! dtaread and not on subsequent calls. dtaread resets grdread
! to 1 after each time a grid/base file is read.
!
!-----------------------------------------------------------------------
!
grdread=iread
RETURN
! 997 CONTINUE
WRITE(6,'(1x,a,a,/1x,i3,a)') &
'Error occured when opening GRIB file ',trim(datafn(1:lendtf)), &
'with the FILE pointer ',nchanl,' Program stopped in DTAREAD.'
CALL arpsstop('arpsstop called from dtaread during Grid read',1)
998 CONTINUE
WRITE(6,'(1x,a,a,/1x,i3,a)') &
'Error occured when opening file ',trim(datafn(1:lendtf)), &
'using FORTRAN unit ',nchanl,' Program stopped in DTAREAD.'
CALL arpsstop('arpsstop called from dtaread during file read',1)
999 CONTINUE
WRITE(6,'(1x,a,a,/1x,i3,a)') &
'Error occured when opening file ',trim(grdbasfn(1:lengbf)), &
'using FORTRAN unit ',nchanl,' Program stopped in DTAREAD.'
CALL arpsstop('arpsstop called from dtaread during file read',1)
END SUBROUTINE dtaread
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE DTADUMP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE dtadump(nx,ny,nz,nstyps, & 24,35
houtfmt,nchout,filnam,grdbas,flcmprs, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv, &
ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the history data dump in various data formats by
! calling the appropriate history dump routine.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 8/23/1992
!
! 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.
!
! 4/22/1995 (M. Xue)
! Modified external boundary condition file naming format, so
! that the files can be directed to directory 'dirname'.
!
! 4/24/1995 (M. Xue)
! Changed boundary file naming convention.
!
! 12/07/95 (Yuhe Liu)
! Added a new data dumping format, GRIB. The parameter definitions
! is mostly in consistance with NMC standard, despite those
! variables which are not in the WMO and NMC's tables (Table 2)
!
! 01/22/1996 (Yuhe Liu)
! Changed the ARPS history dumpings to dump the total values of
! ARPS variables, instead of the perturbated values.
!
! 04/30/1997 (Fanyou Kong -- CMRP)
! Add Vis5D format conversion
!
! 12/09/1998 (Donghai Wang)
! Added the snow cover.
!
!-----------------------------------------------------------------------
!
! 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
!
! nstyps number of soil type
!
! houtfmt The format flag of history data dump
! nchout FORTRAN I/O channel number for history data output.
! filnam The name of history data dump file
! grdbas Flag indicating if this is a call for the data dump
! of grid and base state arrays only. If so, grdbas=1.
! flcmprs Option for file compression (0 or 1), not used in this
! routine. A identical variable, filcmprs, is passed
! through include file globcst.inc
!
! 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 x-velocity component (m/s)
! vbar Base state y-velocity component (m/s)
! wbar Base state vertical 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)
! hterain Terrain height (m)
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
!
! soiltyp Soil type
! vegtyp Vegetation type
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsfc Temperature at ground (K) (in top 1 cm layer)
! tsoil Deep soil temperature (K) (in deep 1 m layer)
! wetsfc Surface soil moisture in the top 1 cm layer
! wetdp Deep soil moisture in the deep 1 m layer
! 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
!
! 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.
! tem2 Temporary work array.
! tem3 Temporary work array.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: houtfmt ! Flag of history data dump format
INTEGER :: nchout ! FORTRAN I/O channel number for output
CHARACTER (LEN=80) :: filnam ! The name of history dump data file
INTEGER :: grdbas ! If this is a grid/base state array dump
INTEGER :: flcmprs ! Option for file compression (0 or 1).
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 x-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state y-velocity (m/s)
REAL :: wbar (nx,ny,nz) ! Base state z-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 :: hterain(nx,ny) ! Terrain height.
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as - d( zp )/d( x )
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as - d( zp )/d( y )
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z )
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 :: tsfc (nx,ny,0:nstyps) ! Temperature at surface (K)
REAL :: tsoil (nx,ny,0:nstyps) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture
REAL :: wetdp (nx,ny,0:nstyps) ! Deep soil moisture
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 rates (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 :: 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*s**2))
REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s))
REAL :: tem1 (nx,ny,nz) ! Temporary work array
REAL :: tem2 (nx,ny,nz) ! Temporary work array
REAL :: tem3 (nx,ny,nz) ! Temporary work array
CHARACTER (LEN=256) :: filnamr
INTEGER :: istat, nchout0, nchout1
INTEGER :: ierr
INTEGER :: packed
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
REAL :: zpmax
!
!-----------------------------------------------------------------------
!
! Declaration for dumping the external boundary variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
INTEGER :: itema, lenstr
INTEGER :: iyr, imon, idy, ihr, imin, isec
CHARACTER (LEN=15) :: ctime
CHARACTER (LEN=80) :: exbcfn
CHARACTER (LEN=80) :: temchar
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'bndry.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
INCLUDE 'exbc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF (sfcout == 1) snowout = 1 ! snowout only independetly zero
! for pre-snow cover versions.
!
! IF( houtfmt.eq.0 ) RETURN ! No data dump
!
!
IF( houtfmt == 1 ) THEN ! Unformatted binary data dump.
!
!-----------------------------------------------------------------------
!
! History dump of unformatted binary data
!
!-----------------------------------------------------------------------
!
CALL getunit( nchout )
!
!-----------------------------------------------------------------------
!
! Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(filnam, '-F f77 -N ieee', ierr)
OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new', &
FORM='unformatted',IOSTAT= istat )
IF( istat /= 0) GO TO 999
CALL bindump(nx,ny,nz,nstyps,nchout, grdbas, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
kmh,kmv, ubar,vbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1)
CLOSE(UNIT=nchout)
CALL retunit( nchout )
ELSE IF( houtfmt == 2 ) THEN ! Formatted ASCII data dump
!
!-----------------------------------------------------------------------
!
! History dump of formatted ASCII data
!
!-----------------------------------------------------------------------
!
CALL getunit( nchout )
OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new', &
FORM='formatted',IOSTAT= istat )
IF( istat /= 0) GO TO 999
CALL ascdump(nx,ny,nz,nstyps, nchout, grdbas, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
kmh,kmv, ubar,vbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1)
CLOSE(UNIT=nchout)
CALL retunit( nchout )
ELSE IF( houtfmt == 3 ) THEN ! HDF4 format dump.
!
!-----------------------------------------------------------------------
!
! History data dump in NCSA HDF4 file format.
!
!-----------------------------------------------------------------------
!
CALL hdfdump(nx,ny,nz,nstyps,filnam, grdbas, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
kmh,kmv,ubar,vbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1)
ELSE IF( houtfmt == 4 ) THEN ! Packed binary data dump.
!
!-----------------------------------------------------------------------
!
! History dump of packed unformatted binary data
!
!-----------------------------------------------------------------------
!
CALL getunit( nchout )
OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new', &
FORM='unformatted',IOSTAT= istat )
IF( istat /= 0) GO TO 999
! CALL pakdump(nx,ny,nz,nstyps,nchout, grdbas, &
! u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
! kmh,kmv, ubar,vbar,ptbar,pbar,rhobar,qvbar, &
! x,y,z,zp,hterain, j1,j2,j3, &
! soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
! raing,rainc,prcrate, &
! radfrc,radsw,rnflx, &
! usflx,vsflx,ptsflx,qvsflx, &
! tem1, tem2)
CALL pakdump
CLOSE(UNIT=nchout)
CALL retunit( nchout )
ELSE IF( houtfmt == 5 ) THEN ! Data dump for Savi3D
!
!-----------------------------------------------------------------------
!
! History dump for Savi3D.
! Data at all time levels are dumped into a single file.
!
!-----------------------------------------------------------------------
!
nchout0 = 79
CALL svidump(nx,ny,nz,nstyps,nchout0,filnam, grdbas, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
kmh,kmv,ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2)
ELSE IF( houtfmt == 6 ) THEN ! Unformatted binary data dump.
!
!-----------------------------------------------------------------------
!
! History dump of unformatted binary data
!
!-----------------------------------------------------------------------
!
CALL getunit( nchout )
OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new', &
FORM='unformatted',IOSTAT= istat )
IF( istat /= 0) GO TO 999
CALL bn2dump(nx,ny,nz,nstyps, nchout, grdbas, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
kmh,kmv, ubar,vbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1)
CLOSE(UNIT=nchout)
CALL retunit( nchout )
!
!-----------------------------------------------------------------------
!
! NetCDF format dump.
!
!-----------------------------------------------------------------------
!
ELSE IF (houtfmt == 7) THEN ! NetCDF format
packed = 0
! CALL netdump(filnam, nx,ny,nz,nstyps,packed, grdbas, &
! u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
! kmh,kmv, ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, &
! x,y,z,zp,hterain, j1,j2,j3, &
! soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
! raing,rainc,prcrate, &
! radfrc,radsw,rnflx, &
! usflx,vsflx,ptsflx,qvsflx, &
! tem1,tem2,tem3)
CALL netdump
!
!-----------------------------------------------------------------------
!
! Packed NetCDF format dump.
!
!-----------------------------------------------------------------------
!
ELSE IF (houtfmt == 8) THEN ! Packed NetCDF format
packed = 1
! CALL netdump(filnam, nx,ny,nz,nstyps,packed, grdbas, &
! u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
! kmh,kmv, ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, &
! x,y,z,zp,hterain, j1,j2,j3, &
! soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
! raing,rainc,prcrate, &
! radfrc,radsw,rnflx, &
! usflx,vsflx,ptsflx,qvsflx, &
! tem1,tem2,tem3)
CALL netdump
ELSE IF( houtfmt == 9 ) THEN ! Unformatted binary data dump.
!
CALL gradsdump(nx,ny,nz,nstyps,nchout, filnam, istager, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
kmh,kmv,ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2)
ELSE IF( houtfmt == 10 ) THEN ! GRIB data dump.
!
!-----------------------------------------------------------------------
!
! History dump of GRIB format data
!
!-----------------------------------------------------------------------
!
lenstr = 80
CALL strlnth( filnam, lenstr )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(filnam, '-F f77 -N ieee', ierr)
CALL getunit( nchout )
OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new', &
FORM='unformatted',IOSTAT= istat )
IF( istat /= 0) GO TO 999
!
CALL gribdump(nx,ny,nz,nstyps, &
nchout, filnam(1:lenstr),grdbas, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv, &
ubar,vbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2, &
tem3(1,1,1),tem3(1,1,2), &
tem3(1,1,3),tem3(1,1,4))
CLOSE(UNIT=nchout)
CALL retunit( nchout )
ELSE IF( houtfmt == 11 ) THEN ! Vis5D data dump.
lenstr = 80
CALL strlnth( filnam, lenstr )
CALL v5ddump(nx,ny,nz, nstyps,filnam(1:lenstr), istager, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
kmh,kmv,ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2)
END IF ! houtfmt switches
!
!-----------------------------------------------------------------------
!
! Compress file except for Savi3D dump.
!
!-----------------------------------------------------------------------
!
IF( houtfmt /= 0 .AND. houtfmt /= 5 .AND. &
houtfmt /= 9 .AND. houtfmt /= 11 .AND. filcmprs == 1 ) &
CALL cmprs( filnam )
!
!-----------------------------------------------------------------------
!
! Create ready file, indicating history dump writing is complete
!
!-----------------------------------------------------------------------
!
IF (readyfl == 1 .AND. houtfmt /= 0) THEN
WRITE (filnamr,'(a)') trim(filnam) // "_ready"
CALL getunit( nchout1 )
OPEN (UNIT=nchout1,FILE=trim(filnamr))
WRITE (nchout1,'(a)') trim(filnam)
CLOSE (nchout1)
CALL retunit ( nchout1 )
END IF
!
!-----------------------------------------------------------------------
!
! Write u,v,w,pt,p, and qv to an additional file for
! external boundary conditions.
!
!-----------------------------------------------------------------------
!
IF ( exbcdmp /= 0 .AND. grdbas /= 1 ) THEN
!
!-----------------------------------------------------------------------
!
! Write the ARPS predicted fields into external boundary format
! files. The purpose for doing so is to produce external boundary
! conditions from ARPS. They can be used to test the external
! BC code, as well as to run the experiments with the ARPS generated
! external external boundary conditions.
!
!-----------------------------------------------------------------------
!
CALL ctim2abss( year,month,day,hour,minute,second, itema )
itema = itema + INT(curtim)
CALL abss2ctim( itema,iyr,imon,idy,ihr,imin,isec )
WRITE (ctime,'(i4.4,2i2.2,a,3i2.2)') &
iyr,imon,idy,'.',ihr,imin,isec
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)
tem2(i,j,k) = pbar (i,j,k) + pprt (i,j,k)
END DO
END DO
END DO
CALL edgfill(u, 1,nx,1,nx, 1,ny,1,ny-1, 1,nz,1,nz-1)
CALL edgfill(v, 1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1)
CALL edgfill(w, 1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz )
CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
CALL edgfill(tem2,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
CALL edgfill(qv, 1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
exbcfn = runname(1:lfnkey)//'.'//ctime
lenstr = lfnkey + 16
IF (mp_opt > 0) THEN
WRITE(exbcfn,'(a,a,a,a,2i2.2)') &
runname(1:lfnkey),'.',ctime,'_',loc_x,loc_y
lenstr = lenstr + 5
END IF
IF( dirname /= ' ' ) THEN
temchar = exbcfn
exbcfn = dirname(1:ldirnam)//'/'//temchar
lenstr = lenstr + ldirnam + 1
END IF
CALL fnversn(exbcfn,lenstr)
WRITE(6,'(1x,a,a)') &
'Dumping to the external boundary format file: ', &
exbcfn(1:lenstr)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! Calculate rayklow, needed for exbcdmp=4 format, if not computed already
IF ( exbcdmp == 4 .and. rayklow <= 0) THEN
rayklow = nz-1
DO k=nz-1,2,-1
zpmax = zp(1,1,k)
DO j=1,ny-1
DO i=1,nx-1
zpmax = MAX( zp(i,j,k), zpmax )
END DO
END DO
IF( zpmax < zbrdmp ) THEN
rayklow = MAX(2, k+1)
EXIT
END IF
END DO
END IF
CALL writexbc(nx,ny,nz,exbcfn,lenstr,ctime, &
1,1,1,1,1,1, &
qcexout,qrexout,qiexout,qsexout,qhexout, &
u,v,w,tem1,tem2,qv,qc,qr,qi,qs,qh)
END IF
RETURN
999 CONTINUE
WRITE(6,'(1x,a,a,a,/1x,i3)') &
'Error occured when opening file ',filnam, &
'using FORTRAN unit ',nchout
CALL arpsstop(' Program stopped in DTADUMP.',1)
END SUBROUTINE dtadump
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GET_DIMS_FROM_DATA ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE get_dims_from_data(hinfmt,grdbasfn, nx,ny,nz,nstyps, ireturn) 14,17
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 7/17/2000
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! hinfmt The format of the history data dump
! grdbasfn Name of the grid/base state array file
!
! OUTPUT:
!
! nx,ny,nz The dimension of data arrays
! nstyps The number of soil types
! ireturn Return status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: hinfmt ! The format of the history data dump
CHARACTER (LEN=* ) :: grdbasfn ! Name of the grid/base state array file
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: nstyps ! Number of soil types
INTEGER :: ireturn,lengbf,nchanl,ierr,istat,packed,i
LOGICAL :: fexist
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ireturn = 0
!
!-----------------------------------------------------------------------
!
! Open and read grid and base state data file depending on the
! values of parameters grdin and basin, which are read in from the
! time dependent data set. If grdin or basin is zero, the grid and
! base state arrays have to be read in from a separate file.
!
!-----------------------------------------------------------------------
!
! IF ( hinfmt.eq.9 ) GOTO 500
INQUIRE(FILE=trim(grdbasfn), EXIST = fexist )
IF( fexist ) GO TO 200
INQUIRE(FILE=trim(grdbasfn)//'.Z', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( trim(grdbasfn)//'.Z' )
GO TO 200
END IF
INQUIRE(FILE=trim(grdbasfn)//'.gz', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( trim(grdbasfn)//'.gz' )
GO TO 200
END IF
WRITE(6,'(/1x,a,/1x,a/)') &
'File '//trim(grdbasfn)// &
' or its compressed version not found.', &
'Program stopped in get_dims_from_data.'
CALL arpsstop('arpsstop called from get_dims_from_data',1)
200 CONTINUE
!
!-----------------------------------------------------------------------
!
! Read grid and base state fields.
!
!-----------------------------------------------------------------------
!
IF( hinfmt == 1 .OR. hinfmt == 6 ) THEN
CALL getunit( nchanl )
!
!-----------------------------------------------------------------------
!
! Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(trim(grdbasfn), '-F f77 -N ieee', ierr)
OPEN(UNIT=nchanl,FILE=trim(grdbasfn), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 999
CALL bin_getdims(nchanl, nx,ny,nz,nstyps, ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF( hinfmt == 2 ) THEN
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(grdbasfn), &
STATUS='old',FORM='formatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 999
CALL asc_getdims(nchanl, nx,ny,nz,nstyps, ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF( hinfmt == 3 ) THEN
CALL hdf_getdims(trim(grdbasfn),nx,ny,nz,nstyps, ireturn)
ELSE IF( hinfmt == 4 ) THEN
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(grdbasfn), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 999
CALL pak_getdims(nchanl, nx,ny,nz,nstyps, ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE IF (hinfmt == 7.OR.hinfmt == 8) THEN ! NetCDF format
PRINT*,'NetCDF format is no longer supported by ARPS.'
PRINT*,'Please check the specified inupt data format.'
STOP
packed = 0
! CALL net_getdims(nchanl, packed, nx,ny,nz,nstyps, ireturn)
ELSE IF( hinfmt == 10 ) THEN
!
!-----------------------------------------------------------------------
!
! Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(trim(grdbasfn), '-F f77 -N ieee', ierr)
CALL getunit( nchanl )
OPEN(UNIT=nchanl,FILE=trim(grdbasfn),STATUS='old', &
FORM='unformatted',IOSTAT= istat )
CALL grib_getdims(nchanl, nx,ny,nz,nstyps, ireturn)
CLOSE(UNIT=nchanl)
CALL retunit( nchanl )
ELSE
WRITE(6,'(a,i3,a)') &
' Data format flag had an invalid value ', &
hinfmt ,' program stopped.'
CALL arpsstop('arpsstop called from get_dims_from_data-wrong flag',1)
END IF
! 500 CONTINUE
RETURN
999 CONTINUE
WRITE(6,'(1x,a,a,/1x,i3,a)') &
'Error occured when opening file ',trim(grdbasfn), &
'using FORTRAN unit ',nchanl,' Program stopped in get_dims_from_data.'
CALL arpsstop('arpsstop called from get_dims_from_data opening file',1)
END SUBROUTINE get_dims_from_data
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE BIN_GETDIMS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE bin_getdims(inch, nx,ny,nz,nstyps, ireturn) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 7/17/2000.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! inch Channel number for binary reading.
!
! OUTPUT:
!
! 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
!
! nstyps Number of soil types
!
! ireturn Return status indicator
! =0, successful read of all data
! =1, error reading data
! =2, end-of-file reached during read attempt
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: nstyps ! Number of soil types
INTEGER :: inch ! Channel number for binary reading
INTEGER :: ireturn ! Return status indicator
CHARACTER (LEN=10) :: tmunit
CHARACTER (LEN=40) :: fmtverin
INTEGER :: i,nocmnt
REAL :: time
CHARACTER (LEN=80) :: runname
INTEGER :: idummy,totin
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
READ(inch,ERR=110,END=120) fmtverin
READ(inch,ERR=110,END=120) runname
WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') runname
READ(inch,ERR=110,END=120) nocmnt
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
READ(inch,ERR=110,END=120)
END DO
END IF
READ(inch,ERR=110,END=120) time,tmunit
READ(inch,ERR=110,END=120) nx, ny, nz
WRITE(6,'(a,3i5,a)') 'nx,ny,nz read in from data are ',nx,ny,nz, &
' respectively.'
READ(inch,ERR=110,END=120) idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,totin, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
IF (totin /= 0) THEN
READ(inch,ERR=110,END=120) ! block of 20 REALs ...
READ(inch,ERR=110,END=120) nstyps,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
WRITE (6,*) "nstyps read in as",nstyps
ELSE
WRITE (6,*) "nstyps not defined in data, set to 4"
nstyps = 4
ENDIF
ireturn = 0
RETURN
110 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in BIN_GETDIMS'
ireturn=1
RETURN
120 CONTINUE
WRITE(6,'(/a/)') ' End of file reached in BIN_GETDIMS'
ireturn=2
RETURN
END SUBROUTINE bin_getdims
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ASC_GETDIMS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE asc_getdims(inch, nx,ny,nz,nstyps, ireturn) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 7/17/2000.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! inch Channel number for ASCII reading.
!
! OUTPUT:
!
! 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
!
! nstyps Number of soil types
!
! ireturn Return status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: nstyps ! Number of soil types
INTEGER :: inch ! Channel number for binary reading
INTEGER :: ireturn ! Return status indicator
CHARACTER (LEN=40) :: fmtverin
CHARACTER (LEN=10) :: tmunit
INTEGER :: i,nocmnt
REAL :: time
CHARACTER (LEN=80) :: runname
INTEGER :: idummy,totin
REAL:: rdummy
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
READ(inch,'(1x,a40)',ERR=110,END=120) fmtverin
READ(inch,'(1x,a80)',ERR=110,END=120) runname
WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') runname
READ(inch,'(1x,i4)',ERR=110,END=120) nocmnt
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
READ(inch,'(1x,a80)',ERR=110,END=120)
END DO
END IF
READ(inch,'(1x,e16.8,1x,a10)',ERR=110,END=120) time,tmunit
READ(inch,'(1x, 3i12)',ERR=110,END=120) nx,ny,nz
WRITE(6,'(a,3i5,a)') 'nx,ny,nz read in from data are ',nx,ny,nz, &
' respectively.'
READ(inch,'(1x,10i8)',ERR=110,END=120) &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,totin, & !should be modified
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
READ(inch,'(1x,8E16.10)',ERR=110,END=120) &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy
IF (totin /= 0) THEN
READ(inch,'(1x,10i8)',ERR=110,END=120) &
nstyps,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
READ(inch,'(1x,8E16.10)',ERR=110,END=120) ! block of 20 REALs ...
WRITE (6,*) "nstyps read in as",nstyps
ELSE
WRITE (6,*) "nstyps not defined in data, set to 4"
nstyps = 4
ENDIF
ireturn = 0
RETURN
110 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in ASC_GETDIMS.'
ireturn=1
RETURN
120 CONTINUE
WRITE(6,'(/a/)') ' End of file reached in ASC_GETDIMS.'
ireturn=2
RETURN
END SUBROUTINE asc_getdims
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HDF_GETDIMS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE hdf_getdims(filename, nx,ny,nz,nstyps, ireturn),6
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 7/17/2000.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! filename Channel number for binary reading.
!
! OUTPUT:
!
! 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
!
! ireturn Return status indicator
! =0, successful read of all data
! =1, error reading data
! =2, end-of-file reached during read attempt
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: stat, sd_id, sds_id
CHARACTER (LEN=*) :: filename
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: nstyps
INTEGER :: ireturn ! Return status indicator
CHARACTER (LEN=10) :: tmunit
CHARACTER (LEN=40) :: fmtverin
INTEGER :: i,nocmnt
REAL :: time
CHARACTER (LEN=80) :: runname
CHARACTER (LEN=80 ) :: cmnt(50) ! String of comments on this model run
INTEGER istat
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "HDFREAD: ERROR opening ", &
trim(filename)," for reading."
GO TO 110
END IF
!-----------------------------------------------------------------------
!
! Get dimensions of data in binary file and check against
! the dimensions passed to HDFREAD
!
!-----------------------------------------------------------------------
CALL hdfrdi(sd_id,"nx",nx,istat)
CALL hdfrdi(sd_id,"ny",ny,istat)
CALL hdfrdi(sd_id,"nz",nz,istat)
CALL hdfrdi(sd_id,"nstyp",nstyps,istat)
WRITE(6,'(a,3i5,a)') 'nx,ny,nz read in from data are ',nx,ny,nz, &
' respectively.'
WRITE (6,*) "nstyps read in as",nstyps
ireturn = 0
GO TO 130
!-----------------------------------------------------------------------
!
! Error during read
!
!-----------------------------------------------------------------------
110 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in HDF_GETDIMS.'
ireturn=1
GO TO 130
!-----------------------------------------------------------------------
!
! End-of-file during read
!
!-----------------------------------------------------------------------
! 120 CONTINUE
! WRITE(6,'(/a/)') ' End of file reached in HDF_GETDIMS.'
! ireturn=2
130 CONTINUE
!tmp stat = sfendacc(sd_id) ! is this necessary?
CALL hdfclose(sd_id,stat)
RETURN
END SUBROUTINE hdf_getdims
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE PAK_GETDIMS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE pak_getdims(inch, nx,ny,nz,nstyps, ireturn) 1,1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 7/17/2000.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! inch Channel number for binary reading.
!
! OUTPUT:
!
! 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
!
! ireturn Return status indicator
! =0, successful read of all data
! =1, error reading data
! =2, end-of-file reached during read attempt
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: nstyps ! Number of soil types
INTEGER :: inch ! Channel number for binary reading
INTEGER :: ireturn ! Return status indicator
CHARACTER (LEN=10) :: tmunit
CHARACTER (LEN=40) :: fmtverin
INTEGER :: i,nocmnt
REAL :: time
CHARACTER (LEN=80) :: runname
INTEGER :: bgrdin,bbasin,bvarin,mstin,brainin,bicein,btkein, &
btrbin,bsfcin,landin,totin,prcin,radin,flxin, &
snowcin,snowin
INTEGER :: ihdlen,ilen
PARAMETER (ihdlen=5000)
INTEGER :: ihead(ihdlen)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
READ(inch,ERR=110,END=120) ihead
! CALL decdhdr(nx,ny,nz,ihdlen,ihead,ilen,time, &
! bgrdin,bbasin,bvarin,mstin,brainin,bicein,btkein, &
! btrbin,bsfcin,landin,totin,prcin,radin,flxin, &
! snowcin,snowin,fmtverin,tmunit)
CALL decdhdr
WRITE (6,*) "nstyps not defined in data, set to 4"
nstyps = 4
ireturn =0
RETURN
110 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in PAK_GETDIMS.'
ireturn=1
RETURN
120 CONTINUE
WRITE(6,'(/a/)') ' End of file reached in PAK_GETDIMS.'
ireturn=2
RETURN
END SUBROUTINE pak_getdims
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GRIB_GETDIMS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE grib_getdims(inch, nx,ny,nz,nstyps, ireturn) 1,1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 7/17/2000.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! inch Channel number for binary reading.
!
! OUTPUT:
!
! 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
!
! nstyps NUmber of soil types.
!
! ireturn Return status indicator
! =0, successful read of all data
! =1, error reading data
! =2, end-of-file reached during read attempt
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: inch
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: nstyps ! Number of soil types
INTEGER :: ireturn ! Return status indicator
CHARACTER (LEN=10) :: tmunit
CHARACTER (LEN=40) :: fmtverin
INTEGER :: i,j,nocmnt
REAL :: time
CHARACTER (LEN=80) :: runname
CHARACTER (LEN=80 ) :: cmnt(50) ! String of comments on this model run
!
!-----------------------------------------------------------------------
!
! Parameters for GRIB packing
!
!-----------------------------------------------------------------------
!
INTEGER :: nbufsz
PARAMETER ( nbufsz = 800000 ) ! Size of GRIB buffer
CHARACTER (LEN=1) :: mgrib(nbufsz) ! Buffer to carry GRIB messages
INTEGER :: npos,hhdlen
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
READ (inch,ERR=110,END=120) (mgrib(i),i=1,3)
npos = 0
hhdlen = ICHAR(mgrib(npos+1))*65536 &
+ ICHAR(mgrib(npos+2))*256 &
+ ICHAR(mgrib(npos+3))
BACKSPACE (inch)
READ (inch,ERR=110,END=120) (mgrib(i),i=1,hhdlen)
npos = npos + 3
DO i=1,40
fmtverin(i:i) = mgrib(npos+i)
END DO
npos = npos + 40
DO i=1,80
runname(i:i) = mgrib(npos+i)
END DO
npos = npos + 80
nocmnt = ICHAR(mgrib(npos+1))
npos = npos + 1
IF( nocmnt > 0 ) THEN
DO j=1,nocmnt
DO i=1,80
cmnt(j)(i:i) = mgrib(npos+i)
END DO
npos = npos + 80
END DO
END IF
CALL ibm2flt( mgrib(npos+1), time )
npos = npos + 4
nx = ICHAR(mgrib(npos+1))*256 + ICHAR(mgrib(npos+2))
npos = npos + 2
ny = ICHAR(mgrib(npos+1))*256 + ICHAR(mgrib(npos+2))
npos = npos + 2
nz = ICHAR(mgrib(npos+1))*256 + ICHAR(mgrib(npos+2))
npos = npos + 2
WRITE (6,*) "nstyps not defined in data, set to 4"
nstyps = 4
ireturn = 0
RETURN
110 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in GRIB_GETDIMS'
ireturn=1
RETURN
120 CONTINUE
WRITE(6,'(/a/)') ' End of file reached in GRIB_GETDIMS'
ireturn=2
RETURN
RETURN
END SUBROUTINE grib_getdims
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EDGFILL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE edgfill(a,nx1,nx2,ibgn,iend,ny1,ny2,jbgn,jend, & 451
nz1,nz2,kbgn,kend)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Fill in the edges of a data array from the valid interior grid
! points so that the arrays are completely filled. This is done primarily
! for the benefit of the compression routine which must have a valid
! value in each array location.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! 7/14/92.
!
! 9/1/94 (Y. Lu)
! Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! a Array to have edges filled in
!
! nx1:nx2 Dimensioned range of the first index of array "a"
! ny1:ny2 Dimensioned range of the second index of array "a"
! nz1:nz2 Dimensioned range of the third index of array "a"
!
! ibgn,iend Valued range of the first index of array "a"
! jgbn,jend Valued range of the second index of array "a"
! kbgn,kend Valued range of the third index of array "a"
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
INTEGER :: nx1,nx2,ny1,ny2,nz1,nz2
REAL :: a(nx1:nx2,ny1:ny2,nz1:nz2)
INTEGER :: ibgn,iend,jbgn,jend,kbgn,kend
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Fill in top and bottom grid levels, as needed.
!
!-----------------------------------------------------------------------
!
IF( kbgn > nz1 ) THEN
DO k=nz1,(kbgn-1)
DO j=jbgn,jend
DO i=ibgn,iend
a(i,j,k)=a(i,j,kbgn)
END DO
END DO
END DO
END IF
IF( kend < nz2 ) THEN
DO k=(kend+1),nz2
DO j=jbgn,jend
DO i=ibgn,iend
a(i,j,k)=a(i,j,kend)
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Fill in north and south edges, as needed.
!
!-----------------------------------------------------------------------
!
IF( jbgn > ny1 ) THEN
DO k=nz1,nz2
DO j=ny1,(jbgn-1)
DO i=ibgn,iend
a(i,j,k)=a(i,jbgn,k)
END DO
END DO
END DO
END IF
IF( jend < ny2 ) THEN
DO k=nz1,nz2
DO j=(jend+1),ny2
DO i=ibgn,iend
a(i,j,k)=a(i,jend,k)
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Fill in east and west edges, as needed.
!
!-----------------------------------------------------------------------
!
IF( ibgn > nx1 ) THEN
DO k=nz1,nz2
DO j=ny1,ny2
DO i=nx1,(ibgn-1)
a(i,j,k)=a(ibgn,j,k)
END DO
END DO
END DO
END IF
IF( iend < nx2 ) THEN
DO k=nz1,nz2
DO j=ny1,ny2
DO i=(iend+1),nx2
a(i,j,k)=a(iend,j,k)
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE edgfill
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE IEDGFILL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE iedgfill(a,nx1,nx2,ibgn,iend,ny1,ny2,jbgn,jend, & 17
nz1,nz2,kbgn,kend)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Fill in the edges of a data array from the valid interior grid
! points so that the arrays are completely filled. This is done primarily
! for the benefit of the compression routine which must have a valid
! value in each array location.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! 7/14/92.
!
! 9/1/94 (Y. Lu)
! Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! a Array to have edges filled in
!
! nx1:nx2 Dimensioned range of the first index of array "a"
! ny1:ny2 Dimensioned range of the second index of array "a"
! nz1:nz2 Dimensioned range of the third index of array "a"
!
! ibgn,iend Valued range of the first index of array "a"
! jgbn,jend Valued range of the second index of array "a"
! kbgn,kend Valued range of the third index of array "a"
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
INTEGER :: nx1,nx2,ny1,ny2,nz1,nz2
INTEGER :: a(nx1:nx2,ny1:ny2,nz1:nz2)
INTEGER :: ibgn,iend,jbgn,jend,kbgn,kend
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Fill in top and bottom grid levels, as needed.
!
!-----------------------------------------------------------------------
!
IF( kbgn > nz1 ) THEN
DO k=nz1,(kbgn-1)
DO j=jbgn,jend
DO i=ibgn,iend
a(i,j,k)=a(i,j,kbgn)
END DO
END DO
END DO
END IF
IF( kend < nz2 ) THEN
DO k=(kend+1),nz2
DO j=jbgn,jend
DO i=ibgn,iend
a(i,j,k)=a(i,j,kend)
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Fill in north and south edges, as needed.
!
!-----------------------------------------------------------------------
!
IF( jbgn > ny1 ) THEN
DO k=nz1,nz2
DO j=ny1,(jbgn-1)
DO i=ibgn,iend
a(i,j,k)=a(i,jbgn,k)
END DO
END DO
END DO
END IF
IF( jend < ny2 ) THEN
DO k=nz1,nz2
DO j=(jend+1),ny2
DO i=ibgn,iend
a(i,j,k)=a(i,jend,k)
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Fill in east and west edges, as needed.
!
!-----------------------------------------------------------------------
!
IF( ibgn > nx1 ) THEN
DO k=nz1,nz2
DO j=ny1,ny2
DO i=nx1,(ibgn-1)
a(i,j,k)=a(ibgn,j,k)
END DO
END DO
END DO
END IF
IF( iend < nx2 ) THEN
DO k=nz1,nz2
DO j=ny1,ny2
DO i=(iend+1),nx2
a(i,j,k)=a(iend,j,k)
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE iedgfill