subroutine route c c + + + PURPOSE + + + c Calculates detachment (when shear stress exceeds critical shear c stress) or deposition at the upper end of the slope segment and c routes sediment through the hillslope profile. c c Called from subroutine CONTIN. c Author(s): G. Foster, M. Nearing, D. Flanagan c Reference in User Guide: c c Changes: 1) Since MSHEAR is a discrete integer variable (flag), c computed GOTO's were substituted for IF's with .LE.'s c to handle cases in "deposition at upper end of segment" c section, and in "detachment at upper end of segment" c section. This seemed to reduce the runtime by about c 3-5% (A 9.9-10.0 second run was reduced to 9.5 seconds.) c 2) Variable ERD was deleted. c 3) KTRATO & QOSTAR were added to the argument lists of c functions DEPEND & DEPC, to eliminate the need for c them to access common blocks. c 4) In the third call to EROD after 155, the 6th argument c was changed from EATA to 0.0 (Dennis's Suggestion). c CRM - 1/22/91. c 5) Re: the values of MSHEAR, since 2 & 3 were lumped c together, and a computed GOTO can't handle a zero, c the original "3" was eliminated, and 0-2 were "pro- c moted" to 1-3. CRM - 1/24/91. c 6) GOTO's were massively changed to nested IF's. c CRM - 2/04/91. c 7) Removed old L2 IF-ELSE logic (v97.3) so that model c always calculates where deposition ends, regardless c of the shape of the segment. This is because with the c adjustments made in PARAM at OFE breaks - slope shape c is not always true predictor of processes on OFE. c dcf (changed L3 to L2 IF/ELSE) c c Version: This module recoded from WEPP version 90.92. c Date recoded: 12/21/90 - 2/08/91. c Recoded by: Charles R. Meyer. c c + + + KEYWORDS + + + c c + + + PARAMETERS + + + include 'pmxelm.inc' include 'pmxhil.inc' include 'pmxnsl.inc' include 'pmxpln.inc' include 'pmxprt.inc' include 'pmxslp.inc' include 'pmxtls.inc' c c + + + COMMON BLOCKS + + + include 'cends1.inc' include 'cends2.inc' include 'cends3.inc' c read: ktrato, qout, qin, strldn include 'cenrpa1.inc' c modify: frcflw(10,mxplan) include 'cerdva1.inc' c write: tc(101),load(101) include 'cinfco1.inc' c read: qostar include 'cinfco2.inc' c read: ainf(mxslp), binf(mxslp), cinf(mxslp) include 'cirriga.inc' c read: noirr include 'cpart1.inc' c read: npart, frac include 'cparva1.inc' c read: eata, tauc, theta, phi include 'cslope1.inc' c read: xl(mxslp,mxplan), nslpts(mxplan), a(mxslp,mxplan) c modify: xu(mxslp,mxplan) include 'cstruc.inc' c read: iplane c c + + + LOCAL VARIABLES + + + real ldlast, loadup, lddend, xdend, xdbeg, xdetst, cdep, dl, du, 1 xc1, xc2 integer i, iendfg, ilast, jj, k, mshear c c ldlast - n.d. sediment load calculated at last point c loadup - n.d. load computed entering deposition region - value c is sent and used in subroutine enrich.f c lddend - n.d. load computed at end of previous deposition c region or overland flow element (i.e., load at xdetst) c xdend - n.d. distance where deposition computed to end c xdbeg - n.d. distance where deposition begins on a segment c xdetst - n.d. distance at start of previous detachment area c ilast - counter index variable at last point c iendfg - flag to subroutine enrich.f indicating last call at c end of an overland flow element c dl - n.d. deposition rate c du - n.d. deposition rate at the top of a segment c cdep - portion of solution to deposition equation c erd - n.d. rill erodibility parameter sent to erod.f c if shear < critical erd=0, otherwise erd=eata c mshear - flag indicating what shear conditions exist on segment c c + + + SUBROUTINES CALLED + + + c xcrit c depos c erod c c + + + FUNCTION DECLARATIONS + + + real depc real depend c c + + + END SPECIFICATIONS + + + c c c********************************************************************* c c ... initialize erosion variables at the top of an OFE. c lddend = strldn ldlast = strldn ilast = 1 ndep = 0 iendfg = 0 xdbeg = 0.0 xdetst = 0.0 c c .... Set sediment loads and transport capacity at each point, c equal to zero. do 10 jj = 1, 101 load(jj) = 0.0 tc(jj) = 0.0 10 continue c c .... Set distance, transport capacity, and sediment load at c first point on OFE. xu(2,iplane) = 0.0 c tc(1)=cinf(2)*ktrato tc(1) = cinftc(2) * ktrato load(1) = strldn c c ... initialize particle fractions in the flow at top of OFE c with no inflow from previous OFE. c if (qout.gt.0.0) then if (qin.le.0.0.or.(noirr.gt.0.and.irsyst.eq.2.and.iplane.eq.1)) 1 then do 20 i = 1, npart frcflw(i,iplane) = frac(i,iplane) 20 continue else c c ........... Initialize particle fractions in the flow for OFE's c with inflow from previous OFE. do 30 i = 1, npart frcflw(i,iplane) = frcflw(i,iplane-1) 30 continue end if else do 40 i = 1, npart frcflw(i,iplane) = 0.0 40 continue end if c c****************************************************************** c upper boundary condition for overland flow elements c c Determine if deposition is occuring at X=0 on OFE. c ...if qostar = 0.0 (no inflow) estimate dl from deposition equation c ...else, estimate dl from the incoming sediment load and flow rate c if (abs(qostar).lt..0011) then dl = phi / (phi+1.0) * (ktrato*binftc(2)-theta) else dl = phi / qostar * (ktrato*cinftc(2)-ldlast) end if c c c ... ************************************************************** c ... *** For Each Slope Segment Within an Overland Flow Element *** c ... ************************************************************** c ... perform detachment and deposition calculations. c c c *** Start of BIG DO-LOOP *** do 170 k = 2, nslpts(iplane) c c For a CASE 4 plane - bypass all calculations if flow c has ended before current segment c *** Start of BIG IF *** if (qout.gt.0.0.or.xu(k,iplane).lt.-qostar) then c c If a Case 2 or 3 plane, or a Case 4 plane segment c on which runoff does not end. if (qout.gt.0.0.or.(qout.le.0.0.and.xl(k,iplane).lt.-qostar)) 1 then c c .......... Calculate shear conditions in segment, then c Find where shear equals critical shear for segment. call xcrit(ainf(k),binf(k),cinf(k),tauc,xu(k,iplane), 1 xl(k,iplane),qostar,xc1,xc2,mshear) c c For a Case 4 plane segment on which the flow ends else call xcrit(ainf(k),binf(k),cinf(k),tauc,xu(k,iplane),- 1 qostar,qostar,xc1,xc2,mshear) end if c c ... determine if there is deposition at the beginning of the c ... segment - if there is - calculate where deposition ends c du = dl c c *** L1 IF *** if (du.lt.0.) then c c deposition at upper end of segment c cdep = depc(xu(k,iplane),ainftc(k),binftc(k),phi,theta,du, 1 ktrato,qostar) c c check for deposition ending within segment c xdend = depend(xu(k,iplane),xl(k,iplane),ainftc(k), 1 binftc(k),cdep,phi,theta,ktrato,qostar) c c deposition does not end c c *** L2 IF *** if (xdend.ge.xl(k,iplane)) then xdend = xl(k,iplane) c loadup = ldlast call depos(xu(k,iplane),xdend,cdep,ainftc(k),binftc(k), 1 cinftc(k),phi,theta,ilast,dl,ldlast) ndep = 0 if (ldlast.gt.0.0.and.qout.gt.0.0) then call enrich(k,xu(k,iplane),xdend,xdetst,loadup,ldlast, 1 lddend,theta,iendfg) lddend = ldlast xdetst = xdend end if c *** L2 ELSE *** else c c Deposition ends in segment c loadup = ldlast call depos(xu(k,iplane),xdend,cdep,ainftc(k),binftc(k), 1 cinftc(k),phi,theta,ilast,dl,ldlast) ndep = 0 if (ldlast.gt.0.0.and.qout.gt.0.0) then call enrich(k,xu(k,iplane),xdend,xdetst,loadup,ldlast, 1 lddend,theta,iendfg) lddend = ldlast xdetst = xdend end if c c c Detachment After Deposition c c c Note: MSHEAR is an integer from 1 to 5. Notice that in the first two c cases, as well as the second two, the only difference is in the c value assigned to ERD. Also, note the games played with the first c two arguments to EROD in the various cases. Except for the first c argument passed to EROD, the 1st case is identical to the ELSE c of the 3rd, and the 2nd to the ELSE of the 4th. There may be c some real opportunities to simplify the code here! -- CRM 1/7/91 c go to (50,60,70,80,90)mshear c c ............. Shear below critical in entire segment. c ** MSHEAR = 1 ** 50 continue call erod(xdend,xl(k,iplane),ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),0.0,tauc,theta,phi, 1 ilast,dl,ldlast,xdbeg,ndep) go to 100 c c ............. Shear exceeds critical in entire segment. c ** MSHEAR = 2 ** 60 continue call erod(xdend,xl(k,iplane),ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),eata,tauc,theta,phi, 1 ilast,dl,ldlast,xdbeg,ndep) go to 100 c c ............. Shear increases down slope and exceeds critical at X=XC1. c ** MSHEAR = 3 ** 70 continue if (xdend.le.xc1) then call erod(xdend,xc1,ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),0.0,tauc,theta,phi, 1 ilast,dl,ldlast,xdbeg,ndep) if (ndep.eq.0) call erod(xc1,xl(k,iplane),ainf(k), 1 binf(k),cinf(k),ainftc(k),binftc(k),cinftc(k), 1 eata,tauc,theta,phi,ilast,dl,ldlast,xdbeg,ndep) else call erod(xdend,xl(k,iplane),ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),eata,tauc,theta,phi, 1 ilast,dl,ldlast,xdbeg,ndep) end if go to 100 c c ............. Shear decreases down slope and drops below critical at X=XC1. c ** MSHEAR = 4 ** 80 continue if (xdend.le.xc1) then call erod(xdend,xc1,ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),eata,tauc,theta,phi, 1 ilast,dl,ldlast,xdbeg,ndep) if (ndep.eq.0) call erod(xc1,xl(k,iplane),ainf(k), 1 binf(k),cinf(k),ainftc(k),binftc(k),cinftc(k),0.0, 1 tauc,theta,phi,ilast,dl,ldlast,xdbeg,ndep) else call erod(xdend,xl(k,iplane),ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),0.0,tauc,theta,phi, 1 ilast,dl,ldlast,xdbeg,ndep) end if go to 100 c c ............ Shear increases down slope -- exceeds critical at X=XCL, c then decreases from XCL to XC2, and drops below critical c at XC2. c ** MSHEAR = 5 ** 90 continue if (xdend.le.xc1) then call erod(xdend,xc1,ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),0.0,tauc,theta,phi, 1 ilast,dl,ldlast,xdbeg,ndep) if (ndep.eq.0) then call erod(xc1,xc2,ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),eata,tauc,theta, 1 phi,ilast,dl,ldlast,xdbeg,ndep) if (ndep.eq.0) call erod(xc2,xl(k,iplane),ainf(k), 1 binf(k),cinf(k),ainftc(k),binftc(k),cinftc(k), 1 0.0,tauc,theta,phi,ilast,dl,ldlast,xdbeg,ndep) end if else if (xdend.gt.xc2) then call erod(xdend,xl(k,iplane),ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),0.0,tauc,theta,phi, 1 ilast,dl,ldlast,xdbeg,ndep) else call erod(xdend,xc2,ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),eata,tauc,theta,phi, 1 ilast,dl,ldlast,xdbeg,ndep) if (ndep.eq.0) call erod(xc2,xl(k,iplane),ainf(k), 1 binf(k),cinf(k),ainftc(k),binftc(k),cinftc(k),0.0, 1 tauc,theta,phi,ilast,dl,ldlast,xdbeg,ndep) end if 100 continue c *** L2 ENDIF *** end if c c detachment at upper end of segment c c *** L1 ELSE *** else c dl = 0.0 du = 0.0 c c Note: MSHEAR is an integer from 1 to 5. Notice that in the first two c cases, as well as the second two, the only difference in in the c value assigned to ERD. Also, note the games played with the first c two arguments to EROD in the various cases. Except for the first c argument passed to EROD, the 1st case is identical to the last half c of the 3rd, and the 2nd to the last half of the 4th. There may be c some real opportunities to simplify the code here! -- CRM 1/7/91 c go to (110,120,130,140,150)mshear c c ............. Shear below critical in entire segment. c ** MSHEAR = 1 ** 110 continue call erod(xu(k,iplane),xl(k,iplane),ainf(k),binf(k), 1 cinf(k),ainftc(k),binftc(k),cinftc(k),0.0,tauc,theta, 1 phi,ilast,dl,ldlast,xdbeg,ndep) go to 160 c c ............. Shear exceeds critical in entire segment. c ** MSHEAR = 2 ** 120 continue call erod(xu(k,iplane),xl(k,iplane),ainf(k),binf(k), 1 cinf(k),ainftc(k),binftc(k),cinftc(k),eata,tauc,theta, 1 phi,ilast,dl,ldlast,xdbeg,ndep) go to 160 c c ............. Shear increases down slope and exceeds critical at X=XC1. c ** MSHEAR = 3 ** 130 continue call erod(xu(k,iplane),xc1,ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),0.0,tauc,theta,phi,ilast, 1 dl,ldlast,xdbeg,ndep) if (ndep.eq.0) call erod(xc1,xl(k,iplane),ainf(k),binf(k), 1 cinf(k),ainftc(k),binftc(k),cinftc(k),eata,tauc,theta, 1 phi,ilast,dl,ldlast,xdbeg,ndep) go to 160 c c ............. Shear decreases down slope and drops below critical at X=XC1. c ** MSHEAR = 4 ** 140 continue call erod(xu(k,iplane),xc1,ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),eata,tauc,theta,phi,ilast, 1 dl,ldlast,xdbeg,ndep) if (ndep.eq.0) call erod(xc1,xl(k,iplane),ainf(k),binf(k), 1 cinf(k),ainftc(k),binftc(k),cinftc(k),0.0,tauc,theta, 1 phi,ilast,dl,ldlast,xdbeg,ndep) go to 160 c c ......... Shear increases down slope -- exceeds critical at X=XCL, c then decreases from XCL to XC2, and drops below critical c at XC2. c ** MSHEAR = 5 ** 150 continue call erod(xu(k,iplane),xc1,ainf(k),binf(k),cinf(k), 1 ainftc(k),binftc(k),cinftc(k),0.0,tauc,theta,phi,ilast, 1 dl,ldlast,xdbeg,ndep) if (ndep.eq.0) then call erod(xc1,xc2,ainf(k),binf(k),cinf(k),ainftc(k), 1 binftc(k),cinftc(k),eata,tauc,theta,phi,ilast,dl, 1 ldlast,xdbeg,ndep) if (ndep.eq.0) call erod(xc2,xl(k,iplane),ainf(k), 1 binf(k),cinf(k),ainftc(k),binftc(k),cinftc(k),0.0, 1 tauc,theta,phi,ilast,dl,ldlast,xdbeg,ndep) end if c c *** L1 ENDIF *** end if c 160 continue c c If this was a detachment section on the segment of the c OFE which went into deposition (NDEP = 1) then call c the deposition routine from the point where load equals c transport capacity to the end of the segment. c if (ndep.ne.0) then if (ilast.lt.102) then dl = 0.0 du = 0.0 cdep = depc(xdbeg,ainftc(k),binftc(k),phi,theta,du,ktrato, 1 qostar) loadup = ldlast if (loadup.lt.lddend) loadup = lddend call depos(xdbeg,xl(k,iplane),cdep,ainftc(k),binftc(k), 1 cinftc(k),phi,theta,ilast,dl,ldlast) ndep = 0 if (ldlast.gt.0.0.and.qout.gt.0.0) then call enrich(k,xdbeg,xl(k,iplane),xdetst,loadup,ldlast, 1 lddend,theta,iendfg) lddend = ldlast xdetst = xl(k,iplane) end if end if end if c c End of BIG IF end if c c End of BIG DO-LOOP 170 continue c c Compute enrichment ratio at the end of each c overland flow element. c iendfg = 1 c c Changed following call to ENRICH 6/30/94 because new code c inserted in ENRICH by Nearing must have actual value for location c of the beginning and end of the deposition region at the last c point on the OFE - thus: xtop = xbot = 1.0 c dcf 6/30/94 c call enrich(k,1.0,1.0,xdetst,ldlast,ldlast,lddend,theta,iendfg) c return end