c$Author: fredfox $ c$Date: 2001-09-27 16:22:31 $ c$Revision: 1.46.2.6 $ 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, 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, 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) c bhzper - accumulated daily drainage (deep percolation) (mm) c bhzrun - accumulated daily runoff (mm) c bhzwid - Water infiltration depth (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: 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 bhzea = 0.0 bhzper = 0.0 bhzrun = 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) evapamp = pi*bhzep*mmtom/lenday/2.0 raindepth = bwzdpt * mmtom !storm total depth (m) rainstart = 1.0 rainend = rainstart + max(2.0, bwdurpt * hrtosec) rainmid = rainstart & + max(rainstart,min(rainend,bwpeaktpt*bwdurpt*hrtosec)) 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 if( daysim .eq. 1 ) then do kindex=1,5 prevvolw(kindex) = 0.0 end do prevvolw(numeq) = 0.0 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(1:layrsn)) * 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) = 2000 !maximum number of steps allowed before error generated if( daysim .eq. 1 ) then istate = 1 beginday = 0.0 t = beginday 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,'2k 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 write(*,*) ' Failed at t = ',t,' with istate = ',istate stop end if goto 40 end if c------------------------------------- c step was sucessful, continue on c------------------------------------- c call subroutine to get flux rate values c neq(2) = 1 c call dvolw(neq,t,volw,netflux) c neq(2) = 0 c if( netflux(1).lt.0.0 ) then c write(*,*)'lay 1:',t,volw(6)/tlay(6),netflux(6), c & netflux(4), netflux(5) c endif 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 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 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 if( beginday .lt. 100*daytosec ) then beginday = beginday + daytosec else istate = 1 !-- initializes solution routines beginday = 0.0 t = beginday do kindex=1,5 prevvolw(kindex) = 0.0 end do prevvolw(numeq) = 0.0 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', c & volw(1)*mtomm,volw(4)*mtomm,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