!
!##################################################################
!###### ######
!###### SUBROUTINE WRTSOIL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE wrtsoil( nx,ny, nstyps, soiloutfl, dx,dy, & 7,24
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
tsfcout,tsoilout,wsfcout,wdpout,wcanpout,snowdout, &
tsfc,tsoil,wetsfc,wetdp,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.
!
!-----------------------------------------------------------------------
!
! 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)
!
! soiloutfl Name of the soil model data output file
!
! tsfcout Flag for output of Temperature at ground surface
! tsoilout Flag for output of Deep soil temperature
! wsfcout Flag for output of Surface soil moisture
! wdpout Flag for output of Deep soil moisture
! wcanpout Flag for output of Canopy water amount
!
! tsfc Temperature at ground surface (K)
! tsoil Deep soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! 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
!
! 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)
!
! tsfcout Flag for output of Temperature at ground surface
! tsoilout Flag for output of Deep soil temperature
! wsfcout Flag for output of Surface soil moisture
! wdpout Flag for output of Deep soil moisture
! wcanpout Flag for output of Canopy water amount
!
! tsfc Temperature at ground surface (K)
! tsoil Deep soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! 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 :: 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 :: tsfcout ! Flag for output of tsfc
INTEGER :: tsoilout ! Flag for output of tsoil
INTEGER :: wsfcout ! Flag for output of wetsfc
INTEGER :: wdpout ! Flag for output of wetdp
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 :: tsfc (nx,ny,0:nstyps) ! Temperature at ground surface (K)
REAL :: tsoil (nx,ny,0:nstyps) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture
REAL :: wetdp (nx,ny,0:nstyps) ! Deep soil moisture
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
REAL :: snowdpth(nx,ny) ! Snow depth (m)
INTEGER :: soiltyp (nx,ny,nstyps) ! Soil type in model domain
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,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 :: stat, sd_id
!
!-----------------------------------------------------------------------
!
! 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 ( tsfcout+tsoilout+wsfcout+wdpout+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)
OPEN (UNIT=flunit, FILE=trim(soiloutfl), &
STATUS='unknown', FORM='unformatted', ACCESS='sequential')
WRITE (flunit) nx, ny
WRITE (flunit) mapproj, tsfcout, tsoilout, wsfcout, wdpout, &
wcanpout, 0, snowdout, stypout, idummy, &
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 ( tsfcout /= 0 ) THEN
IF ( nstyp <= 1 ) THEN
WRITE (flunit) ((tsfc(i,j,0),i=1,nx),j=1,ny)
ELSE
!--------------------- A bug fixed by Yunheng ------------
! DO is=0,nstyp
! WRITE (flunit) ((tsfc(i,j,is),i=1,nx),j=1,ny)
! END DO
!---------------------------------------------------------
WRITE (flunit) tsfc
END IF
END IF
IF ( tsoilout /= 0 ) THEN
IF ( nstyp <= 1 ) THEN
WRITE (flunit) ((tsoil(i,j,0),i=1,nx),j=1,ny)
ELSE
!--------------------- A bug fixed by Yunheng ------------
! DO is=0,nstyp
! WRITE (flunit) ((tsoil(i,j,is),i=1,nx),j=1,ny)
! END DO
!---------------------------------------------------------
WRITE (flunit) tsoil
END IF
END IF
IF ( wsfcout /= 0 ) THEN
IF ( nstyp <= 1 ) THEN
WRITE (flunit) ((wetsfc(i,j,0),i=1,nx),j=1,ny)
ELSE
!--------------------- A bug fixed by Yunheng ------------
! DO is=0,nstyp
! WRITE (flunit) ((wetsfc(i,j,is),i=1,nx),j=1,ny)
! END DO
!---------------------------------------------------------
WRITE (flunit) wetsfc
END IF
END IF
IF ( wdpout /= 0 ) THEN
IF ( nstyp <= 1 ) THEN
WRITE (flunit) ((wetdp(i,j,0),i=1,nx),j=1,ny)
ELSE
!--------------------- A bug fixed by Yunheng ------------
! DO is=0,nstyp
! WRITE (flunit) ((wetdp(i,j,is),i=1,nx),j=1,ny)
! END DO
!---------------------------------------------------------
WRITE (flunit) wetdp
END IF
END IF
IF ( wcanpout /= 0 ) THEN
IF ( nstyp <= 1 ) THEN
WRITE (flunit) ((wetcanp(i,j,0),i=1,nx),j=1,ny)
ELSE
!--------------------- A bug fixed by Yunheng ------------
! DO is=0,nstyp
! WRITE (flunit) ((wetcanp(i,j,is),i=1,nx),j=1,ny)
! END DO
!---------------------------------------------------------
WRITE (flunit) wetcanp
END IF
END IF
IF ( snowdout /= 0 ) THEN
WRITE (flunit) ((snowdpth(i,j),i=1,nx),j=1,ny)
END IF
IF ( stypout /= 0 ) THEN
WRITE (flunit) soiltyp
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.
!
!-----------------------------------------------------------------------
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 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 ( nstyp == 0 ) nstyp = 1
IF ( tsfcout /= 0 ) THEN
CALL hdfwrt3d(tsfc,nx,ny,nstyp+1,sd_id,0,0, &
'tsfc','Surface ground temperature','K', &
itmp,atmp1,atmp2)
END IF
IF ( tsoilout /= 0 ) THEN
CALL hdfwrt3d(tsoil,nx,ny,nstyp+1,sd_id,0,0, &
'tsoil','Deep soil temperature','K', &
itmp,atmp1,atmp2)
END IF
IF ( wsfcout /= 0 ) THEN
CALL hdfwrt3d(wetsfc,nx,ny,nstyp+1,sd_id,0,0, &
'wetsfc','Surface soil moisture','fraction', &
itmp,atmp1,atmp2)
END IF
IF ( wdpout /= 0 ) THEN
CALL hdfwrt3d(wetdp,nx,ny,nstyp+1,sd_id,0,0, &
'wetdp','Deep soil moisture','fraction', &
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,nstyps,soilfile,dx,dy, & 6,30
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin, &
tsfc,tsoil,wetsfc,wetdp,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.
!
!-----------------------------------------------------------------------
!
! 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:
!
! tsfc Temperature at ground surface (K)
! tsoil Deep soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! 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 :: 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 :: tsfcin
INTEGER :: tsoilin
INTEGER :: wsfcin
INTEGER :: wdpin
INTEGER :: wcanpin
INTEGER :: snowdin
INTEGER :: snowcin
INTEGER :: stypin
REAL :: tsfc (nx,ny,0:nstyps) ! Temperature at ground surface (K)
REAL :: tsoil (nx,ny,0:nstyps) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture
REAL :: wetdp (nx,ny,0:nstyps) ! Deep soil moisture
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
INTEGER :: soiltyp (nx,ny,nstyps) ! Soil type in model domain
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL, ALLOCATABLE :: tsfc_in (:,:,:) ! Temperature at ground surface (K)
REAL, ALLOCATABLE :: tsoil_in (:,:,:) ! Deep soil temperature (K)
REAL, ALLOCATABLE :: wetsfc_in (:,:,:) ! Surface soil moisture
REAL, ALLOCATABLE :: wetdp_in (:,:,:) ! Deep soil moisture
REAL, ALLOCATABLE :: wetcanp_in(:,:,:) ! Canopy water amount
INTEGER, ALLOCATABLE :: soiltyp_in (:,:,:) ! Soil type in model domain
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: nxin,nyin
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,is
INTEGER :: istat, ierr
INTEGER :: ireturn
INTEGER :: itmp1,itmp2
REAL :: rtmp1,rtmp2,rtmp3,rtmp4,rtmp5,rtmp6,rtmp7
CHARACTER (LEN=128) :: savename !Message passing code.
!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 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=998) nxin,nyin
READ (flunit,ERR=998) mprojin,tsfcin, tsoilin,wsfcin,wdpin, &
wcanpin,snowcin,snowdin,stypin,idummy, &
idummy, idummy, idummy, idummy,idummy, &
idummy, idummy, idummy, idummy,nstyp1
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
!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 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 ...
END IF
nstyp1 = MAX( nstyp1, 1 )
ALLOCATE (tsfc_in(nx,ny,0:nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "READSOIL: ERROR allocating tsfc_in, returning"
RETURN
END IF
ALLOCATE (tsoil_in(nx,ny,0:nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "READSOIL: ERROR allocating tsoil_in, returning"
RETURN
END IF
ALLOCATE (wetsfc_in(nx,ny,0:nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "READSOIL: ERROR allocating wetsfc_in, returning"
RETURN
END IF
ALLOCATE (wetdp_in(nx,ny,0:nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "READSOIL: ERROR allocating wetdp_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.
!
!-----------------------------------------------------------------------
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
!
!-----------------------------------------------------------------------
!
! Read in the soil data from the soil data file.
!
!-----------------------------------------------------------------------
!
IF (soilfmt == 0) THEN ! Fortran unformatted
WRITE (6, '(a/a,i2/a,i2/a,i2/a,i2/a,i2/a,i2/a,i2)') &
' Surface data flags for: ', &
' tsfcin --', tsfcin, &
' tsoilin --', tsoilin, &
' wsfcin --', wsfcin, &
' wdpin --', wdpin, &
' wcanpin --', wcanpin, &
' snowdin --', snowdin, &
' stypin --', stypin
IF ( tsfcin /= 0 ) THEN
WRITE(6,'(a)') 'Read in the surface skin temperature data'
!--------------------- A bug fixed by Yunheng ------------
IF ( nstyp <= 1 ) THEN
READ (flunit,ERR=998) ((tsfc_in(i,j,0),i=1,nx),j=1,ny)
ELSE
READ (flunit,ERR=998) tsfc_in
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) ((tsoil_in(i,j,0),i=1,nx),j=1,ny)
ELSE
READ (flunit,ERR=998) tsoil_in
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) ((wetsfc_in(i,j,0),i=1,nx),j=1,ny)
ELSE
READ (flunit,ERR=998) wetsfc_in
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) ((wetdp_in(i,j,0),i=1,nx),j=1,ny)
ELSE
READ (flunit,ERR=998) wetdp_in
END IF
ELSE
WRITE(6,'(a)') &
'Variable wetdp is not in the data set.'
END IF
IF ( wcanpin /= 0 ) THEN
WRITE (6, '(a)') 'Read in the canopy water amount data'
IF ( nstyp1 == 1 ) THEN
READ (flunit,ERR=998) ((wetcanp_in(i,j,0),i=1,nx),j=1,ny)
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(i,j),i=1,nx),j=1,ny)
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
CLOSE ( flunit )
CALL retunit ( flunit )
ELSE
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block
CALL hdfrd3d(sd_id,"tsfc",nx,ny,nstyp1+1,tsfc_in,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
CALL hdfrd3d(sd_id,"tsoil",nx,ny,nstyp1+1,tsoil_in,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
CALL hdfrd3d(sd_id,"wetsfc",nx,ny,nstyp1+1,wetsfc_in,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
CALL hdfrd3d(sd_id,"wetdp",nx,ny,nstyp1+1,wetdp_in,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
CALL hdfrd3d(sd_id,"wetcanp",nx,ny,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",nx,ny,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",nx,ny,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 ...
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt begin block:
! nstypin = MIN(nstyp1, nstyps)
!
! IF (tsfcin == 1) tsfc(:,:,0:nstypin) = tsfc_in(:,:,0:nstypin)
! IF (tsoilin == 1) tsoil(:,:,0:nstypin) = tsoil_in(:,:,0:nstypin)
! IF (wsfcin == 1) wetsfc(:,:,0:nstypin) = wetsfc_in(:,:,0:nstypin)
! IF (wdpin == 1) wetdp(:,:,0:nstypin) = wetdp_in(:,:,0:nstypin)
! IF (wcanpin == 1) wetcanp(:,:,0:nstypin) = wetcanp_in(:,:,0:nstypin)
!
! CALL fix_soil_nstyp(nx,ny,nstyp1,nstyp,tsfc,tsoil, &
! wetsfc,wetdp,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 (tsfcin == 1) tsfc(:,:,0:nstypin) = tsfc_in(:,:,0:nstypin)
IF (tsoilin == 1) tsoil(:,:,0:nstypin) = tsoil_in(:,:,0:nstypin)
IF (wsfcin == 1) wetsfc(:,:,0:nstypin) = wetsfc_in(:,:,0:nstypin)
IF (wdpin == 1) wetdp(:,:,0:nstypin) = wetdp_in(:,:,0:nstypin)
IF (wcanpin == 1) wetcanp(:,:,0:nstypin) = wetcanp_in(:,:,0:nstypin)
CALL fix_soil_nstyp(nx,ny,nstyp1,nstyp,tsfc,tsoil, &
wetsfc,wetdp,wetcanp)
ELSE
IF (nstyp1 == 1) THEN
tsfc_in (:,:,1) = tsfc_in (:,:,0)
tsoil_in (:,:,1) = tsoil_in (:,:,0)
wetsfc_in (:,:,1) = wetsfc_in (:,:,0)
wetdp_in (:,:,1) = wetdp_in (:,:,0)
wetcanp_in(:,:,1) = wetcanp_in(:,:,0)
ENDIF
CALL remap_soil_vars(nx,ny,nstyp1,nstyp, &
tsfc_in,tsoil_in,wetsfc_in,wetdp_in,wetcanp_in,soiltyp_in, &
tsfcin,tsoilin,wsfcin,wdpin,wcanpin, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,soiltyp)
ENDIF
!wdt end block
IF (mp_opt > 0) soilfile(1:128) = savename(1:128)
DEALLOCATE (tsfc_in,stat=istat)
DEALLOCATE (tsoil_in,stat=istat)
DEALLOCATE (wetsfc_in,stat=istat)
DEALLOCATE (wetdp_in,stat=istat)
DEALLOCATE (wetcanp_in,stat=istat)
DEALLOCATE (soiltyp_in,stat=istat)
! 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,29
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
INTEGER :: itmp1,itmp2
REAL :: rtmp1,rtmp2,rtmp3,rtmp4,rtmp5,rtmp6,rtmp7
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...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! 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
!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 ...
END IF
nstyp1 = MAX( nstyp1, 1 )
!-----------------------------------------------------------------------
!
! Check the data file for consistent grid parameters.
!
!-----------------------------------------------------------------------
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
!-----------------------------------------------------------------------
!
! 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
!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 ...
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
INTEGER :: itmp1,itmp2
REAL :: rtmp1,rtmp2,rtmp3,rtmp4,rtmp5,rtmp6,rtmp7
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/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,'(a6,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, soiloutfl, & 5,7
tsfcout,tsoilout,wsfcout,wdpout,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)
!
! 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=*) :: soiloutfl ! Surface data file name
INTEGER :: tsfcout,tsoilout,wsfcout,wdpout,wcanpout,snowdout
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) :: soilctlfl
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
!
!-----------------------------------------------------------------------
!
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 (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'
DO i=1,varnumax
varlev(i) = 0
varparam(i) = '99'
END DO
varnum = 0
IF ( tsfcout == 1 ) THEN
IF ( nstyp <= 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'ts'
vartit(varnum) = 'Surface skin temperature (K)'
ELSE
DO is=0,nstyp
WRITE (chr1,'(i1)') is
varnum = varnum + 1
varnam(varnum) = 'ts'//chr1
vartit(varnum) = 'Surface skin temperature (K) ' &
//'for soil type '//chr1
END DO
END IF
END IF
IF ( tsoilout == 1 ) THEN
IF ( nstyp <= 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 't2'
vartit(varnum) = 'Deep soil temperature (K)'
ELSE
DO is=0,nstyp
WRITE (chr1,'(i1)') is
varnum = varnum + 1
varnam(varnum) = 't2'//chr1
vartit(varnum) = 'Deep soil temperature (K)' &
//' for soil type '//chr1
END DO
END IF
END IF
IF ( wsfcout == 1 ) THEN
IF ( nstyp <= 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'wg'
vartit(varnum) = 'Surface skin soil moisture (m**3/m**3)'
ELSE
DO is=0,nstyp
WRITE (chr1,'(i1)') is
varnum = varnum + 1
varnam(varnum) = 'wg'//chr1
vartit(varnum) = 'Surface skin soil moisture (m**3/m**3)' &
//' for soil type '//chr1
END DO
END IF
END IF
IF ( wdpout == 1 ) THEN
IF ( nstyp <= 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'w2'
vartit(varnum) = 'Deep soil moisture (m**3/m**3)'
ELSE
DO is=0,nstyp
WRITE (chr1,'(i1)') is
varnum = varnum + 1
varnam(varnum) = 'w2'//chr1
vartit(varnum) = 'Deep soil moisture (m**3/m**3) ' &
//' for soil type '//chr1
END DO
END IF
END IF
IF ( wcanpout == 1 ) THEN
IF ( nstyp <= 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'wr'
vartit(varnum) = 'Canopy water amount (m)'
ELSE
DO is=0,nstyp
WRITE (chr1,'(i1)') is
varnum = varnum + 1
varnam(varnum) = 'wr'//chr1
vartit(varnum) = 'Canopy water amount (m) ' &
//' for soil type '//chr1
END DO
END IF
END IF
IF ( snowdout == 1 ) THEN
varnum = varnum + 1
varnam(varnum) = 'snowd'
vartit(varnum) = 'Snow depth (m)'
END IF
WRITE (nchout0,'(a/a)') &
'TITLE ARPS soil variable data control for ' &
//runname(1:lfnkey),'*'
WRITE (nchout0,'(a,a)') &
'DSET ',soiloutfl(1:fnlen)
WRITE (nchout0,'(a)') &
'OPTIONS sequential'
WRITE (nchout0,'(a,i10)') &
'FILEHEADER ',192 !!! 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,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,'(a6,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, temxy1,temxy2) 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