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