subroutine infpar(nowcrp,jdplt,plant,irsyst) c c + + + PURPOSE + + + c Calculates sat. hydraulic conductivity and effective matric c potential from : c 1) bare soil hyd. cond. c 2) avg. potential across wetting front c 3) effective porosity c 4) percent ground cover c 5) percent canopy cover c 6) relative effective saturation c c Called from SOIL c Author(s): Savabi,Risse,Zhang c Reference in User Guide: Chapter 4 c c NOTE: Computation of the fraction of soil surface covered by both c canopy and ground cover (COVU) assumes location of ground c cover and canopy cover are independent. Reza Savabi says c this is correct. I would expect the location of the residue c to be somewhat dependent on the location of the plants that c generated it. -- CRM (9/14/92 conversation with R. Savabi) c c Changes: c 1) Common block SOLVAR was not used. It was de-referenced. c 2) The generic "SAVE", which saves ALL local variables, c was eliminated. c 3) Eliminated local variables TOTADJ, IPLUG, & NCOUNT. c 4) Introduced intermediate local variables TMPVR1 c to TMPVR4 to make calculations more efficient. c 5) The statement: c if (wetfrt.eq.tc) wetfrt=tc+.00001 c was changed to: c if(abs(wetfrt-tc) .lt. 0.00001) wetfrt=tc+.00001 c 6) Local variable SAT11 was computed, but the result was c never used. It was eliminated. c 7) Added local variable RFCUMX so INFPAR could tell if c rainfall had occurred. In my test data sets, this c eliminated about 80% of the executions of the code for c soil crusting adjustments. c 8) Changed statement at end of subroutine to prevent a c divide by zero occurring in SR ROCHEK - prevents the c value of sm(iplane) from becoming zero. dcf 8/16/93 c 9) Moved RFCUMX to common block cifpar.inc jca2 8/31/93 c 10) Added new Ksat adjustments from Risse and Zhang c dcf 1/11/94 c 11) Changes made to Ksat adjustments of Risse and c Zhang - dcf 2/4/94 c 12) Changes made to Ksat adjustment equations again c by Nearing - dcf 3/8/94 c 13) Change to Ksat adjustment equations for established c perennials from John Zhang - dcf 5/26/94 c 14) Change to Ksat adjustment equations for surface c cover(both residue and canopy) adjustments to c hydraulic conductivity and also for first year c of perennial growth. dcf 12/14/94 c 15) Change to exclude canopy cover factor in the c adjustment for hydraulic conductivity for the case c of furrow irrigation water addition. dcf 12/14/94 c 16) Change statement which checks for water content to c be above upper limit so that saturation can occur for c a single storm simulation. savabi and dcf 4/95 c FROM: if (st(i,iplane).ge.ul(i)) then c TO: if (st(i,iplane).ge.0.95*ul(i)) then c c Version: This module recoded from WEPP Version 92.25. c Date recoded: 09/03/92 & 10/20/92. c Recoded by: Charles R. Meyer. c c + + + KEYWORDS + + + c c + + + PARAMETERS + + + include 'pmxcrp.inc' include 'pmxelm.inc' include 'pmxhil.inc' include 'pmxtil.inc' include 'pmxtls.inc' include 'pmxnsl.inc' include 'pmxpln.inc' include 'pmxpnd.inc' include 'pmxres.inc' include 'pmxsrg.inc' include 'pntype.inc' include 'ptilty.inc' include 'pxstep.inc' c c + + + ARGUMENT DECLARATIONS + + + integer nowcrp, jdplt, plant, irsyst c c c + + + ARGUMENT DEFINITIONS + + + c nowcrp - index of current crop on the OFE c jdplt - julian planting date of current crop on OFE c plant - index of current plant type on the OFE c irsyst - the type of irrigation system: 0 - none c 1 - sprinkler c 2 - furrow c c + + + COMMON BLOCKS + + + c include 'cavepar.inc' c modify: avsat(mxplan), avbd(mxplan) c include 'cclim.inc' c read: tave c include 'ccons.inc' c read: cpm(mxnsl,mxplan) c modify: avpor(mxplan) c include 'ccover.inc' c read: cancov(mxplan), lanuse(mxplan) c modify: gcover(mxplan) c include 'ccrpout.inc' c read: rescov(mxplan), bd(mxnsl,mxplan), rtmass(mxplan) c include 'ccrpvr1.inc' c read: rtm(mxres,mxplan) c include 'cgcovr.inc' c read: gcvplt(mxplan),mantyp(mxplan) c modify: gcvplt(mxplan) c include 'chydrol.inc' c read: rain(mxplan),prcp c include 'cirfurr.inc' c modify: dtheta c include 'cke.inc' c read: rkecum(mxplan) c include 'cparame.inc' c write: ks(mxplan), sm(mxplan) c include 'cperen1.inc' c read: imngmt(mxcrop,mxplan),rtmmax(ntype) c include 'cstmflg.inc' c read: norain(mxplan) c include 'cstruc.inc' c read: iplane c include 'ctemp.inc' c read: avclay(mxplan), avsand(mxplan) c include 'ctillge.inc' c read: tillay(2,mxplan) c include 'cwater.inc' c read: st(mxnsl,mxplan), ul(mxnsl,mxplan), thetdr(mxnsl,mxplan), c nsl(mxplan), dg(mxnsl,mxplan) c c include 'cifpar.inc' c modify: rfcumx c include 'cwint.inc' c cd Added by S. Dun, Jan 08, 2008 include 'cupdate.inc' cd End adding c c + + + LOCAL VARIABLES + + + real avsm15, avks, sf, wetfrt, a, rra, bbbb, tc, crust(mxplan), 1 eke, sc, cke, crstad, kbare, ktmp, ccovef, scovef, kcov 1 ,fzul,ffi,tmpvr1,tmpvr2,tmpvr3,tmpvr4,avcpm integer i,tpfzfg c c 1 eke, sc, cke, crstad, kbare, ktmp, rtmt, rtmtef, ccovef, c 1 scovef c c + + + LOCAL DEFINITIONS + + + c avsm15 - average 1500 KPa (15 bar) soil water content c avks - average saturated hydralic conductivity for the tillage c layer c sf - matric potential across wetting front (m) c wetfrt - average depth of wetting front (m) c cf - canopy cover adjustment for saturated hydraulic c conductivity c a - macroporosity adjustment for saturated hydraulic c conductivity c bareu - bare area under canopy (fraction) c bareo - bare area outside canopy (fraction) c covu - ground cover under canopy (fraction) c covo - ground cover outside canopy (fraction) c tc - crust thickness c crust - crust adjustment for saturated hydraulic conductivity c eke - effective hydraulic conductivity in fill layer c (m/sec) c sc - reduction factor for subcrust hydraulic conductivity c avcpm - average rock fragment correction factor for the tillage c layer c grdcov - ground cover value used in macroporosity calculations c for cropland annuals - assigned a value on date of c last planting. For perennials and range - use actual c cover values c cke - coefficient relating amount of kinetic energy since last c tillage to the speed of crust formation c rtmt - total dead and live root mass in top 15 cm of soil c (kg/m**2) (NOT USED) cc rtmtef - transformed (live plus dead) root mass c ccovef - effective canopy cover corrected for the effect of c canopy height c scovef - effective total surface cover c ktmp - same as AVKS, only in units of (mm/hr) c kbare - effective AVKS after adjustment for crusting/tillage c kcov - portion of equation to compute effective AVKS c (surface cover adjustment) c tpfzfg - top layer frozen flag c 0: no influence of top frozen layer on infiltration c 1: for there is influence of top frozen layer c c + + + SAVES + + + save crust c c + + + DATA INITIALIZATIONS + + + c c + + + END SPECIFICATIONS + + + c c thickness. The layers assumed to affect infiltration are the c primary (deepest), and [the average of] the secondary tillage c layers.) c avbd(iplane) = 0.0 avpor(iplane) = 0.0 avsm15 = 0.0 avcpm = 0.0 tpfzfg = 0 c do 10 i = 1, 2 c ------ ratio of tillage layer thick. to primary (deepest) til. depth tmpvr1 = dg(i,iplane) / tillay(2,iplane) c ------ bulk density avbd(iplane) = avbd(iplane) + bd(i,iplane) * tmpvr1 c ------ porosity for infiltration avpor(iplane) = avpor(iplane) + por(i,iplane) * tmpvr1 c ------ 15-bar soil water content avsm15 = avsm15 + thetdr(i,iplane) * tmpvr1 c ------ correction factor for rock fragments avcpm = avcpm + cpm(i,iplane) * tmpvr1 10 continue c c Range checks: if (avpor(iplane).gt.1.0) avpor(iplane) = 1.0 if (avpor(iplane).lt.0.0) avpor(iplane) = 0.0 if (avbd(iplane).lt.800.0) avbd(iplane) = 800.0 if (avbd(iplane).gt.2200.0) avbd(iplane) = 2200.0 c c ---- Calculate the harmonic mean of Ks in the tillage layer (AVKS) c for the infiltration routine. c cd Modified by S. Dun, April 17, 2008 c for infiltration c if ((frozen(1,iplane) .gt. 0.0).and. (sscv(1,iplane).gt.0.0) 1 .and.(sscv(2,iplane).gt.0.0) ) then c avks = solthk(2,iplane) / (solthk(1,iplane)/sscv(1,iplane) 1 + (solthk(2,iplane)-solthk(1,iplane))/sscv(2,iplane)) c eke = avks tpfzfg = 1 c else avks = solthk(2,iplane) / (solthk(1,iplane)/ssc(1,iplane)+( 1 solthk(2,iplane)-solthk(1,iplane))/ssc(2,iplane)) endif cd End Modifying c c ---- Calculate average water content in tillage layer (AVSAT) c for the infiltration routine. (ST is constant > 15 bars) avsat(iplane) = (st(1,iplane)+st(2,iplane)) / tillay(2,iplane) + 1 avsm15 if (avsat(iplane).gt.avpor(iplane)) avsat(iplane) = 1 avpor(iplane) * .98 c c ---- Compute the matric potential of the infiltration zone (SF) c (WEPP Equation 4.3.2 ff) c -- XXX -- This equation needs a *number* in the User Document! c (See top of p 4.5) -- CRM -- 9/14/92. tmpvr2 = avclay(iplane) ** 2 tmpvr3 = avsand(iplane) ** 2 tmpvr4 = avpor(iplane) ** 2 c sf = 0.01 * exp(6.531-7.33*avpor(iplane)+15.8*tmpvr2+3.81*tmpvr4+ 1 avsand(iplane)*(3.4*avclay(iplane)-4.98*avpor(iplane))+tmpvr4* 1 (16.1*tmpvr3+16.0*tmpvr2)-14.0*tmpvr3*avclay(iplane)- 1 avpor(iplane)*(34.8*tmpvr2+8.0*tmpvr3)) if (sf.gt.0.5) sf = 0.5 c if (tpfzfg .ne. 1) then c no influence of top frozen layer on infitration. c when top frozen layer influence exists, sealing adjudtment is skipped c c *** L0 IF *** c CROPLAND if (lanuse(iplane).eq.1) then c c ------ ground cover (fraction) gcover(iplane) = rescov(iplane) c c ------ Compute average depth of the wetting front (meters) c (WEPP Equation 4.3.11) wetfrt = 0.147 - 0.15 * tmpvr3 - (0.0003*avclay(iplane)* 1 avbd(iplane)) if (wetfrt.lt.0.005) wetfrt = 0.005 c c ------ crust thickness tc = 0.005 c c *** L0 ELSE-IF *** c RANGELAND else if (lanuse(iplane).eq.2) then c c ------ average wetting front depth (meters) c (WEPP Equation 4.3.11) wetfrt = 0.147 - 0.15 * tmpvr3 - (0.0003*avclay(iplane)* 1 avbd(iplane)) if (wetfrt.lt.0.01) wetfrt = 0.01 c c ------ crust thickness tc = 0.01 c c ------ canopy cover adjustment factor for sat. hydraulic conductivity c (WEPP Equation 4.3.12) c cf = 1.0 + cancov(iplane) c c ------ macroporosity adjustment factor for sat. hydraulic conductivity a = exp(6.10-(10.3*avsand(iplane))-(3.7*avclay(iplane))) if (a.lt.1.0) a = 1.0 if (a.gt.10.0) a = 10.0 c c *** L0 ELSE *** c FOREST else c c (This branch left intentionally blank.) c c *** L0 ENDIF *** end if c c c *** M0 IF *** c If rainfall or tillage has occurred, compute saturated c hydraulic conductivity adjustment for crusted soil surface. c if (rfcum.ne.rfcumx) then c c ------ correction factor for partial saturation of the sub-crust layer c (WEPP Equation 4.3.8) sc = 0.736 + (0.19*avsand(iplane)) c c ------ compute maximum potential crust adjustment fraction c (WEPP Equation 4.3.7) if (abs(wetfrt-tc).lt.0.00001) wetfrt = tc + .00001 c c Use new equation for forming crust - see Rawls et al. 1989 c Note: This is now maximum adjustment ffi = 45.19 - (46.68*sc) crust(iplane) = sc / (1.0+(ffi/(wetfrt*100))) if (crust(iplane).lt.0.20) crust(iplane) = 0.20 c c ************************************************************ c * Note: CRUST is multiplied times Ksat. When there is no * c * crust, Ksat is not adjusted; ie, CRUST = 1. CRUST * c * is computed in 2 parts. The second part is the * c * adjustment for cumulative rainfall since tillage. * c * It is performed ONLY if RF since tillage is less * c * than 0.1 m * c ************************************************************ c c-----for CROPLAND, if cumulative RF since tillage is less than c 1/10 meter, update crust reduction factor for rainfall. c New equation inserted Risse c crust adjustment=maxadj+(1-maxadj) exp(-C kecum (1-rr/4)) c where maxadj is the same as it was in previous version with c correction for units(crust), C is calculated based on analysis c by Risse, and instead of using a linear relationship between 0 c and 100mm of rfcum we used the exponential relationship of c Brakensiek and Rawls, 1983 c (WEPP Equation 4.3.10) c if(lanuse(iplane).eq.1) then c Original Code: c 1 crust=1.-(((1.-crust)/0.1)*rfcum) c 1 crust = 1.0 - (1.0-crust)*10.0*rfcum c cke is coeffient relating to speed with which crust forms c based on analysis of natural runoff plot data by Risse c c-----for CROPLAND, if cumulative RF since tillage is less than c 1/10 meter, update crust reduction factor for rainfall. c if (lanuse(iplane).eq.1) then cke = -0.0028 + 0.0113 * avsand(iplane) + 0.125 * 1 avclay(iplane) / cec1(1,iplane) if (cke.gt.0.01) cke = 0.01 if (cke.lt.0.0001) cke = 0.0001 c make sure random roughness does not cause positive exponent if (rrc(iplane).le.0.04) then rra = rrc(iplane) else rra = 0.04 end if c bbbb = -cke * rkecum(iplane) * (1-rra/0.04) if (bbbb.lt.-25.0) bbbb = -25.0 c Calculate crusting/tillage adjustment crstad = crust(iplane) + (1-crust(iplane)) * exp(bbbb) else crstad = 1.0 end if c c *** M0 ENDIF *** c endif c c Adjust AVKS for soil surface characteristics (canopy & cover) c and for dead and live roots. (from zhang 1/94) dcf c kbare = avks * crstad ktmp = avks * 3.6e6 c c EXCLUDE the adjustment for canopy cover for the case of c snow melt or furrow irrigation. dcf 12/15/94 c if(tave.gt.0.0)then if(prcp.gt.0.001 .or. irsyst.eq.1)then c c Adjust the effectiveness of canopy cover by canopy height c ccovef = cancov(iplane) * exp(-0.3358*canhgt(iplane)/2.0) c c Calculate the total effective surface cover c scovef = ccovef + rescov(iplane) - ccovef * rescov(iplane) c else scovef = rescov(iplane) endif else scovef = rescov(iplane) endif c c Calculate the final effective conductivity c c IF the user has indicated that he/she wants the internal c Ksat adjustments used in the SOIL input file - then c adjust final effective conductivity for crusting/tillage/ c crop/rainfall dcf 1/11/94 c if (ksflag.eq.1) then c if(tave.gt.0.0)then c c XXX NOTE - We should really add in the amount of sprinkler c XXX irrigation water here (if any exists for the day) c XXX but unfortunately at this point in the program we c XXX do not yet know what this amount will be (have not c XXX call subroutine IRRIG yet). dcf 12/15/94 c kcov = (0.0534 + 0.01179*ktmp) * prcp * 1000.0 / 3.6e6 else kcov = 0.0 endif c c Zhang change 12/9/94 c if(kcov .lt. avks)kcov = avks if(kcov .lt. 0.5*avks) kcov = 0.5*avks c eke = kbare*(1.0 - scovef) + kcov*scovef c c If crop adjusted eke is smaller than that of crust adjusted c set it back to crust adjusted value c if (eke.lt.kbare) eke = kbare c c c ADJUST FOR ESTABLISHED PERENNIAL CROP (meadows, etc.) c c Note - Changed coefficient in equation below from 1.7965 to 1.81 c when given change from Zhang. 7/1/94 dcf c if(lanuse(iplane).eq.1)then c c Changes to include perennial adjustment for first c year of perennial growth when plant is sufficiently c developed. dcf 12/9/94 c if(imngmt(nowcrp,iplane).eq.2 .and. jdplt.eq.0) c 1 eke = 1.81 * eke c if(imngmt(nowcrp,iplane).eq.2)then if(jdplt .eq. 0)then eke = 1.81 * eke else c c If the root mass present for the perennial exceeds c half of the input maximum possible, then adjust the c EKE value by the full amount. The 0.5 value is c arbitrary and could be changed here and below. dcf c if(rtmass(iplane).gt.0.5*rtmmax(plant))then eke = 1.81*eke else c c Linearly change the eke adjustment based upon the c ratio of root mass present to half the maximum c root mass possible (input) for the perennial. The c use of the 0.5 value is arbitrary and could be c changed here and above. dcf 12/9/94 c eke = eke + (1.62*eke/rtmmax(plant))*rtmass(iplane) endif endif endif endif else eke = avks end if c Added by S. Dun, May 17, 2008 c for considering the influence of top frozen layer on infiltration endif c end adding cd Modified by S. Dun 06/20/2002 c LIMIT MINIMUM KSAT TO 1.94E-08 m/s (0.07 mm/h) for pre c 2006 soil files. c if (solwpv.eq.2006) then if (eke.le.1.0e-14) eke = 3.0e-14 else if (eke.le.1.94e-08) eke = 1.94e-08 endif c cd adjust the lower limit to e-14 m/s cd (the reference we are using is "Physical and Chemical Hydrogeology" cd by P.A. Domenico and F.W. Schwartz) cd end modifying. c c c *** BEGIN N0 LOOP *** c c In case a restricted soil layer controls percolation and c infiltration.... c i = 0 20 continue i = i + 1 c ---- If the water content is above the upper limit for this layer.... c cd Modified by S. Dun, April 17, 2008 c for frozen soil effect cd if (st(i,iplane).ge.0.95*ul(i,iplane)) then c fzul =ul(i,iplane) - frzw(i,iplane) if (st(i,iplane).ge.0.95*fzul) then c ------ If this layer's Ksat is less than the average Ksat for the c plow layer .... if (ssc(i,iplane).le.eke) then eke = ssc(i,iplane) endif c if ((frdp(iplane) .gt. 0.0).and. (sscv(i,iplane).gt.0.0) 1 .and. (sscv(i,iplane).le.eke) ) then c eke = sscv(i,iplane) endif cd End modifying else c ------ (force exit from loop) i = nsl(iplane) end if c *** END N0 LOOP *** if (i.lt.nsl(iplane)) go to 20 c ks(iplane) = eke c c Write(60, 1505) day, mon, year, eke, kbare, scovef,crstad c 1 ,crust(iplane),bbbb 1505 format(1x, 3I6, 6E12.3) c c Compute effective matric potential (SM), correcting for rock c fragments (using AVCPM). c if (avsat(iplane).ge.(avpor(iplane)*avcpm)) avsat(iplane) = ( 1 avpor(iplane)*avcpm) * 0.99 c ---- compute water above field capacity dtheta(iplane) = avpor(iplane) * avcpm - avsat(iplane) c ---- compute effective matric potential (SM) sm(iplane) = dtheta(iplane) * sf c c rfcumx = rfcum c return end