!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WRTSOIL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE wrtsoil( nx,ny,nzsoil,nstyps, soiloutfl, dx,dy, zpsoil, & 7,27
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
zpsoilout,tsoilout,qsoilout,wcanpout,snowdout, &
tsoil,qsoil,wetcanp,snowdpth,soiltyp )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Write the soil model variables into file soiloutfl.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 04/20/1995
!
! MODIFICATION HISTORY:
!
! 08/07/1995 (M. Xue)
! Added file name to the argument list.
!
! 2000/03/29 (Gene Bassett)
! Removed the globcst.inc include.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2000/03/29 (Gene Bassett)
! Added HDF4 format.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2000/10/27 (Gene Bassett)
! Added soiltyp to soil file to track consistency of soil data.
!
! 05/13/2002 (J. Brotzge)
! Modified arrays, call statements to allow for multiple soil schemes.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of model grid points in the x-dir. (east/west)
! ny Number of model grid points in the y-dir. (north/south)
! nzsoil Number of model grid points in the soil.
!
! soiloutfl Name of the soil model data output file
!
! zpsoilout Flag for output of soil level depth
! tsoilout Flag for output of soil temperature
! qsoilout Flag for output of soil moisture
! wcanpout Flag for output of Canopy water amount
!
! tsoil Soil temperature (K)
! qsoil Soil moisture (m3/m3)
! wetcanp Canopy water amount
! snowdpth Snow depth (m)
! soiltyp Soil type in model domain
!
! OUTPUT:
!
!
! Output written to the surface data file, soilinfl:
!
! nx Number of model grid points in the x-direction
! ny Number of model grid points in the y-direction
! nzsoil Number of model grid points in the soil
!
! mapproj Type of map projection used to setup the analysis grid.
! trulat1 1st real true latitude of map projection.
! trulat2 2nd real true latitude of map projection.
! trulon Real true longitude of map projection.
! sclfct Map scale factor. At latitude = trulat1 and trulat2
!
! dx Model grid spacing in the x-direction east-west (meters)
! dy Model grid spacing in the y-direction east-west (meters)
! ctrlat Lat. at the origin of the model grid (deg. N)
! ctrlon Lon. at the origin of the model grid (deg. E)
!
! tsoilout Flag for output of soil temperature
! qsoilout Flag for output of soil moisture
! wcanpout Flag for output of Canopy water amount
!
! tsoil Soil temperature (K)
! qsoil Soil moisture (m3/m3)
! wetcanp Canopy water amount
! snowdpth Snow depth (m)
! soiltyp Soil type in model domain
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx ! Number of grid points in the x-direction
INTEGER :: ny ! Number of grid points in the y-direction
INTEGER :: nzsoil ! Number of grid points in the soil
INTEGER :: nstyps ! Number of soil types at each grid point
CHARACTER (LEN=*) :: soiloutfl ! Name of the output file
REAL :: dx
REAL :: dy
INTEGER :: mapproj ! Map projection
REAL :: trulat1 ! 1st real true latitude of map projection
REAL :: trulat2 ! 2nd real true latitude of map projection
REAL :: trulon ! Real true longitude of map projection
REAL :: sclfct ! Map scale factor
REAL :: ctrlat ! Center latitude of the model domain (deg. N)
REAL :: ctrlon ! Center longitude of the model domain (deg. E)
INTEGER :: zpsoilout ! Flag for output of zpsoil
INTEGER :: tsoilout ! Flag for output of tsoil
INTEGER :: qsoilout ! Flag for output of qsoil
INTEGER :: wcanpout ! Flag for output of wetcanp
INTEGER :: stypout ! Flag for output of soiltyp (set to 1 if any of
! the above flags are set)
INTEGER :: snowdout ! Flag for output of snowdpth
REAL :: zpsoil (nx,ny,nzsoil) ! Soil depth (m)
REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K)
REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m3/m3)
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
REAL :: snowdpth(nx,ny) ! Snow depth (m)
INTEGER :: soiltyp (nx,ny,nstyps) ! Soil type in model domain
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,is
!wdt remove:
! INTEGER :: lenfl
INTEGER :: flunit
INTEGER :: idummy
REAL :: rdummy
INTEGER :: ierr
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! INTEGER(KIND=selected_int_kind(4)) :: itmp(1)
! ! unused array in hdf routines since NO COMPRESSION
! REAL :: atmp1(1),atmp2(1)
! ! unused arrays in hdf routines since NO COMPRESSION
INTEGER(KIND=selected_int_kind(4)) :: itmp(nx,ny,nzsoil,0:nstyps)
REAL :: atmp1(nzsoil),atmp2(nzsoil)
INTEGER :: stat, sd_id
!
! 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,fmtver500
INTEGER :: intver,intver500
PARAMETER (fmtver500='* 005.00 GrADS Soilvar Data',intver500=500)
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
idummy = 0
rdummy = 0.0
IF (soildmp == 0) RETURN
WRITE (6,'(a,a)') 'WRTSOIL: Opening file ',trim(soiloutfl)
IF ( tsoilout+qsoilout+wcanpout /= 0 ) THEN
stypout = 1
ELSE
stypout = 0
END IF
!-----------------------------------------------------------------------
!
! Write out in Fortran unformatted.
!
!-----------------------------------------------------------------------
IF (soildmp == 1) THEN
CALL getunit
( flunit )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(soiloutfl, '-F f77 -N ieee', ierr)
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 == intver500) THEN
fmtver=fmtver500
ELSE
WRITE(6,'(/1x,a,i10,a/)') &
'Data format, intver=',intver,', not found. The Job stopped.'
CALL arpsstop('arpstop called from wrtsoil.',1)
END IF
OPEN (UNIT=flunit, FILE=trim(soiloutfl), &
STATUS='unknown', FORM='unformatted', ACCESS='sequential')
WRITE (flunit) fmtver
WRITE (flunit) nx, ny, nzsoil
WRITE (flunit) mapproj,tsoilout, qsoilout, &
wcanpout, 0,snowdout, stypout, zpsoilout, &
idummy, idummy, idummy, idummy, idummy, &
idummy, idummy, idummy, idummy, nstyp
WRITE (flunit) dx, dy, ctrlat, ctrlon,trulat1, &
trulat2, trulon, sclfct, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy
IF ( zpsoilout /= 0) THEN
DO k=1,nzsoil
WRITE (flunit) ((zpsoil(i,j,k),i=1,nx),j=1,ny)
END DO
END IF
IF ( tsoilout /= 0 ) THEN
DO is=0,nstyp
DO k=1,nzsoil
WRITE (flunit) ((tsoil(i,j,k,is),i=1,nx),j=1,ny)
END DO
END DO
END IF
IF ( qsoilout /= 0 ) THEN
DO is=0,nstyp
DO k=1,nzsoil
WRITE (flunit) ((qsoil(i,j,k,is),i=1,nx),j=1,ny)
END DO
END DO
END IF
IF ( wcanpout /= 0 ) THEN
DO is=0,nstyp
WRITE (flunit) ((wetcanp(i,j,is),i=1,nx),j=1,ny)
END DO
END IF
IF ( snowdout /= 0 ) THEN
WRITE (flunit) ((snowdpth(i,j),i=1,nx),j=1,ny)
END IF
IF ( stypout /= 0 ) THEN
DO is=1,nstyp
WRITE (flunit) ((soiltyp(i,j,is),i=1,nx),j=1,ny)
END DO
ENDIF
CLOSE ( flunit )
CALL retunit ( flunit )
ELSE IF (soildmp == 2) THEN
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block
!-----------------------------------------------------------------------
!
! Write out in HDF4.
!
!-----------------------------------------------------------------------
!EMK 11 July 2002
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 == intver500) THEN
fmtver=fmtver500
ELSE
WRITE(6,'(/1x,a,i10,a/)') &
'Data format, intver=',intver,', not found. The Job stopped.'
CALL arpsstop('arpstop called from wrtsoil.',1)
END IF
!EMK END 11 July 2002
CALL hdfopen(trim(soiloutfl), 2, sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "WRTSOIL: ERROR creating HDF4 file: ", &
trim(soiloutfl)
RETURN
END IF
CALL hdfwrtc(sd_id,40,"fmtver",fmtver,stat)
CALL hdfwrti(sd_id, 'nstyp', nstyp, stat)
CALL hdfwrti(sd_id, 'nx', nx, stat)
CALL hdfwrti(sd_id, 'ny', ny, stat)
CALL hdfwrti(sd_id, 'nzsoil', nzsoil, stat)
CALL hdfwrtr(sd_id, 'dx', dx, stat)
CALL hdfwrtr(sd_id, 'dy', dy, stat)
CALL hdfwrti(sd_id, 'mapproj', mapproj, stat)
CALL hdfwrtr(sd_id, 'trulat1', trulat1, stat)
CALL hdfwrtr(sd_id, 'trulat2', trulat2, stat)
CALL hdfwrtr(sd_id, 'trulon', trulon, stat)
CALL hdfwrtr(sd_id, 'sclfct', sclfct, stat)
CALL hdfwrtr(sd_id, 'ctrlat', ctrlat, stat)
CALL hdfwrtr(sd_id, 'ctrlon', ctrlon, stat)
nstyp=max(1,nstyp)
IF ( zpsoilout /= 0) THEN
CALL hdfwrt3d(zpsoil,nx,ny,nzsoil,sd_id,0,0, &
'zpsoil','Soil level depth','m', &
itmp,atmp1,atmp2)
END IF
IF ( tsoilout /= 0 ) THEN
CALL hdfwrt4d(tsoil,nx,ny,nzsoil,nstyp+1,sd_id,0,0, &
'tsoil','Soil temperature','K', &
itmp,atmp1,atmp2)
END IF
IF ( qsoilout /= 0 ) THEN
CALL hdfwrt4d(qsoil,nx,ny,nzsoil,nstyp+1,sd_id,0,0, &
'qsoil','Soil moisture','m3/m3', &
itmp,atmp1,atmp2)
END IF
IF ( wcanpout /= 0 ) THEN
CALL hdfwrt3d(wetcanp,nx,ny,nstyp+1,sd_id,0,0, &
'wetcanp','Canopy water amount','fraction', &
itmp,atmp1,atmp2)
END IF
IF ( snowdout /= 0 ) THEN
CALL hdfwrt2d(snowdpth,nx,ny,sd_id,0,0, &
'snowdpth','Snow depth','m',itmp)
END IF
IF ( stypout /= 0 ) THEN
CALL hdfwrt3di(soiltyp,nx,ny,nstyp,sd_id,0,0, &
'soiltyp','Soil type','index')
ENDIF
! 200 CONTINUE
CALL hdfclose(sd_id,stat)
IF (stat /= 0) THEN
WRITE (6,*) "WRTSOIL: ERROR on closing file ",trim(soiloutfl), &
" (status",stat,")"
END IF
!wdt end block
! alternate dump format ...
END IF
RETURN
END SUBROUTINE wrtsoil
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE READSOIL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE readsoil( nx,ny,nzsoil,nstyps,soilfile,dx,dy,zpsoil, & 6,38
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
zpsoilin,tsoilin,qsoilin,wcanpin,snowdin, &
tsoil,qsoil,wetcanp,snowdpth,soiltyp )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read the soil variables from file soilfile.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 04/20/95
!
! MODIFICATION HISTORY:
!
! 2000/03/29 (Gene Bassett)
! Removed the globcst.inc include.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2000/03/29 (Gene Bassett)
! Added HDF4 format.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2000/10/27 (Gene Bassett)
! Added soiltyp to soil file to track consistency of soil data.
!
! 05/14/2002 (J. Brotzge)
! Added additional arrays, changed call statements to allow for
! multiple soil schemes.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of model grid points in the x-dir. (east/west)
! ny Number of model grid points in the y-dir. (north/south)
! nzsoil Number of model grid points in the soil.
!
! mapproj Type of map projection used to setup the analysis grid.
! trulat1 1st real true latitude of map projection.
! trulat2 2nd real true latitude of map projection.
! trulon Real true longitude of map projection.
! sclfct Map scale factor. At latitude = trulat1 and trulat2
!
! dx Model grid spacing in the x-direction east-west (meters)
! dy Model grid spacing in the y-direction east-west (meters)
! ctrlat Lat. at the origin of the model grid (deg. N)
! ctrlon Lon. at the origin of the model grid (deg. E)
!
! OUTPUT:
!
! tsoil Soil temperature (K)
! qsoil Soil moisture (m3/m3)
! wetcanp Canopy water amount
! snowdpth Snow depth (m)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx ! Number of grid points in the x-direction
INTEGER :: ny ! Number of grid points in the y-direction
INTEGER :: nzsoil ! Number of grid points in the soil.
INTEGER :: nstyps ! Number of soil types for each grid point
CHARACTER (LEN=*) :: soilfile ! Name of the soil file
REAL :: dx
REAL :: dy
INTEGER :: mapproj ! Map projection
REAL :: trulat1 ! 1st real true latitude of map projection
REAL :: trulat2 ! 2nd real true latitude of map projection
REAL :: trulon ! Real true longitude of map projection
REAL :: sclfct ! Map scale factor
REAL :: ctrlat ! Center latitude of the model domain (deg. N)
REAL :: ctrlon ! Center longitude of the model domain (deg. E)
INTEGER :: zpsoilin
INTEGER :: tsoilin
INTEGER :: qsoilin
INTEGER :: tsfcin ! for backward compatibility Zuwen He, 07/01/02
INTEGER :: wsfcin ! for backward compatibility Zuwen He, 07/01/02
INTEGER :: wdpin ! for backward compatibility Zuwen He, 07/01/02
INTEGER :: wcanpin
INTEGER :: snowdin
INTEGER :: snowcin
INTEGER :: stypin
REAL :: zpsoil (nx,ny,nzsoil) ! Soil depths (m)
REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K)
REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m3/m3)
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
INTEGER :: soiltyp (nx,ny,nstyps) ! Soil type in model domain
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL, ALLOCATABLE :: zpsoil_in (:,:,:) ! Soil level depth (m)
REAL, ALLOCATABLE :: tsoil_in (:,:,:,:) ! Soil temperature (K)
REAL, ALLOCATABLE :: qsoil_in (:,:,:,:) ! Soil moisture (m3/m3)
REAL, ALLOCATABLE :: wetcanp_in(:,:,:) ! Canopy water amount
INTEGER, ALLOCATABLE :: soiltyp_in (:,:,:) ! Soil type in model domain
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: nxin,nyin,nzsoilin
REAL :: dxin,dyin
INTEGER :: mprojin
INTEGER :: nstyp1,nstypin
REAL :: trlat1in, trlat2in, trlonin, sclfctin
REAL :: ctrlonin, ctrlatin
INTEGER :: flunit
INTEGER :: idummy
REAL :: rdummy
INTEGER :: i,j,k,is
INTEGER :: istat, ierr
INTEGER :: ireturn
CHARACTER (LEN=128) :: savename !Message passing code.
CHARACTER :: amm*2, ayear*4, aday*2
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! INTEGER(KIND=selected_int_kind(4)) :: itmp(1)
! ! unused array in hdf routines since NO COMPRESSION
! REAL :: atmp1(1),atmp2(1)
! ! unused arrays in hdf routines since NO COMPRESSION
INTEGER(KIND=selected_int_kind(4)) :: itmp(nx,ny,nzsoil,0:nstyps)
! unused array in hdf routines since NO COMPRESSION
REAL :: atmp1(nzsoil),atmp2(nzsoil)
! unused arrays in hdf routines since NO COMPRESSION
INTEGER :: stat, sd_id
!
! 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 GrADS Soilvar Data',intver410=410)
PARAMETER (fmtver500='* 005.00 GrADS Soilvar Data',intver500=500)
CHARACTER (LEN=40) :: fmtverin
REAL, allocatable :: tem1(:,:,:) ! Temporary array
!
!-----------------------------------------------------------------------
!
! Include file:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Open the surface data file. Read the parameters first, check if
! the data are consistant with the model. If everything is OK, then
! read the surface data, soiltyp, vegtyp, lai, and roufns.
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0) THEN
savename(1:128) = soilfile(1:128)
WRITE(soilfile, '(a,a,2i2.2)') trim(savename),'_',loc_x,loc_y
END IF
WRITE (6,*) "READSOIL: reading in external supplied soil data ", &
"from file ",trim(soilfile)
!-----------------------------------------------------------------------
!
! Read in header information.
!
!-----------------------------------------------------------------------
IF (soilfmt == 0) THEN
!-----------------------------------------------------------------------
!
! Fortran unformatted dump.
!
!-----------------------------------------------------------------------
CALL getunit( flunit )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(soilfile, '-F f77 -N ieee', ierr)
OPEN (UNIT=flunit,FILE=trim(soilfile),FORM='unformatted', &
STATUS='old',IOSTAT=istat)
IF( istat /= 0 ) THEN
WRITE(6,'(/1x,a,i2,/1x,a/)') &
'Error occured when opening the surface data file ' &
//soilfile//' using FORTRAN unit ',flunit, &
' Program stopped in READSOIL.'
CALL arpsstop("arpsstop called from READSOIL while opening file",1)
END IF
WRITE(6,'(/1x,a,/1x,a,i2/)') &
'This run will start from an external supplied soil ', &
'data file '//soilfile//' using FORTRAN unit ',flunit
READ (flunit,ERR=997) fmtverin
!
! 07/01/2002 Zuwen He
!
! The following code is not a safe practice.
!
! However, this may be the only way to distinguish versions
! prior to 500.
!
IF (fmtverin == fmtver500) THEN
intver=intver500
ELSE
WRITE(6,'(/1x,a/)') &
'WARNING: Incoming data format are older than version 5.00!!! '
END IF
WRITE(6,'(/1x,a,a/)') 'Incoming data format, fmtverin=',fmtverin
READ (flunit,ERR=998) nxin,nyin,nzsoilin
GOTO 996
997 WRITE(6,'(/1x,a,a/)') &
'Incoming data format: fmtver=fmtver410. Data read-in may be wrong.'
CLOSE (flunit)
OPEN (UNIT=flunit,FILE=trim(soilfile),FORM='unformatted', &
STATUS='old',IOSTAT=istat)
READ (flunit,ERR=998) nxin,nyin
nzsoilin=2
intver=intver410 ! there is no fmtverin prior to version 500
fmtver=fmtver410
WRITE(6,'(/1x,a/,a/)') &
'WARNING: Incoming data format are to read as version 4.10'
996 CONTINUE
IF (intver == intver410) THEN
READ (flunit,ERR=998) mprojin,tsfcin,tsoilin,wsfcin,wdpin, &
wcanpin,snowcin,snowdin,stypin,zpsoilin, &
idummy, idummy, idummy, idummy,idummy, &
idummy, idummy, idummy, idummy,nstypin
ELSE IF (intver >= intver500) THEN
READ (flunit,ERR=998) mprojin,tsoilin,qsoilin, &
wcanpin,snowcin,snowdin,stypin,zpsoilin, &
idummy, idummy, idummy, idummy,idummy, &
idummy, idummy, idummy, idummy,nstypin
END IF
READ (flunit,ERR=998) dxin,dyin, ctrlatin,ctrlonin,trlat1in, &
trlat2in,trlonin,sclfctin,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy
ELSE IF (soilfmt == 1) THEN !HDF4
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block
!-----------------------------------------------------------------------
!
! HDF4 format.
!
!-----------------------------------------------------------------------
CALL hdfopen(trim(soilfile), 1, sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "READSOIL: ERROR opening ", &
trim(soilfile)," for reading."
GO TO 998
END IF
CALL hdfrdc(sd_id,40,"fmtver",fmtverin,istat)
!
! 07/01/2002 Zuwen He
!
! The following code is a dangerous practice
! but it may be the only way to distinguish
! versions prior to 500.
!
IF (fmtverin == fmtver500) THEN
intver=intver500
ELSE
intver=intver410 ! prior to 500, there is no fmtver variable
istat=0
WRITE(6,'(/1x,a/,a/)') &
'WARNING: Incoming data format are older than version 5.00!!! ', &
'It is to be read as if it version 4.10!!! '
END IF
CALL hdfrdi(sd_id,"nstyp",nstypin,istat)
CALL hdfrdi(sd_id,"nx",nxin,istat)
CALL hdfrdi(sd_id,"ny",nyin,istat)
IF (intver >= intver500) THEN
CALL hdfrdi(sd_id,"nzsoil",nzsoilin,istat)
ELSE
nzsoilin = 2 ! prior to version 500, it is 2 layer soil
END IF
CALL hdfrdr(sd_id,"dx",dxin,istat)
CALL hdfrdr(sd_id,"dy",dyin,istat)
CALL hdfrdi(sd_id,"mapproj",mprojin,istat)
CALL hdfrdr(sd_id,"trulat1",trlat1in,istat)
CALL hdfrdr(sd_id,"trulat2",trlat2in,istat)
CALL hdfrdr(sd_id,"trulon",trlonin,istat)
CALL hdfrdr(sd_id,"sclfct",sclfctin,istat)
CALL hdfrdr(sd_id,"ctrlat",ctrlatin,istat)
CALL hdfrdr(sd_id,"ctrlon",ctrlonin,istat)
!wdt end block
! alternate dump format ...
ELSE IF (sfcfmt == 2) THEN !Read data directly
nstypin = 1
END IF
nstyp1 = MAX( nstypin, 1 )
ALLOCATE (zpsoil_in(nx,ny,nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "READSOIL: ERROR allocating zpsoil_in, returning"
RETURN
END IF
ALLOCATE (tsoil_in(nx,ny,nzsoil,0:nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "READSOIL: ERROR allocating tsoil_in, returning"
RETURN
END IF
ALLOCATE (qsoil_in(nx,ny,nzsoil,0:nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "READSOIL: ERROR allocating qsoil_in, returning"
RETURN
END IF
ALLOCATE (wetcanp_in(nx,ny,0:nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "READSOIL: ERROR allocating wetcanp_in, returning"
RETURN
END IF
ALLOCATE (soiltyp_in(nx,ny,nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "READSOIL: ERROR allocating soiltyp_in, returning"
RETURN
END IF
!-----------------------------------------------------------------------
!
! Check the data file for consistent grid parameters.
!
!-----------------------------------------------------------------------
IF (soilfmt /= 2) THEN !(OASIS testing, sfcfmt = 2)
CALL checkgrid2d(nx,ny,nxin,nyin, &
dx,dy,ctrlat,ctrlon, &
mapproj,trulat1,trulat2,trulon,sclfct, &
dxin,dyin,ctrlatin,ctrlonin, &
mprojin,trlat1in,trlat2in,trlonin,sclfctin,ireturn)
IF (ireturn /= 0) THEN
WRITE (6,*) "READSOIL: ERROR, grid parameter mismatch"
CALL arpsstop("arpsstop called from READSOIL parameter mismatch",1)
END IF
WRITE (6,'(a,a//a,i1/6(a,e15.8/))') &
' The map projection and griding information for the ', &
' surface data: ', &
' Projection: ', mprojin, &
' The 1st real true latitude: ', trlat1in, &
' The 2nd real true latitude: ', trlat2in, &
' The real true longitude: ', trlonin, &
' Map scale factor: ', sclfctin, &
' Latitude at the origin: ', ctrlatin, &
' Longitude at the origin: ', ctrlonin
! ENDIF
IF (nzsoil /= nzsoilin) THEN
WRITE(6,*) &
'ERROR -- nzsoilin in soil file not equal to expected nzsoil.'
WRITE(6,*)'nzsoil = ',nzsoil,' nzsoilin = ',nzsoilin
CALL arpsstop("arpsstop called from READSOIL parameter mismatch",1)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Read in the soil data from the soil data file.
!
!-----------------------------------------------------------------------
!
IF (intver == intver410) THEN
ALLOCATE (tem1(nx,ny,0:nstyps),stat=istat) ! for reading old version
! tsfc,tsoil,wetsfc,wetdp
END IF
IF (soilfmt == 0) THEN ! Fortran unformatted
IF (intver == intver410) THEN
WRITE (6, '(a/8(a,i2/))') &
' Surface data flags for: ', &
' zpsoilin --', zpsoilin, &
' tsfcin --', tsfcin, &
' tsoilin --', tsoilin, &
' wsfcin --', wsfcin, &
' wdpin --', wdpin, &
' wcanpin --', wcanpin, &
' snowdin --', snowdin, &
' stypin --', stypin
ELSE IF (intver == intver500) THEN
WRITE (6, '(a/6(a,i2/))') &
' Surface data flags for: ', &
' zpsoilin --', zpsoilin, &
' tsoilin --', tsoilin, &
' qsoilin --', qsoilin, &
' wcanpin --', wcanpin, &
' snowdin --', snowdin, &
' stypin --', stypin
END IF
IF (intver == intver410) THEN
zpsoil_in(:,:,1)=0.
zpsoil_in(:,:,2)=-1.
ELSE IF (intver >= intver500) THEN
IF (zpsoilin /= 0) THEN
WRITE(6,'(a)') 'Read in the soil depth data'
DO k=1,nzsoilin
READ (flunit,ERR=998) ((zpsoil_in(i,j,k),i=1,nxin),j=1,nyin)
END DO
END IF
END IF ! intver
IF (intver == intver410) THEN
IF ( tsfcin /= 0 ) THEN
WRITE(6,'(a)') 'Read in the surface skin temperature data'
IF ( nstyp1 == 1 ) THEN
READ (flunit,ERR=998) ((tem1(i,j,0),i=1,nxin),j=1,nyin)
tsoil_in(:,:,1,0)=tem1(:,:,0)
ELSE
READ (flunit,ERR=998) tem1
tsoil_in(:,:,1,:)=tem1(:,:,:)
END IF
ELSE
WRITE(6,'(a)') 'Variable tsfc is not in the data set.'
END IF
IF ( tsoilin /= 0 ) THEN
WRITE(6,'(a)') 'Read in the deep soil temperature data'
IF ( nstyp1 == 1 ) THEN
READ (flunit,ERR=998) ((tem1(i,j,0),i=1,nxin),j=1,nyin)
tsoil_in(:,:,2,0)=tem1(:,:,0)
ELSE
READ (flunit,ERR=998) tem1
tsoil_in(:,:,2,:)=tem1(:,:,:)
END IF
ELSE
WRITE(6,'(a)') 'Variable tsoil is not in the data set.'
END IF
IF ( wsfcin /= 0 ) THEN
WRITE(6,'(a)') 'Read in the skin soil moisture data'
IF ( nstyp1 == 1 ) THEN
READ (flunit,ERR=998) ((tem1(i,j,0),i=1,nxin),j=1,nyin)
qsoil_in(:,:,1,0)=tem1(:,:,0)
ELSE
READ (flunit,ERR=998) tem1
qsoil_in(:,:,1,:)=tem1(:,:,:)
END IF
ELSE
WRITE(6,'(a)') 'Variable wetsfc is not in the data set.'
END IF
IF ( wdpin /= 0 ) THEN
WRITE(6,'(a)') 'Read in the deep soil moisture data'
IF ( nstyp1 == 1 ) THEN
READ (flunit,ERR=998) ((tem1(i,j,0),i=1,nxin),j=1,nyin)
qsoil_in(:,:,2,0)=tem1(:,:,0)
ELSE
READ (flunit,ERR=998) tem1
qsoil_in(:,:,2,:)=tem1(:,:,:)
END IF
ELSE
WRITE(6,'(a)') 'Variable wetdp is not in the data set.'
END IF
IF ( wcanpin /= 0 ) THEN
IF ( nstyp1 == 1 ) THEN
READ (flunit,ERR=998) ((wetcanp_in(i,j,0),i=1,nxin),j=1,nyin)
ELSE
READ (flunit,ERR=998) wetcanp_in
END IF
ELSE
WRITE (6, '(a)') 'Variable wetcanp is not in the data set.'
END IF
IF ( snowcin /= 0 ) THEN
WRITE (6, '(a)') 'File contains snowcvr -- discarding'
READ (flunit,ERR=998)
END IF
IF ( snowdin /= 0 ) THEN
WRITE (6, '(a)') 'Read in the snow depth data'
READ (flunit,ERR=998) snowdpth
ELSE
WRITE (6, '(a)') 'Variable snowdpth is not in the data set.'
END IF
IF ( stypin /= 0 ) THEN
WRITE (6, '(a)') 'Read soil type of soil data.'
READ (flunit,ERR=998) soiltyp_in
END IF
ELSE IF (intver >= intver500) THEN
IF ( tsoilin /= 0 ) THEN
DO is=0,nstypin
WRITE(6,'(a,i4)') 'Read in the soil temperature for soil type ',is
DO k=1,nzsoilin
READ (flunit,ERR=998) ((tsoil_in(i,j,k,is),i=1,nxin),j=1,nyin)
END DO
END DO
ELSE
WRITE(6,'(a)') 'Variable tsoil is not in the data set.'
END IF
IF ( qsoilin /= 0 ) THEN
DO is=0,nstypin
WRITE(6,'(a,i4)') 'Read in the soil moisture data for soil type ',is
DO k=1,nzsoilin
READ (flunit,ERR=998) ((qsoil_in(i,j,k,is),i=1,nxin),j=1,nyin)
END DO
END DO
ELSE
WRITE(6,'(a)') 'Variable qsoil is not in the data set.'
END IF
IF ( wcanpin /= 0 ) THEN
DO is=0,nstypin
WRITE (6, '(a,i4)') 'Read in the canopy water amount data for soil type ',is
READ (flunit,ERR=998) ((wetcanp_in(i,j,is),i=1,nxin),j=1,nyin)
END DO
ELSE
WRITE (6, '(a)') 'Variable wetcanp is not in the data set.'
END IF
IF ( snowcin /= 0 ) THEN
WRITE (6, '(a)') 'File contains snowcvr -- discarding'
READ (flunit,ERR=998)
END IF
IF ( snowdin /= 0 ) THEN
WRITE (6, '(a)') 'Read in the snow depth data'
READ (flunit,ERR=998) ((snowdpth(i,j),i=1,nxin),j=1,nyin)
ELSE
WRITE (6, '(a)') 'Variable snowdpth is not in the data set.'
END IF
IF ( stypin /= 0 ) THEN
DO is=1,nstypin
WRITE (6, '(a,i4)') 'Read soil type of soil data for soil type ',is
READ (flunit,ERR=998) ((soiltyp_in(i,j,is),i=1,nxin),j=1,nyin)
END DO
END IF
END IF
CLOSE ( flunit )
CALL retunit ( flunit )
ELSE IF (soilfmt == 1) THEN
IF (intver <= intver410) THEN
WRITE(6,'(a)') 'WARNING: No zpsoil is defined in this version. '
WRITE(6,'(a)') 'Assume zpsoil_in(,,1)=0 and zpsoil(,,2)=-1.'
zpsoil_in(:,:,1)=0.
zpsoil_in(:,:,2)=-1.
CALL hdfrd3d(sd_id,"tsfc",nxin,nyin,nstyp1+1,tem1,stat, &
itmp,atmp1,atmp2)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE(6,'(a)') 'Read in the surface skin temperature data'
tsfcin = 1
ELSE
WRITE(6,'(a)') 'Variable tsfc is not in the data set.'
tsfcin = 0
END IF
tsoil_in(:,:,1,:)=tem1(:,:,:)
CALL hdfrd3d(sd_id,"tsoil",nxin,nyin,nstyp1+1,tem1,stat, &
itmp,atmp1,atmp2)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE(6,'(a)') 'Read in the deep soil temperature data'
tsoilin = 1
ELSE
WRITE(6,'(a)') 'Variable tsoil is not in the data set.'
tsoilin = 0
END IF
tsoil_in(:,:,2,:)=tem1(:,:,:)
CALL hdfrd3d(sd_id,"wetsfc",nxin,nyin,nstyp1+1,tem1,stat, &
itmp,atmp1,atmp2)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE(6,'(a)') 'Read in the skin soil moisture data'
wsfcin = 1
ELSE
WRITE(6,'(a)') 'Variable wetsfc is not in the data set.'
wsfcin = 0
END IF
qsoil_in(:,:,1,:)=tem1(:,:,:)
CALL hdfrd3d(sd_id,"wetdp",nxin,nyin,nstyp1+1,tem1,stat, &
itmp,atmp1,atmp2)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE(6,'(a)') 'Read in the deep soil moisture data'
wdpin = 1
ELSE
WRITE(6,'(a)') 'Variable wetdp is not in the data set.'
wdpin = 0
END IF
qsoil_in(:,:,2,:)=tem1(:,:,:)
ELSE IF (intver >= intver500) THEN
CALL hdfrd3d(sd_id,"zpsoil",nxin,nyin,nzsoil,zpsoil_in,stat, &
itmp,atmp1,atmp2)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE(6,'(a)') 'Read in the soil layer depth data'
zpsoilin = 1
ELSE
WRITE(6,'(a)') 'Variable zpsoil is not in the data set.'
zpsoilin = 0
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block
CALL hdfrd4d(sd_id,"tsoil",nxin,nyin,nzsoil,nstyp1+1,tsoil_in,stat, &
itmp,atmp1,atmp2)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE(6,'(a)') 'Read in the soil temperature data'
tsoilin = 1
ELSE
WRITE(6,'(a)') 'Variable tsoil is not in the data set.'
tsoilin = 0
END IF
CALL hdfrd4d(sd_id,"qsoil",nxin,nyin,nzsoil,nstyp1+1,qsoil_in,stat, &
itmp,atmp1,atmp2)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE(6,'(a)') 'Read in the soil moisture data'
qsoilin = 1
ELSE
WRITE(6,'(a)') 'Variable qsoil is not in the data set.'
qsoilin = 0
END IF
END IF
CALL hdfrd3d(sd_id,"wetcanp",nxin,nyin,nstyp1+1,wetcanp_in,stat, &
itmp,atmp1,atmp2)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read in the canopy water amount data'
wcanpin = 1
ELSE
WRITE (6, '(a)') 'Variable wetcanp is not in the data set.'
wcanpin = 0
END IF
CALL hdfrd2d(sd_id,"snowdpth",nxin,nyin,snowdpth,stat,itmp)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read in the snow depth data'
snowdin = 1
ELSE
WRITE (6, '(a)') 'Variable snowdpth is not in the data set.'
snowdin = 0
END IF
CALL hdfrd3di(sd_id,"soiltyp",nxin,nyin,nstyp1,soiltyp_in,stat)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read soil type of soil data'
stypin = 1
ELSE
WRITE (6, '(a)') 'Soil type of soil data is not in the data set.'
stypin = 0
END IF
CALL hdfclose(sd_id, stat)
!wdt end block
! alternate dump format ...
!OASIS code
ELSE IF (soilfmt == 2) THEN !Data read directly (JAB)
CALL initztime(ayear,amm,aday)
CALL readjsoil(nx,ny,nzsoil, nstyp,ayear,amm,aday,ztime, &
zpsoil,tsoil,qsoil,wetcanp,snowdpth)
END IF
!-------------------------------------------------------------------------
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block:
! nstypin = MIN(nstyp1, nstyps)
!
! IF (tsoilin == 1) tsoil(:,:,:,0:nstypin) = tsoil_in(:,:,:,0:nstypin)
! IF (qsoilin == 1) qsoil(:,:,:,0:nstypin) = qsoil_in(:,:,:,0:nstypin)
! IF (wcanpin == 1) wetcanp(:,:,0:nstypin) = wetcanp_in(:,:,0:nstypin)
!
! CALL fix_soil_nstyp(nx,ny,nzsoil,nstyp1,nstyp,tsoil, &
! qsoil,wetcanp)
IF (stypin == 0) THEN
WRITE (6,*) "READSOIL: WARNING, no check made for consistency ", &
"between soil types in surface and soil data sets."
nstypin = MIN(nstyp1, nstyps)
IF (zpsoilin == 1) zpsoil(:,:,:) = zpsoil_in(:,:,:)
IF (tsoilin == 1) tsoil(:,:,:,0:nstypin) = tsoil_in(:,:,:,0:nstypin)
IF (qsoilin == 1) qsoil(:,:,:,0:nstypin) = qsoil_in(:,:,:,0:nstypin)
IF (wcanpin == 1) wetcanp(:,:,0:nstypin) = wetcanp_in(:,:,0:nstypin)
CALL fix_soil_nstyp(nx,ny,nzsoil,nstyp1,nstyp,tsoil,qsoil, &
wetcanp)
ELSE
IF (nstyp1 == 1) THEN
tsoil_in (:,:,:,1) = tsoil_in (:,:,:,0)
qsoil_in (:,:,:,1) = qsoil_in (:,:,:,0)
wetcanp_in(:,:,1) = wetcanp_in(:,:,0)
ENDIF
IF (soiltestmeso /= 1) THEN
CALL remap_soil_vars(nx,ny,nzsoil,nstyp1,nstyp, &
tsoil_in,qsoil_in,wetcanp_in,soiltyp_in, &
tsfcin,tsoilin,wsfcin,wdpin,qsoilin,wcanpin, &
intver, &
tsoil,qsoil,wetcanp,soiltyp)
ENDIF
ENDIF
IF (mp_opt > 0) soilfile(1:128) = savename(1:128)
IF (intver == intver410) THEN
DEALLOCATE (tem1) ! for reading old version tsfc,tsoil,wetsfc,wetdp
END IF
DEALLOCATE (tsoil_in,stat=istat)
DEALLOCATE (qsoil_in,stat=istat)
DEALLOCATE (wetcanp_in,stat=istat)
DEALLOCATE (soiltyp_in,stat=istat)
IF (intver == intver410) THEN
IF (tsfcin /= tsoilin .OR. wsfcin /= wdpin) THEN
WRITE (6,'(a,a/,a/,a/)') &
'READSOIL: WARNING: The soilvar data is of version ', fmtver410, &
'. The inconsistency flag between tsfcin and tsoilin, ', &
' or between wsfin and wdpin, may cause some problems. '
END IF
tsoilin = max(tsfcin,tsoilin)
qsoilin = max(wsfcin,wdpin)
END IF
! Correct only for flint, otherwise flint will never end
! RETURN
GO TO 999
998 WRITE (6,'(/a,i2/a)') &
' Read error in surface data file ' &
//soilfile//' with the I/O unit ',flunit, &
'The model will STOP in subroutine READSOIL.'
CALL arpsstop("arpsstop called from READSOIL reading surface data",1)
999 CONTINUE
RETURN
END SUBROUTINE readsoil
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WRTSFCDT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE wrtsfcdt( nx,ny,nstyps, sfcoutfl, dx,dy, & 4,24
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
stypout,vtypout,laiout,rfnsout,vegout,ndviout, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg,ndvi )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Write a surface property data .
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 2/20/94
!
! MODIFICATION HISTORY:
!
! 2000/03/29 (Gene Bassett)
! Removed the globcst.inc include.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2000/03/29 (Gene Bassett)
! Added HDF4 format.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of model grid points in the x-dir. (east/west)
! ny Number of model grid points in the y-dir. (north/south)
! sfcoutfl Name of the output surface property data file
!
! soiltyp Soil type in model domain
! vegtyp Vegetation type in model domain
! lai Leaf Area Index in model domain
! roufns Surface roughness
! veg Vegetation fraction
! ndvi NDVI
!
! OUTPUT:
!
! Output written to the surface data file, sfcdtfl:
!
! nx Number of model grid points in the x-direction
! ny Number of model grid points in the y-direction
!
! mapproj Type of map projection used to setup the analysis grid.
! trulat1 1st real true latitude of map projection.
! trulat2 2nd real true latitude of map projection.
! trulon Real true longitude of map projection.
! sclfct Map scale factor. At latitude = trulat1 and trulat2
!
! dx Model grid spacing in the x-direction east-west (meters)
! dy Model grid spacing in the y-direction east-west (meters)
! ctrlat Lat. at the origin of the model grid (deg. N)
! ctrlon Lon. at the origin of the model grid (deg. E)
!
! stypout Flag for output of soil type
! vtypout Flag for output of vegetation type
! laiout Flag for output of Leaf Area Index
! rfnsout Flag for output of surface roughness
! vegout Flag for output of vegetation fraction
! ndviout Flag for output of NDVI
!
! soiltyp Soil type in model domain
! stypfrct Fraction of each soil type in each grid box
! vegtyp Vegetation type in model domain
! lai Leaf Area Index in model domain
! roufns Surface roughness
! veg Vegetation fraction
! ndvi NDVI
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx ! Number of grid points in the x-direction
INTEGER :: ny ! Number of grid points in the y-direction
INTEGER :: nstyps ! Max number of soil types in a grid box
CHARACTER (LEN=* ) :: sfcoutfl ! Surface property data file name
REAL :: dx
REAL :: dy
INTEGER :: mapproj ! Map projection
REAL :: trulat1 ! 1st real true latitude of map projection
REAL :: trulat2 ! 2nd real true latitude of map projection
REAL :: trulon ! Real true longitude of map projection
REAL :: sclfct ! Map scale factor
REAL :: ctrlat ! Center latitude of the model domain (deg. N)
REAL :: ctrlon ! Center longitude of the model domain (deg. E)
INTEGER :: stypout ! Flag for output of soiltyp
INTEGER :: vtypout ! Flag for output of vegtyp
INTEGER :: laiout ! Flag for output of lai
INTEGER :: rfnsout ! Flag for output of roufns
INTEGER :: vegout ! Flag for output of veg
INTEGER :: ndviout ! Flag for output of ndvi
INTEGER :: soiltyp (nx,ny,nstyps) ! Soil type in model domain
REAL :: stypfrct(nx,ny,nstyps) ! Fraction of soil types
INTEGER :: vegtyp (nx,ny) ! Vegetation type in model domain
REAL :: lai (nx,ny) ! Leaf Area Index in model domain
REAL :: roufns (nx,ny) ! Surface roughness
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: ndvi (nx,ny) ! NDVI
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
!wdt remove:
! INTEGER :: lenfl
INTEGER :: flunit
INTEGER :: idummy
REAL :: rdummy
INTEGER :: ierr
INTEGER :: i,j,is
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
INTEGER(KIND=selected_int_kind(4)) :: itmp(1)
! unused array in hdf routines since NO COMPRESSION
REAL :: atmp1(1),atmp2(1)
! unused arrays in hdf routines since NO COMPRESSION
INTEGER :: stat, sd_id
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
idummy = 0
rdummy = 0.0
IF (sfcdmp == 0) RETURN
WRITE (6,'(a,a)') 'WRTSFCDT: Opening file ',trim(sfcoutfl)
!-----------------------------------------------------------------------
!
! Write out in Fortran unformatted.
!
!-----------------------------------------------------------------------
IF (sfcdmp == 1) THEN
CALL getunit( flunit )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(sfcoutfl, '-F f77 -N ieee', ierr)
OPEN (UNIT = flunit, FILE = trim(sfcoutfl), &
STATUS = 'unknown', &
FORM = 'unformatted', ACCESS = 'sequential')
WRITE (flunit) nx, ny
WRITE (flunit) mapproj, stypout, vtypout, laiout, rfnsout, &
vegout, ndviout, nstyp, idummy, idummy, &
idummy, idummy, idummy, idummy, idummy, &
idummy, idummy, idummy, idummy, idummy
WRITE (flunit) dx, dy, ctrlat, ctrlon, trulat1, &
trulat2, trulon, sclfct, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy
IF ( stypout /= 0 ) THEN
IF ( nstyp <= 1 ) THEN
WRITE (flunit) ((soiltyp(i,j,1),i=1,nx),j=1,ny)
ELSE
DO is=1,nstyp
WRITE (flunit) ((soiltyp (i,j,is),i=1,nx),j=1,ny)
WRITE (flunit) ((stypfrct(i,j,is),i=1,nx),j=1,ny)
END DO
END IF
END IF
IF ( vtypout /= 0 ) WRITE (flunit) vegtyp
IF ( laiout /= 0 ) WRITE (flunit) lai
IF ( rfnsout /= 0 ) WRITE (flunit) roufns
IF ( vegout /= 0 ) WRITE (flunit) veg
IF ( ndviout /= 0 ) WRITE (flunit) ndvi
CLOSE ( flunit )
CALL retunit ( flunit )
ELSE IF (sfcdmp == 2) THEN
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block
!-----------------------------------------------------------------------
!
! Write out in HDF4.
!
!-----------------------------------------------------------------------
CALL hdfopen(trim(sfcoutfl), 2, sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "WRTSFCDT: ERROR creating HDF4 file: ", &
trim(sfcoutfl)
RETURN
END IF
! Make sure nstyp & stypfrct are valid for one soil type
IF (nstyp == 0) nstyp = 1
IF (nstyp == 1) stypfrct = 1
CALL hdfwrti(sd_id, 'nstyp', nstyp, stat)
CALL hdfwrti(sd_id, 'nx', nx, stat)
CALL hdfwrti(sd_id, 'ny', ny, stat)
CALL hdfwrtr(sd_id, 'dx', dx, stat)
CALL hdfwrtr(sd_id, 'dy', dy, stat)
CALL hdfwrti(sd_id, 'mapproj', mapproj, stat)
CALL hdfwrtr(sd_id, 'trulat1', trulat1, stat)
CALL hdfwrtr(sd_id, 'trulat2', trulat2, stat)
CALL hdfwrtr(sd_id, 'trulon', trulon, stat)
CALL hdfwrtr(sd_id, 'sclfct', sclfct, stat)
CALL hdfwrtr(sd_id, 'ctrlat', ctrlat, stat)
CALL hdfwrtr(sd_id, 'ctrlon', ctrlon, stat)
IF ( stypout /= 0 ) THEN
CALL hdfwrt3di(soiltyp,nx,ny,nstyp,sd_id,0,0, &
'soiltyp','Soil type','index')
CALL hdfwrt3d(stypfrct,nx,ny,nstyp,sd_id,0,0, &
'stypfrct','Soil type fractional coverage','fraction', &
itmp,atmp1,atmp2)
END IF
IF ( vtypout /= 0 ) CALL hdfwrt2di(vegtyp,nx,ny,sd_id,0,0, &
'vegtyp','Vegetation type','index')
IF ( laiout /= 0 ) CALL hdfwrt2d(lai,nx,ny,sd_id,0,0, &
'lai','Leaf Area Index','index',itmp)
IF ( rfnsout /= 0 ) CALL hdfwrt2d(roufns,nx,ny,sd_id,0,0, &
'roufns','Surface roughness','0-1',itmp)
IF ( vegout /= 0 ) CALL hdfwrt2d(veg,nx,ny,sd_id,0,0, &
'veg','Vegetation fraction','fraction',itmp)
!wdt update
IF ( ndviout /= 0 ) CALL hdfwrt2d(ndvi,nx,ny,sd_id,0,0, &
'ndvi','Normalized differential vegetation index','index',itmp)
CALL hdfclose(sd_id,stat)
IF (stat /= 0) THEN
WRITE (6,*) "WRTSFCDT: ERROR on closing file ",trim(sfcoutfl), &
" (status",stat,")"
END IF
!wdt end block
! alternate dump format ...
END IF
RETURN
END SUBROUTINE wrtsfcdt
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE READSFCDT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE readsfcdt( nx,ny,nstyps,sfcfile,dx,dy, & 3,33
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg,ndvi )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read the surface data sets from file sfcfile.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 5/3/94
!
! MODIFICATION HISTORY:
!
! 02/07/1995 (Yuhe Liu)
! Added a new 2-D array, veg(nx,ny), to the soil and vegetation data
! set file.
!
! 2000/03/29 (Gene Bassett)
! Removed the globcst.inc include.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2000/03/29 (Gene Bassett)
! Added HDF4 format.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of model grid points in the x-dir. (east/west)
! ny Number of model grid points in the y-dir. (north/south)
!
! mapproj Type of map projection used to setup the analysis grid.
! trulat1 1st real true latitude of map projection.
! trulat2 2nd real true latitude of map projection.
! trulon Real true longitude of map projection.
! sclfct Map scale factor. At latitude = trulat1 and trulat2
!
! dx Model grid spacing in the x-direction east-west (meters)
! dy Model grid spacing in the y-direction east-west (meters)
! ctrlat Lat. at the origin of the model grid (deg. N)
! ctrlon Lon. at the origin of the model grid (deg. E)
!
! OUTPUT:
!
! soiltyp Soil type in model domain
! vegtyp Vegetation type in model domain
! lai Leaf Area Index in model domain
! roufns Surface roughness
! veg Vegetation fraction
! ndvi NDVI
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx ! Number of grid points in the x-direction
INTEGER :: ny ! Number of grid points in the y-direction
INTEGER :: nstyps ! Max number of soil types in a grid box
CHARACTER (LEN=*) :: sfcfile ! Name of the surface data file
REAL :: dx
REAL :: dy
INTEGER :: mapproj ! Map projection
REAL :: trulat1 ! 1st real true latitude of map projection
REAL :: trulat2 ! 2nd real true latitude of map projection
REAL :: trulon ! Real true longitude of map projection
REAL :: sclfct ! Map scale factor
REAL :: ctrlat ! Center latitude of the model domain (deg. N)
REAL :: ctrlon ! Center longitude of the model domain (deg. E)
INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type in model domain
REAL :: stypfrct(nx,ny,nstyps) ! Fraction of soil types
INTEGER :: vegtyp (nx,ny) ! Vegetation type in model domain
REAL :: lai (nx,ny) ! Leaf Area Index in model domain
REAL :: roufns (nx,ny) ! NDVI in model domain
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: ndvi (nx,ny) ! NDVI
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: nxin,nyin
REAL :: dxin,dyin
INTEGER :: mprojin
INTEGER :: nstyp1, nstypin
REAL :: trlat1in, trlat2in, trlonin, sclfctin
REAL :: ctrlonin, ctrlatin
INTEGER :: stypin,vtypin,laiin,roufin,vegin,ndviin
INTEGER :: idummy
REAL :: rdummy
INTEGER :: istat, ierr,i,j,is
INTEGER :: ireturn
CHARACTER (LEN=128) :: savename
CHARACTER :: ayear*4, amm*2, aday*2
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
INTEGER(KIND=selected_int_kind(4)) :: itmp(1)
! unused array in hdf routines since NO COMPRESSION
REAL :: atmp1(1),atmp2(1)
! unused arrays in hdf routines since NO COMPRESSION
INTEGER :: stat, sd_id
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Open the surface data file. Read the parameters first, check if
! the data are consistant with the model. If everything is OK, then
! read the surface data, soiltyp, vegtyp, lai, roufns, veg, and
! ndvi.
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0) THEN
savename(1:128) = sfcfile(1:128)
WRITE(sfcfile, '(a,a,2i2.2)') trim(savename),'_',loc_x,loc_y
END IF
WRITE (6,*) "READSFCDT: reading in external supplied surface", &
"data from file ",trim(sfcfile)
!-----------------------------------------------------------------------
!
! Read in header information.
!
!-----------------------------------------------------------------------
IF (sfcfmt == 0) THEN
!-----------------------------------------------------------------------
!
! Fortran unformatted dump.
!
!-----------------------------------------------------------------------
CALL getunit( sfcunit )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(sfcfile, '-F f77 -N ieee', ierr)
OPEN(UNIT=sfcunit,FILE=trim(sfcfile),FORM='unformatted', &
STATUS='old',IOSTAT=istat)
IF( istat /= 0 ) THEN
WRITE(6,'(/1x,a,i2,/1x,a/)') &
'Error occured when opening the surface data file ' &
//sfcfile//' using FORTRAN unit ',sfcunit, &
' Program stopped in READSFCDT.'
CALL arpsstop("arpsstop called from READSFCDT opening file",1)
END IF
WRITE(6,'(/1x,a,/1x,a,i2/)') &
'This run will start from an external supplied surface ', &
'data file '//sfcfile//' using FORTRAN unit ',sfcunit
READ (sfcunit,ERR=998) nxin,nyin
READ (sfcunit,ERR=998) mprojin,stypin,vtypin,laiin,roufin, &
vegin, ndviin,nstyp1,idummy,idummy, &
idummy, idummy,idummy,idummy,idummy, &
idummy, idummy,idummy,idummy,idummy
READ (sfcunit,ERR=998) dxin,dyin, ctrlatin,ctrlonin,trlat1in, &
trlat2in,trlonin,sclfctin,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy
ELSE IF (sfcfmt == 1) THEN
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block
!-----------------------------------------------------------------------
!
! HDF4 format.
!
!-----------------------------------------------------------------------
CALL hdfopen(trim(sfcfile), 1, sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "READSFCDT: ERROR opening ", &
trim(sfcfile)," for reading."
GO TO 998
END IF
CALL hdfrdi(sd_id,"nstyp",nstyp1,istat)
CALL hdfrdi(sd_id,"nx",nxin,istat)
CALL hdfrdi(sd_id,"ny",nyin,istat)
CALL hdfrdr(sd_id,"dx",dxin,istat)
CALL hdfrdr(sd_id,"dy",dyin,istat)
CALL hdfrdi(sd_id,"mapproj",mprojin,istat)
CALL hdfrdr(sd_id,"trulat1",trlat1in,istat)
CALL hdfrdr(sd_id,"trulat2",trlat2in,istat)
CALL hdfrdr(sd_id,"trulon",trlonin,istat)
CALL hdfrdr(sd_id,"sclfct",sclfctin,istat)
CALL hdfrdr(sd_id,"ctrlat",ctrlatin,istat)
CALL hdfrdr(sd_id,"ctrlon",ctrlonin,istat)
!wdt end block
! alternate dump format ...
ELSE IF (sfcfmt == 2) THEN !Direct output (JAB)
CALL initztime(ayear,amm,aday)
CALL readjveg(nx,ny,nstyp,ayear,amm,aday,ztime, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg,ndvi)
END IF !sfcfmt loop
nstyp1 = MAX( nstyp1, 1 )
!-----------------------------------------------------------------------
!
! Check the data file for consistent grid parameters.
!
!-----------------------------------------------------------------------
IF (sfcfmt /= 2) THEN
CALL checkgrid2d(nx,ny,nxin,nyin, &
dx,dy,ctrlat,ctrlon, &
mapproj,trulat1,trulat2,trulon,sclfct, &
dxin,dyin,ctrlatin,ctrlonin, &
mprojin,trlat1in,trlat2in,trlonin,sclfctin,ireturn)
IF (ireturn /= 0) THEN
WRITE (6,*) "READSFCDT: ERROR, grid parameter mismatch"
CALL arpsstop("arpsstop called from READSFCDT grid param. mismatch",1)
END IF
WRITE (6,'(a,a//a,i1/6(a,e15.8/))') &
' The map projection and griding information for the ', &
' surface data: ', &
' Projection: ', mprojin, &
' The 1st real true latitude: ', trlat1in, &
' The 2nd real true latitude: ', trlat2in, &
' The real true longitude: ', trlonin, &
' Map scale factor: ', sclfctin, &
' Latitude at the origin: ', ctrlatin, &
' Longitude at the origin: ', ctrlonin
ENDIF
!-----------------------------------------------------------------------
!
! Read in the surface data from the surface data file.
!
!-----------------------------------------------------------------------
IF (sfcfmt == 0) THEN ! Fortran unformatted
WRITE (6, '(a/a,i2/a,i2/a,i2/a,i2/a,i2/a,i2)') &
' Surface data flags for: ', &
' soiltyp --', stypin, &
' vegtyp --', vtypin, &
' lai --', laiin, &
' roufns --', roufin, &
' veg --', vegin, &
' ndvi --', ndviin
WRITE (6, '(a/a,i2)') &
' Number of soil types in each grid box:', &
' nstyp --', nstyp
IF(stypin == 1) THEN
IF ( nstyp1 == 1 ) THEN
WRITE (6, '(a)') 'Read in the soil type data'
READ (sfcunit,ERR=998) ((soiltyp(i,j,1),i=1,nx),j=1,ny)
DO j=1,ny
DO i=1,nx
stypfrct(i,j,1) = 1.0
END DO
END DO
ELSE
DO is=1,nstyp1
IF (is <= nstyps) THEN
WRITE (6, '(a)') 'Read in the soil type data'
READ (sfcunit,ERR=998) ((soiltyp(i,j,is),i=1,nx),j=1,ny)
WRITE (6, '(a)') 'Read in the fraction of soil types'
READ (sfcunit,ERR=998) ((stypfrct(i,j,is),i=1,nx),j=1,ny)
ELSE
READ (sfcunit,ERR=998)
READ (sfcunit,ERR=998)
ENDIF
END DO
END IF
END IF
IF(vtypin == 1) THEN
WRITE (6, '(a)') 'Read in the vegetation type data'
READ (sfcunit,ERR=998) vegtyp
END IF
IF(laiin == 1) THEN
WRITE (6, '(a)') 'Read in the Leaf Area Index data'
READ (sfcunit,ERR=998) lai
END IF
IF(roufin == 1) THEN
WRITE (6, '(a)') 'Read in the surface roughness data'
READ (sfcunit,ERR=998) roufns
END IF
IF(vegin == 1) THEN
WRITE (6, '(a)') 'Read in the vegetatin fraction data'
READ (sfcunit,ERR=998) veg
END IF
IF (ndviin == 1) THEN
WRITE (6, '(a)') 'Read in the NDVI data'
READ (sfcunit,ERR=998) ndvi
END IF
CLOSE ( sfcunit )
CALL retunit ( sfcunit )
ELSE IF (sfcfmt == 1) THEN
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block
nstypin = MIN(nstyp1, nstyps)
CALL hdfrd3di(sd_id,"soiltyp",nx,ny,nstypin,soiltyp,stat)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read in soiltyp'
stypin = 1
ELSE
WRITE (6, '(a)') 'Variable soiltyp is not in the data set.'
stypin = 0
END IF
CALL hdfrd3d(sd_id,"stypfrct",nx,ny,nstypin, &
stypfrct,stat,itmp,atmp1,atmp2)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read in stypfrct'
ELSE
WRITE (6, '(a)') 'Variable stypfrct is not in the data set.'
stypfrct(:,:,1) = 1.
stypfrct(:,:,2:nstypin) = 0.
END IF
CALL hdfrd2di(sd_id,"vegtyp",nx,ny,vegtyp,stat)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read in vegtyp'
vtypin = 1
ELSE
WRITE (6, '(a)') 'Variable vegtyp is not in the data set.'
vtypin = 0
END IF
CALL hdfrd2d(sd_id,"lai",nx,ny,lai,stat,itmp)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read in lai'
laiin = 1
ELSE
WRITE (6, '(a)') 'Variable lai is not in the data set.'
laiin = 0
END IF
CALL hdfrd2d(sd_id,"roufns",nx,ny,roufns,stat,itmp)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read in roufns'
roufin = 1
ELSE
WRITE (6, '(a)') 'Variable roufns is not in the data set.'
roufin = 0
END IF
CALL hdfrd2d(sd_id,"veg",nx,ny,veg,stat,itmp)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read in veg'
vegin = 1
ELSE
WRITE (6, '(a)') 'Variable veg is not in the data set.'
vegin = 0
END IF
CALL hdfrd2d(sd_id,"ndvi",nx,ny,ndvi,stat,itmp)
IF (stat > 1) GO TO 998
IF (stat == 0) THEN
WRITE (6, '(a)') 'Read in ndvi'
ndviin = 1
ELSE
WRITE (6, '(a)') 'Variable ndvi is not in the data set.'
ndviin = 0
END IF
CALL hdfclose(sd_id, stat)
!wdt end block
! alternate dump format ...
ELSE IF (sfcfmt == 2) THEN !Direct output from OASIS data
CALL initztime(ayear,amm,aday)
CALL readjveg(nx,ny,nstyp,ayear,amm,aday,ztime, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg,ndvi)
END IF
CALL fix_stypfrct_nstyp(nx,ny,nstyp1,nstyp,stypfrct)
IF(stypin == 0) THEN
WRITE(6,'(/1x,a,a)') &
'Soil type data not present in the surface propery ', &
'data file, set it to styp.'
DO i=1,nx
DO j=1,ny
soiltyp(i,j,1) = styp
END DO
END DO
END IF
IF(vtypin == 0) THEN
WRITE(6,'(/1x,a,a)') &
'Vegetation type data not present in the surface propery ', &
'data file, set it to vtyp.'
DO i=1,nx-1
DO j=1,ny-1
vegtyp(i,j) = vtyp
END DO
END DO
END IF
IF(laiin == 0) THEN
WRITE(6,'(/1x,a,a)') &
'Leaf Area Index data not present in the surface propery ', &
'data file, set it to lai0.'
DO i=1,nx-1
DO j=1,ny-1
lai(i,j) = lai0
END DO
END DO
END IF
IF(roufin == 0) THEN
WRITE(6,'(/1x,a,a)') &
'Roughness data not present in the surface propery ', &
'data file, set it to roufns0.'
DO i=1,nx-1
DO j=1,ny-1
roufns(i,j) = roufns0
END DO
END DO
END IF
IF(vegin == 0) THEN
WRITE(6,'(/1x,a,a)') &
'Vegetation fraction not present in the surface propery ', &
'data file, set it to veg0.'
DO i=1,nx-1
DO j=1,ny-1
veg(i,j) = veg0
END DO
END DO
END IF
IF (mp_opt > 0) sfcfile(1:128) = savename(1:128)
RETURN
998 WRITE (6,'(/a,i2/a)') &
'READSFCDT: Read error in surface data file ' &
//sfcfile//' with the I/O unit ',sfcunit, &
'The model will STOP in subroutine READSFCDT.'
CALL arpsstop("arpsstop called from READSFCDT reading sfc file",1)
END SUBROUTINE readsfcdt
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WRITTRN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE writtrn(nx,ny,ternfn,dx,dy, & 5,17
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
hterain)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Write out a terrain data file to be used by ARPS.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 5/1/1995.
!
! MODIFICATION HISTORY:
!
! 2000/03/29 (Gene Bassett)
! Removed the globcst.inc include.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2000/03/29 (Gene Bassett)
! Added HDF4 format.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
!
!
! mapproj Type of map projection used to setup the analysis grid.
! trulat1 1st real true latitude of map projection.
! trulat2 2nd real true latitude of map projection.
! trulon Real true longitude of map projection.
! sclfct Map scale factor. At latitude = trulat1 and trulat2
!
! dx Model grid spacing in the x-direction east-west (meters)
! dy Model grid spacing in the y-direction east-west (meters)
! ctrlat Lat. at the origin of the model grid (deg. N)
! ctrlon Lon. at the origin of the model grid (deg. E)
!
! ternfn Terrain data file name
! hterain Terrain height (m)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny
CHARACTER (LEN=*) :: ternfn
REAL :: dx
REAL :: dy
INTEGER :: mapproj ! Map projection
REAL :: trulat1 ! 1st real true latitude of map projection
REAL :: trulat2 ! 2nd real true latitude of map projection
REAL :: trulon ! Real true longitude of map projection
REAL :: sclfct ! Map scale factor
REAL :: ctrlat ! Center latitude of the model domain (deg. N)
REAL :: ctrlon ! Center longitude of the model domain (deg. E)
REAL :: hterain(nx,ny) ! The height of terrain.
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
INTEGER :: idummy, ierr, nunit
REAL :: rdummy
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
INTEGER(KIND=selected_int_kind(4)) :: itmp(1)
! unused array in hdf routines since NO COMPRESSION
REAL :: atmp1(1),atmp2(1)
! unused arrays in hdf routines since NO COMPRESSION
INTEGER :: stat, sd_id
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF (terndmp == 0) RETURN
WRITE (6,*) 'WRITTRN: Opening Terrain file ',ternfn
!-----------------------------------------------------------------------
!
! Write out in Fortran unformatted.
!
!-----------------------------------------------------------------------
IF (terndmp == 1) THEN
CALL getunit( nunit )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(ternfn, '-F f77 -N ieee', ierr)
OPEN(nunit,FILE=trim(ternfn),FORM='unformatted', &
STATUS='unknown')
WRITE(nunit) nx,ny
idummy = 0
rdummy = 0.0
WRITE(nunit) idummy,mapproj,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
WRITE(nunit) dx ,dy ,ctrlat,ctrlon,rdummy , &
rdummy,trulat1,trulat2,trulon,sclfct, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy
WRITE(nunit) hterain
CLOSE(nunit)
ELSE IF (terndmp == 2) THEN
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block
!-----------------------------------------------------------------------
!
! Write out in HDF4.
!
!-----------------------------------------------------------------------
CALL hdfopen(trim(ternfn), 2, sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "WRTTRN: ERROR creating HDF4 file: ", &
trim(ternfn)
RETURN
END IF
CALL hdfwrti(sd_id, 'nx', nx, stat)
CALL hdfwrti(sd_id, 'ny', ny, stat)
CALL hdfwrtr(sd_id, 'dx', dx, stat)
CALL hdfwrtr(sd_id, 'dy', dy, stat)
CALL hdfwrti(sd_id, 'mapproj', mapproj, stat)
CALL hdfwrtr(sd_id, 'trulat1', trulat1, stat)
CALL hdfwrtr(sd_id, 'trulat2', trulat2, stat)
CALL hdfwrtr(sd_id, 'trulon', trulon, stat)
CALL hdfwrtr(sd_id, 'sclfct', sclfct, stat)
CALL hdfwrtr(sd_id, 'ctrlat', ctrlat, stat)
CALL hdfwrtr(sd_id, 'ctrlon', ctrlon, stat)
CALL hdfwrt2d(hterain,nx,ny,sd_id,0,0, &
'hterain','Terrain Height','m',itmp)
! 200 CONTINUE
CALL hdfclose(sd_id,stat)
IF (stat /= 0) THEN
WRITE (6,*) "WRITTRN: ERROR on closing file ",trim(ternfn), &
" (status",stat,")"
END IF
!wdt end block
! alternate dump format ...
END IF
RETURN
END SUBROUTINE writtrn
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE READTRN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE readtrn(nx,ny, dx,dy, ternfile, & 2,22
mapproj,trulat1,trulat2,trulon,sclfct, &
ctrlat,ctrlon,hterain )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read the terrain data into model array hterain from a specified
! terrain data file.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 2/27/1994.
!
! MODIFICATION HISTORY:
!
! 9/10/94 (Weygandt & Y. Lu)
! Cleaned up documentation.
!
! 2000/03/29 (Gene Bassett)
! Removed the globcst.inc include.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2000/03/29 (Gene Bassett)
! Added HDF4 format.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
!
! dx Grid interval in x-direction
! dy Grid interval in y-direction
!
! ternfile Terrain data file name
!
! mapproj Type of map projection used to setup the analysis grid.
! trulat1 1st real true latitude of map projection.
! trulat2 2nd real true latitude of map projection.
! trulon Real true longitude of map projection.
! sclfct Map scale factor. At latitude = trulat1 and trulat2
!
! ctrlat Lat. at the origin of the model grid (deg. N)
! ctrlon Lon. at the origin of the model grid (deg. E)
!
! OUTPUT:
!
! hterain Terrain height (m)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx ! Number of grid points in the x-direction
INTEGER :: ny ! Number of grid points in the y-direction
REAL :: dx ! Grid interval in x-direction
REAL :: dy ! Grid interval in y-direction
CHARACTER (LEN=*) :: ternfile ! Terrain data file name
INTEGER :: mapproj ! Map projection
REAL :: trulat1 ! 1st real true latitude of map projection
REAL :: trulat2 ! 2nd real true latitude of map projection
REAL :: trulon ! Real true longitude of map projection
REAL :: sclfct ! Map scale factor
REAL :: ctrlat ! Center latitude of the model domain (deg. N)
REAL :: ctrlon ! Center longitude of the model domain (deg. E)
REAL :: hterain(nx,ny) ! Terrain height.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: inunit,istat
INTEGER :: nxin,nyin,idummy,ierr
REAL :: dxin,dyin,rdummy,amin,amax
INTEGER :: ireturn
REAL :: ctrlatin,ctrlonin, &
trulat1in,trulat2in,trulonin,sclfctin
INTEGER :: mapprojin
CHARACTER (LEN=128) :: savename
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
INTEGER(KIND=selected_int_kind(4)) :: itmp(1)
! unused array in hdf routines since NO COMPRESSION
REAL :: atmp1(1),atmp2(1)
! unused arrays in hdf routines since NO COMPRESSION
INTEGER :: stat, sd_id
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Read in the terrain data.
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0) THEN
savename(1:128) = ternfile(1:128)
WRITE(ternfile, '(a,a,2i2.2)') trim(savename),'_',loc_x,loc_y
END IF
WRITE (6,*) "READTRN: reading in external supplied terrain ", &
"data from file ",trim(ternfile)
!-----------------------------------------------------------------------
!
! Read in header information.
!
!-----------------------------------------------------------------------
IF (ternfmt == 0) THEN
CALL getunit( inunit )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(ternfile, '-F f77 -N ieee', ierr)
OPEN(UNIT=inunit,FILE=trim(ternfile),FORM='unformatted', &
STATUS='old',IOSTAT=istat)
IF( istat /= 0) THEN
WRITE(6,'(/1x,a,a,/1x,a/)') &
'Error occured when opening terrain data file ', &
ternfile,' Job stopped in READTRN.'
CALL arpsstop("arpsstop called from READTRN reading terrain file",1)
END IF
READ(inunit,ERR=999) nxin,nyin
READ(inunit,ERR=999) idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
READ(inunit,ERR=999) dxin ,dyin ,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy
ELSE
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block
CALL hdfopen(trim(ternfile), 1, sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "READTRN: ERROR opening ", &
trim(ternfile)," for reading."
GO TO 999
END IF
CALL hdfrdi(sd_id,"nx",nxin,istat)
CALL hdfrdi(sd_id,"ny",nyin,istat)
CALL hdfrdr(sd_id,"dx",dxin,istat)
CALL hdfrdr(sd_id,"dy",dyin,istat)
CALL hdfrdi(sd_id,"mapproj",mapprojin,istat)
CALL hdfrdr(sd_id,"trulat1",trulat1in,istat)
CALL hdfrdr(sd_id,"trulat2",trulat2in,istat)
CALL hdfrdr(sd_id,"trulon",trulonin,istat)
CALL hdfrdr(sd_id,"sclfct",sclfctin,istat)
CALL hdfrdr(sd_id,"ctrlat",ctrlatin,istat)
CALL hdfrdr(sd_id,"ctrlon",ctrlonin,istat)
!wdt end block
! alternate dump format ...
END IF
IF (ternfmt == 0) THEN
WRITE (6,*) &
"READTRN: WARNING, not checking map projection parameters"
CALL checkgrid2d(nx,ny,nxin,nyin, &
dx,dy,ctrlat,ctrlon, &
mapproj,trulat1,trulat2,trulon,sclfct, &
dxin,dyin,ctrlat,ctrlon, &
mapproj,trulat1,trulat2,trulon,sclfct,ireturn)
ELSE
CALL checkgrid2d(nx,ny,nxin,nyin, &
dx,dy,ctrlat,ctrlon, &
mapproj,trulat1,trulat2,trulon,sclfct, &
dxin,dyin,ctrlatin,ctrlonin, &
mapprojin,trulat1in,trulat2in,trulonin,sclfctin,ireturn)
END IF
IF (ireturn /= 0) THEN
WRITE (6,*) "READTRN: ERROR, grid parameter mismatch"
CALL arpsstop("arpsstop called from READTRN parameter mismatch",1)
END IF
IF (ternfmt == 0) THEN
READ(inunit,ERR=999) hterain
CALL retunit( inunit )
CLOSE (UNIT=inunit)
ELSE
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
CALL hdfrd2d(sd_id,"hterain",nx,ny,hterain,stat,itmp)
END IF
WRITE(6,'(1x,a/)') 'Minimum and maximum terrain height:'
CALL a3dmax0(hterain,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'htrnmin = ', amin,', htrnmax=',amax
IF (mp_opt > 0) ternfile(1:128) = savename(1:128)
RETURN
999 WRITE(6,'(1x,a)') &
'Error in reading terrain data. Job stopped in READTRN.'
CALL arpsstop("arpsstop called from READTRN reading file",1)
END SUBROUTINE readtrn
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE WRITESND ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE writesnd(nx,ny,nz, & 1,2
ubar,vbar,ptbar,pbar,qvbar, &
zp, rhobar, zs)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Print the base state sounding profile at the center of the domain.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 3/17/1991.
!
! MODIFICATION HISTORY:
!
! 5/1/98 (G. Bassett, D. Weber)
! Moved from inibase.
!
!-----------------------------------------------------------------------
!
! 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
!
! OUTPUT:
!
! 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)
! qvbar Base state water vapor specific humidity (kg/kg)
!
! zp Vertical coordinate of grid points in physical space(m)
!
! rhobar Base state air density (kg/m**3)
! zs Physical coordinate height at scalar point (m)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number of grid points in 3
! directions
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 :: qvbar (nx,ny,nz) ! Base state water specific humidity
! (kg/kg).
REAL :: zp (nx,ny,nz) ! Physical height coordinate defined at
! w-point of staggered grid.
REAL :: rhobar(nx,ny,nz) ! Base state air density (kg/m**3)
REAL :: zs (nx,ny,nz) ! Physical coordinate height at scalar
! point (m)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k, istat, nunit, itrnmin, jtrnmin
REAL :: ternmin
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF(myproc == 0) THEN
ternmin = zp(1,1,2)
itrnmin = 1
jtrnmin = 1
DO j=1,ny-1
DO i=1,nx-1
IF(zp(i,j,2) < ternmin) THEN
ternmin = zp(i,j,2)
itrnmin = i
jtrnmin = j
END IF
END DO
END DO
i = itrnmin
j = jtrnmin
WRITE(6,'(/a,i4,a,i4/)') &
' Base state z-profiles at I=',i,', J=',j
CALL getunit( nunit )
OPEN(UNIT=nunit, FILE=runname(1:lfnkey)//'.sound', &
STATUS='unknown',IOSTAT=istat)
IF(istat /= 0) THEN
WRITE (6,'(1x,a,a)') 'Error opening file ', &
runname(1:lfnkey)//'.sound'
CALL arpsstop("arpsstop called from WRITESND opening file",1)
END IF
WRITE(6,'(1x,2a)') &
' k z(m) pbar(Pascal) ptbar(K) rhobar(kg/m**3)', &
' Qvbar(kg/kg) U(m/s) V(m/s)'
WRITE(nunit,'(1x,2a)') &
' k z(m) pbar(Pascal) ptbar(K) rhobar(kg/m**3)', &
' Qvbar(kg/kg) U(m/s) V(m/s)'
DO k=nz-1,1,-1
WRITE(6, '(1x,i4,2f12.3,f12.5,2f12.8,2f12.5)') &
k,zs(i,j,k),pbar(i,j,k),ptbar(i,j,k),rhobar(i,j,k) &
,qvbar(i,j,k),ubar(i,j,k), vbar(i,j,k)
WRITE(nunit,'(1x,i4,2f12.3,f12.5,2f12.8,2f12.5)') &
k,zs(i,j,k),pbar(i,j,k),ptbar(i,j,k),rhobar(i,j,k) &
,qvbar(i,j,k),ubar(i,j,k), vbar(i,j,k)
END DO
CLOSE( UNIT = nunit )
CALL retunit( nunit )
END IF !Message Passing
RETURN
END SUBROUTINE writesnd
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SFCCNTL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE sfccntl(nx,ny, sfcoutfl, & 3,7
stypout,vtypout,laiout,rfnsout,vegout,ndviout, &
x,y, temxy1,temxy2)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! ARPS surface property data description file for GrADS display
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 10/05/1998
!
! 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)
!
! x The x-coord. of the physical and computational grid.
! Defined at u-point.
! y The y-coord. of the physical and computational grid.
! Defined at v-point.
!
! WORK ARRAY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny ! Number of grid points in 3 directions
CHARACTER (LEN=* ) :: sfcoutfl ! Surface data file name
INTEGER :: stypout ! Flag for output of soiltyp
INTEGER :: vtypout ! Flag for output of vegtyp
INTEGER :: laiout ! Flag for output of lai
INTEGER :: rfnsout ! Flag for output of roufns
INTEGER :: vegout ! Flag for output of veg
INTEGER :: ndviout ! Flag for output of ndvi
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 :: temxy1(nx,ny) ! Temporary array
REAL :: temxy2(nx,ny) ! Temporary array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: varnumax
PARAMETER ( varnumax = 100 )
CHARACTER (LEN=256) :: sfcctlfl
CHARACTER (LEN=256) :: sfcdatfl
CHARACTER (LEN=15) :: chrstr, chr1
CHARACTER (LEN=8) :: varnam(varnumax)
CHARACTER (LEN=60) :: vartit(varnumax)
CHARACTER (LEN=10) :: varparam(varnumax)
INTEGER :: varlev(varnumax)
INTEGER :: varnum
INTEGER :: nchout0
CHARACTER (LEN=3) :: monnam(12) ! Name of months
DATA monnam/'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', &
'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/
CHARACTER (LEN=2) :: dtunit
INTEGER :: ntm,tinc
REAL :: lonmin,latmin,lonmax,latmax
REAL :: xbgn,ybgn
REAL :: xinc,yinc
REAL :: lat11,lat12,lon11,lon12,lat21,lat22,lon21,lon22
REAL :: latinc,loninc
INTEGER :: ctllen,fnlen
INTEGER :: i,k, is
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Open the GrADS data control file: runname.ctl
!
!-----------------------------------------------------------------------
!
fnlen = LEN( sfcoutfl )
CALL strlnth( sfcoutfl, fnlen )
ctllen = fnlen + 4
sfcctlfl(1:ctllen) = sfcoutfl(1:fnlen)//'.ctl'
CALL fnversn( sfcctlfl, ctllen )
CALL getunit (nchout0)
WRITE (6,'(a)') 'The GrADS data control file is ' &
//sfcctlfl(1:ctllen)
OPEN (nchout0, FILE = sfcctlfl(1:ctllen), STATUS = 'unknown')
xbgn = 0.5 * (x(1) + x(2))
ybgn = 0.5 * (y(1) + y(2))
xinc = (x(2) - x(1))
yinc = (y(2) - y(1))
CALL xytoll(nx,ny,x,y,temxy1,temxy2)
CALL xytoll(1,1,xbgn,ybgn,lat11,lon11)
CALL a3dmax0(temxy1,1,nx,1,nx,1,ny,1,ny-1,1,1,1,1, &
latmax,latmin)
CALL a3dmax0(temxy2,1,nx,1,nx,1,ny,1,ny-1,1,1,1,1, &
lonmax,lonmin)
latinc = (latmax-latmin)/(ny-1)
loninc = (lonmax-lonmin)/(nx-1)
WRITE (6,'(a,f10.4,a,f10.4,a,f10.4)') &
'latmin:latmax:latinc = ', &
latmin,':',latmax,':',latinc
WRITE (6,'(a,f10.4,a,f10.4,a,f10.4)') &
'lonmin:lonmax:loninc = ', &
lonmin,':',lonmax,':',loninc
WRITE (chrstr,'(i2.2,a,i2.2,a,i2.2,a3,i4.4)') &
hour,':',minute,'Z',day,monnam(month),year
ntm = 1
tinc = 1
dtunit = 'MO'
DO i=1,varnumax
varlev(i) = 0
varparam(i) = '99'
END DO
varnum = 0
IF ( stypout == 1 ) THEN
IF ( nstyp == 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'styp'
vartit(varnum) = 'Soil type'
varparam(varnum) = '-1,40,4'
ELSE
DO is=1,nstyp
WRITE (chr1,'(i1)') is
varnum = varnum + 1
varnam(varnum) = 'styp'//chr1
vartit(varnum) = 'Soil type '//chr1
varparam(varnum) = '-1,40,4'
varnum = varnum + 1
varnam(varnum) = 'sfct'//chr1
vartit(varnum) = 'Fraction of soil type '//chr1
END DO
END IF
END IF
IF ( vtypout == 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'vtyp'
vartit(varnum) = 'Vegetation type'
varparam(varnum) = '-1,40,4'
END IF
IF ( laiout == 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'lai'
vartit(varnum) = 'Leaf Area Index'
END IF
IF ( rfnsout == 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'roufns'
vartit(varnum) = 'Surface roughness (m)'
END IF
IF ( vegout == 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'veg'
vartit(varnum) = 'Vegetation fraction'
END IF
IF ( ndviout == 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'ndvi'
vartit(varnum) = 'NDVI'
END IF
WRITE (nchout0,'(a/a)') &
'TITLE ARPS surface property data control for ' &
//runname(1:lfnkey),'*'
WRITE (nchout0,'(a,a)') &
'DSET ',sfcoutfl(1:fnlen)
WRITE (nchout0,'(a)') &
'OPTIONS sequential'
! WRITE (nchout0,'(a,i10)') &
! 'FILEHEADER ',192 !!! The number depends on the file
! !!! structure of sfcoutfl. See iolib3d.f
WRITE (nchout0,'(a,i10)') &
'FILEHEADER ',192 !!! The number depends on the file
!!! structure of sfcoutfl. See iolib3d.f
WRITE (nchout0,'(a/a)') &
'UNDEF -9.e+33','*'
IF ( mapproj == 2 ) THEN
WRITE (nchout0,'(a)') &
'* For lat-lon-lev display, umcomment the following 4 lines.'
WRITE (nchout0,'(a,1x,i8,1x,i3,a,2f12.6,2i3,3f12.6,2f12.2)') &
'PDEF',nx,ny,' LCC',lat11,lon11,1,1, &
trulat1,trulat2,trulon,xinc,yinc
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'XDEF',nx, ' LINEAR ',lonmin,loninc
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'YDEF',ny, ' LINEAR ',latmin,latinc
ELSE
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'XDEF',nx, ' LINEAR ',xbgn,xinc
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'YDEF',ny, ' LINEAR ',ybgn,yinc
END IF
WRITE (nchout0,'(a,1x,i8,a,i1)') &
'ZDEF',1,' LEVELS ',0
WRITE (nchout0,'(a,1x,i8,a,a,3x,i2.2,a/a)') &
'TDEF',ntm,' LINEAR ',chrstr,tinc,dtunit,'*'
WRITE (nchout0,'(a,1x,i3)') &
'VARS',varnum
DO i = 1, varnum
WRITE (nchout0,'(a8,1x,i3,2x,a,2x,a)') &
varnam(i),varlev(i),varparam(i),vartit(i)
END DO
WRITE (nchout0,'(a)') &
'ENDVARS'
CLOSE (nchout0)
CALL retunit(nchout0)
RETURN
END SUBROUTINE sfccntl
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SOILCNTL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE soilcntl(nx,ny,nzsoil,zpsoil,soiloutfl, & 5,7
zpsoilout,tsoilout,qsoilout,wcanpout,snowdout,x,y)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! ARPS soil data description file for GrADS display
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 10/05/1998
!
! 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)
! nzsoil Number of grid points in the soil
!
! x The x-coord. of the physical and computational grid.
! Defined at u-point.
! y The y-coord. of the physical and computational grid.
! Defined at v-point.
!
! WORK ARRAY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny ! Number of grid points in 3 directions
INTEGER :: nzsoil ! Number of grid points in the soil
CHARACTER (LEN=*) :: soiloutfl ! Surface data file name
INTEGER :: zpsoilout,tsoilout,qsoilout,wcanpout,snowdout
INTEGER :: stypout
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 :: zpsoil(nx,ny,nzsoil) ! Depth of soil layers.
REAL :: temxy1(nx,ny) ! Temporary array
REAL :: temxy2(nx,ny) ! Temporary array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: varnumax
PARAMETER ( varnumax = 100 )
CHARACTER (LEN=256) :: soilctlfl
CHARACTER (LEN=15) :: timestr,chrstr, chr1
CHARACTER (LEN=8) :: varnam(varnumax)
CHARACTER (LEN=60) :: vartit(varnumax)
CHARACTER (LEN=10) :: varparam(varnumax)
INTEGER :: varlev(varnumax)
INTEGER :: varnum
INTEGER :: nchout0
CHARACTER (LEN=3) :: monnam(12) ! Name of months
DATA monnam/'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', &
'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/
CHARACTER (LEN=2) :: dtunit
INTEGER :: ntm,tinc
REAL :: lonmin,latmin,lonmax,latmax
REAL :: xbgn,ybgn
REAL :: xinc,yinc
REAL :: lat11,lat12,lon11,lon12,lat21,lat22,lon21,lon22
REAL :: latinc,loninc
INTEGER :: ctllen, fnlen
INTEGER :: i,j,k, is
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Open the GrADS data control file
!
!-----------------------------------------------------------------------
!
fnlen = LEN( soiloutfl )
CALL strlnth( soiloutfl, fnlen )
ctllen = fnlen + 4
soilctlfl(1:ctllen) = soiloutfl(1:fnlen)//'.ctl'
CALL fnversn( soilctlfl, ctllen )
CALL getunit (nchout0)
WRITE (6,'(a)') 'The GrADS data control file is ' &
//soilctlfl(1:ctllen)
OPEN (nchout0, FILE = soilctlfl(1:ctllen), STATUS = 'unknown')
xbgn = 0.5 * (x(1) + x(2))
ybgn = 0.5 * (y(1) + y(2))
xinc = (x(2) - x(1))
yinc = (y(2) - y(1))
CALL xytoll(nx,ny,x,y,temxy1,temxy2)
CALL xytoll(1,1,xbgn,ybgn,lat11,lon11)
CALL a3dmax0(temxy1,1,nx,1,nx,1,ny,1,ny-1,1,1,1,1, &
latmax,latmin)
CALL a3dmax0(temxy2,1,nx,1,nx,1,ny,1,ny-1,1,1,1,1, &
lonmax,lonmin)
latinc = (latmax-latmin)/(ny-1)
loninc = (lonmax-lonmin)/(nx-1)
WRITE (6,'(a,f10.4,a,f10.4,a,f10.4)') &
'latmin:latmax:latinc = ', &
latmin,':',latmax,':',latinc
WRITE (6,'(a,f10.4,a,f10.4,a,f10.4)') &
'lonmin:lonmax:loninc = ', &
lonmin,':',lonmax,':',loninc
WRITE (timestr,'(i2.2,a,i2.2,a,i2.2,a3,i4.4)') &
hour,':',minute,'Z',day,monnam(month),year
ntm = 1
tinc = 1
dtunit = 'MN'
DO i=1,varnumax
varlev(i) = 0
varparam(i) = '99'
END DO
varnum = 0
!
! 06/19/2002 Zuwen He
!
! Code changed due to new ou soil model.
!
! IF ( tsoilout+qsoilout+wcanpout /= 0 ) THEN
! stypout = 1
! ELSE
! stypout = 0
! END IF
IF ( zpsoilout == 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'zpsoil'
vartit(varnum) = 'Soil layer depth (m) '
varlev(varnum) = nzsoil
END IF
IF ( tsoilout == 1 ) THEN
DO is=0,nstyp
WRITE (chrstr,'(i2.2)') is
varnum = varnum + 1
varnam(varnum) = 'tsoil'//chrstr(1:2)
vartit(varnum) = 'Soil temperature (K) for soil type '//chrstr(1:2)
if (is == 0) vartit(varnum) = 'Soil temperature (K)'
varlev(varnum) = nzsoil
END DO
END IF
IF ( qsoilout == 1 ) THEN
DO is=0,nstyp
WRITE (chrstr,'(i2.2)') is
varnum = varnum + 1
varnam(varnum) = 'qsoil'//chrstr(1:2)
vartit(varnum) = 'Soil moisture (m**3/m**3) for soil type '//chrstr(1:2)
if (is == 0) vartit(varnum) = 'Soil moisture (m**3/m**3)'
varlev(varnum) = nzsoil
END DO
END IF
IF ( wcanpout == 1 ) THEN
DO is=0,nstyp
WRITE (chrstr,'(i2.2)') is
varnum = varnum + 1
varnam(varnum) = 'wr'//chrstr(1:2)
vartit(varnum) = 'Canopy water amount (m) for soil type '//chrstr(1:2)
if (is == 0) vartit(varnum) = 'Canopy water amount (m)'
varlev(varnum) = 0
END DO
END IF
IF ( snowdout == 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'snowd'
vartit(varnum) = 'Snow depth (m)'
varlev(varnum) = 0
END IF
! IF ( stypout == 1 ) THEN
! DO is=1,nstyp
! WRITE (chrstr,'(i2.2)') is
! varnum = varnum + 1
! varnam(varnum) = 'styp'//chrstr(1:2)
! vartit(varnum) = 'Soil type '//chrstr(1:2)
! varlev(varnum) = 0
! END DO
! END IF
WRITE (nchout0,'(a/a)') &
'TITLE ARPS soil variable data control for ' &
//runname(1:lfnkey),'*'
WRITE (nchout0,'(a,a)') &
'DSET ',soiloutfl(1:fnlen)
! PRINT *,'TEST #1 = fnlen = ',fnlen
WRITE (nchout0,'(a)') &
'OPTIONS sequential'
WRITE (nchout0,'(a,i10)') &
'FILEHEADER ',236 !!! The number depends on the file
!!! structure of soiloutfl. See iolib3d.f
WRITE (nchout0,'(a/a)') &
'UNDEF -9.e+33','*'
IF ( mapproj == 2 ) THEN
WRITE (nchout0,'(a)') &
'* For lat-lon-lev display, umcomment the following 4 lines.'
WRITE (nchout0,'(a,1x,i8,1x,i3,a,2f12.6,2i3,3f12.6,2f12.2)') &
'PDEF',nx,ny,' LCC',lat11,lon11,1,1, &
trulat1,trulat2,trulon,xinc,yinc
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'XDEF',nx, ' LINEAR ',lonmin,loninc
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'YDEF',ny, ' LINEAR ',latmin,latinc
ELSE
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'XDEF',nx, ' LINEAR ',xbgn,xinc
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'YDEF',ny, ' LINEAR ',ybgn,yinc
END IF
WRITE (nchout0,'(a,1x,i8,a)') &
'ZDEF',nzsoil,' LEVELS '
!
! 06/19/2002 Zuwen He
!
! The soil model is upside down.
!
WRITE (nchout0,'(8f10.2)') &
(zpsoil(1,1,k)-zpsoil(1,1,1),k=1,nzsoil)
WRITE (nchout0,'(a,1x,i8,a,a,3x,i2.2,a/a)') &
'TDEF',ntm,' LINEAR ',timestr,tinc,dtunit,'*'
WRITE (nchout0,'(a,1x,i3)') &
'VARS',varnum
DO i = 1, varnum
WRITE (nchout0,'(a8,1x,i3,2x,a,2x,a)') &
varnam(i),varlev(i),varparam(i),vartit(i)
END DO
WRITE (nchout0,'(a)') &
'ENDVARS'
CLOSE (nchout0)
CALL retunit(nchout0)
RETURN
END SUBROUTINE soilcntl
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TRNCNTL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE trncntl(nx,ny, ternfn, x,y) 1,7
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! ARPS terrain data description file for GrADS display
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 10/05/1998
!
! 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)
!
! x The x-coord. of the physical and computational grid.
! Defined at u-point.
! y The y-coord. of the physical and computational grid.
! Defined at v-point.
!
! WORK ARRAY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny ! Number of grid points in 3 directions
CHARACTER (LEN=* ) :: ternfn ! Terrain data file name
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 :: temxy1(nx,ny) ! Temporary array
REAL :: temxy2(nx,ny) ! Temporary array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=256) :: trnctlfl
CHARACTER (LEN=15) :: chrstr, chr1
INTEGER :: nchout0
CHARACTER (LEN=3) :: monnam(12) ! Name of months
DATA monnam/'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', &
'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/
CHARACTER (LEN=2) :: dtunit
INTEGER :: ntm,tinc
REAL :: lonmin,latmin,lonmax,latmax
REAL :: xbgn,ybgn
REAL :: xinc,yinc
REAL :: lat11,lat12,lon11,lon12,lat21,lat22,lon21,lon22
REAL :: latinc,loninc
INTEGER :: ctllen, fnlen
INTEGER :: i,k, is
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Open the GrADS data control file
!
!-----------------------------------------------------------------------
!
fnlen = LEN( ternfn )
CALL strlnth( ternfn, fnlen )
ctllen = fnlen + 4
trnctlfl(1:ctllen) = ternfn(1:fnlen)//'.ctl'
CALL fnversn( trnctlfl, ctllen )
CALL getunit (nchout0)
WRITE (6,'(a)') 'The GrADS data control file is ' &
//trnctlfl(1:ctllen)
OPEN (nchout0, FILE = trnctlfl(1:ctllen), STATUS = 'unknown')
xbgn = 0.5 * (x(1) + x(2))
ybgn = 0.5 * (y(1) + y(2))
xinc = (x(2) - x(1))
yinc = (y(2) - y(1))
CALL xytoll(nx,ny,x,y,temxy1,temxy2)
CALL xytoll(1,1,xbgn,ybgn,lat11,lon11)
CALL a3dmax0(temxy1,1,nx,1,nx,1,ny,1,ny-1,1,1,1,1, &
latmax,latmin)
CALL a3dmax0(temxy2,1,nx,1,nx,1,ny,1,ny-1,1,1,1,1, &
lonmax,lonmin)
latinc = (latmax-latmin)/(ny-1)
loninc = (lonmax-lonmin)/(nx-1)
WRITE (6,'(a,f10.4,a,f10.4,a,f10.4)') &
'latmin:latmax:latinc = ', &
latmin,':',latmax,':',latinc
WRITE (6,'(a,f10.4,a,f10.4,a,f10.4)') &
'lonmin:lonmax:loninc = ', &
lonmin,':',lonmax,':',loninc
WRITE (chrstr,'(i2.2,a,i2.2,a,i2.2,a3,i4.4)') &
hour,':',minute,'Z',day,monnam(month),year
ntm = 1
tinc = 1
dtunit = 'MN'
WRITE (nchout0,'(a/a)') &
'TITLE ARPS terrain data control for ' &
//runname(1:lfnkey),'*'
WRITE (nchout0,'(a,a)') &
'DSET ',ternfn(1:fnlen)
WRITE (nchout0,'(a)') &
'OPTIONS sequential'
WRITE (nchout0,'(a,i10)') &
'FILEHEADER ',192 !!! The number depends on the file
!!! structure of ternfn. See iolib3d.f
WRITE (nchout0,'(a/a)') &
'UNDEF -9.e+33','*'
IF ( mapproj == 2 ) THEN
WRITE (nchout0,'(a)') &
'* For lat-lon-lev display, umcomment the following 4 lines.'
WRITE (nchout0,'(a,1x,i8,1x,i3,a,2f12.6,2i3,3f12.6,2f12.2)') &
'PDEF',nx,ny,' LCC',lat11,lon11,1,1, &
trulat1,trulat2,trulon,xinc,yinc
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'XDEF',nx, ' LINEAR ',lonmin,loninc
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'YDEF',ny, ' LINEAR ',latmin,latinc
ELSE
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'XDEF',nx, ' LINEAR ',xbgn,xinc
WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)') &
'YDEF',ny, ' LINEAR ',ybgn,yinc
END IF
WRITE (nchout0,'(a,1x,i8,a,i1)') &
'ZDEF',1,' LEVELS ',0
WRITE (nchout0,'(a,1x,i8,a,a,3x,i2.2,a/a)') &
'TDEF',ntm,' LINEAR ',chrstr,tinc,dtunit,'*'
WRITE (nchout0,'(a,1x,i3)') &
'VARS',1
WRITE (nchout0,'(a)') &
'trn 0 99 ARPS terrain (m)'
WRITE (nchout0,'(a)') &
'ENDVARS'
CLOSE (nchout0)
CALL retunit(nchout0)
RETURN
END SUBROUTINE trncntl
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CHECKGRID3D ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE checkgrid3d(nx,ny,nz,nxin,nyin,nzin, & 2
dx,dy,dz,dzmin,ctrlat,ctrlon, &
strhopt,zrefsfc,dlayer1,dlayer2,zflat,strhtune, &
mapproj,trulat1,trulat2,trulon,sclfct, &
dxin,dyin,dzin,dzminin,ctrlatin,ctrlonin, &
strhoptin,zrefsfcin,dlayer1in,dlayer2in,zflatin,strhtunein, &
mapprojin,trulat1in,trulat2in,trulonin,sclfctin,ireturn)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Check grid information to see if input data is consistent with
! ARPS grid.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/03/24
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx
! ny
! nz
! nxin
! nyin
! nzin
! dx
! dy
! dz
! dzmin
! ctrlat
! ctrlon
! strhopt
! zrefsfc
! dlayer1
! dlayer2
! zflat
! strhtune
! mapproj
! trulat1
! trulat2
! trulon
! sclfct
! dxin
! dyin
! dzin
! dzminin
! ctrlatin
! ctrlonin
! strhoptin
! zrefsfcin
! dlayer1in
! dlayer2in
! zflatin
! strhtunein
! mapprojin
! trulat1in
! trulat2in
! trulonin
! sclfctin
!
! OUTPUT:
!
! ireturn Flag indicating if the grids are the same (0-okay,
! 1-significant differences)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nxin
INTEGER :: nyin
INTEGER :: nzin
INTEGER :: nx,ny,nz
REAL :: dx
REAL :: dy
REAL :: dz
REAL :: dzmin
REAL :: ctrlat
REAL :: ctrlon
INTEGER :: strhopt
REAL :: zrefsfc
REAL :: dlayer1
REAL :: dlayer2
REAL :: zflat
! INTEGER :: strhtune
REAL :: strhtune
INTEGER :: mapproj
REAL :: trulat1
REAL :: trulat2
REAL :: trulon
REAL :: sclfct
REAL :: dxin
REAL :: dyin
REAL :: dzin
REAL :: dzminin
REAL :: ctrlatin
REAL :: ctrlonin
INTEGER :: strhoptin
REAL :: zrefsfcin
REAL :: dlayer1in
REAL :: dlayer2in
REAL :: zflatin
! INTEGER :: strhtunein
REAL :: strhtunein
INTEGER :: mapprojin
REAL :: trulat1in
REAL :: trulat2in
REAL :: trulonin
REAL :: sclfctin
INTEGER :: ireturn
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
REAL :: eps
PARAMETER (eps= 0.01)
!-----------------------------------------------------------------------
!
! Include file:
!
!-----------------------------------------------------------------------
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!-----------------------------------------------------------------------
!
! Compare all the variables. (Note that for 3-D data sets one ought
! to also check hterain, which is not done here.)
!
!-----------------------------------------------------------------------
ireturn = 0
IF ((nxin /= nx) .OR. (nyin /= ny) .OR. (nzin /= nz)) THEN
WRITE (6,'(/2(/5x,a),2(/5x,3(a,i9))/)') &
'ERROR: nx, ny or nz in file does not match', &
'the values in the program.', &
'In file, nx=', nxin, ', ny=', nyin, ', nz=',nzin, &
'In executable, nx=', nx, ', ny=', ny, ', nz=',nz
ireturn = 1
END IF
IF(ABS(dxin-dx) > eps .OR. ABS(dyin-dy) > eps .OR. &
ABS(dzin-dz) > eps .OR. &
ABS(dzminin-dzmin) > eps ) THEN
WRITE(6,'(/2(/5x,a),4(/5x,2(a,f13.3))/)') &
'ERROR: dx, dy, dz or dzmin in the file ', &
'do not match those specified in the input file.', &
'In file, dx=', dxin, ', dy=', dyin, &
' dz=', dzin, ', dzmin=', dzminin, &
'In input, dx=', dx, ', dy=', dy, &
' dz=', dz, ', dzmin=', dzmin
ireturn = 1
END IF
IF(ABS(ctrlatin-ctrlat) > eps .OR. ABS(ctrlonin-ctrlon) > eps ) THEN
WRITE(6,'(/2(/5x,a),2(/5x,2(a,f13.3))/)') &
'ERROR: Central latitude and/or longitude of the grid ', &
'in the file do not match those specified in input file.', &
'In file, ctrlat=',ctrlatin,', ctrlon=',ctrlonin, &
'In input, ctrlat=',ctrlat ,', ctrlon=',ctrlon
ireturn = 1
END IF
IF(strhoptin /= strhopt .OR. ABS(zrefsfcin-zrefsfc) > eps .OR. &
ABS(dlayer1in-dlayer1) > eps .OR. &
ABS(dlayer2in-dlayer2) > eps .OR. &
ABS(zflatin-zflat) > eps .OR. &
ABS(strhtunein-strhtune) > eps .OR. &
mapprojin /= mapproj .OR. &
ABS(trulat1in-trulat1) > eps .OR. &
ABS(trulat2in-trulat2) > eps .OR. &
ABS(trulonin-trulon ) > eps .OR. &
ABS(sclfctin-sclfct ) > eps) THEN
WRITE(6,'(/2(/5x,a),2(/5x,3(a,f13.3),2(a,i3),5(a,f13.3))/)') &
'ERROR: Map projection or other grid parameters do not ', &
'match those specified in input file.', &
'In file, trulat1=',trulat1in,', trulat2=',trulat2in, &
', trulon=',trulonin,', mapproj=',mapprojin, &
', strhopt=',strhoptin,', zrefsfc=',zrefsfcin, &
', dlayer1=',dlayer1in,', dlayer2=',dlayer2in, &
', zflat=',zflatin,', strhtune=',strhtunein, &
'In input, trulat1=',trulat1 ,', trulat2=',trulat2, &
', trulon=',trulon, ', mapproj=',mapproj, &
', strhopt=',strhopt,', zrefsfc=',zrefsfc, &
', dlayer1=',dlayer1,', dlayer2=',dlayer2, &
', zflat=',zflat,', strhtune=',strhtune
ireturn = 1
END IF
RETURN
END SUBROUTINE checkgrid3d
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CHECKGRID2D ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE checkgrid2d(nx,ny,nxin,nyin, & 6
dx,dy,ctrlat,ctrlon, &
mapproj,trulat1,trulat2,trulon,sclfct, &
dxin,dyin,ctrlatin,ctrlonin, &
mapprojin,trulat1in,trulat2in,trulonin,sclfctin,ireturn)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Check grid information to see if input data is consistent with
! ARPS grid.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Gene Bassett
! 2000/03/24
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx
! ny
! nxin
! nyin
! dx
! dy
! ctrlat
! ctrlon
! mapproj
! trulat1
! trulat2
! trulon
! sclfct
! dxin
! dyin
! ctrlatin
! ctrlonin
! mapprojin
! trulat1in
! trulat2in
! trulonin
! sclfctin
!
! OUTPUT:
!
! ireturn Flag indicating if the grids are the same (0-okay,
! 1-significant differences)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nxin
INTEGER :: nyin
INTEGER :: nx,ny
REAL :: dx
REAL :: dy
REAL :: ctrlat
REAL :: ctrlon
INTEGER :: mapproj
REAL :: trulat1
REAL :: trulat2
REAL :: trulon
REAL :: sclfct
REAL :: dxin
REAL :: dyin
REAL :: ctrlatin
REAL :: ctrlonin
INTEGER :: mapprojin
REAL :: trulat1in
REAL :: trulat2in
REAL :: trulonin
REAL :: sclfctin
INTEGER :: ireturn
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
REAL :: eps
PARAMETER (eps= 0.1)
!-----------------------------------------------------------------------
!
! Include file:
!
!-----------------------------------------------------------------------
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!-----------------------------------------------------------------------
!
! Compare all the variables.
!
!-----------------------------------------------------------------------
ireturn = 0
IF ((nxin /= nx) .OR. (nyin /= ny)) THEN
WRITE (6,'(/2(/5x,a),2(/5x,(a,i9))/)') &
'ERROR: nx and/or ny in file does not match', &
'the values in the program.', &
'In file, nx=', nxin, ', ny=', nyin, &
'In executable, nx=', nx, ', ny=', ny
ireturn = 1
END IF
IF(ABS(dxin-dx) > eps .OR. ABS(dyin-dy) > eps) THEN
WRITE(6,'(/2(/5x,a),2(/5x,2(a,f13.3))/)') &
'ERROR: dx or dy in the file ', &
'do not match those specified in the input file.', &
'In file, dx=', dxin, ', dy=', dyin, &
'In input, dx=', dx, ', dy=', dy
ireturn = 1
END IF
IF(ABS(ctrlatin-ctrlat) > eps .OR. ABS(ctrlonin-ctrlon) > eps ) THEN
WRITE(6,'(/2(/5x,a),2(/5x,2(a,f13.3))/)') &
'ERROR: Central latitude and/or longitude of the grid ', &
'in the file do not match those specified in input file.', &
'In file, ctrlat=',ctrlatin,', ctrlon=',ctrlonin, &
'In input, ctrlat=',ctrlat ,', ctrlon=',ctrlon
ireturn = 1
END IF
IF(mapprojin /= mapproj .OR. ABS(trulat1in-trulat1) > eps .OR. &
ABS(trulat2in-trulat2) > eps .OR. &
ABS(trulonin-trulon ) > eps .OR. &
ABS(sclfctin-sclfct ) > eps) THEN
WRITE(6,'(/2(/5x,a),2(/5x,3(a,f13.3),(a,i3,a,f13.3))/)') &
'ERROR: Map projection or other grid parameters do not ', &
'match those specified in input file.', &
'In file, trulat1=',trulat1in,', trulat2=',trulat2in, &
', trulon=',trulonin,', mapproj=',mapprojin,' sclfct=',sclfctin, &
'In input, trulat1=',trulat1 ,', trulat2=',trulat2, &
', trulon=',trulon, ', mapproj=',mapproj,' sclfct=',sclfct
ireturn = 1
END IF
RETURN
END SUBROUTINE checkgrid2d
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INITZTIME ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE initztime(ayear,amm,aday) 7
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To calculate Zulu time based upon Mesonet observations read directly
! for any given day.
!
!-----------------------------------------------------------------------
!
! AUTHOR: J.A. Brotzge
! 2001/08/23
!
!-----------------------------------------------------------------------
!
!
! INPUT:
!
! curtim Computational time (seconds)
!
! OUTPUT:
!
! ztime Zulu time, calculated (seconds)
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
!
! INPUT
!
! INTEGER :: hour, minute, second ! Computational time
! REAL :: curtim ! Computational time (seconds)
!
! OUTPUT
!
! REAL :: ztime !Zulu time (seconds)
!
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
INTEGER :: dayt, initday
CHARACTER :: ayear*4, amm*2, aday*2
REAL :: modelinitime, day1time, day2time, day3time, day4time
REAL :: day5time, day6time, day7time, day8time
REAL :: day9time, day10time
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
modelinitime = real(hour*3600 + minute*60 + second)
day1time = 86400.0 - modelinitime
day2time = day1time + 86400.0
day3time = day2time + 86400.0
day4time = day3time + 86400.0
day5time = day4time + 86400.0
day6time = day5time + 86400.0
day7time = day6time + 86400.0
day8time = day7time + 86400.0
day9time = day8time + 86400.0
day10time= day9time + 86400.0
IF (curtim.lt.day1time) THEN
dayt = day
ztime = modelinitime + curtim
ENDIF
IF (curtim.ge.day1time.and.curtim.lt.day2time) THEN
dayt = day + 1
ztime = curtim - day1time
ENDIF
IF (curtim.ge.day2time.and.curtim.lt.day3time) THEN
dayt = day + 2
ztime = curtim - day2time
ENDIF
IF (curtim.ge.day3time.and.curtim.lt.day4time) THEN
dayt = day + 3
ztime = curtim - day3time
ENDIF
IF (curtim.ge.day4time.and.curtim.lt.day5time) THEN
dayt = day + 4
ztime = curtim - day4time
ENDIF
IF (curtim.ge.day5time.and.curtim.lt.day6time) THEN
dayt = day + 5
ztime = curtim - day5time
ENDIF
IF (curtim.ge.day6time.and.curtim.lt.day7time) THEN
dayt = day + 6
ztime = curtim - day6time
ENDIF
IF (curtim.ge.day7time.and.curtim.lt.day8time) THEN
dayt = day + 7
ztime = curtim - day7time
ENDIF
IF (curtim.ge.day8time.and.curtim.lt.day9time) THEN
dayt = day + 8
ztime = curtim - day8time
ENDIF
IF (curtim.ge.day9time.and.curtim.lt.day10time) THEN
dayt = day + 9
ztime = curtim - day9time
ENDIF
IF (curtim.ge.day10time) THEN
dayt = day + 10
ztime = curtim - day10time
ENDIF
IF (month.eq.2.and.day.gt.29) THEN
month=month+1
dayt=1
ENDIF
IF (month.eq.2.and.day.eq.29) THEN
IF (year.ne.2000) THEN
month=month+1
dayt=1
ENDIF
ENDIF
IF (dayt.eq.31) THEN
IF (month.eq.4.or.month.eq.6) THEN
month=month+1
dayt=1
ENDIF
IF (month.eq.9.or.month.eq.11) THEN
month=month+1
dayt=1
ENDIF
ENDIF
IF (dayt.eq.32) THEN
month=month+1
dayt=1
ENDIF
write(ayear,'(i4)') year
write(amm,'(i2)') month
write(aday,'(i2)') dayt
IF (month.lt.10) amm(1:1)='0'
IF (dayt.lt.10) aday(1:1)='0'
RETURN
END SUBROUTINE initztime
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE READJSOIL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE readjsoil(nx,ny,nzsoil,numbsty,ayear,amm,aday,arpstime, & 1,2
zpsoil,tsoil,qsoil,wetcanp,snowdpth)
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read in the mesonet radt'n data qc'd and processed by Jerry Brotzge.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber and Jerry Brotzge
! 8/28/2001.
!
! 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)
! nzsoil Number of soil layers
! numbsty Number of soil types within grid box
!
! OUTPUT:
!
! tsoil Soil temperature (K)
! qsoil Soil moisture (m^3/m^3)
! wetcanp Canopy soil moisture (unitless)
! snowdpth Snow depth (m)
! soiltyp Soil type (unitless)
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
!
! INPUT
!
INTEGER :: nx, ny ! Number of grid points in 3 directions
INTEGER :: nzsoil ! Number of soil layers
INTEGER :: numbsty ! Number of soil types within grid box
INTEGER, PARAMETER :: nmeso=48, mmeso=288
INTEGER, PARAMETER :: nzsoilext=5
INTEGER, PARAMETER :: nzsoil_ext=5
REAL :: t105 (nmeso)
REAL :: t125 (nmeso)
REAL :: t160 (nmeso)
REAL :: t175 (nmeso)
REAL :: smo05 (nmeso)
REAL :: smo25 (nmeso)
REAL :: smo60 (nmeso)
REAL :: smo75 (nmeso)
REAL :: wrxx (nmeso)
REAL :: qvsf (nmeso)
REAL :: snow (nmeso)
REAL :: rnet (mmeso) ! Mesonet CNR1 rnet (W/m**2)
REAL :: hflx (mmeso) ! Mesonet sonic sensible heat flux (W/m**2)
REAL :: lflx (mmeso) ! Mesonet sonic latent heat flux (W/m**2)
REAL :: gflx (mmeso) ! Mesonet ground heat flux (W/m**2)
REAL :: l2fx (mmeso) ! Mesonet residual latent heat flux (W/m**2)
REAL :: irtx (mmeso) ! Mesonet rainfall rate (kg/m^2/s)
REAL :: irthour (mmeso)
REAL :: prof (mmeso)
CHARACTER*4 :: temp1 (nmeso) ! Mesonet station name
INTEGER :: stnm1 (nmeso) ! Mesonet station number
INTEGER :: mesotime(nmeso) ! Mesonet time (minutes then seconds)
CHARACTER*4 :: ftemp1 (mmeso) ! Mesonet station name
INTEGER :: fstnm1 (mmeso) ! Mesonet station number
INTEGER :: fmesotime(mmeso) ! Mesonet time (minutes then seconds)
REAL :: mesot(nmeso) ! Mesonet time (seconds)
!
! OUTPUT
!
REAL :: zpsoil(nx,ny,nzsoil) ! Depth of soil (m)
REAL :: tsoil(nx,ny,nzsoil,numbsty) ! Soil temperature (K)
REAL :: qsoil(nx,ny,nzsoil,numbsty) ! Soil moisture (m3/m3)
REAL :: wetcanp(nx,ny,numbsty) ! Canopy moisture
REAL :: snowdpth(nx,ny) ! Snow depth
INTEGER :: soiltyp(nx,ny,numbsty) !Soil type
!
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,n,nt
CHARACTER :: ayear*4, amm*2, aday*2
CHARACTER :: atemp*4, btemp*2, ctemp*2, dtemp*3
CHARACTER :: astid*4, astnm*4, aime*4
CHARACTER :: a125*4, a160*4, a175*4, atm05*4, atm25*4
CHARACTER :: atm60*4, atm75*4, awrxx*4, aqvsf*4, asnow*4
CHARACTER :: a105*4
CHARACTER :: arnet*4, ahflx*4, alflx*4, al2fx*4
CHARACTER :: agflx*4, airtx*4
CHARACTER :: bstid*4, bstnm*4, bime*4
CHARACTER :: etemp*4, ftemp*4, gtemp*4, htemp*4
INTEGER :: stnm, cmm, cdd
INTEGER :: fmm, fdd, filelength
REAL :: tema, temb, convertp
REAL :: arpstime
REAL :: zpsoil_ext(nx,ny,nzsoilext)
REAL :: tsoil_ext(nx,ny,nzsoilext), qsoil_ext(nx,ny,nzsoilext)
! DATA zpsoil_ext/0.0,-0.05,-0.25,-0.60,-0.75/
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'indtflg.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
filelength = len_trim(sitesoil)
! nzsoil_ext = 5
OPEN(unit=60,file=sitesoil(1:filelength)// &
ayear//amm//aday//'.mts',status='old',form='formatted')
READ(60,311) dtemp
READ(60,313) btemp, atemp, cmm, cdd, ctemp, ctemp, ctemp
READ(60,301) astid, astnm, aime, a105, atm05, &
& a125, atm25, a160, atm60, a175, atm75, &
& awrxx, aqvsf, asnow
301 FORMAT (a4,13(1x,a4))
302 FORMAT (1x,a4,4x,i3,4x,i5,11(1x,f9.4))
311 FORMAT (6x,a3)
313 FORMAT (7x,a2,1x,a4,1x,i2,1x,i2,1x,a2,1x,a2,1x,a2)
filelength = len_trim(siteflux)
OPEN(unit=20,file=siteflux(1:filelength)// &
& ayear//amm//aday//'.mts',status='old',form='formatted')
read(20,511) htemp
read(20,513) ftemp, etemp, fmm, fdd, gtemp, gtemp, gtemp
read(20,501) bstid, bstnm, bime, arnet, ahflx, alflx, &
agflx, al2fx, airtx
501 FORMAT (a4,8(1x,a4))
502 format (1x,a4,4x,i3,5x,i4,7(1x,f9.4))
511 FORMAT (6x,a3)
513 FORMAT (7x,a2,1x,a4,1x,i2,1x,i2,1x,a2,1x,a2,1x,a2)
nt = 0
DO n=1,288 ! perform the file read....
read(20,502) ftemp1(n), fstnm1(n), fmesotime(n), &
& rnet(n), hflx(n), lflx(n), &
& gflx(n), l2fx(n), irtx(n), prof(n)
IF (n.eq.1.or.mod(n,6).eq.0.0) THEN
nt = nt + 1
irthour(nt) = irtx(n) + 273.15 !Convert C to K
ENDIF
ENDDO !end of file read for Flux files
DO j=1,ny
DO i=1,nx
zpsoil_ext(i,j,1) = 0.00
zpsoil_ext(i,j,2) = -0.05
zpsoil_ext(i,j,3) = -0.25
zpsoil_ext(i,j,4) = -0.60
zpsoil_ext(i,j,5) = -0.75
END DO
END DO
DO n=1,48 ! perform the file read ....
READ(60,302) temp1(n),stnm1(n),mesotime(n), &
t105(n),smo05(n),t125(n),smo25(n), &
t160(n),smo60(n),t175(n),smo75(n), &
wrxx(n),qvsf(n),snow(n)
t105(n) = t105(n) + 273.15
t125(n) = t125(n) + 273.15
t160(n) = t160(n) + 273.15
t175(n) = t175(n) + 273.15
mesotime(n) = mesotime(n) * 60
mesot(n) = real(mesotime(n))
IF (arpstime.le.mesot(n)) THEN
IF(arpstime.eq.mesot(n)) THEN ! we have an exact match
DO j=1,ny
DO i=1,nx
tsoil_ext(i,j,1) = irthour(n)
tsoil_ext(i,j,2) = t105(n)
tsoil_ext(i,j,3) = t125(n)
tsoil_ext(i,j,4) = t160(n)
tsoil_ext(i,j,5) = t175(n)
qsoil_ext(i,j,1) = smo05(n)
qsoil_ext(i,j,2) = smo05(n)
qsoil_ext(i,j,3) = smo25(n)
qsoil_ext(i,j,4) = smo60(n)
qsoil_ext(i,j,5) = smo75(n)
END DO
END DO
CALL initsoil(nx,ny,nzsoil,nzsoil_ext,numbsty,zpsoil, &
zpsoil_ext,tsoil,tsoil_ext,qsoil,qsoil_ext)
DO j = 1,ny
DO i = 1,nx
wetcanp(i,j,1) = wrxx(n)
snowdpth(i,j) = snow(n)
END DO
END DO
ELSE IF(arpstime.gt.mesot(n-1).and.arpstime.le.mesot(n)) THEN
tema = (arpstime-mesot(n-1))/1800.0
temb = 1.0-tema
DO j = 1,ny
DO i = 1,nx
tsoil_ext(i,j,1) = temb*irthour(n-1)+tema*irthour(n)
tsoil_ext(i,j,2) = temb*t105(n-1)+tema*t105(n)
tsoil_ext(i,j,3) = temb*t125(n-1)+tema*t125(n)
tsoil_ext(i,j,4) = temb*t160(n-1)+tema*t160(n)
tsoil_ext(i,j,5) = temb*t175(n-1)+tema*t175(n)
qsoil_ext(i,j,1) = temb*smo05(n-1)+tema*smo05(n)
qsoil_ext(i,j,2) = temb*smo05(n-1)+tema*smo05(n)
qsoil_ext(i,j,3) = temb*smo25(n-1)+tema*smo25(n)
qsoil_ext(i,j,4) = temb*smo60(n-1)+tema*smo60(n)
qsoil_ext(i,j,5) = temb*smo75(n-1)+tema*smo75(n)
wetcanp(i,j,1) = temb*wrxx(n-1)+tema*wrxx(n)
snowdpth(i,j) = temb*snow(n-1)+tema*snow(n)
END DO
END DO
CALL initsoil(nx,ny,nzsoil,nzsoil_ext,numbsty,zpsoil, &
zpsoil_ext,tsoil,tsoil_ext,qsoil,qsoil_ext)
END IF ! end of the time if loop
END IF ! end of the time if loop for matching arpstime and mesot
ENDDO !Time loop
close (20)
close (60)
sfcin = 1
return
RETURN
END SUBROUTINE readjsoil
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE READJVEG ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE readjveg(nx,ny,numbsty,ayear,amm,aday,arpstime, & 2
soiltyp,stypfrct,vegtyp,lai,roufns,veg,ndvi)
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read in the mesonet radt'n data qc'd and processed by Jerry Brotzge.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dan Weber and Jerry Brotzge
! 8/28/2001.
!
! 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)
! numbsty Number of soil types within grid box
!
! OUTPUT:
!
! soiltyp Soil type
! stypfrct Surface type fraction
! vegtyp Vegetation type (unitless)
! lai Leaf area index ()
! roufns Roughness (unitless)
! veg Vegetation fraction (unitless)
! ndvi Normalized Diff Veg Index (unitless)
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!-----------------------------------------------------------------------
!
!
IMPLICIT NONE ! Force explicit declarations
!
! INPUT
!
INTEGER :: nx, ny ! Number of grid points in 3 directions
INTEGER :: numbsty ! Number of soil types within grid box
INTEGER :: nst ! Number of soil types
INTEGER, PARAMETER :: nmeso = 48
REAL :: veg2 (nmeso)
REAL :: veg4 (nmeso)
REAL :: veg5 (nmeso)
REAL :: veg6 (nmeso)
REAL :: veg7 (nmeso)
INTEGER :: veg1(nmeso)
INTEGER :: veg3(nmeso)
CHARACTER*4 :: temp1 (nmeso) ! Mesonet station name
INTEGER :: stnm1 (nmeso) ! Mesonet station number
INTEGER :: mesotime(nmeso) !Mesonet time (minutes then seconds)
REAL :: mesot(nmeso) ! Mesonet time (seconds)
!
! OUTPUT
!
REAL :: stypfrct(nx,ny,numbsty) !Surface type fraction
REAL :: lai(nx,ny) !Leaf Area Index
REAL :: roufns(nx,ny) !Surface roughness
REAL :: veg(nx,ny) !Vegetation
REAL :: ndvi(nx,ny) !Normalized Diff Veg Index
INTEGER :: soiltyp(nx,ny,numbsty) !Soil type
INTEGER :: vegtyp(nx,ny) !Veg type
!
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,n
CHARACTER :: ayear*4, amm*2, aday*2
CHARACTER :: atemp*4, btemp*2, ctemp*2, dtemp*3
CHARACTER :: astid*4, astnm*4, aime*4
CHARACTER :: astyp*4, asfct*4, avtyp*4, alaix*4, arfns*4
CHARACTER :: avfct*4, andvi*4
INTEGER :: stnm, cmm, cdd, filelength
REAL :: tema, temb, convertp
REAL :: arpstime
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
filelength = len_trim(siteveg)
OPEN(unit=70,file=siteveg(1:filelength)// &
& ayear//amm//aday//'.mts',status='old',form='formatted')
! loop that reads Jerry's mesonet data file (qc'd)
! note that the ARPS time must match the mesonet datafile time.
READ(70,711) dtemp
READ(70,713) btemp, atemp, cmm, cdd, ctemp, ctemp, ctemp
READ(70,701) astid, astnm, aime, astyp, asfct, &
avtyp, alaix, arfns, avfct, andvi
701 FORMAT (a4,9(1x,a4))
702 FORMAT (1x,a4,4x,i3,5x,i4,1x,i2,1x,f4.2,1x,i2, &
4(1x,f4.2))
711 FORMAT (6x,a3)
713 FORMAT (7x,a2,1x,a4,1x,a2,1x,a2,1x,a2,1x,a2,1x,a2)
do n=1,48 ! perform the file read....
READ(70,702) temp1(n),stnm1(n),mesotime(n), &
veg1(n),veg2(n),veg3(n),veg4(n), &
veg5(n),veg6(n),veg7(n)
mesotime(n) = mesotime(n) * 60
mesot(n) = real(mesotime(n))
IF (arpstime.le.mesot(n)) THEN
IF(arpstime.eq.mesot(n)) THEN ! we have an exact match
DO j = 1,ny
DO i = 1,nx
soiltyp(i,j,1) = veg1(n)
stypfrct(i,j,1) = veg2(n)
vegtyp(i,j) = veg3(n)
lai(i,j) = veg4(n)
roufns(i,j) = veg5(n)
veg(i,j) = veg6(n)
ndvi(i,j) = veg7(n)
END DO
END DO
ELSE IF(arpstime.gt.mesot(n-1).and.arpstime.le.mesot(n))THEN
tema = (arpstime-mesot(n-1))/1800.0
temb = 1.0-tema
DO j = 1,ny
DO i = 1,nx
soiltyp(i,j,1) = veg1(n)
stypfrct(i,j,1) = veg2(n)
vegtyp(i,j) = veg3(n)
lai(i,j) = veg4(n)
roufns(i,j) = veg5(n)
veg(i,j) = veg6(n)
ndvi(i,j) = veg7(n)
END DO
END DO
END IF ! end of the time if loop
END IF ! end of the time if loop for matching arpstime and mesot
ENDDO !Time loop
close (70)
return
RETURN
END SUBROUTINE readjveg
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INITSOIL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE initsoil (nx,ny,nzsoil,nzsoil_ext, nstyp,zpsoil, & 2
zpsoil_ext,tsoil,tsoil_ext,qsoil,qsoil_ext)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Take soil data from external file and interpolate in
! the vertical to form the ARPS 2D data set profile
!
!-----------------------------------------------------------------------
!
! AUTHOR: Jerry Brotzge
! 05/15/2002
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! for the ARPS grid
! ny Number of grid points in the y-direction (north/south)
! for the ARPS grid
! nzsoil Number of grid points in the soil
!
! nzsoil_ext Number of grid points in the soil external file
!
! nstyp Number of soil types
!
! zpsoil Physical depth of soil layers (m)
!
! zpsoil_ext Physical depth of external model soil layers (m)
!
! tsoil Soil temperature (K)
!
! qsoil Soil moisture (m**3/m**3)
!
! OUTPUT:
!
! tsoil Soil temperature profile (K)
!
! qsoil Soil moisture profile (m**3/m**3)
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!
!-----------------------------------------------------------------------
!
! Input/output variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,kk
INTEGER :: nx,ny ! ARPS grid dimensions
INTEGER :: nzsoil ! ARPS soil levels
INTEGER :: nzsoil_ext ! External file soil levels
INTEGER :: nstyp ! Number of soil types
!
REAL :: zpsoil(nx,ny,nzsoil) ! Physical depth of ARPS soil layers
REAL :: zpsoil_ext(nx,ny,nzsoil_ext) ! Physical depth of ext. file soil layers
REAL :: tsoil(nx,ny,nzsoil) ! Soil temperature (K)
REAL :: qsoil(nx,ny,nzsoil) ! Soil moisture (m**3/m**3)
REAL :: tsoil_ext(nx,ny,nzsoil_ext) ! External file soil temperature (K)
REAL :: qsoil_ext(nx,ny,nzsoil_ext) ! External file soil moisture (m**3/m**3)
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
INTEGER :: index
REAL :: dampdepth
REAL :: tsoil15
REAL :: qsoil15
REAL :: dampdz
REAL :: depthdz(nx,ny,nzsoil)
REAL :: depthdzext
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Vertical interpolation
!
!-----------------------------------------------------------------------
!
! NOTE: Index INCREASES WITH DEPTH. k=1 at the top; nzsoil at the bottom BC.
!
!----------------------------------------------------------------
! Set top boundary condition at level = 1
!-------------------------------------------------------------
DO i=1,nx
DO j=1,ny
tsoil (i,j,1) = tsoil_ext(i,j,1)
qsoil (i,j,1) = qsoil_ext(i,j,1)
DO kk=1,nzsoil_ext
zpsoil_ext(i,j,kk) = zpsoil_ext(i,j,kk) + zrefsfc
END DO
END DO
END DO
!----------------------------------------------------------------------
! Initialize soil temp and moisture profiles
!----------------------------------------------------------------------
! Note: All soil depths are negative (downward)
DO k=2,nzsoil
DO j=1,ny
DO i=1,nx
dampdepth = zpsoil(i,j,1) - 0.15 !Typical damping depth = 15 cm
dampdz = zpsoil(i,j,k) - dampdepth
depthdz(i,j,k) = zrefsfc - zpsoil(i,j,k)
tsoil15 = 0.0
qsoil15 = 0.0
index = 0
DO kk=2,nzsoil_ext
depthdzext = zrefsfc - zpsoil_ext(i,j,kk)
IF (((zpsoil_ext(i,j,kk-1) >= dampdepth).AND.(zpsoil_ext(i,j,kk) <= &
dampdepth)).AND.(index < 1)) THEN
tsoil15 = tsoil_ext(i,j,kk-1)+((tsoil_ext(i,j,kk)- &
tsoil_ext(i,j,kk-1))* ((dampdepth - zpsoil_ext(i,j,kk))/ &
(zpsoil_ext(i,j,kk-1) - zpsoil_ext(i,j,kk))) )
qsoil15 = qsoil_ext(i,j,kk-1)+((qsoil_ext(i,j,kk)- &
qsoil_ext(i,j,kk-1))* ((dampdepth - zpsoil_ext(i,j,kk))/ &
(zpsoil_ext(i,j,kk-1) - zpsoil_ext(i,j,kk))) )
index = index + 1
ENDIF
ENDDO
DO kk=2,nzsoil_ext
!----------------------------------------------------------------------
! Initialize levels equal to obs when levels are equal.
!----------------------------------------------------------------------
IF (zpsoil(i,j,k) == zpsoil_ext(i,j,kk)) THEN
tsoil(i,j,k) = tsoil_ext(i,j,kk)
qsoil(i,j,k) = qsoil_ext(i,j,kk)
!----------------------------------------------------------------------
! Exponential fit to initialization above damping depth (15 cm)
!----------------------------------------------------------------------
ELSE IF (zpsoil(i,j,k) >= dampdepth) THEN !Soil level within top 15 cm
IF ((zpsoil(i,j,k) > zpsoil_ext(i,j,kk)).AND.(zpsoil_ext(i,j,kk) >= &
dampdepth)) THEN
tsoil(i,j,k) = tsoil_ext(i,j,kk)+((tsoil(i,j,k-1)- &
tsoil_ext(i,j,kk))*EXP(- (depthdz(i,j,k)/depthdzext)) )
qsoil(i,j,k) = qsoil_ext(i,j,kk)+((qsoil(i,j,k-1)- &
qsoil_ext(i,j,kk))*EXP(- (depthdz(i,j,k)/depthdzext)) )
ELSE IF ((zpsoil(i,j,k) > zpsoil_ext(i,j,kk)).AND.(zpsoil_ext(i,j,kk) < &
dampdepth)) THEN !(Est's linear fit upward to 15cm from tsoil_ext)
tsoil(i,j,k)=tsoil15 + (tsoil(i,j,k-1)-tsoil15)* &
EXP( - (depthdz(i,j,k)/dampdz) )
qsoil(i,j,k)=qsoil15 + (qsoil(i,j,k-1)-qsoil15)* &
EXP( - (depthdz(i,j,k)/dampdz) )
ELSE IF ((zpsoil(i,j,k) < zpsoil_ext(i,j,kk-1)).AND.(zpsoil(i,j,k) > &
zpsoil_ext(i,j,kk))) THEN
tsoil(i,j,k)=tsoil_ext(i,j,kk)+(tsoil_ext(i,j,kk-1) - &
tsoil_ext(i,j,kk)) * EXP( - (depthdz(i,j,k)/depthdzext) )
qsoil(i,j,k)=qsoil_ext(i,j,kk)+(qsoil_ext(i,j,kk-1) - &
qsoil_ext(i,j,kk)) * EXP( - (depthdz(i,j,k)/depthdzext) )
END IF
!----------------------------------------------------------------------
! Linear fit between initialized levels below damping depth (15 cm)
!----------------------------------------------------------------------
ELSE IF (zpsoil(i,j,k) < dampdepth) THEN !Soil level below top 15 cm
IF (zpsoil(i,j,k) < zpsoil_ext(i,j,kk-1) .AND. &
zpsoil(i,j,k) > zpsoil_ext(i,j,kk)) THEN
tsoil(i,j,k) = tsoil_ext(i,j,kk)+((tsoil_ext(i,j,kk-1)- &
tsoil_ext(i,j,kk)) * (zpsoil(i,j,k) - zpsoil_ext(i,j,kk))/ &
(zpsoil_ext(i,j,kk-1) - zpsoil_ext(i,j,kk)) )
qsoil(i,j,k) = qsoil_ext(i,j,kk)+((qsoil_ext(i,j,kk-1)- &
qsoil_ext(i,j,kk)) * (zpsoil(i,j,k) - zpsoil_ext(i,j,kk))/ &
(zpsoil_ext(i,j,kk-1) - zpsoil_ext(i,j,kk)) )
ELSE IF (zpsoil(i,j,k) > zpsoil_ext(i,j,2)) THEN
tsoil(i,j,k) = tsoil_ext(i,j,kk)+((tsoil15 - tsoil_ext(i,j,kk)) * &
((zpsoil(i,j,k) - zpsoil_ext(i,j,kk))/(dampdepth - &
zpsoil_ext(i,j,kk))) )
qsoil(i,j,k) = qsoil_ext(i,j,kk)+((qsoil15 - qsoil_ext(i,j,kk)) * &
((zpsoil(i,j,k) - zpsoil_ext(i,j,kk))/(dampdepth - &
zpsoil_ext(i,j,kk))) )
ELSE IF (zpsoil(i,j,k) < zpsoil_ext(i,j,nzsoil_ext)) THEN
tsoil(i,j,k) = tsoil_ext(i,j,kk)
qsoil(i,j,k) = qsoil_ext(i,j,kk)
END IF
END IF !Damping depth if/then
END DO
END DO
END DO
END DO
RETURN
END SUBROUTINE initsoil