c file name: darcy.for C23456789*123456789*123456789*123456789*123456789*123456789*123456789*12 subroutine darcy (layrsn, bszlyd, bheaep, bhzper, sunris, & bhrsk, bh0cb, netrad, bhrwc0, bszlyt, daysim, * bsdblk1) integer daysim ! temp fix c + + + PURPOSE + + + c This subroutine predicts on an hourly basis soil water profile, c soil water content at the soil-air interface, potential and c actual soil evaporation, and deep percolation. c DATE: 09/16/93 c MODIFIED: 12/10/93 c MODIFIED: 08/25/95 c + + + KEYWORDS + + + c soil water redistribution, evaporation, deep percolation c + + + ARGUMENT DECLARATIONS + + + real sunris integer layrsn real bhrwc0(*) real bh0cb(*) real bszlyd(*) real bhzper real netrad real bheaep(*) real bhrsk(*) real bszlyt(*) real bsdblk1 c + + + ARGUMENT DEFINITIONS + + + c bhrwc0 - Hourly surface water content - volumetric c bh0cb - Power of Campbell's water release curve model c bszlyd - Depth to bottom of soil layer from surface (mm) c bhzper - Output variable for amount of daily deep percolation c sunris - Time of sunrise c laysrn - Number of soil layers used in simulation c netrad - Net radiation (mj/m^2/day) c bheaep - Soil air entry potential (j/kg) c bhrsk - Saturated soil hydraulic conductivity (m/s) c (mm/day) c bszlyt - Layer thickness (mm) c + + + PARAMETERS + + + c + + + COMMON BLOCKS + + + *$noereference include 'p1werm.inc' include 'p1const.inc' include 'p1unconv.inc' include 'wpath.inc' C *** include 'm1sim.inc' include 'm1subr.inc' include 'm1flag.inc' include 'h1et.inc' c + + + LOCAL COMMON BLOCKS + + + include 'hydro/htheta.inc' include 'hydro/hh2o.inc' C *** include 'hydro/hlayrs.inc' *$reference c + + + LOCAL VARIABLES + + + integer day,mo,yr integer hourstep integer k, l integer lrx real*8 amep real*8 aswc real*8 aswcr real*8 avep real*8 cond(mnsz), conda(mnsz+1) real*8 dph(24) real*8 eah(24) real*8 eph(24) real*8 ephc real*8 remtim , deltim, olddeltim, newdeltim real*8 minstep, maxchange real*8 wfluxr(mnsz+1) real*8 wfluxn(mnsz) real*8 thetad(0:10) real*8 thetah(0:10,24) real*8 wct(mnsz) c$ timestep test real*8 thetatemp(0:mnsz) c$ end timestep test logical dbgflg real*8 dltm real*8 wcad(mnsz) real*8 dist(mnsz+1) real*8 lenday real*8 sunset real*8 perbeg,perend integer flipflop real*8 div_deltim c temp variable to show time step reductions c integer timecount c + + + LOCAL DEFINITIONS + + + c amep - Daily wave amplitude of potential bare soil c evaporation (mm/hr) c aswc - Max available surface soil water content (m^3/m^3) c aswcr - Ratio of actual to maximum available surface soil c water content c avep - Time-average value of daily potential bare soil c evaporation (mm/hr) c cond - Unsaturated soil hydraulic conductivity (m/s) c dph - Hourly deep percolation (mm/hr) c eah - Hourly actual bare soil evaporation (mm/hr) c eph - Hourly potential bare soil evaporation (mm/hr) c ephc - Time slice potential bare soil evaporation (mm) c swh - Soil water hydraulic head (mm) c wfluxr - Rate of soil water flow (mm/sec) at top of each layer c wfluxn - Net flow rate of soil water flow (mm/sec) into each layer c thetad - Average daily soil wetness at the soil-atmosphere c interface and at each of the simulation layers, m^3/m^3 c thetah - Hourly simulation soil wetness at the soil-atmosphere c interface and at each of the simulation layers, m^3/m^3 c wct - temporary layer water content (mm) c dist - Distance between midpoints of layers. Used in both flux c equations and gravitational heads. c lenday - length of day (seconds) c sunset - time of sunset (seconds) c perbeg - beginning time of timeslice being computed (seconds) c perend - ending time of timeslice being computed (seconds) c flipflop - indicator that adjacent layers went from wet/dry to dry/wet c div_deltim - time step devisor calculated from magnitude of change c + + + SUBROUTINES CALLED + + + c caldatw c + + + FUNCTION DECLARATIONS + + + real calctht0 real availwc c + + + DATA INITIALIZATIONS + + + save deltim if (daysim.eq.1) deltim = 3600 olddeltim = deltim data dbgflg /.false./ call dbgwrt(dbgflg,'darcy 1') data cond /mnsz*0.0/ c *** write(*,*) 'in darcy' c parameters regulating automatic time step adjustments c minstep - in seconds. This small value allows small changes to occur with large gradients c maxchange - 0.02 is a good value. larger values allow getting into unstable excursions. c Smaller values slightly improve hourly output values. minstep = 0.0001 maxchange = 0.02 flipflop = 0 ahzea = 0.0 swc = 0.0 bhzper = 0.0 sunris = sunris*3600 lenday = 24.*3600-2.*sunris sunset = sunris + lenday c + + + OUTPUT FORMATS + + + 2000 format('*','Time of Sunrise: ',f10.2,/,'*','Length of Day: ' &,f10.2,/,'*','Time of Sunset: ',f10.2) 2009 format ('*',' Hourly HYDROLOGY output ') 2010 format('*','date: ',i2,'/',i2,'/',i4,' Subregion # ' & ,i4) 2020 format ('*',22x,3('*'),' hourly soil water data ',3('*')) 2030 format ('*',79('-')/'*','hr',t7,'eph',t13,'eah',t17,'dph',t24, & 'swc',t28,'theta(0)',1x,11('-'),'water content of soil layers', & 10('-')/,'*',t7,'--------mm--------',2x,25('-'), & 'm^3/m^3',24('-')) 2040 format('*',18x,f8.3,1x,f6.3,1x,10(1x,f6.3)) 2050 format(1x,i2,2f6.3,f4.1,f8.3,1x,f6.3,1x,10(1x,f6.3)) 2060 format('*',79('-')) 2065 format('*','average daily wetness',6x,f6.3,1x,10(1x,f6.3)) 2070 format('*','ep=',f7.2,' mm',' (cumulative amount of daily', & ' potential soil evaporation)') 2080 format('*','ea=',f7.2,' mm',' (cumulative amount of daily', & ' actual soil evaporation)') 2090 format('*','bhzper=',f6.2,' mm',' (cumulative amount of daily', & ' deep percolation)') 2100 format('darcy: lrx, wc, wct, wfluxn, deltim, wfluxr(top), wflu &xr(bottom)',/,i2,7f11.3,/,7f11.3) c c + + + END SPECIFICATIONS + + + c c output initial daily conditions if ((am0ifl .eqv. .true.).and.((am0hfl .eq. 2).or.(am0hfl .eq. 6) & .or. (am0hfl .eq. 3) .or. (am0hfl .eq. 7))) then open(12,file= 'water.out', & access='sequential', status='unknown') end if C calc the solar energy delivered to surface for evaporation loss C this does not correspond to the docs but is a better method C initialize layer vars that don't change during iterations and theta related C variables that are not used until now dist(1) = bszlyd(1) / 2. do lrx = 2, layrsn dist(lrx) = (bszlyt(lrx) + bszlyt(lrx-1))/2. end do dist(layrsn+1) = bszlyt(layrsn)/2. do lrx = 1, layrsn C *** theta(lrx) = wc(lrx) / bszlyt(lrx) cm(lrx)= 2. * bh0cb(lrx) + 3. !H-76 cond(lrx) = bhrsk(lrx)*(theta(lrx)/thetas(lrx))**cm(lrx) * mtomm !H-74,75 potm(lrx) = bheaep(lrx)*((theta(lrx)/thetas(lrx))** * (-bh0cb(lrx))) !H-73 swm(lrx) = (potm(lrx) / 9.8) * mtomm !H-47 Convert from J/kg to mm of H20 wc(lrx) = theta(lrx) * bszlyt(lrx) wcad(lrx) = 0.0 end do swci = sum(wc(1:layrsn)) c write(12,*) 'darcy: tim',daysim,(theta(k), k=1,layrsn) c write(12,*) 'darcy: tim',daysim,(swm(k), k=1,layrsn), c &(cond(k), k=1,layrsn),(theta(k), k=1,layrsn) C *** write(*,*) ' swm ', swm(1:layrsn) C *** write(*,*) ' cond ', cond(1:layrsn) C *** write(*,*) ' wc ', wc(1:layrsn) C eph & eah are both declared 0 here; it is the middle of the night and C they are, by previous calculations, 0 theta(0) = calctht0(bszlyd, thetes, thetar, theta, * thetaw, thetaf(1) - thetaw(1), 0.0_8, 0.0_8) !H-64,65,66 c write(*,*) 'calctht0: darcy - theta(0) after',theta(0) c print out zero hour initialization values if ((am0hfl .eq. 2) .or. (am0hfl .eq. 3) .or. (am0hfl .eq. 6) .or. & (am0hfl .eq. 7)) then call caldatw(day,mo,yr) write(12,2009) write(12,2010) day,mo,yr,am0csr c write(12,2000) sunris,lenday,sunset write(12,2020) write(12,2030) write(12,2040) swci,theta(0),(theta(k), k=1,layrsn) end if c temp variable to show time step reductions c timecount = 0 c start the hourly cycle here do 300 hourstep=1,24 c temp variable to show time step reductions c timecount = 0 C initialize various parameters, deltim == time increment, remtim == time remaining in hour remtim = 3600.0 eah(hourstep) = 0.0 eph(hourstep) = 0.0 dph(hourstep) = 0.0 202 continue !start of timeslice wct = wc thetatemp = theta C calc the relative amount of available surface water from 0.7 field capacity to reduced wilting point aswcr = availwc( theta(1),thetar(1),0.7*thetaf(1) ) !H-60,62 if (thetatemp(layrsn).gt.thetaf(layrsn)) then wfluxr(layrsn+1) = cond(layrsn) else wfluxr(layrsn+1) = 0.0 endif C calc the conductivity and rate of flow in for each layer do lrx = 2,layrsn conda(lrx) = (cond(lrx-1)*bszlyt(lrx-1)+cond(lrx)*bszlyt(lrx)) ! H-48 & / (bszlyt(lrx)+bszlyt(lrx-1)) wfluxr(lrx) = (swm(lrx-1)-swm(lrx)) * conda(lrx) / dist(lrx) & + conda(lrx) ! H-44,H-45 end do C calc the net flow into (+) or out of (-) each layer do lrx = 2, layrsn wfluxn(lrx) = wfluxr(lrx) - wfluxr(lrx+1) end do 203 continue c temp variable to show time step reductions c timecount = timecount+1 perbeg = (hourstep-1)*3600+(3600-remtim) perend = perbeg+deltim if ((perend.lt.sunris).or.(perbeg.gt.sunset)) then ephc = 0.0 ! the entire timeslice is at night else if (perbeg.lt.sunris) perbeg = sunris ! period starts before sunrise if (perend.gt.sunset) perend = sunset ! period ends after sunset ephc = (ahzep / 2) *(-cos(pi*(perend-sunris)/lenday) + * cos(pi*(perbeg-sunris)/lenday)) endif C adjust soil evap down for drier layer 1 wfluxr(1) = (-ephc / deltim) * aswcr wfluxn(1) = wfluxr(1) - wfluxr(2) C add water moving down into layer, sub water moving down, out of layer C compute new values for wc & theta C that is water content per layer (wc), flow rate (theta) 204 continue do lrx = 1,layrsn wc(lrx)= wct(lrx) + wfluxn(lrx) * deltim !H-70&71 theta(lrx) = wc(lrx) / bszlyt(lrx) !H-72 C *** debugging fix C drop least significant digits c theta(lrx) = int(theta(lrx)*100000.)/100000. C *** eodf c if theta changes more than maxchange fraction of diff between saturated and really dry C adjust time step and recalculate div_deltim = maxchange*(thetas(lrx)-thetar(lrx)) if (abs(theta(lrx)-thetatemp(lrx)) .gt. div_deltim) then c write(*,*) 'darcy: reducing time step, change', lrx, c & theta(lrx),thetatemp(lrx),remtim,hourstep,daysim if( deltim .gt. minstep ) then div_deltim = abs(theta(lrx)-thetatemp(lrx))/div_deltim div_deltim = max( div_deltim, 2.0 ) deltim = deltim / div_deltim go to 203 else write(*,*) 'darcy: min time step, change',lrx, & theta(lrx),thetatemp(lrx),remtim,hourstep,daysim endif endif C if water content is going from within to out of bounds c ie. to greater than saturation or to less than really dry, C adjust time step and recalculate if ((thetatemp(lrx) .lt. thetas(lrx)) .and. & (theta(lrx) .gt. thetas(lrx))) then c write(*,*) 'darcy: reducing time step, too wet', lrx, c & theta(lrx),thetas(lrx),remtim,hourstep,daysim if( deltim .gt. minstep ) then deltim = deltim / 2.0 go to 203 else write(*,*) 'darcy: min time step, too wet' endif else if ((thetatemp(lrx) .gt. thetar(lrx)) .and. & (theta(lrx) .lt. thetar(lrx))) then c write(*,*) 'darcy: reducing time step, too dry', lrx, c & theta(lrx),thetaf(lrx),remtim,hourstep,daysim if( deltim .gt. minstep ) then deltim = deltim / 2.0 go to 203 else c desperation step to keep water content from going negative write(*,*) & 'darcy: min time step, too dry, amount(mm)', & (thetar(lrx)-theta(lrx))*bszlyt(lrx) theta(lrx) = thetar(lrx) endif endif end do c check for layer content flip flops newdeltim = deltim do lrx = 2,layrsn if (((thetatemp(lrx-1) .gt. thetatemp(lrx)).and. & (theta(lrx-1) .lt. theta(lrx))) .or. & ((thetatemp(lrx-1) .lt. thetatemp(lrx)).and. & (theta(lrx-1) .gt. theta(lrx)))) then if ( flipflop .gt. 0 ) then c write(*,*) 'darcy: flipflop',lrx-1,lrx, c & theta(lrx-1),theta(lrx), c & thetatemp(lrx-1),thetatemp(lrx), c & remtim,hourstep,daysim c Here we want to accept this time slice, and cut the c next time step in half to avoid over correction. In c accepting the time step we need to retain DELTIM until c the accumulators have been updated, then do the adjustment. if( deltim .gt. minstep ) then newdeltim = deltim / 4.0 else write(*,*) 'darcy: min time step, flipflop' endif else flipflop = 1 endif else flipflop = 0 endif end do C compute new values for potm, swm, and cond C potential head (ptom), soil/water hydr head (swm) C and conductivity (cond) do lrx = 1,layrsn potm(lrx) = bheaep(lrx) * ((theta(lrx)/thetas(lrx)) * **(-bh0cb(lrx))) swm(lrx) = (potm(lrx) / 9.8) * mtomm !H-47 Convert from J/kg to mm of H20 cond(lrx)= (bhrsk(lrx)*(theta(lrx)/thetas(lrx))**cm(lrx))* * mtomm end do C completed timeslice, sum up, etc. c write(*,*) 'darcy: timeslice:',deltim,remtim,hourstep,daysim C *** eah(hourstep) = eah(hourstep) + ephc * deltim eah(hourstep) = eah(hourstep) + (-wfluxr(1)) * deltim eph(hourstep) = eph(hourstep) + ephc dph(hourstep) = dph(hourstep) + wfluxr(layrsn+1) * deltim c write(12,*) 'darcy: tim', c &daysim+hourstep/24.0-remtim/86400.0,(theta(k), k=1,layrsn) c write(12,*) 'darcy: tim', c &daysim+hourstep/24.0-remtim/86400.0,(swm(k), k=1,layrsn), c &(cond(k), k=1,layrsn),(theta(k), k=1,layrsn) c temp variable to show time step reductions c timecount = timecount+1 c update time variables for end of time slice remtim = remtim - deltim deltim = newdeltim * 2.0 if ( remtim .gt. 0.0 ) then if (deltim .gt. remtim) then olddeltim = deltim deltim = remtim endif go to 202 else deltim = olddeltim if (deltim .gt. 3600.0) deltim = 3600.0 endif c temp variable to show time step reductions c if ( timecount .gt. 1 ) then c write(*,*) 'darcy: steps-hr ',timecount,hourstep,daysim c endif C completed the hour, sum up theta(0) = calctht0(bszlyd, thetes, thetar, theta, & thetaw, thetaf(1) - thetaw(1), eph(hourstep), eah(hourstep)) !H-64,65,66 swc = sum(wc(1:layrsn)) bhrwc0(hourstep) = theta(0) / bsdblk1 c output the hourly info here if ((am0hfl .eq. 2) .or. (am0hfl .eq. 3) .or. (am0hfl .eq. 6) .or. & (am0hfl .eq. 7)) then write(12,2050) hourstep,eph(hourstep),eah(hourstep), & dph(hourstep),swc,theta(0),(theta(k), k=1,layrsn) end if call dbgwrt(dbgflg,'darcy 12') c Storing hourly soil wetness value thetah(0,hourstep) = theta(0) do k=1,layrsn thetah(k,hourstep) = theta(k) end do 300 continue C completed the day C *** do lrx=1,layrsn C *** theta(lrx) = wcad(lrx)/ (3600*24*bszlyt(lrx)) C *** end do C *** write(*,*) ' darcy 2: wc ', (wc(k), k=1,layrsn) ahzea = sum(eah(1:24)) bhzper = sum(dph(1:24)) bhzper = bhzper + dflout c calculating average daily soil wetness do k=0,layrsn thetad(k) = sum(thetah(k,1:24))/24 end do call dbgwrt(dbgflg,'darcy 13') c output the daily info here if ((am0hfl .eq. 2) .or. (am0hfl .eq. 3) .or. (am0hfl .eq. 6) .or. & (am0hfl .eq. 7)) then write(12,2065) (thetad(k), k=0, layrsn) write(12,2060) write(12,2070) ahzep write(12,2080) ahzea write(12,2090) bhzper end if call dbgwrt(dbgflg,'darcy 14') c write(12,*) 'darcytrack out ',(swm(k), k=1,layrsn), c &(cond(k), k=1,layrsn),(theta(k), k=1,layrsn) c c close (unit = 12) C *** write(*,9905) sum(wcad(1:layrsn)/(3600*24)), C *** * wcad(1:layrsn)/(3600*24) 9905 format(' darcy: wcad ', 10f10.2) return end