!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKSOILVAR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mksoilvar(nx,ny,nz,nstyps, & 1,24
mxdays,mxstns, &
xs,ys, &
hterain,zsc, &
ptbar,ptprt, &
pbar,pprt, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
soiltyp,stypfrct, &
precfile,staid, &
iptstn,jptstn,iptapi0,jptapi0, &
api0,obsprec,difprec,totpreca,prec, &
itema,initapi,api1,api2,kk2dep, &
totwt,dprec,tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Set the surface variables for ARPS model.
! Includes calculation of API from precip data to initialize
! soil moisture.
!
!-----------------------------------------------------------------------
!
! AUTHOR: John Mewes (precip and API) and Keith Brewster (temps)
! March, 1997
!
! MODIFICATION HISTORY:
!
! 04/11/1997 (Keith Brewster)
! Added processing of NCEP real-time precip data to API calculation.
! at the same time delete the table data array veg(14).
!
! 04/15/1997 (Keith Brewster)
! Added option of reading in pre-calculated soil moisture data,
! since the real-time precip data are typically only once per day.
! Minor reorganization of parameters. Bells and whistles added.
! Renamed to version 2.0
!
! 01/05/1998 (Donghai Wang)
! Fixed a problem in the case of k1< 0.0.
!
! 12/09/1998 (Donghai Wang)
! Added the snow cover.
!
! 2000/01/03 (Gene Bassett)
! Renamed mstinit to soilinit2 and added capability to update soil
! data over water and snow cover with data in a soil data file.
!
!-----------------------------------------------------------------------
!
! Set various parameters used in API calculation
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations....
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz,nstyps
INTEGER :: mxdays,mxstns
!
REAL :: xs (nx) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: ys (ny) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal)
! humidity (kg/kg)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal)
INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type
REAL :: stypfrct (nx,ny,nstyps)
INTEGER :: soilinit2 ! Moisture initialization option
INTEGER :: prdata ! Precip data option
CHARACTER (LEN=80) :: apifile ! Name of file containing names of the
! precip data files.
CHARACTER (LEN=80) :: apiinit ! File containing initial API's and
! their corresponding locations
CHARACTER (LEN=80) :: prcpdir ! Disk directory containing NCEP precip files
CHARACTER (LEN=80) :: prcplst ! File containing list of NCEP precip stns
CHARACTER (LEN=10) :: inidate ! Day for start of precip files
CHARACTER (LEN=10) :: fnldate ! Day for stop of precip files
REAL :: respapi ! Desired response of the initial API
! analysis to wavelengths equal to
! 2 x (mean initAPI stn sep).
INTEGER :: ndays ! Number of days between initial date
! and history file date
REAL :: k1,k2 ! API depletion minimum coefficients for
! the ground surface layer and the root
! zone layer
REAL :: range ! Range parameter in Barnes' wt. fcn.
! denominator
REAL :: respprec ! Desired response of the final daily
! precip analyses at wavelengths equal
! to 2 x (mean stn sep).
REAL :: gamma(2) ! Range multipliers for 1st and 2nd
! passes
!
!-----------------------------------------------------------------------
!
! API generation arrays:
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: precfile(mxstns)
CHARACTER (LEN=7) :: staid(mxstns)
INTEGER :: iptstn(mxstns)
INTEGER :: jptstn(mxstns)
INTEGER :: iptapi0(mxstns)
INTEGER :: jptapi0(mxstns)
REAL :: api0(mxstns)
REAL :: obsprec(mxdays,mxstns)
REAL :: difprec(mxdays,mxstns)
REAL :: totpreca(mxstns)
REAL :: prec(nx,ny)
!
!-----------------------------------------------------------------------
!
! Computed variables
!
!-----------------------------------------------------------------------
!
REAL :: tsfc (nx,ny,0:nstyps) ! Temperature at surface (K)
REAL :: tsoil (nx,ny,0:nstyps) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture
REAL :: wetdp (nx,ny,0:nstyps) ! Deep soil moisture
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL :: hterain(nx,ny)
REAL :: zsc(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
! Work Arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: itema(nx,ny)
REAL :: initapi(nx,ny,nstyps)
REAL :: api1(nx,ny,nstyps)
REAL :: api2(nx,ny,nstyps)
REAL :: kk2dep(nx,ny)
REAL :: totwt(nx,ny)
REAL :: dprec(nx,ny)
REAL :: tem1(nx,ny,nz)
! temparary arrays for soilinit2 option 5:
REAL, ALLOCATABLE :: tsfc2 (:,:,:)
REAL, ALLOCATABLE :: tsoil2 (:,:,:)
REAL, ALLOCATABLE :: wetsfc2 (:,:,:)
REAL, ALLOCATABLE :: wetdp2 (:,:,:)
REAL, ALLOCATABLE :: wetcanp2(:,:,:)
INTEGER, ALLOCATABLE :: soiltyp2(:,:,:)
!
!-----------------------------------------------------------------------
!
! API parameters
!
!-----------------------------------------------------------------------
!
COMMON /apicom1/ apiinit,apifile,prcpdir,prcplst,inidate,fnldate
COMMON /apicom2/ soilinit2,prdata
COMMON /apicom3/ respapi,k1,k2,range,respprec,gamma
!
!-----------------------------------------------------------------------
!
! Map projection variables
!
!-----------------------------------------------------------------------
!
REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
INCLUDE 'indtflg.inc'
!
!-----------------------------------------------------------------------
!
! Miscellaneous local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: sec24hr
PARAMETER( sec24hr = 86400 )
!
INTEGER :: i,j,ii,ij,is,is2,istn,jstn,iday,ix,iy,ireturn
INTEGER :: ipt,jpt,out
INTEGER :: tstns,totapi0,minday,maxday,length
INTEGER :: iniyr,inimo,inidy,initi,nowti,jlday
INTEGER :: itime,iyr,imon,idy,ihr,imin,isec
INTEGER :: idel,jdel,istr,iend,jstr,jend
INTEGER :: i1,i2,i3,i4,i5,j1,j2,j3,j4,j5
!
REAL :: tmp2,tmp3,prs2,prs3,w2,w3,ptk
REAL :: p0inv,tmpk,pres
REAL :: xpt,ypt,xctr,yctr,xll,yll
REAL :: lat,lon,ddx,ddy
REAL :: wmin,wmax,kk1,kk2,adj1,adj2
REAL :: apirng,apir2,r2,dist,const,thrdst
REAL :: sumwt,sumapi,wgt,precip,maxprec
REAL :: bias,rms,wsfcscl,wdpscl
REAL :: mindist,totdist,total
REAL :: d_0,d_1,d_1star,depth1
!
INTEGER :: tsfcout,tsoilout,wsfcout,wdpout,wcanpout,snowdout
INTEGER :: tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin
COMMON /intgcom/tsfcout,tsoilout,wsfcout,wdpout,wcanpout,snowdout
REAL :: twater
COMMON /realcom/twater
REAL :: amax, amin
INTEGER :: nbwater,nbw_max
PARAMETER (nbw_max=128)
REAL :: tbwater(nbw_max),blat1(nbw_max),blat2(nbw_max)
REAL :: blon1(nbw_max),blon2(nbw_max)
COMMON /bwcomi/ nbwater
COMMON /bwcomr/ tbwater,blat1,blat2,blon1,blon2
INTEGER :: ib
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF (soilinit2 == 5) THEN
ALLOCATE(tsfc2 (nx,ny,0:nstyps))
ALLOCATE(tsoil2 (nx,ny,0:nstyps))
ALLOCATE(wetsfc2 (nx,ny,0:nstyps))
ALLOCATE(wetdp2 (nx,ny,0:nstyps))
ALLOCATE(wetcanp2(nx,ny,0:nstyps))
ALLOCATE(soiltyp2(nx,ny,nstyps))
END IF
READ(inidate,'(i4,1x,i2,1x,i2)') iniyr,inimo,inidy
CALL ctim2abss
(iniyr,inimo,inidy,12,0,0, initi)
READ(fnldate,'(i4,1x,i2,1x,i2)') year,month,day
hour = 12
minute = 0
second = 0
curtim = 0
CALL ctim2abss(year, month, day, hour, minute, second, nowti)
nowti=nowti+curtim
ndays=((nowti-initi)/sec24hr)+1
WRITE(6,'(a,i6)') ' Number of days of data for API : ',ndays
IF(ndays > mxdays.AND.(soilinit2 == 2 .OR. soilinit2 == 3)) THEN
WRITE(6,'(a,i6,/,a)') &
' Number of days of data for API exceeds mxdays dimension: ', &
mxdays,' Increase mxdays or change inidate.'
STOP
END IF
!
!-----------------------------------------------------------------------
!
! Find the center lat,lon and the location (in map coord.) of the
! lower left corner (xll,yll).
!
!-----------------------------------------------------------------------
!
latnot(1)=trulat1
latnot(2)=trulat2
CALL setmapr(mapproj,sclfct,latnot,trulon)
CALL lltoxy(1,1,ctrlat,ctrlon,xctr,yctr)
WRITE(6,*) ' dx= ',dx,' dy= ',dy
xll=xctr-(0.5*(nx-3)*dx)
yll=yctr-(0.5*(ny-3)*dy)
!
WRITE(6,'(2x,a,f10.2,a,f10.2//)') &
' Grid center: lat= ',ctrlat,' lon= ',ctrlon
!
IF( sfcdat == 1 ) THEN
DO j=1,ny
DO i=1,nx
soiltyp(i,j,1) = styp
stypfrct(i,j,1) = 1.0
END DO
END DO
DO is=2,nstyps
DO j=1,ny
DO i=1,nx
soiltyp (i,j,is) = 0
stypfrct(i,j,is) = 0.0
END DO
END DO
END DO
ELSE IF (sfcdat == 2 .OR. (sfcdat == 3.AND.landin /= 1) ) THEN
DO is=1,nstyps
DO j=1,ny
DO i=1,nx
soiltyp (i,j,is) = 0
stypfrct(i,j,is) = 0.0
END DO
END DO
END DO
CALL readsfcdt(nx,ny,nstyps,sfcdtfl,dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
soiltyp,stypfrct,itema,tem1,tem1,tem1,tem1)
WRITE(6,*) ' nstyp = ',nstyp
IF( nstyp == 1 ) THEN
DO j=1,ny-1
DO i=1,nx-1
stypfrct(i,j,1) = 1.0
END DO
END DO
END IF
ELSE IF (sfcdat == 3 .AND. landin == 1) THEN
WRITE(6,'(1x,a/)') &
'Surface property data in the history file was used.'
ELSE
WRITE(6,'(1x,a,i3,a/)') &
'Invalid surface data input option. sfcdat =',sfcdat, &
'. Program stopped in MKSOILVAR.'
STOP
END IF
DO is=0,nstyp
PRINT*,'In MKSOILVAR for is =', is
CALL a3dmax0(tsoil(1,1,is), &
1,nx,1,nx-1,1,ny,1,ny-1,1,1 ,1,1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'tsoilmin= ', amin,', tsoilmax =',amax
CALL a3dmax0(tsfc(1,1,is), &
1,nx,1,nx-1,1,ny,1,ny-1,1,1 ,1,1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'tsfcmin= ', amin,', tsfcmax =',amax
END DO
!
IF(soilinit2 /= 5) THEN
p0inv=1./p0
DO j=1,ny-1
DO i=1,nx-1
tmp3 = ptbar(i,j,3) + ptprt(i,j,3)
prs3 = ALOG(pbar (i,j,3) + pprt (i,j,3))
tmp2 = ptbar(i,j,2) + ptprt(i,j,2)
prs2 = ALOG(pbar (i,j,2) + pprt (i,j,2))
w2=1.+(zsc(i,j,2)-hterain(i,j))/(zsc(i,j,3)-zsc(i,j,2))
w3=1.-w2
ptk=w2*tmp2 + w3*tmp3
pres=EXP((w2*prs2 + w3*prs3))
tmpk = ptk * (p0inv*pres)**rddcp
IF (nbwater > 0) CALL xytoll(1,1,xll+dx*(i+0.5),yll+dy*(j+0.5),lat,lon)
!
DO is=1,nstyp
IF(soiltyp(i,j,is) /= 0) THEN
IF( soiltyp(i,j,is) == 13 ) THEN
tsfc (i,j,is) = twater + 273.15
tsoil (i,j,is) = twater + 273.15
DO ib = 1,nbwater
IF ((lat >= blat1(ib)) .AND. (lat <= blat2(ib)) &
.AND. (lon >= blon1(ib)) .AND. (lon <= blon2(ib))) THEN
tsfc (i,j,is) = tbwater(ib) + 273.15
tsoil (i,j,is) = tbwater(ib) + 273.15
END IF
END DO
ELSE
tsfc (i,j,is) = tmpk + tsprt
tsoil (i,j,is) = tmpk + t2prt
END IF
END IF
END DO
END DO
END DO
END IF
IF ( snowin /= 1 ) THEN
DO j=1,ny
DO i=1,nx
snowdpth(i,j)=snowdpth0
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Modify or initialize soil properties using soilinit2 option
!
!-----------------------------------------------------------------------
!
IF( soilinit2 == 1 ) THEN
!
DO is=1,nstyp
DO j=1,ny-1
DO i=1,nx-1
!
IF( soiltyp(i,j,is) /= 0 ) THEN
wetsfc (i,j,is) = wwlt(soiltyp(i,j,is)) + wgrat &
*( wsat(soiltyp(i,j,is)) - wwlt(soiltyp(i,j,is)) )
wetdp (i,j,is) = wwlt(soiltyp(i,j,is)) + w2rat &
*( wsat(soiltyp(i,j,is)) - wwlt(soiltyp(i,j,is)) )
wetcanp(i,j,is) = wetcanp0
END IF
END DO
END DO
END DO
!
ELSE IF( soilinit2 == 2 .OR. soilinit2 == 3 ) THEN
!
!-----------------------------------------------------------------------
!
! Read the initial API data and find its location in the grid
!
!-----------------------------------------------------------------------
!
totapi0=0
OPEN(UNIT=3,ERR=205,FILE=apiinit,STATUS='unknown')
!
out=0
WRITE(6,*)
DO istn=1,mxstns
jstn=istn-out
READ(3,*,END=201) lat,lon,api0(jstn)
CALL lltoxy(1,1,lat,lon,xpt,ypt)
xpt=xpt-xll
ypt=ypt-yll
CALL findlc(nx,ny,xs,ys,xpt,ypt,ipt,jpt,ireturn)
IF( ireturn == 0 ) THEN
iptapi0(jstn)=ipt
jptapi0(jstn)=jpt
WRITE(6,*) 'API obs:',jstn,' found near ',ipt,jpt, &
' with initial API of ',api0(jstn)
WRITE(6,*) ' --> main soil type:',soiltyp(ipt,jpt,1), &
' wilting:',wwlt(soiltyp(ipt,jpt,1)), &
' saturation:',wsat(soiltyp(ipt,jpt,1))
WRITE(6,*)
ELSE
WRITE(6,*)'Initial API station #',jstn, &
' is outside the grid..'
WRITE(6,*)
out=out+1
END IF
END DO
201 CONTINUE
totapi0=istn-1-out
CLOSE(3)
205 CONTINUE
!
!-----------------------------------------------------------------------
!
! Read the observed precip data and find its location in the grid
!
!-----------------------------------------------------------------------
!
IF( prdata == 1 ) THEN
out=0
WRITE(6,*)
OPEN(UNIT=1,ERR=301,FILE=apifile,STATUS='old')
DO istn=1,mxstns
jstn=istn-out
READ(1,1010,END=301) precfile(jstn)
1010 FORMAT(a80)
length=LEN(precfile(jstn))
CALL strlnth(precfile(jstn),length)
WRITE(6,*) 'File:',precfile(jstn)(1:length)
OPEN(UNIT=2,FILE=precfile(jstn),STATUS='old')
READ(2,*) lat,lon
CALL lltoxy(1,1,lat,lon,xpt,ypt)
xpt=xpt-xll
ypt=ypt-yll
CALL findlc(nx,ny,xs,ys,xpt,ypt,ipt,jpt,ireturn)
IF( ireturn == 0 ) THEN
!
iptstn(jstn)=ipt
jptstn(jstn)=jpt
totpreca(jstn)=0.0
!
WRITE(6,*)
WRITE(6,*) 'Precip stn # ',jstn,' was found near i,j:', &
ipt,jpt
IF( jpt >= nint(ny*0.5) ) THEN
WRITE(6,*) 'North of Center'
ELSE
WRITE(6,*) 'South of Center'
END IF
IF( ipt >= nint(nx*0.5)) THEN
WRITE(6,*) 'East of Center'
ELSE
WRITE(6,*) 'West of Center'
END IF
!
WRITE(6,*)
WRITE(6,*)
DO iday=1,ndays
READ(2,*) obsprec(iday,jstn) ! in inches...
obsprec(iday,jstn)=obsprec(iday,jstn)*2.54
!
! Assume (arbitrarily) that the ground can only absorb 4.0 cm of
! rainfall per day, with the rest being runoff.
!
obsprec(iday,jstn)=AMIN1(obsprec(iday,jstn),4.0 )
END DO
WRITE(6,*)
ELSE
PRINT *, &
'Precip station #',jstn,' is outside of the grid..'
PRINT *
out=out+1
END IF
CLOSE(2)
END DO
301 CONTINUE
CLOSE(1)
CLOSE(2)
tstns=istn-1-out
ELSE IF (prdata == 2) THEN
CALL rdnceppr(nx,ny,mxstns,mxdays,xs,ys,xll,yll, &
year,month,day,ndays, &
prcpdir,prcplst, &
staid,iptstn,jptstn,obsprec,totpreca, &
tstns)
ELSE
WRITE(6,'(i6,a)') prdata,' is not a valid prdata option!'
STOP
END IF
!
!---------------------------------------------------------------------
!
! Analyze the initial root layer API to the grid. Determine the
! range parameter using respapi.
!
!---------------------------------------------------------------------
!
i1=1
i2=nint(nx*.25)
i3=nint(nx*.50)
i4=nint(nx*.75)
i5=nx-1
!
j1=1
j2=nint(ny*.25)
j3=nint(ny*.50)
j4=nint(ny*.75)
j5=ny-1
!
WRITE(6,*)
WRITE(6,*) 'Analyzing initial API to grid...'
WRITE(6,*)
!
IF( totapi0 > 1 ) THEN
!
totdist=SQRT((nx*dx*0.001)*(ny*dy*0.001))* &
(1.+SQRT(FLOAT(totapi0)))/FLOAT(totapi0-1) ! using separation that
! would result if obs were randomly scattered
apirng= &
SQRT(-((2.*totdist/3.14159)**2)*LOG(respapi))
apir2=apirng*apirng
const=1.0E-06/apir2
!
WRITE(6,*)
WRITE(6,*) 'Using range:',apirng,' km, for initializing API'
WRITE(6,*)
DO iy=1,ny-1
DO ix=1,nx-1
sumwt=0.
sumapi=0.
DO istn=1,totapi0
ddx=xs(iptapi0(istn))-xs(ix)
ddy=ys(jptapi0(istn))-ys(iy)
dist=(ddx*ddx+ddy*ddy)
wgt=EXP(-dist*const)
sumwt=sumwt+wgt
sumapi=sumapi+wgt*api0(istn)
END DO
IF( sumwt > 0. ) sumapi=sumapi/sumwt
DO is=1,nstyp
IF( soiltyp(ix,iy,is) /= 0 ) initapi(ix,iy,is)= &
AMIN1(sumapi,(wsat(soiltyp(ix,iy,is))*d2*100.))
END DO
END DO
END DO
ELSE ! if no initial values available, set to 0.5*wsat
DO is=1,nstyp
DO iy=1,ny-1
DO ix=1,nx-1
IF( soiltyp(ix,iy,is) /= 0) &
initapi(ix,iy,is)=0.5*wsat(soiltyp(ix,iy,is))*d2*100.0
END DO
END DO
END DO
END IF
!
WRITE(6,*)
WRITE(6,*) 'Initial API analysis:'
WRITE(6,*) '---------------------'
WRITE(6,*)
!
WRITE(6,1040) j5,initapi(i1,j5,1),initapi(i2,j5,1), &
initapi(i3,j5,1),initapi(i4,j5,1), &
initapi(i5,j5,1)
WRITE(6,1040) j4,initapi(i1,j4,1),initapi(i2,j4,1), &
initapi(i3,j4,1),initapi(i4,j4,1), &
initapi(i5,j4,1)
WRITE(6,1040) j3,initapi(i1,j3,1),initapi(i2,j3,1), &
initapi(i3,j3,1),initapi(i4,j3,1), &
initapi(i5,j3,1)
WRITE(6,1040) j2,initapi(i1,j2,1),initapi(i2,j2,1), &
initapi(i3,j2,1),initapi(i4,j2,1), &
initapi(i5,j2,1)
WRITE(6,1040) j1,initapi(i1,j1,1),initapi(i2,j1,1), &
initapi(i3,j1,1),initapi(i4,j1,1), &
initapi(i5,j1,1)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
minday=15
depth1=0.01
!
DO iy=1,ny-1
DO ix=1,nx-1
kk2dep(ix,iy)=1.0
END DO
END DO
DO is=1,nstyp
DO iy=1,ny-1
DO ix=1,nx-1
api1(ix,iy,is)=initapi(ix,iy,is)*depth1/d2
api2(ix,iy,is)=initapi(ix,iy,is)
END DO
END DO
END DO
!
maxprec=0.0
DO iday=1,ndays
precip=0.0
DO istn=1,tstns
IF( obsprec(iday,istn) >= 0. ) precip=precip+obsprec(iday,istn)
END DO
precip=precip/tstns
IF( precip > maxprec ) THEN
maxprec=precip
maxday=iday
END IF
END DO
!
!----------------------------------------------------------------------
!
! If the range is not specified in the input file (i.e. = 0.0)
! determine the own range based on the input desired response.
!
!----------------------------------------------------------------------
!
totdist=0.0
DO ii=1,tstns
mindist=1.0E30
DO ij=1,tstns
IF( ii /= ij ) THEN
ddx=(xs(iptstn(ii))-xs(iptstn(ij)))
ddy=(ys(jptstn(ii))-ys(jptstn(ij)))
mindist=AMIN1(mindist,(ddx*ddx+ddy*ddy))
END IF
END DO
totdist=totdist+0.001*SQRT(mindist)
END DO
!
WRITE(6,*) 'Average station separation (km):',totdist/tstns
WRITE(6,*)
!
IF( range == 0. ) THEN
DO ii=1,1000000
d_0=.000001*ii
d_1=d_0**(gamma(2)**1.0)
d_1star=d_0+(1-d_0)*d_1
IF( d_1star >= respprec ) EXIT
END DO
396 CONTINUE
!
range=SQRT(-((2*(totdist/(tstns*1.0))/3.14159)**2)*LOG(d_0))
WRITE(6,*)
WRITE(6,*) 'Range not specified, using range:',range
WRITE(6,1026) ' -- 1st pass response at 2 x (stn sep):',d_0
WRITE(6,1026) ' -- 2nd pass response at 2 x (stn sep):', &
d_1star
1026 FORMAT(a40,f10.8)
WRITE(6,*)
END IF
!
r2=range*range
!
!----------------------------------------------------------------------
!
! Begin main day loop
! Objectively analyze the precipitation data on 'ndays' different
! days using a 2-pass analysis.
!
!----------------------------------------------------------------------
!
DO iday=1,ndays
!
DO iy=1,ny-1
DO ix=1,nx-1
prec(ix,iy)=0.
totwt(ix,iy)=0.
END DO
END DO
itime=initi+(iday-1)*sec24hr
CALL abss2ctim( itime,iyr,imon,idy,ihr,imin,isec )
!
! 1st pass....
!
const=1.0E-06/(r2*gamma(1))
thrdst=1.0E06*(6.*totdist/tstns)*(6.*totdist/tstns)
idel=INT(SQRT(thrdst)/dx)+1
jdel=INT(SQRT(thrdst)/dy)+1
!
DO istn=1,tstns
IF( obsprec(iday,istn) >= 0. ) THEN
istr =MAX(1,(iptstn(istn)-idel))
iend =MIN((nx-1),(iptstn(istn)+idel))
jstr =MAX(1,(jptstn(istn)-jdel))
jend =MIN((ny-1),(jptstn(istn)+jdel))
DO iy=jstr,jend
DO ix=istr,iend
ddx=(xs(iptstn(istn))-xs(ix))
ddy=(ys(jptstn(istn))-ys(iy))
dist=(ddx*ddx+ddy*ddy)
IF( dist < thrdst ) THEN
wgt=EXP(-dist*const)
prec(ix,iy)=prec(ix,iy)+wgt*obsprec(iday,istn)
totwt(ix,iy)=totwt(ix,iy)+wgt
END IF
END DO
END DO
END IF
END DO
DO iy=1,ny-1
DO ix=1,nx-1
IF( totwt(ix,iy) > 0. ) prec(ix,iy)=prec(ix,iy)/totwt(ix,iy)
END DO
END DO
!
IF( iday == maxday ) THEN
!
WRITE(6,'(a,i4,a,i2.2,a,i2.2)') &
'Example precip analysis on ', &
iyr,'-',imon,'-',idy
WRITE(6,'(a)') '-------------------------------------'
WRITE(6,*)
!
1035 FORMAT(7X,i3,5X,i3,5X,i3,5X,i3,5X,i3)
1040 FORMAT(i3,3X,f5.2,3X,f5.2,3X,f5.2,3X,f5.2,3X,f5.2)
WRITE(6,1040) j5,prec(i1,j5),prec(i2,j5), &
prec(i3,j5),prec(i4,j5), &
prec(i5,j5)
WRITE(6,1040) j4,prec(i1,j4),prec(i2,j4), &
prec(i3,j4),prec(i4,j4), &
prec(i5,j4)
WRITE(6,1040) j3,prec(i1,j3),prec(i2,j3), &
prec(i3,j3),prec(i4,j3), &
prec(i5,j3)
WRITE(6,1040) j2,prec(i1,j2),prec(i2,j2), &
prec(i3,j2),prec(i4,j2), &
prec(i5,j2)
WRITE(6,1040) j1,prec(i1,j1),prec(i2,j1), &
prec(i3,j1),prec(i4,j1), &
prec(i5,j1)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
rms=0.0
bias=0.0
total=0.0
DO istn=1,tstns
IF( obsprec(iday,istn) >= 0. ) THEN
1045 FORMAT(a4,1X,f6.3,1X,a6,1X,f6.3)
WRITE(6,1045) 'obs:',obsprec(iday,istn),' anal:', &
prec(iptstn(istn),jptstn(istn))
rms=rms+(obsprec(iday,istn)- &
prec(iptstn(istn),jptstn(istn)))**2
bias=bias+obsprec(iday,istn)- &
prec(iptstn(istn),jptstn(istn))
total=total+1.
END IF
END DO
IF( total > 0. ) THEN
rms=SQRT(rms/total)
bias=-bias/total
WRITE(6,*)
WRITE(6,*) 'RMS Error:',rms,' Bias:',bias
WRITE(6,*)
END IF
END IF
!
! 2nd pass....
!
DO istn=1,tstns
IF( obsprec(iday,istn) >= 0. ) THEN
difprec(iday,istn)=obsprec(iday,istn)- &
prec(iptstn(istn),jptstn(istn))
ELSE
difprec(iday,istn)=-900.
END IF
END DO
DO iy=1,ny-1
DO ix=1,nx-1
dprec(ix,iy)=0.
totwt(ix,iy)=0.
END DO
END DO
!
const=1.0E-06/(r2*gamma(2))
thrdst=1.0E06*(2.5*totdist/tstns)*(2.5*totdist/tstns)
idel=INT(SQRT(thrdst)/dx)+1
jdel=INT(SQRT(thrdst)/dy)+1
DO istn=1,tstns
IF( difprec(iday,istn) > -100. ) THEN
istr =MAX(1,(iptstn(istn)-idel))
iend =MIN((nx-1),(iptstn(istn)+idel))
jstr =MAX(1,(jptstn(istn)-jdel))
jend =MIN((ny-1),(jptstn(istn)+jdel))
DO iy=jstr,jend
DO ix=istr,iend
ddx=xs(iptstn(istn))-xs(ix)
ddy=ys(jptstn(istn))-ys(iy)
dist=(ddx*ddx+ddy*ddy)
IF( dist < thrdst ) THEN
wgt=EXP(-dist*const)
dprec(ix,iy)=dprec(ix,iy)+wgt*difprec(iday,istn)
totwt(ix,iy)=totwt(ix,iy)+wgt
END IF
END DO
END DO
END IF
END DO
!
!
! Sum difference and impose limits on daily precip values...
!
DO iy=1,ny-1
DO ix=1,nx-1
IF( totwt(ix,iy) > 0. ) &
prec(ix,iy)=prec(ix,iy)+(dprec(ix,iy)/totwt(ix,iy))
prec(ix,iy)=AMAX1(prec(ix,iy),0.0)
prec(ix,iy)=AMIN1(prec(ix,iy),4.0)
END DO
END DO
!
IF( iday == maxday) THEN
!
itime=initi+(iday-1)*sec24hr
CALL abss2ctim( itime,iyr,imon,idy,ihr,imin,isec )
WRITE(6,'(a,i4,a,i2.2,a,i2.2)') &
'Example precip analysis on ', &
iyr,'-',imon,'-',idy
WRITE(6,*) '-------------------------------------'
WRITE(6,*)
!
WRITE(6,1040) j5,prec(i1,j5),prec(i2,j5), &
prec(i3,j5),prec(i4,j5), &
prec(i5,j5)
WRITE(6,1040) j4,prec(i1,j4),prec(i2,j4), &
prec(i3,j4),prec(i4,j4), &
prec(i5,j4)
WRITE(6,1040) j3,prec(i1,j3),prec(i2,j3), &
prec(i3,j3),prec(i4,j3), &
prec(i5,j3)
WRITE(6,1040) j2,prec(i1,j2),prec(i2,j2), &
prec(i3,j2),prec(i4,j2), &
prec(i5,j2)
WRITE(6,1040) j1,prec(i1,j1),prec(i2,j1), &
prec(i3,j1),prec(i4,j1), &
prec(i5,j1)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
rms=0.0
bias=0.0
total=0.0
DO istn=1,tstns
IF(obsprec(iday,istn) >= 0.) THEN
WRITE(6,1045) 'obs:',obsprec(iday,istn),' anal:', &
prec(iptstn(istn),jptstn(istn))
rms=rms+(obsprec(iday,istn)- &
prec(iptstn(istn),jptstn(istn)))**2
bias=bias+obsprec(iday,istn)- &
prec(iptstn(istn),jptstn(istn))
total=total+1.
END IF
END DO
IF( total > 0. ) THEN
rms=SQRT(rms/total)
bias=-bias/total
WRITE(6,*)
WRITE(6,*) 'RMS Error:',rms,' Bias:',bias
WRITE(6,*)
END IF
END IF
!
!---------------------------------------------------------------------
!
! Calculate the API - assume that API will be relatively
! insensitive to the initial API after about 90 days.
! Recall d2=100.0cm and depth1 (m) is set below. Note that
! depth1 is not necessarilly equal to d1 in soilcst.inc, which
! is only a normalizing depth. In Noilhan-Planton (1989) they
! suggest that depth1 is only a few millimeters.
!
! k1 and k2 are minimum depletion coefficients for the topsoil
! and root layer. kk1 and kk2 are the time-dependent depletion
! coefficients where kk? is at a maximum on minday (for minimum
! moisture depletion - set to January 15, or Julian day 15) and
! kk? is at a minimum 6 months later (July 15)
!
! The time-dependent depletion coefficients can THEN be modified
! up or down according to the apparent volumetric soil water
! content (i.e. the API) relative to wilting and saturation using
! the adj1 and adj2 terms, the values of which can be derived
! empirically using ARM data. With the limited data available
! at this time, the adj1 and adj2 terms were simply set to 1.0
! due to inability to show value in adding this term. However,
! intuitively there should be some relationship.
!
! The API depletion coefficient multiplies the difference
! between the calculated API and the wilting point (or the
! initial API, whichever is lower) instead of the multiplying
! the prior API as in the common formulation. The new method
! was derived empirically and appeared to give better correlation
! with observed data as well as a lower dependence on the initial
! API using the depletion coefficient formula below.
!
!---------------------------------------------------------------------
!
!
adj1=1.0
adj2=1.0
CALL julday( iyr, imon, idy, jlday )
kk1=1.0-adj1* &
(1.0-k1)*SIN((3.14159*(jlday-minday)) &
/365.0)
kk2=1.0-adj2* &
(1.0-k2)*SIN((3.14159*(jlday-minday)) &
/365.0)
DO iy=1,ny-1
DO ix=1,nx-1
kk2dep(ix,iy)=kk2dep(ix,iy)*kk2
END DO
END DO
DO is=1,nstyp
DO iy=1,ny-1
DO ix=1,nx-1
!
IF(soiltyp(ix,iy,is) > 0) THEN
IF( k1 > 0. ) THEN
api1(ix,iy,is)=prec(ix,iy)+kk1*api1(ix,iy,is)
ELSE
api1(ix,iy,is)=0.0
END IF
api1(ix,iy,is)= &
AMIN1(api1(ix,iy,is), &
(wsat(soiltyp(ix,iy,is))*depth1*100.))
!
api2(ix,iy,is)=prec(ix,iy)+ &
kk2*(api2(ix,iy,is)-AMIN1(initapi(ix,iy,is), &
wwlt(soiltyp(ix,iy,is))*(d2*100.)))+ &
AMIN1(initapi(ix,iy,is), &
wwlt(soiltyp(ix,iy,is))*(d2*100.))
api2(ix,iy,is)= &
AMIN1(api2(ix,iy,is),(wsat(soiltyp(ix,iy,is))*d2*100.))
!
! Carry out 1/2 of the daily depletion on the final day (the
! day of the run) without adding in a precip term.
!
IF( iday == ndays) THEN
IF( k1 > 0. ) &
api1(ix,iy,is)=(api1(ix,iy,is)+api1(ix,iy,is)*kk1)/2.0
api2(ix,iy,is)=(api2(ix,iy,is)+api2(ix,iy,is)*kk2)/2.0
END IF
END IF
!
END DO
END DO
END DO
DO istn=1,tstns
! write(6,1020) 'Stn #:',istn,'Day:',iyr,'-',imon,'-',idy,
! : 'Obs. precip (cm):',obsprec(iday,istn),
! : 'Anal. prec (cm):',prec(iptstn(istn),jptstn(istn))
totpreca(istn)=totpreca(istn)+ &
prec(iptstn(istn),jptstn(istn))
END DO
1020 FORMAT(a6,i4,3X,a5,i4,a,i2.2,a,i2.2,3X,a17,f5.2,3X,a16,f5.2)
END DO
!
WRITE(6,*)
WRITE(6,*) 'Top layer API analysis:'
WRITE(6,*) '-----------------------'
WRITE(6,*)
!
WRITE(6,1040) j5,api1(i1,j5,1),api1(i2,j5,1), &
api1(i3,j5,1),api1(i4,j5,1), &
api1(i5,j5,1)
WRITE(6,1040) j4,api1(i1,j4,1),api1(i2,j4,1), &
api1(i3,j4,1),api1(i4,j4,1), &
api1(i5,j4,1)
WRITE(6,1040) j3,api1(i1,j3,1),api1(i2,j3,1), &
api1(i3,j3,1),api1(i4,j3,1), &
api1(i5,j3,1)
WRITE(6,1040) j2,api1(i1,j2,1),api1(i2,j2,1), &
api1(i3,j2,1),api1(i4,j2,1), &
api1(i5,j2,1)
WRITE(6,1040) j1,api1(i1,j1,1),api1(i2,j1,1), &
api1(i3,j1,1),api1(i4,j1,1), &
api1(i5,j1,1)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
WRITE(6,*)
WRITE(6,*) 'Root zone API analysis:'
WRITE(6,*) '-----------------------'
WRITE(6,*)
!
WRITE(6,1040) j5,api2(i1,j5,1),api2(i2,j5,1), &
api2(i3,j5,1),api2(i4,j5,1), &
api2(i5,j5,1)
WRITE(6,1040) j4,api2(i1,j4,1),api2(i2,j4,1), &
api2(i3,j4,1),api2(i4,j4,1), &
api2(i5,j4,1)
WRITE(6,1040) j3,api2(i1,j3,1),api2(i2,j3,1), &
api2(i3,j3,1),api2(i4,j3,1), &
api2(i5,j3,1)
WRITE(6,1040) j2,api2(i1,j2,1),api2(i2,j2,1), &
api2(i3,j2,1),api2(i4,j2,1), &
api2(i5,j2,1)
WRITE(6,1040) j1,api2(i1,j1,1),api2(i2,j1,1), &
api2(i3,j1,1),api2(i4,j1,1), &
api2(i5,j1,1)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
WRITE(6,*)
WRITE(6,*) 'Dependence of root zone API on initial API: '
WRITE(6,*) '(most useful when adj2 .ne. 1.0 -- see code)'
WRITE(6,*) '--------------------------------------------'
WRITE(6,*)
!
WRITE(6,1040) j5,kk2dep(i1,j5),kk2dep(i2,j5), &
kk2dep(i3,j5),kk2dep(i4,j5), &
kk2dep(i5,j5)
WRITE(6,1040) j4,kk2dep(i1,j4),kk2dep(i2,j4), &
kk2dep(i3,j4),kk2dep(i4,j4), &
kk2dep(i5,j4)
WRITE(6,1040) j3,kk2dep(i1,j3),kk2dep(i2,j3), &
kk2dep(i3,j3),kk2dep(i4,j3), &
kk2dep(i5,j3)
WRITE(6,1040) j2,kk2dep(i1,j2),kk2dep(i2,j2), &
kk2dep(i3,j2),kk2dep(i4,j2), &
kk2dep(i5,j2)
WRITE(6,1040) j1,kk2dep(i1,j1),kk2dep(i2,j1), &
kk2dep(i3,j1),kk2dep(i4,j1), &
kk2dep(i5,j1)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
!--------------------------------------------------------------------
!
! Convert the API's to wetsfc and wetdp. Limit both 'wetsfc' and
! 'wetdp' to ~'wsat' as a maximum. Limit 'wetsfc' to 'wgeq' as
! a minimum and limit 'wetdp' to ~'wwlt' as a minimum. If k1=0.0
! set 'wetsfc' equal to 'wgeq' and if k1<0.0 set 'wetsfc' ~ equal
! to 0.8*'wgeq'. See input file for more details.
!
! Also, convert API from scalar points back to x,y grid.
!
!--------------------------------------------------------------------
!
IF(soilinit2 == 3) THEN
wsfcscl=wgrat
wdpscl=w2rat
ELSE
wsfcscl=1.0
wdpscl=1.0
END IF
!
DO is=1,nstyp
DO iy=1,ny-1
DO ix=1,nx-1
IF( soiltyp(ix,iy,is) /= 0 ) THEN
wetsfc(ix,iy,is) =wsfcscl*api1(ix,iy,is)/(depth1*100.0)
wetdp(ix,iy,is) =wdpscl*api2(ix,iy,is)/(d2*100.0)
wetcanp(ix,iy,is)=wetcanp0
!
IF( soiltyp(ix,iy,is) /= 13 ) THEN
wmin=0.95*wwlt(soiltyp(ix,iy,is)) ! arbitrary min,max values
wmax=0.95*wsat(soiltyp(ix,iy,is)) ! for wetdp...
ELSE ! if soiltyp is water
wmin=wwlt(soiltyp(ix,iy,is))
wmax=wsat(soiltyp(ix,iy,is))
END IF
IF( wetdp(ix,iy,is) < wmin ) THEN
wetdp(ix,iy,is)=wmin ! keep wetdp close to wwlt
ELSE IF( wetdp(ix,iy,is) > wmax ) THEN
wetdp(ix,iy,is)=wmax ! keep wetdp slightly less than saturated
END IF
wmin=wetdp(ix,iy,is)/wsat(soiltyp(ix,iy,is))
wmin=wmin - ( awgeq(soiltyp(ix,iy,is)) * &
wmin**pwgeq(soiltyp(ix,iy,is)) ) * &
( 1.0 - wmin**(8.0*pwgeq(soiltyp(ix,iy,is))) )
wmin=wmin*wsat(soiltyp(ix,iy,is))
IF( wetsfc(ix,iy,is) > wmax ) THEN
wetsfc(ix,iy,is)=wsat(soiltyp(ix,iy,is))
ELSE IF( wetsfc(ix,iy,is) < wmin ) THEN
wetsfc(ix,iy,is)=wmin
END IF
!
IF( k1 == 0. ) THEN
wetsfc(ix,iy,is)=wmin
ELSE IF( k1 < 0. .AND. soiltyp(ix,iy,is) /= 13) THEN
wetsfc(ix,iy,is)=wmin-AMIN1(0.2*wmin,0.05)
END IF
END IF
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
tsfc (i,j,0) = 0.0
tsoil (i,j,0) = 0.0
wetsfc (i,j,0) = 0.0
wetdp (i,j,0) = 0.0
wetcanp(i,j,0) = 0.0
END DO
END DO
DO is = 1,nstyp
DO j=1,ny-1
DO i=1,nx-1
tsfc (i,j,0) = tsfc (i,j,0) &
+ tsfc (i,j,is) * stypfrct(i,j,is)
tsoil (i,j,0) = tsoil (i,j,0) &
+ tsoil (i,j,is) * stypfrct(i,j,is)
wetsfc (i,j,0) = wetsfc (i,j,0) &
+ wetsfc (i,j,is) * stypfrct(i,j,is)
wetdp (i,j,0) = wetdp (i,j,0) &
+ wetdp (i,j,is) * stypfrct(i,j,is)
wetcanp(i,j,0) = wetcanp(i,j,0) &
+ wetcanp(i,j,is) * stypfrct(i,j,is)
END DO
END DO
END DO
DO is=0,nstyp
CALL edgfill(tsfc (1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL edgfill(tsoil (1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL edgfill(wetsfc (1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL edgfill(wetdp (1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
CALL edgfill(wetcanp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
END DO
!
!-------------------------------------------------------------------
!
! Print out the calculated values of wetsfc, wetdp, and wetcanp.
!
!-------------------------------------------------------------------
!
i5=nx
j5=ny
!
WRITE(6,*)
WRITE(6,*) 'wetsfc analysis:'
WRITE(6,*) '-----------------------'
WRITE(6,*)
!
WRITE(6,1040) j5,wetsfc(i1,j5,0),wetsfc(i2,j5,0), &
wetsfc(i3,j5,0),wetsfc(i4,j5,0), &
wetsfc(i5,j5,0)
WRITE(6,1040) j4,wetsfc(i1,j4,0),wetsfc(i2,j4,0), &
wetsfc(i3,j4,0),wetsfc(i4,j4,0), &
wetsfc(i5,j4,0)
WRITE(6,1040) j3,wetsfc(i1,j3,0),wetsfc(i2,j3,0), &
wetsfc(i3,j3,0),wetsfc(i4,j3,0), &
wetsfc(i5,j3,0)
WRITE(6,1040) j2,wetsfc(i1,j2,0),wetsfc(i2,j2,0), &
wetsfc(i3,j2,0),wetsfc(i4,j2,0), &
wetsfc(i5,j2,0)
WRITE(6,1040) j1,wetsfc(i1,j1,0),wetsfc(i2,j1,0), &
wetsfc(i3,j1,0),wetsfc(i4,j1,0), &
wetsfc(i5,j1,0)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
WRITE(6,*)
WRITE(6,*) 'wetdp analysis:'
WRITE(6,*) '-----------------------'
WRITE(6,*)
!
WRITE(6,1040) j5,wetdp(i1,j5,0),wetdp(i2,j5,0), &
wetdp(i3,j5,0),wetdp(i4,j5,0), &
wetdp(i5,j5,0)
WRITE(6,1040) j4,wetdp(i1,j4,0),wetdp(i2,j4,0), &
wetdp(i3,j4,0),wetdp(i4,j4,0), &
wetdp(i5,j4,0)
WRITE(6,1040) j3,wetdp(i1,j3,0),wetdp(i2,j3,0), &
wetdp(i3,j3,0),wetdp(i4,j3,0), &
wetdp(i5,j3,0)
WRITE(6,1040) j2,wetdp(i1,j2,0),wetdp(i2,j2,0), &
wetdp(i3,j2,0),wetdp(i4,j2,0), &
wetdp(i5,j2,0)
WRITE(6,1040) j1,wetdp(i1,j1,0),wetdp(i2,j1,0), &
wetdp(i3,j1,0),wetdp(i4,j1,0), &
wetdp(i5,j1,0)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
! DO 720 istn=1,tstns
! totpreco=0.0
! DO 710 iday=1,ndays
! IF(obsprec(iday,istn).gt.0.)
! : totpreco=totpreco+obsprec(iday,istn)
! 710 CONTINUE
! write(6,*)
! write(6,1021) 'Total recorded precip:',totpreco
! write(6,1021) 'Total analyzed precip:',totpreca(istn)
!1021 format(a22,f6.2)
! write(6,1022) 'API:',api2(iptstn(istn),jptstn(istn),1)
!1022 format(a4,18x,f6.2)
! write(6,1023) 'Wetsfc:',wetsfc(iptstn(istn),jptstn(istn),0)
!1023 format(a7,15x,f6.4)
! write(6,1024) 'Wetdp:',wetdp(iptstn(istn),jptstn(istn),0)
!1024 format(a6,16x,f6.4)
! write(6,*)
! write(6,*)
! 720 CONTINUE
!
IF( totapi0 > 0 ) THEN
WRITE(6,*)
WRITE(6,*) 'Root zone volumetric water contents' &
//' at initial API stations:'
WRITE(6,*)
DO i=1,totapi0
WRITE(6,1042) 'Stn #:',i,' has wetdp of:', &
wetdp(iptapi0(i),jptapi0(i),0)
1042 FORMAT(a6,1X,i2,a15,1X,f5.3)
END DO
WRITE(6,*)
END IF
!
WRITE(6,*)
WRITE(6,*) 'wetcanp analysis:'
WRITE(6,*) '-----------------------'
WRITE(6,*)
!
WRITE(6,1040) j5,wetcanp(i1,j5,0),wetcanp(i2,j5,0), &
wetcanp(i3,j5,0),wetcanp(i4,j5,0), &
wetcanp(i5,j5,0)
WRITE(6,1040) j4,wetcanp(i1,j4,0),wetcanp(i2,j4,0), &
wetcanp(i3,j4,0),wetcanp(i4,j4,0), &
wetcanp(i5,j4,0)
WRITE(6,1040) j3,wetcanp(i1,j3,0),wetcanp(i2,j3,0), &
wetcanp(i3,j3,0),wetcanp(i4,j3,0), &
wetcanp(i5,j3,0)
WRITE(6,1040) j2,wetcanp(i1,j2,0),wetcanp(i2,j2,0), &
wetcanp(i3,j2,0),wetcanp(i4,j2,0), &
wetcanp(i5,j2,0)
WRITE(6,1040) j1,wetcanp(i1,j1,0),wetcanp(i2,j1,0), &
wetcanp(i3,j1,0),wetcanp(i4,j1,0), &
wetcanp(i5,j1,0)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
1041 FORMAT(i3,3X,f5.1,3X,f5.1,3X,f5.1,3X,f5.1,3X,f5.1)
!
WRITE(6,*)
WRITE(6,*) 'tsfc analysis:'
WRITE(6,*) '-----------------------'
WRITE(6,*)
!
WRITE(6,1041) j5,tsfc(i1,j5,0),tsfc(i2,j5,0), &
tsfc(i3,j5,0),tsfc(i4,j5,0), &
tsfc(i5,j5,0)
WRITE(6,1041) j4,tsfc(i1,j4,0),tsfc(i2,j4,0), &
tsfc(i3,j4,0),tsfc(i4,j4,0), &
tsfc(i5,j4,0)
WRITE(6,1041) j3,tsfc(i1,j3,0),tsfc(i2,j3,0), &
tsfc(i3,j3,0),tsfc(i4,j3,0), &
tsfc(i5,j3,0)
WRITE(6,1041) j2,tsfc(i1,j2,0),tsfc(i2,j2,0), &
tsfc(i3,j2,1),tsfc(i4,j2,1), &
tsfc(i5,j2,1)
WRITE(6,1041) j1,tsfc(i1,j1,0),tsfc(i2,j1,0), &
tsfc(i3,j1,0),tsfc(i4,j1,0), &
tsfc(i5,j1,0)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
WRITE(6,*)
WRITE(6,*) 'tsoil analysis:'
WRITE(6,*) '-----------------------'
WRITE(6,*)
!
WRITE(6,1041) j5,tsoil(i1,j5,0),tsoil(i2,j5,0), &
tsoil(i3,j5,0),tsoil(i4,j5,0), &
tsoil(i5,j5,0)
WRITE(6,1041) j4,tsoil(i1,j4,0),tsoil(i2,j4,0), &
tsoil(i3,j4,0),tsoil(i4,j4,0), &
tsoil(i5,j4,0)
WRITE(6,1041) j3,tsoil(i1,j3,0),tsoil(i2,j3,0), &
tsoil(i3,j3,0),tsoil(i4,j3,0), &
tsoil(i5,j3,0)
WRITE(6,1041) j2,tsoil(i1,j2,0),tsoil(i2,j2,0), &
tsoil(i3,j2,0),tsoil(i4,j2,0), &
tsoil(i5,j2,0)
WRITE(6,1041) j1,tsoil(i1,j1,0),tsoil(i2,j1,0), &
tsoil(i3,j1,0),tsoil(i4,j1,0), &
tsoil(i5,j1,0)
WRITE(6,*)
WRITE(6,1035) i1,i2,i3,i4,i5
WRITE(6,*)
!
!--------------------------------------------------------------------
!
! Read in soil moisture data from a file.
!
!--------------------------------------------------------------------
!
ELSE IF( soilinit2 == 4 ) THEN
WRITE(6,'(a,a)') 'Reading soil moisture from ',soilinfl
CALL readsoil(nx,ny,nstyps,soilinfl,dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin, &
tem1,tem1,wetsfc,wetdp,wetcanp,snowdpth,soiltyp)
DO is = 1,nstyp
DO j=1,ny-1
DO i=1,nx-1
wetsfc (i,j,is) = wetsfc (i,j,0)
wetdp (i,j,is) = wetdp (i,j,0)
wetcanp(i,j,is) = wetcanp(i,j,0)
tsfc (i,j,0) = tsfc (i,j,0) &
+ tsfc (i,j,is) * stypfrct(i,j,is)
tsoil (i,j,0) = tsoil (i,j,0) &
+ tsoil (i,j,is) * stypfrct(i,j,is)
END DO
END DO
END DO
!
!--------------------------------------------------------------------
!
! Merge in soil data from soilinfl into soil data provided by history
! dump where ever there is water, also use snowdpth in soilinfl
! if present.
!
!--------------------------------------------------------------------
!
ELSE IF( soilinit2 == 5 ) THEN
WRITE(6,'(a,a)') 'Merging water soil data from ',soilinfl
! Note that nstyp=1 is assumed (or at least only the composite
! type is used) for data in soilinfl
soiltyp2 = soiltyp
CALL readsoil(nx,ny,nstyp,soilinfl,dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin, &
tsfc2,tsoil2,wetsfc2,wetdp2,wetcanp2,snowdpth,soiltyp2)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! DO is=1,nstyp
! DO j=1,ny-1
! DO i=1,nx-1
! IF (soiltyp(i,j,is) == 13) THEN
! IF (tsfcin /= 0) tsfc(i,j,is) = tsfc2(i,j,0)
! IF (tsoilin /= 0) tsoil(i,j,is) = tsoil2(i,j,0)
! IF (wsfcin /= 0) wetsfc(i,j,is) = wetsfc2(i,j,0)
! IF (wdpin /= 0) wetdp(i,j,is) = wetdp2(i,j,0)
! IF (wcanpin /= 0) wetcanp(i,j,is) = wetcanp2(i,j,0)
! END IF
! END DO
! END DO
! END DO
DO is=1,nstyp
DO j=1,ny-1
DO i=1,nx-1
IF (soiltyp(i,j,is) == 13) THEN
DO is2=1,nstyp
IF (soiltyp2(i,j,is2) == 13) THEN
IF (tsfcin /= 0) tsfc(i,j,is) = tsfc2(i,j,is2)
IF (tsoilin /= 0) tsoil(i,j,is) = tsoil2(i,j,is2)
IF (wsfcin /= 0) wetsfc(i,j,is) = wetsfc2(i,j,is2)
IF (wdpin /= 0) wetdp(i,j,is) = wetdp2(i,j,is2)
IF (wcanpin /= 0) wetcanp(i,j,is) = wetcanp2(i,j,is2)
EXIT
END IF
END DO
END IF
END DO
END DO
END DO
! is=0 filled in when program returns to main program
ELSE
WRITE(6,'(i6,a)') soilinit2,' is not a valid soilinit2 option!'
STOP
END IF
!wdt update
IF (soilinit2 == 5) THEN
DEALLOCATE(tsfc2)
DEALLOCATE(tsoil2)
DEALLOCATE(wetsfc2)
DEALLOCATE(wetdp2)
DEALLOCATE(wetcanp2)
DEALLOCATE(soiltyp2)
END IF
!
!-------------------------------------------------------------------
!
! End of subroutine
!
!-------------------------------------------------------------------
!
PRINT *,'Finishes mksoilvar OK!!'
!
RETURN
!
END SUBROUTINE mksoilvar
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FINDLC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE findlc(nx,ny,xs,ys,xpt,ypt,ipt,jpt,ireturn) 7
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Searches in x and y to find i,j location of xpt, ypt.
!
! X and Y do not have to be on a regular grid, however it is
! assumed that x and y are monotonically increasing as i and j
! indices increase.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! April 1992.
!
! MODIFICATION HISTORY:
!
! February, 1993 (K. Brewster)
! Additional documentation for ARPS 3.1 release
!
! October, 1994 (K. Brewster)
! Changed to reference scalar points.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! xs x coordinate of scalar points in physical/comp. space (m)
! ys y coordinate of scalar points in physical/comp. space (m)
!
! xpt location to find in x coordinate (m)
! ypt location to find in y coordinate (m)
!
! OUTPUT:
!
! ipt i index to the west of desired location
! jpt j index to the south of desired location
!
!-----------------------------------------------------------------------
!
! Arguments
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny ! Dimensions of ARPS grids
REAL :: xs(nx) ! x coordinate of scalar grid points in
! physical/comp. space (m)
REAL :: ys(ny) ! y coordinate of grid points in
! physical/comp. space (m)
REAL :: xpt ! location to find in x coordinate
REAL :: ypt ! location to find in y coordinate
INTEGER :: ipt ! i index to the west of desired
! location
INTEGER :: jpt ! j index to the south of desired
! location
INTEGER :: ireturn
!
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ireturn=0
!
DO i=2,nx
IF(xpt < xs(i)) EXIT
END DO
101 CONTINUE
ipt=i-1
IF( xpt > xs(nx-1) .OR. xpt < xs(1) ) ireturn=-1
DO j=2,ny
IF( ypt < ys(j)) EXIT
END DO
201 CONTINUE
jpt=j-1
IF( ypt > ys(ny-1) .OR. ypt < ys(1) ) ireturn=-2
RETURN
END SUBROUTINE findlc
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RDNCEPPR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rdnceppr(nx,ny,mxstns,mxdays,xs,ys,xll,yll, & 1,5
iendyr,iendmo,ienddy,ndays, &
prcpdir,prcplst, &
staid,iptstn,jptstn,obsprec,totpreca, &
tstns)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Program to use history data dump and surface property data to
! create the soil variables
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! 04/11/1995
!
! MODIFICATION HISTORY:
!
! 05/07/1998 (Yuhe Liu)
! Hub-CAPS's modification to the precipitation format change in
! extracting the information from each line.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Arguments
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,mxstns,mxdays
!
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: xll
REAL :: yll
INTEGER :: iendyr,iendmo,ienddy,ndays
CHARACTER (LEN=80) :: prcpdir
CHARACTER (LEN=80) :: prcplst
CHARACTER (LEN=7) :: staid(mxstns)
INTEGER :: iptstn(mxstns)
INTEGER :: jptstn(mxstns)
REAL :: obsprec(mxdays,mxstns)
REAL :: totpreca(mxstns)
INTEGER :: tstns
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: mxstafile
PARAMETER( mxstafile = 99000)
INTEGER :: sec24hr
PARAMETER( sec24hr = 86400 )
REAL :: in2cm,maxprcm
PARAMETER (in2cm=2.54,maxprcm=4.00)
!
CHARACTER (LEN=7) :: rdstaid
CHARACTER (LEN=7) :: rdst(3)
CHARACTER (LEN=7) :: dummy
CHARACTER (LEN=80) :: rfcname
CHARACTER (LEN=80) :: line
INTEGER :: iunit,istn,jstn,kstn,iday,out
INTEGER :: ilat(3),ilon(3)
INTEGER :: ipt,jpt,ios,ireturn,lendir
INTEGER :: istart,itime,iyr,imon,idy,ihr,imin,isec
INTEGER :: nread,nmatch
REAL :: lat,lon,xpt,ypt
REAL :: rdlat,rdlon,rdobpr
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-------------------------------------------------------------------
!
! Initialize precipitation observations to missing
!
!-------------------------------------------------------------------
!
DO istn=1,mxstns
staid(istn)=' '
DO iday=1,mxdays
obsprec(iday,istn)=-9.0
END DO
END DO
!
!-------------------------------------------------------------------
!
! Get list of precip stations
! and keep those that are within the domain.
!
!-------------------------------------------------------------------
!
iunit=38
OPEN(iunit,ERR=910,FILE=prcplst,IOSTAT=ios, &
STATUS='old',FORM='formatted')
READ(iunit,'(a5)') dummy
out=0
nread=0
DO istn=1,mxstns,3
READ(iunit,'(a7,i5,i6,8x,a7,i5,i6,8x,a7,i5,i6)', &
ERR=101,END=101) &
rdst(1),ilat(1),ilon(1), &
rdst(2),ilat(2),ilon(2), &
rdst(3),ilat(3),ilon(3)
nread=nread+3
DO kstn=1,3
jstn=istn-out+(kstn-1)
lat = 0.01*FLOAT(ilat(kstn))
lon =-0.01*FLOAT(ilon(kstn))
CALL lltoxy(1,1,lat,lon,xpt,ypt)
xpt=xpt-xll
ypt=ypt-yll
CALL findlc(nx,ny,xs,ys,xpt,ypt,ipt,jpt,ireturn)
IF( ireturn == 0 ) THEN
staid(jstn)=rdst(kstn)
iptstn(jstn)=ipt
jptstn(jstn)=jpt
totpreca(jstn)=0.0
WRITE(6,*) ' Precip stn ',rdst(kstn),' was found near i,j:', &
ipt,jpt
ELSE
! write(6,*)
! : ' Precip stn ',rdst(kstn),' is outside of the grid.'
out=out+1
END IF
END DO
END DO
101 CONTINUE
tstns=istn-out-1
CLOSE(iunit)
WRITE(6,'(a,i6,a)') ' Read ',nread,' precip station locations.'
WRITE(6,'(i12,a)') tstns,' inside grid.'
!
!-------------------------------------------------------------------
!
! Begin main day loop, reading precip file for each day.
!
!-------------------------------------------------------------------
!
lendir = LEN( prcpdir )
CALL strlnth( prcpdir, lendir )
CALL ctim2abss(iendyr,iendmo,ienddy,12,0,0, istart)
istart=istart-((ndays-1)*sec24hr)
!
DO iday=1,ndays
itime = istart + ((iday-1)*sec24hr)
CALL abss2ctim( itime, iyr, imon, idy, ihr, imin, isec )
WRITE(rfcname,'(a,a,i4.4,i2.2,i2.2)') &
prcpdir(1:lendir),'/usa-dlyprcp-',iyr,imon,idy
WRITE(6,'(/a,a)') ' Opening: ',rfcname
OPEN(iunit,ERR=400,FILE=rfcname, &
STATUS='old',FORM='formatted')
nread=0
nmatch=0
READ(iunit,'(a5)',ERR=351,END=351) dummy
DO kstn=1,mxstafile
READ(iunit,'(a80)',ERR=350,END=351) line
READ(line,*,ERR=350,END=350) rdlat,rdlon,rdobpr
IF (line(23:23) == ' ') THEN
rdstaid = line(24:31)
ELSE
rdstaid = line(23:30)
END IF
nread=nread+1
DO istn=1,tstns
IF( rdstaid == staid(istn) ) THEN
nmatch=nmatch+1
obsprec(iday,istn)=in2cm*rdobpr
obsprec(iday,istn)=AMIN1(maxprcm,obsprec(iday,istn))
EXIT
END IF
END DO
301 CONTINUE
END DO
350 CONTINUE
351 CONTINUE
CLOSE(iunit)
WRITE(6,'(a,i6,a)') ' Read ',nread,' precip obs.'
WRITE(6,'(i12,a)') nmatch,' obs matched inside grid.'
END DO
400 CONTINUE
401 CONTINUE
!
!-------------------------------------------------------------------
!
! Normal end
!
!-------------------------------------------------------------------
!
RETURN
!
!-------------------------------------------------------------------
!
! Error destinations
!
!-------------------------------------------------------------------
!
910 CONTINUE
WRITE(6,'(a,a,/a,i6)') &
'Error opening list of precip stations,',prcplst,' iostat = ',ios
STOP
END SUBROUTINE rdnceppr