c$Author: fredfox $ c$Date: 2003-02-27 20:27:04 $ c$Revision: 1.52.2.1 $ c$Source: /weru/cvs/weps/weps.src/hydro/darcy.for,v $ subroutine darcy(daysim, numeq, bszlyt, bszlyd, theta, bulkden, & bthetas, bthetar, bhrsk, bheaep, bh0cb, & bsfcla, bsfom, bhtsav, & bwtdmxprev, bwtdmn, bwtdmx, bwtdmnnext, bwtdpt, & rise, daylength, bhzep, bwzdpt, bwdurpt, bwpeaktpt, & bbdstm, bbffcv, bszrro, bszrr, bmzele, bhrwc0, & bhzea, bhzper, bhzrun, bhzinf, bhzwid ) 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, runoff, ponding and deep percolation. c + + + KEYWORDS + + + c soil water redistribution, evaporation, runoff, deep percolation c + + + ARGUMENT DECLARATIONS + + + integer daysim, numeq real theta(0:*), bszlyt(*), bulkden(*), bszlyd(*), & bthetas(*), bthetar(*), bhrsk(*), bheaep(*), bh0cb(*), & bsfcla(*), bsfom(*), bhtsav(*), & bwtdmxprev, bwtdmn, bwtdmx, bwtdmnnext, bwtdpt, & rise, daylength, bhzep, bwzdpt, bwdurpt, bwpeaktpt, & bbdstm, bbffcv, bszrro, bszrr, bmzele, bhrwc0(*), & bhzea, bhzper, bhzrun, bhzinf, bhzwid c + + + ARGUMENT DEFINITIONS + + + c daysim - day of the simulation (very useful for debugging, not necessary otherwise) c numeq - number of equations to be solved = layrsn + 6 c bszlyt(*) - thickness of the soil layer (mm) c bszlyd(*) - thickness of the soil layer (mm) c theta(*) - volumetric water content (m^3/m^3) c bulkden(*) - soil bulk density Mg/m^3) c bthetas(*) - saturated volumetric water content (m^3/m^3) c bthetar(*) - residual (conductivity) volumetric water content (m^3/m^3) c bhrsk(*) - saturated hydraulic conductivity (m/s) c bheaep(*) - air entry potential (J/kg) c bh0cb(*) - exponent of Campbell soil water release curve (unitless) c bsfcla(*) - fraction of soil mineral content which is clay (unitless) c bsfom(*) - fraction of total soil which is organic (unitless) c bhtsav(*) - daily average soil temperature (C) c bwtdmxprev - maximum air temperature of previous day (C) c bwtdmn - minimum air temperature (C) c bwtdmx - maximum air temperature (C) c bwtdmnnext - minimum air temperature of next day (C) c bwtdpt - dew point temperature (C) c rise - time of the sunrise (hour of day) c daylength - time of daylight (hours) c bhzep - potential soil evaporation (mm/day) c bwzdpt - rainfall depth (mm) c bwdurpt - duration of precipitation (hours) c bwpeaktpt - normalized time to peak of precipitation (time to peak/duration) c bbdstm - total number of stems (#/m^2) c bszrro - original random roughness height, after tillage, mm c bszrr - Allmaras random roughness parameter (mm) c bmzele - Average site elevation (m) c bhrwc0(*) - Hourly values of surface soil water content (not as in soil) c bhzea - accumulated daily evaporation (mm) (comes in with a value set from snow evap) c bhzper - accumulated daily drainage (deep percolation) (mm) c bhzrun - accumulated daily runoff (mm) c bhzinf - depth of water infiltrated (mm) c bhzwid - Water infiltration depth into soil profile (mm) c + + + PARAMETERS + + + real pi parameter( pi = 3.1415927 ) c + + + COMMON BLOCKS + + + include 'p1werm.inc' include 'p1unconv.inc' include 'm1flag.inc' include 'm1subr.inc' c + + + LOCAL COMMON BLOCKS + + + include 'hydro/lsoda.inc' include 'hydro/dvolwparam.inc' include 'hydro/vapprop.inc' include 'hydro/clayomprop.inc' c + + + LOCAL VARIABLES + + + integer itol, itask, jt, iopt, neq(2) real tday,tout,relerr(1),abserr(numeq), & volw(numeq) integer kindex, hourstep, day, mo, yr real swc, swci, ref_ranrough, ranrough, thetad(0:numeq) integer lstep,lfunc,ljac real temp, netflux(numeq) real evap, runoff, drain, infil, tdaystart integer iminlay, imaxlay c + + + LOCAL DEFINITIONS + + + c itol - lsoda: setting to make relerr and/or abserr scalar or array c itask - lsoda: stepping mode c jt - lsoda: flag for selection of jacobian matrix c iopt - lsoda: flag for use of extra inputs c neq(1) - lsoda: number of equations to be solved c tday - time of day (seconds) c tout - time of day for end of integration step (seconds) c relerr(1) - lsoda: relative error value to control integration accuracy c abserr(numeq) - lsoda: absolute error value to control integration accuracy c volw(1-5) - accumulated values integrate in time c (1) accumulated rainfall depth (total surface water application)(m) c (2) runoff depth (m) c (3) evaporation depth (m) c (4) depth of water infiltrated (m) c (5) ponded water depth (m) c volw(numeq) - (iminlay:imaxlay) volume of water in a soil layer (m) c (imaxlay+1) depth of water drained (m) c kindex - array index for loops c hourstep - step counter for 24 hourly steps c day - day of month c mo - month of year c yr - year of simulation c swc - total depth of water in soil profile (mm) c swci - total depth of water in soil profile a start of day (mm) c ref_ranrough - random roughness after last tillage (m) c ranrough - present random roughness (m) c thetad(numeq) - daily average volumetric water content for layer c lstep - preserve number of steps from previous lsoda call c lfunc - preserve number of function evaluations from previous lsoda call c ljac - preserve number of jacobian evaluations from previous lsoda call c temp - temporary value c prevevap - previous hour accumulated soil surface evaporation depth (mm) c evap - accumulated soil surface evaporation depth (mm) c prevrunoff - previous hour accumulated soil surface runoff (mm) c runoff - accumulated soil surface runoff (mm) c prevdrain - previous hour accumulated soil drainage (mm) c drain - accumulated soil drainage (mm) c tdaystart - time in seconds at the beginning of the day (reset to zero when solver reinitialized) c iminlay - minimum index for placement of soil layers in volw array c imaxlay - maximum index for placement of soil layers in volw array c + + + SUBROUTINES CALLED + + + c slsoda - livermore solver for ordinary differential equations c + + + FUNCTION DECLARATIONS + + + external dvolw, jac real*4 soilrelhum, volwatadsorb, & volwat_matpot_bc, unsatcond_bc, & atmpreselev, depstore, fricfact, store c + + + DATA INITIALIZATIONS + + + 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,'evap',t13,'run',t17,'dper',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 *** write(*,*) 'in darcy' c initialize step reporting counters lstep = 0 lfunc = 0 ljac = 0 c bhzea not zeroed since snow evap must be included bhzper = 0.0 bhzrun = 0.0 bhzinf = 0.0 c initialize LSODA options call xsetf(0) !do not print internal error messages c initialize values for the DVOLWPARAM common block c values are place in volw array position, which is structured with c surface variables first, then the soil surface to depth and c then drainage at depth. layrsn = numeq - 6 iminlay = 6 imaxlay = layrsn + 5 do kindex = 1, layrsn tlay(kindex) = bszlyt(kindex) * mmtom thetas(kindex) = bthetas(kindex) thetar(kindex) = bthetar(kindex) ksat(kindex) = bhrsk(kindex) airentry(kindex) = bheaep(kindex) / gravconst lambda(kindex) = 1.0 / bh0cb(kindex) temp = bulkden(kindex)*1000.0 !convert Mg/m^3 to kg/m^3 theta80rh(kindex) = volwatadsorb( temp, & bsfcla(kindex), bsfom(kindex), & claygrav80rh, orggrav80rh ) thetaw(kindex) = volwat_matpot_bc(potwilt, thetar(kindex), & thetas(kindex), airentry(kindex), lambda(kindex)) c temporary field capacity c if( daysim.eq.1 ) then c temp = -3.33 c temp = volwat_matpot_bc(temp, thetar(kindex), c & thetas(kindex), airentry(kindex), lambda(kindex)) c write(*,*) 'darcy:layer:',kindex, temp, thetaw(kindex) c end if c temporary field capacity soiltemp(kindex) = bhtsav(kindex) end do dist(1) = tlay(1) / 2. depth(1) = dist(1) do kindex = 2, layrsn dist(kindex) = (tlay(kindex) + tlay(kindex-1))/2. depth(kindex) = depth(kindex-1) + dist(kindex) end do airtmaxprev = bwtdmxprev airtmin = bwtdmn airtmax = bwtdmx airtminnext = bwtdmnnext tdew = bwtdpt lenday = daylength * hrtosec sunrise = rise * hrtosec sunset = sunrise + lenday c amplitude of evaporation (m) if( lenday.gt.0.0 ) then evapamp = pi*bhzep*mmtom/lenday/2.0 else evapamp = 0.0 end if raindepth = bwzdpt * mmtom !storm total depth (m) rainstart = 1.0 if( (bwdurpt.lt.0.08333333).and.(bwzdpt.gt.0.001) ) then !0.0833333 = 5 minutes = 300sec !enforce minimum time so solver does not miss rainfall. The !minimum is really between 60 and 120 for the example run, but to make sure used 300 !add sanity check for duration amount relationship from !Linsley,R.K., M.A. Kohler, J.L.H. Paulhus. 1982. Hydrology for Engineers. p80. rainend = rainstart & + max(300.0,10.0**(log10(bwzdpt)/0.48-1.90231428)) rainmid = rainstart & + max(rainstart,min(rainend,bwpeaktpt*300.0)) else rainend = rainstart + max(300.0, bwdurpt * hrtosec) rainmid = rainstart & + max(rainstart,min(rainend,bwpeaktpt*bwdurpt*hrtosec)) end if ref_ranrough = bszrro * mmtom ranrough = bszrr * mmtom soilslope = 0.01 pondmax = depstore( ranrough, soilslope )!maximum ponding depth in meters dw_friction = fricfact( ref_ranrough, ranrough, bbdstm, bbffcv ) slopelength = 50.0 !temporarily set since value is not supplied atmpres = atmpreselev( bmzele ) c initialize state of soil do kindex=1,layrsn if(theta(kindex).le.0.0) then write(*,*)'darcy:begin theta<0',kindex-5,theta(kindex) stop endif volw(kindex+5) = theta(kindex)*tlay(kindex) lastvolw(kindex+5) = -1.0e15 end do do kindex=iminlay,imaxlay prevvolw(kindex) = volw(kindex) end do bhzwid = 0.0 c initialize auxiliary variables for integration c if( daysim .eq. 1 ) then !commented out to initialized every day anew do kindex=1,5 prevvolw(kindex) = 0.0 end do prevvolw(numeq) = 0.0 c endif do kindex=1,5 volw(kindex) = prevvolw(kindex) end do volw(numeq) = prevvolw(numeq) 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) sunrise,lenday,sunset write(12,2020) write(12,2030) swci = sum(volw(6:layrsn+5)) * mtomm c theta(0) = calctht0(bszlyd, thetes, thetar, theta, c * thetaw, thetaf(1) - thetaw(1), 0.0, 0.0) !H-64,65,66 write(12,2040) swci,theta(0),(theta(kindex), kindex=1,layrsn) end if neq(1) = numeq !can't declare parameter as array, but must be passed as an array neq(2) = 0 !this flag is used to activate printing internal to dvolw c if(daysim.eq.48) then c neq(2) = 1 c else c neq(2) = 0 !this flag is used to activate printing internal to dvolw c endif itol = 2 !lsoda relerr is scalar and abserr is an array c itol = 1 !lsoda relerr is scalar and abserr is scalar relerr(1) = 1.0e-4 ! relative error tolerance thetad(0) = 0.0 ! layer daily average do kindex=1,numeq abserr(kindex) = 1.0e-6 ! absolute error tolerance thetad(kindex) = 0.0 ! layer daily average end do itask = 1 !lsoda returns value at tout c iopt = 0 !lsoda no extra inputs iopt = 1 !lsoda using optional inputs on rwork and iwork(5-10) do kindex=5,10 rwork(kindex) = 0.0 iwork(kindex) = 0 end do rwork(6) = 7200 !maximum allowed step size (seconds) c jt = 2 !lsoda internally generated Jacobian matrix jt = 5 !lsoda internally generated banded Jacobian matrix iwork(1) = 3 !ml, lower half band width of jacobian matrix iwork(2) = 1 !mu, upper half band width of jacobian matrix iwork(6) = 500 !maximum number of steps allowed before error generated c if( daysim .eq. 1 ) then !commented out to initialized every day anew istate = 1 beginday = 0.0 t = beginday c endif hourstep = 1 c start loop in time 30 continue tday = t - beginday tout = beginday + hourstep*hrtosec temp = (rainstart+rainend)/2.0 if( (raindepth.gt.0.0) & .and.(tday.le.rainstart) & .and.(tout.gt.rainstart) & .and.(tout.gt.temp)) then istate = 1 !-- initializes solution routines tout = min(temp, hourstep*hrtosec) !guarantees that integration will find rain event beginday = 0.0 t = tday c lstep = 0 !step reporting counters c lfunc = 0 c ljac = 0 c itask = 2 !put in single step mode for more error reporting end if c settings for single step mode, retaining values on the hour c itask = 5 c rwork(1) = tout 40 continue call slsoda(dvolw,neq,volw,t,tout,itol,relerr,abserr,itask, & istate,iopt,rwork,lrw,iwork,liw,jac,jt) c if( daysim.eq.6321 ) then c write(*,*) 'steps, functions, jacobians, order:', c & iwork(11)-lstep, c & iwork(12)-lfunc,iwork(13)-ljac,iwork(14) c end if c lstep = iwork(11) c lfunc = iwork(12) c ljac = iwork(13) c------------------------------------------------------------- c Was the step successful? If not, quit with an explanation. c------------------------------------------------------------- if(istate .lt. 0) then if(istate.eq.-1) then istate=2 write(*,*) 'day',daysim,'time',t,'3k steps ', & 'infil=',( volw(4) - prevvolw(4) ) * mtomm, & 'drain=',(volw(numeq)-prevvolw(numeq))*mtomm c do kindex=1,layrsn c write(*,*) 'darcy:s,r,e,b',bthetas(kindex), c & bthetar(kindex),bheaep(kindex),bh0cb(kindex) c theta(kindex) = volw(kindex+5)/tlay(kindex) c end do c write(*,*) 'darcy:theta',(theta(kindex),kindex=1,layrsn) c call dvolw(neq,t,volw,netflux) c write(*,*) 'darcy:netflux', c & (netflux(kindex+5),kindex=1,layrsn) else if(istate.eq.-5) then if( jt .eq. 5) then write(*,*)'Darcy: Switching to full matrix: day:', & & daysim,' time:',t,' istate:',istate istate = 2 jt = 2 else write(*,*)'Darcy: Band to full matrix switch failed' stop end if else write(*,*)'Darcy: Failed day:',daysim,' time:',t,' istate:',& & istate stop end if goto 40 end if c------------------------------------- c step was sucessful, continue on c------------------------------------- c if( daysim.eq.112 ) then c call subroutine to get flux rate values c neq(2) = 1 c call dvolw(neq,t,volw,netflux) c neq(2) = 0 c write(*,*)'r,o,e,i,p',t, c & volw(1),netflux(1), c & volw(2),netflux(2), c & volw(3),netflux(3), c & volw(4),netflux(4), c & volw(5),netflux(5) c end if if( tout.lt.(hourstep*hrtosec)) goto 30 c------------------------------------- C completed the hour, sum up c------------------------------------- do kindex=1,layrsn theta(kindex) = volw(kindex+5)/tlay(kindex) end do swci = sum(volw(6:layrsn+5)) * mtomm theta(0) = theta(1) c theta(0) = calctht0(bszlyd, thetes, thetar, theta, c & thetaw, thetaf(1) - thetaw(1), eph(hourstep), eah(hourstep)) !H-64,65,66 bhrwc0(hourstep) = theta(0)/bulkden(1) c create hourly and daily output values runoff = ( volw(2) - prevvolw(2) ) * mtomm bhzrun = bhzrun + runoff evap = ( volw(3) - prevvolw(3) ) * mtomm bhzea = bhzea + evap infil = ( volw(4) - prevvolw(4) ) * mtomm bhzinf = bhzinf + infil drain = ( volw(numeq) - prevvolw(numeq) ) * mtomm bhzper = bhzper + drain c update prevvolw to carry over to next hour or day do kindex=1,5 prevvolw(kindex) = volw(kindex) end do prevvolw(numeq) = volw(numeq) c output the hourly info here if ((am0hfl .eq. 2) .or. (am0hfl .eq. 3) .or. (am0hfl .eq. 6) .or. & (am0hfl .eq. 7)) then swc = sum(volw(iminlay:imaxlay)) * mtomm write(12,2050) hourstep,evap, runoff, drain, & swc,theta(0),(theta(kindex), kindex=1,layrsn) end if c Accumulating hourly soil wetness values do kindex=0,layrsn thetad(kindex) = thetad(kindex) + theta(kindex) if( (theta(kindex).le.0.001) .and. (kindex.gt.0) ) then write(*,*)'darcy:end theta<0.001',kindex,theta(kindex) call dvolw(neq,t,volw,netflux) write(*,*) 'lay',kindex,':',t,netflux(kindex+5), & netflux(kindex+6), netflux(3) write(*,*) 'tcur,hu:',rwork(13), rwork(11) c stop endif end do if( (raindepth.gt.0.0) & .and.(tout.ge.rainend) & .and.(bhzwid.le.0.0) ) then bhzwid = store( iminlay, imaxlay, prevvolw, volw, bszlyd ) endif hourstep = hourstep + 1 c------------------------------------- c If not done yet, take another step. c------------------------------------- if(hourstep.le.24) goto 30 C completed the day c this section should be enabled to extend solution over c multiple days when no outside process changes water contents c if( beginday .lt. 100*daytosec ) then c beginday = beginday + daytosec c else c istate = 1 !-- initializes solution routines c beginday = 0.0 c t = beginday c do kindex=1,5 c prevvolw(kindex) = 0.0 c end do c prevvolw(numeq) = 0.0 c end if c calculating average daily soil wetness do kindex=0,layrsn thetad(kindex) = thetad(kindex)/24 end do 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(kindex), kindex=0, layrsn) write(12,2060) write(12,2070) bhzep write(12,2080) bhzea write(12,2090) bhzper end if c if(raindepth.gt.0.0) write(*,*)'rain:p,i,r,p-i-r', c & volw(1)*mtomm,volw(4)*mtomm,volw(2)*mtomm, c & (volw(1)-volw(4)-volw(2))*mtomm c if( daysim.eq.6321 ) then c if(raindepth.gt.0.0) write(*,*)'rain:p,i,r', c & volw(1)*mtomm,volw(4)*mtomm,volw(2)*mtomm c write(*,*) 'daysim,steps,func,jac,ord:',daysim, c & iwork(11),iwork(12),iwork(13),iwork(14) c end if return end