!****************************************************************************** !Reference in User Guide: Chapters 4, 6, 7, and 8. ! Also see "A Compendium of Soil Erodibility Data from ! WEPP Cropland Soil Field Erodibility Experiments 1987 ! & 88", NSERL Report No. 3. ! Also see "EPIC -- Erosion Productivity Impact Calcu- ! lator (EPIC Model Documentation) printed Sept. 1990. ! Also see "Microclimate, the Biological Environment" ! by Norman J. Rosenberg, et. al. 1983. ! Also see "Soil Cover Prediction with Various Amounts ! and Types of Crop Residue". 1982. James M. Gregory. ! Trans. ASAE. pp 1333-1337. !****************************************************************************** subroutine watbal(fwidthPrev, slplenPrev, runoffPrev,efflenPrev, sbrunfPrev, & & nowcrp,norun,iflag,nsl,cntflg, & & iflget, idrain, sdate, & & itype, elevm, & & s1, s2, surdra, resint, plaint, & & sbrunf,ul,fc,hk,soilw, & & sep,pintlv, st, ep, unsdep, & & watstr, satdep, drainq, runoff, drawat, & & es, por, thetdr, dg, thetfc, & & solthk, rtd, lai, fin, rain, & & wmelt, irdept, iraplo, roffon, fwidth, & & slplen, efflen, tillay, cv, rmagt, & & rmogt,coca, ssc, radinc, ddrain, & & ytn, daylen, canhgt, & & iwind, mon, j1, & & j2, tave, salb, radly, obmaxt, & & obmint, eo, tdpt, radpot, vwind, & & tu, su, tv, et, tmax, & & tmin, kcb, kcbcon, rawp, drseq, & & sdrain, drdiam, drainc, pltol, ul4) implicit none include 'constants.inc' real :: watcon, fcdep !, totlqc - probably need to be saved/passed !************************************************************************************************************* ! fwidthPrev - width of field (m) for previous OFE ! slplenPrev - slope length (m) for previous OFE ! runoffPrev - daily runoff amount (m) for previous OFE ! efflenPrev - effective flow length of a previous overland flow element for continuous flow planes - efflen is the sum of ! the plane lengths having flow - for a Case 4 plane on which runoff ends - efflen is the length of the ! portion of the Case 4 OFE on which runoff is present. ! sbrunfPrev - daily subsurface lateral flow (m) for previous OFE ! nowcrp - ? ! norun - flag for runoff occurrence ! iflag - flag for hillslope initialization calls ! nsl - number of soil layers ! cntflg - flag indicating contour farming is used. ! iflget - type of ET method (1 - WEPP original Penman, 2 - Penman_Monteith) ! idrain - flag for tile drainage (0 - no tile 1 - tile drainage) ! sdate - date of year in Julian date ! itype - plant type, plant growth scenario index(?) ! elevm - elevation of the climate station in meters ! s1 - stage 1, soil evaporation ! s2 - stage 2, soil evaporation ! surdra - surface drainage water (m) ! resint - interception of rain water by flat residue(m) ! plaint - interception of rain water by live plants (m) ! sbrunf - daily subsurface lateral flow (m) ! ul - upper limit of water content per soil layer ! fc - soil field capacity (m) ! hk - a parameter that causes SC approach zero as soil water approaches FC ! soilw - soil water content per layer ! sep - seepage ! pintlv - plant interception left over from the day before ! st - current available water content per soil layer (m) ! ep - plant transpiration (m/day) ! unsdep - unsaturated depth from surface to water table (m) ! watstr - water stress parameter for plant growth (0-1) ! satdep - ? ! drainq - drainage spacing, m ! runoff - daily runoff amount (m) ! drawat - potential water which can be drained from a soil layer ! es - soil evaporation (m/day) ! por - ? ! thetdr - 15-bar soil water content (wilting point) ! dg - depth of each soil layer, meters ! thetfc - 1/3-bar soil water content (field capacity) ! solthk - cumulative thickness of soil layer (m) ! rtd - root depth (m) ! lai - leaf area index ! fin - infiltrated water amount ! rain - daily rainfall amount (m) ! wmelt - amount of snow melt occurring in a day on an OFE ! irdept - depth of water applied by irrigation (m) ! iraplo - irrigation water volume supplied to an overland flow element per unit area (m) ! roffon - ? ! fwidth - width of field (m) ! slplen - slope length (m) ! efflen - effective flow length of an overland flow element for continuous flow planes - efflen is the sum of ! the plane lengths having flow - for a Case 4 plane on which runoff ends - efflen is the length of the ! portion of the Case 4 OFE on which runoff is present ! tillay - depth of secondary (tillay (1)) and primary (tillay (2)) tillage for the current crop on the OFE ! cv - residue amount ! rmagt - surface residue mass above ground today (kg/m^2) ! rmogt - surface residue mass on the ground today (kg/m^2) ! coca - entrapped air correction factor ! ssc - saturated hydraulic conductivity ! radinc - average slope in degree ! ddrain - depth of drain tiles from surface, m ! ytn - funtion of latitude - used in day length calculation (EPIC 2.194) ! daylen - day length (hours) ! canhgt - canopy height (m) ! iwind - wind data flag in climate input file ( 0 - wind data, 1 - no wind data ) ! mon - month of year ! j1 - soil layers subjected to soil evaporation ! j2 - soil layers(?) ! tave - average daily temperature (degrees C) ! salb - soil albedo (0-1) ! radly - daily solar radiation (Langleys/day) ! obmaxt - observed average monthly maximum temperature (degrees C) ! obmint - observed average monthly minimum temperature (degrees C) ! eo - potential daily evapotranspiration (m/day) ! tdpt - mean daily dew point temperature (degrees C) ! radpot - potential radiation at slope (cal/cm^2/day) ! vwind - mean daily wind speed (m/s) ! tu - upper limit soil evaporation, stage 1 ! su - soil water available for evaporation ! tv - ? ! et - ? ! tmax - maximum daily temperature (degrees C) ! tmin - minimum daily temperature (degrees C) ! kcb - mid-season crop coefficient for Penman-Monteith ! kcbcon - basal crop coefficient. ! rawp - coefficient p for readily available root zone ! drseq - drainage sequence ! sdrain - drainage spacing, m ! drdiam - drain tile diameter, m ! drainc - drainage coefficient, m/day ! ul4 - parameter to adjust potential water use by plants ! pltol - plant drought resistance factor (0 to 1), unitless ! nowcrp - current crop ! !************************************************************************************************ integer, intent(in) :: nowcrp,iflag,nsl,cntflg integer, intent(in) :: iflget, idrain, sdate integer, intent(in) :: itype(mxcrop) integer, intent(in) :: iwind, mon, j1, j2, drseq(mxnsl) real, intent(in) :: elevm,fwidthPrev,slplenPrev,runoffPrev,efflenPrev,sbrunfPrev real, intent(in) :: por(mxnsl), thetdr(mxnsl), dg(mxnsl) real, intent(in) :: thetfc(mxnsl), solthk(mxnsl), rain, wmelt real, intent(in) :: irdept, iraplo, roffon, fwidth, slplen, efflen real, intent(in) :: tillay(2), coca(mxnsl) real, intent(in) :: radinc, ddrain, ytn real, intent(in) :: tave, salb, radly, obmaxt(12), obmint(12), eo, tdpt real, intent(in) :: radpot, vwind, tu, su, tv, et real, intent(in) :: tmax, tmin, kcb(mxcrop), kcbcon, rawp(mxcrop) real, intent(in) :: sdrain, drdiam, drainc, pltol, ul4 integer, intent(inout) :: norun real, intent(inout) :: s1, s2, surdra, resint, plaint, sbrunf real, intent(inout) :: ul(mxnsl), hk(mxnsl), soilw(mxnsl), sep, fin, pintlv real, intent(inout) :: st(mxnsl), cv, ep, es, unsdep, drawat(mxnsl), satdep real, intent(inout) :: drainq, runoff, daylen, watstr, fc(mxnsl) ! the following are inout because of calls to other subs real, intent(inout) :: canhgt, lai, rmagt, rtd real, intent(inout) :: rmogt(mxres),ssc(mxnsl+1) ! ! + + + LOCAL VARIABLES + + + ! real xfin, soldep, sd, ch, h, tmpvr1, & & drfc(mxnsl), latqcc, latk, avcoca, avfca, avhk, avpora, & & avstt, avul, fffx, fslope, rm, sstz, subq, totdg, totk, & & totlqc, watyld,runoffin integer i, mn, jj, lflag ! ! + + + LOCAL DEFINITIONS + + + ! watcon - water content of the root zone ! xfin - water available to infiltrate into the current soil layer ! soldep - soil profile depth ! sd - sun's declination angle ! ch - ! h - ! tmpvr1 - ! drawat - ! drfc - ! latqc - ! latqcc - ! latk = hydraulic conductivity for the subsurface lateral flow calculations, m/day ! fcdep - ! runoffin - runoff flow into the segment. ! If this an initialization call, initialize constants. if (iflag.eq.0) then watcon = 0.0 s1 = 0.0 s2 = 0.0 surdra = 0.0 resint = 0.0 plaint = 0.0 totlqc = 0.0 latqcc = 0.0 subq = 0.0 sbrunf =0.0 do 10 i = 1, nsl ! upper limit of water content for current layer ul(i) = (por(i)-thetdr(i)) * dg(i) ! field capacity for current layer fc(i) = dg(i) * (thetfc(i)-thetdr(i)) ! Used in PERC to adjust sat. hyd. cond. on non-saturated soils. hk(i) = -2.655 / alog10(fc(i)/ul(i)) ! calculate total soil water in the root zone (m) soilw(i) = st(i) + (thetdr(i)*dg(i)) 10 continue end if sep = 0.0 soldep = solthk(nsl) ! If this is NOT an initialization call, proceed; otherwise, RETURN if (iflag.ne.0) then ! Note: Intercepted rainfall is subtracted from rainfall Reza 9/93 if((plaint.lt.1.0E-8).and.(pintlv.gt.1.0e-8)) then if((rtd.gt.0.0).and.(lai.gt.0.0)) then plaint = pintlv else resint = resint + pintlv endif endif fin = rain-(plaint-pintlv+resint)+ & & wmelt + irdept + iraplo pintlv = 0.0 ! For channel ! ??? Check impondment part later if (cntflg.eq.1) then fin = fin + roffon - runoff runoffin = roffon else fin = fin + sbrunfPrev & & *(fwidthPrev*slplenPrev) & & /(fwidth*slplen) & & + (runoffPrev*efflenPrev & & - runoff*efflen)/slplen runoffin = runoffPrev*efflenPrev/slplen endif if (fin.lt.0.0) then if (cntflg.eq.1) then fin = roffon - runoff else fin = (runoffPrev*efflenPrev & & - runoff*efflen)/slplen endif endif ! If there is infiltration, prorate water into the plow layer. if (fin.gt.0.0) then ! available water xfin = fin ! *** Begin L2-Loop *** ! Starting at top, infiltrate water into each tilled layer. i = 0 20 continue i = i + 1 if (solthk(i).lt.tillay(2)) then ! plow layer does not end in current soil layer, add water tmpvr1 = xfin * dg(i) / tillay(2) st(i) = st(i) + tmpvr1 xfin = xfin - tmpvr1 else ! plow layer ends in current soil layer, dump remaining ! water into this layer. st(i) = st(i) + xfin xfin = 0.0 end if ! *** End L2-Loop *** if ((i.lt.nsl).and.xfin.gt.0.0) go to 20 end if ! combined standing & fallen residue (used in EVAP to compute % bare soil) cv = rmagt do 30 i = 1, mxres cv = rmogt(i) + cv 30 continue ! percolation of water through soil layers call purk(nsl, st, fc, ul, hk, ssc, sep) ep = 0.0 es = 0.0 ! compute evapotranspiration (ET). ! S. Dun switched the evpotranspiration method to Penman-Monteith if (iflget.eq.1) then call evap(iwind,mon,j2,j1,nsl,elevm,cv,tave,lai, & & salb,radly,obmaxt,obmint,eo,tdpt,radpot,vwind, & & canhgt,resint,s1,tu,s2,fin,su,es,tv,ep,et,st,solthk) else call evappm(nowcrp,itype,nsl,j2,j1, & & elevm,radly,tave,tdpt,tmax,tmin,radpot,vwind,eo, & & lai,rtd, kcb, canhgt, kcbcon, solthk, thetfc, thetdr, & & st, fin, rawp, resint, cv, es, ep, et) endif ! prevent the soil water content of any layer from exceeding the ! upper limit for that layer. Add water to next higher layer. ! Moved to after lateral flow calculations and et - reza 7/25/93 do 40 i = nsl, 2, -1 if (st(i).gt.ul(i)) then st(i-1) = st(i-1) + (st(i)-ul(i)) st(i) = ul(i) end if 40 continue ! subsurface lateral flow calculation, Reza Savabi, 7/93 fcdep = 0.0 unsdep = 0.0 do 50 i = 1, nsl drfc(i) = fc(i) + ((1-coca(i))*dg(i)) if (st(i).ge.drfc(i).and.unsdep.eq.0.0) then fcdep = fcdep + dg(i) else unsdep = unsdep + dg(i) end if 50 continue latqcc = 0.0 subq = 0.0 sbrunf = 0.0 ! the big subsurface lateral flow if (fcdep.gt.0.0) then ! Calculate average saturated soil hydraulic conductivity in ! saturated zone, (below water table). avpora = 0.0 avfca = 0.0 avcoca = 0.0 totk = 0.0 totdg = 0.0 avstt = 0.0 avul = 0.0 avhk = 0.0 fffx = 1.0 lflag = 1 do 60 mn = 1, nsl if (st(mn).ge.drfc(mn).and.lflag.eq.1) then totdg = totdg + dg(mn) totk = totk + (ssc(mn)*dg(mn)) avpora = avpora + (por(mn)*(dg(mn)/fcdep)) avfca = avfca + (thetfc(mn)*(dg(mn)/fcdep)) avcoca = avcoca + (coca(mn)*(dg(mn)/fcdep)) avstt = avstt + (st(mn)*(dg(mn)/fcdep)) avul = avul + (ul(mn)*(dg(mn)/fcdep)) avhk = avhk + (hk(mn)*(dg(mn)/fcdep)) else if (totdg.gt.0.0) lflag = 0 end if 60 continue if (avul.gt.0.0) then sstz = avstt / avul if (sstz.lt.0.95) then fffx = sstz ** avhk if (fffx.lt.0.002) fffx = 0.002 else fffx = 1. end if end if fslope = sin(radinc) latk = 86400. * (totk/totdg) * fffx subq = fcdep * 1 * latk * fslope latqcc = subq /slplen sbrunf = latqcc if (latqcc.lt.0.0) latqcc = 0.0 ! calculate new hoo2 (depth of saturated depth at the bottom ! of hillslope) for next day (m) watyld = avpora - (avfca+(1.0-avcoca)) fcdep = fcdep - (latqcc/watyld) if (fcdep.lt.0.0) fcdep = 0.0 unsdep = soldep - fcdep ! Convert subsurface lateral flux to soil water depth using average ! drainable porosity (water yield). ! If there is lateral flow, latqcc, update the ! available water content (ST) in each soil layer to reflect ! the lateral flow from each layer, starting at the top layer. if (latqcc.gt.0.0) then do 70 jj = 1, nsl ! If layer is above field capacity, drain it. if (st(jj).gt.drfc(jj).and.latqcc.ge.0.0) then ! potential water which can be flow laterally from this layer drawat(jj) = st(jj) - drfc(jj) if (latqcc.gt.0.0) then ! If water excess in this layer exceeds or equals what has ! run from latqcc. if (drawat(jj).ge.latqcc) then ! subtract amount that has run out the drain from this layer. st(jj) = st(jj) - latqcc if (st(jj).lt.1e-10) st(jj) = 1e-10 latqcc = 0.0 ! If water excess in this layer is less than what has run ! from profile ... else ! subtract from latqcc what this layer can contribute. latqcc = latqcc - drawat(jj) ! adjust soil layer water content by same amount of water st(jj) = drfc(jj) end if end if end if 70 continue end if end if ! end subsurface flow calculations ! If unsaturated depth is greater than or equal to the depth ! of the drains than DRAIN will not be called and DRAINQ ! will be 0.0 satdep = 0.0 unsdep = 0.0 do 80 i = 1, nsl drfc(i) = fc(i) + ((1-coca(i))*dg(i)) if (st(i).ge.drfc(i).and.satdep.eq.0.0) then satdep = soldep - (solthk(i)-dg(i)) end if 80 continue unsdep = soldep - satdep drainq = 0.0 if (idrain.eq.1.and.unsdep.lt.ddrain) & & call drain(solthk(nsl),nowcrp, nsl, drseq, drainq, & & solthk, unsdep, satdep, drawat, st, ssc, dg, ddrain, sdrain, & & drdiam, drainc, por, thetfc, coca, fc) ! Compute surface drainage water (SURDRA). if (st(1).gt.ul(1)) then ! Problem currently exists of how to handle this surface ! drainage water computed - it should really be added back ! into the runoff for the current day - HOWEVER - this ! runoff has already been routed and erosion computations ! made at this point in the model. POSSIBLE SOLUTION: ! move the call to WATBAL to within SR IRS or CONTIN - ! check to see if SURDRA(iplane) > 0.0 for any OFE - if it ! is then somehow add this water as an input and go back ! and redo the IRS-GRNA computations. ANOTHER SOLUTION - ! discard this water amount as done currently. DCF - 5/21/93 surdra = st(1) - ul(1) ! Added by S. Dun. May 22, 2003, Modified Jan 28, 2004 ! Looks like it is added to runoff and then flag for rerun set. FAF if (cntflg.eq.1) then runoff = runoff + surdra else runoff = runoff + surdra * slplen/efflen endif if (surdra.gt. 1.0e-6) norun = 1 st(1) = ul(1) unsdep = unsdep-dg(1) else surdra = 0.0 end if ! QUESTION?? Is the following equation applicable to all ! world hemispheres????? (will equation work in ! Brazil or Australia???) dcf 5/21/93 ! Sun's declination angle (EPIC Equation 2.195) sd = 0.4102 * sin((sdate-80.25)/58.09) ! Code changed because of bug in Lahey compiler using tangent ! function on AT&T 6300 machine. dcf -- 10/1/91 ! Original code: ! ch = -ytn * tan(sd) ch = -ytn * sin(sd) / cos(sd) if (ch.ge.1.0) then h = 0.0 else if (ch.le.-1.0) then h = 3.1416 else h = acos(ch) end if ! day length (EPIC Equation 2.194) daylen = 7.72 * h ! If evapotranspiration and root depth are positive.... if (ep.gt.0.0.and.rtd.gt.0.0) then ! estimate evaporation and plant water use call swu(nowcrp,nsl,itype,rtd,solthk,plaint,ep,pintlv, & & pltol, ul4, ul, st, watstr) ! Added by S. Dun July 10,2003 for the havest day's water balance. ! Probably we should have a better way to solve this problem. elseif(plaint.gt.0.0) then if (plaint.gt.ep) then pintlv = plaint-ep else ep = plaint pintlv = 0.0 endif plaint = 0.0 else ! no evapotranspiration ep = 0.0 ! no water stress QUESTION - do we want to do this ? dcf watstr = 1.0 end if ! variable waty is used below in water balance checking ! code which is currently commented out dcf 5/21/93 ! waty=watcon watcon = 0.0 do 90 i = 1, nsl ! soil water content by layer soilw(i) = st(i) + (thetdr(i)*dg(i)) ! water content of the root zone watcon = watcon + soilw(i) 90 continue ! check the water balance (equation 7.1) ! watsm = abs( rain + wmelt + irdept + irapld & ! & - (ep + es + sep + runoff)) ! watdel = watcon - waty ! if(watdel.gt.0.0)then ! watdif = watdel - watsm ! else ! watdif = watdel + watsm ! endif ! write (6,*)'watdif',watdif ! if(abs(watdif).ge.0.001.and.sdate.gt.1)then ! write (6,*)' ' ! write (6,*)'IN SR WATBAL - watdif = ',watdif ! write (6,*)'SDATE = ',sdate,' rain= ',rain ! write (6,*)'iplane = ',iplane,' wmelt=',wmelt ! write (6,*)'irdept,irapld = ',irdept,irapld ! write (6,*)'ep = ',ep,'es = ',es ! write (6,*)'sep = ',sep,'runoff = ',runoff ! write (6,*)' ' ! endif rm=(rain+wmelt+irdept+iraplo)*1000. end if return end subroutine