!********************************************************************** ! subroutine sbqout ver 9/2/98 file name sbqout.for !********************************************************************** subroutine sbqout & & (flg, wus, wust, wusp, sf1, sf10, sf84, & & sf200, szcr, sfcr, sflos, smlos, & & szrgh, sxrgs, sxprg, slrr, & & seags, secr, sfcla, sfsan,sfvfs, & & svroc, brsai, bzht, bffcv, time, & & lx, qi, qssi, q10i, i, j, imax, jmax, & & smaglos,dmlos, sf84mn, sf84ic,svrocic, & & qo, qsso, q10o ) ! ! +++ PURPOSE +++ ! calculate the saltation/creep, suspension, and PM-10 discharge ! from each control volume ! ! +++ ARGUMENT DECLARATIONS +++ integer flg !Surface update flag (1=on, 0=off) real wus, wust, wusp, sf1, sf10, sf84 real sf200, szcr, sfcr, sflos, smlos real szrgh, sxrgs, sxprg, slrr real seags, secr, sfcla, sfsan,sfvfs real svroc, brsai, bzht, bffcv, time real lx, qi, qssi, q10i, qo, qsso, q10o real smaglos, dmlos, sf84mn, sf84ic, svrocic integer i, j, imax, jmax include 'erosion/p1erode.inc' !Parameters defining scope of erosion submodel ! +++ ARGUMENT DEFINITIONS +++ ! wus = friction velocity (m/s) ! wust = friction velocity threhold at emission (m/s) ! wusp = friction velocity threshold at transport cap.(m/s) ! sf1, sf10 soil fractions less than 0.010 and 0.10 mm ! sf84, sf200 soil fractions less than 0.84 and 2.0 mm ! sf84ic = soil surface fraction less than 0.84 at start of event ! svrocic= soil surface volume rock at start of event ! sfcr = soil fraction area crusted ! sflos = soil fraction area of loose soil (only on crust) ! sarrc = soil angle r. roughness weibull c parm ! szrgh = soil ridge height (mm) ! sxprg = soil ridge spacing parallel wind direction (mm) ! sxrg = soil ridge spacing (mm) ! seags = soil aggregate stability (ln(J/kg)) ! secr = soil crust stability (ln(j/kg)) ! sfcla = soil fraction clay by mass ! sfsan = soil fraction sand by mass ! sfvfs = soil fraction very fine sand by mass (0.05-0.1 mm) ! svroc = soil fraction rock >2.0 mm by volume ! brsai = biomass stem area index ! bzht = biomass height (m) ! bffcv = biomass fraction flat cover ! lx = node spacing in x-direction (m) ! qi, qssi, q10i=input to C.V. of saltion, creep, pm-10 (kg/m*2) ! qo, qsso, q10o=output from C.V. of saltion, creep, pm-10 (kg/m*s) ! qav = average saltation/creep discharge in cell (kg/m*s) ! szc = tmp variable for prior crust thickness (mm) ! sz = tmp variable for roughness (mm) ! crlos = factor that decreases loose cover with roughness ! dmlos = change in loose mass on aggregated sfc. (kg/m^2) ! dmcld = change in clod mass on aggregated sfc. (kg/m^2) ! szv = change in height based on volume change (mm) ! sacd = specific area of clods per unit volume (mm^2/mm^3) ! ! +++ PARAMETERS +++ real cs, cmp, ctf, cdp, c10dp parameter (cs=0.3, cmp=0.0001, ctf=1.2, cdp=0.02, & & c10dp=0.001) ! +++ PARAMETER DEFINITIONS +++ ! cs = saltation transport coef. (kg*s^2/m^4) ! cmp = mixing parameter coef. ! ct = ridge trapping coef. ! ctf = ridge trapping fraction ! cdp = deposition coef. for saltation ! c10dp = deposition coef. for pm10 ! ! +++ LOCAL VARIABLES +++ ! real cen, ceno, renb, renv, qen, qcp, ci real sfa12, sfcv, sargc, sarrc, sfsn real fan, fanag, fancr, canag, cancr, fancan, sfssen real sfssan, sf10an, sf10en, cbk, sf10bk,a,b,c,f,g,h,k real s, p, t1,t2, ct real cm, tmp, srrg, dmtlos, fdm, fsm real qav, szc, sz, crlos real szv, qitmp, sf84old integer kk, m,n ! ! +++ LOCAL VARIABLE DEFINITIONS +++ ! cen = coef. of emission (1/m) ! ceno = Coef. of emission for bare, loose, erodible sfc. ! renb = reduction in emission on bare soil by cover and rough. ! sargc = soil angle of ridge shelter weibull 'c' (deg.) ! sarrc = soil angle of random roughness weibull 'c'(deg.) ! sfa12 = soil fraction area with shelter angle > 12 deg. ! sfcv = soil fraction clod & crust cover ! srrg = ratio of ridge spacing to ridge spacing parallel wind direction ! renv = reductin in emission on bare soil by flat biomass ! qen = transport capacity using emission threshold (kg/m*s) ! qcp = transport capacity using trans. cap. threshold (kg/m*s) ! ci = coef. of saltation interception by biomass stems (1/m) ! a2 = intermediate calc. parameter ! fan = fraction saltation impacting clods & crust ! fanag = fraction saltation impacting clods ! fancr = fraction saltation impacting crust ! canag = coef. of abrasion of clods (1/m) ! cancr = coef. of abrasion of crust (1/m) ! fancan = product of fan & can ie. effective abrasion coef. (1/m) ! sfssen = soil fraction to suspension from emission ! sfssan = soil fraction from abrasion to emission ! sf10en = soil fraction of emitted suspension that is pm-10 ! sf10an = soil fraction of abraded suspension that is pm-10 ! a,b,c = composite parameters for saltation discharge soln eq. ! f, g = composite parameters for suspension discharge soln eq. ! h, k = composite parameter for pm-10 discharge eq. ! s, p = intermediate calc. param. for saltation discharge soln. ! t1,t2 = intermediate calc. param. for salt. and susp. soln eq. ! tmp = intermdediate calc. param. for transport cap. soln eq. ! dmtlos = soil loss/dep. for each time-step (+ = deposition) ! fdm = split of dmtlos to agg. & rock reservoir ! fsm = split of dmtlos to crust reservoir ! ! +++ END SPECIFICATIONS +++ ! ! calc. fraction of area with shelter > 12 degrees sarrc = 2.3 * sqrt (slrr) if ((sxprg .gt. 10.0) .and. (szrgh .gt. 1.0)) then sargc = 65.4*(szrgh/sxprg)**0.65 sfa12 =(1.- exp(-(12./sargc)**0.77))*(exp(-(12./sarrc)**0.77)) & & + exp(-(12./sargc)**0.77) else sfa12 = exp(-(12./sarrc)**0.77) endif ! ! update SF84 when net loose soil loss with smaglos for current wus if (dmlos < 0.0) then ! calc loose reservoir in area with roughness shelter < 12 degrees smaglos = smaglos*(1-sfa12) if (dmlos > -0.989*smaglos) then sf84 = sf84ic + alog((smaglos - abs(dmlos))/smaglos) & *(sf84ic - sf84mn)/4.6 else sf84 = sf84mn endif ! set limits on sf84 sf84 = amax1(0.001, sf84) sf84 = amin1(1.0, sf84) endif ! ! calc. transport cap. for emission edit 5/30/01 LH qen = cs * wus * wus * (wus - wust) ! calc. transport cap. for trapping qcp = cs * wus * wus * (wus - wusp) if (qcp .lt. 0.0) then qcp = 0.001 endif ! ! store tmp qi qitmp = qi ! test for trap region with both saltation & suspension deposition if (qen .le. 0.) then qo = 0.0 qsso = qssi*exp(-cdp*lx) q10o = q10i*exp(-c10dp*lx) go to 100 ! test for saltation deposition region with suspension emission ! (a more exact procedure would be to find new q(x) for deposition ! and then put that solution into the suspension eqn and develop ! another soln for the suspension eqn. For now we have set ! qi to qen to calc. suspension and then reset qi to the real qi ! value to calc. the deposition in sberod. This explanation added ! 12-5-2000 LH) ! elseif (qen .le. qi) then qi = qen endif ! ! Calc. emission params for saltation/creep ! emission coef. for bare, loose, erodible soil (1/m) ceno = 0.06 ! fraction emission coef. reduction by flat biomass renv = 0.075 + 0.934*exp(-bffcv/0.149) ! fraction emission coef. reduction by roughness and area not emitting sfcv = ((1.-sfcr)*(1.-sf84) + sfcr - sflos*sfcr)*(1.-svroc) & & + svroc ! ! renb = exp(-4.0*sfcv)*exp(-2.5*sfa12) edit 3-29-01 LH ! renb = exp(-3.0*sfcv)*(1.0 - sfa12) ! edit 7-17-01 LH renb = (-0.051 + 1.051*exp(-sfcv/0.33050512))*(1 - sfa12) cen = ceno*renv*renb ! if (cen .lt. 0.0001) cen = 0.0 ! soil fraction suspension in emitted soil; lower for loose on cr sfssen = (sf10/(sf200+0.001))*(1. - sfcr*0.8) ! soil mixing parameter for suspension cm = cmp*sfssen ! ! calc. trap params. for saltation/creep ! calc. coef of stem interception ! for plants above 5 mm height edited 3-29-01 LH if (bzht .lt. 0.005) then ci = 0. else if (sxrgs .gt. 10) then srrg = sxrgs/sxprg else srrg = 1.0 endif ci = 0.05*srrg*(1 - exp(-brsai/bzht))*(1 - exp(-bzht/0.15)) endif ! ! Calc. abrasion params. for saltation/creep ! soil fraction saltation/creep in saltation ! changed estimate edit 5/30/01 ! sfsn = sfa12**0.09 ! sfsn = amax1(sfsn,0.7) sfsn = 1 ! ! fraction of saltation abrading clods, crust and loose fan = (exp(- 4.0*bffcv))*(exp(-3*svroc))*sfsn fan = amax1(fan,0.) ! fraction abrading aggregates fanag = (1. - sf84)*(1. - sfcr)*fan ! fraction abrading crust fancr = (sfcr - sflos*sfcr)*fan ! calc. coef. of abrasion canag = exp(-2.07-0.077*seags**2.5-0.119*alog(seags)) cancr = exp(-2.07-0.077*secr**2.5-0.119*alog(secr)) ! sum fancan = fancr*cancr + fanag*canag if (fancan .lt. 0.0001) fancan = 0.0001 ! calc fraction of abraded of suspension size sfssan = (0.4 - 4.83*sfcla + 27.18*sfcla**2 & & - 53.7*sfcla**3 + 42.25*sfcla**4 - 10.7*sfcla**5) ! ! set upper limit on sfssan for high sand content sfssan = min(sfssan, (1.0 - (sfsan-sfvfs))) ! ! calc suspension breakage coef. for saltation creep ! edit 6-9-01 LH cbk = 0.11*canag*(1. - (sfsan-sfvfs))*sfsn ! ! calc total trap coef. if ((qcp .lt. qen) & & .and. ((slrr .gt. 10.) .or. (szrgh .gt. 50.))) then ct = (1. - sfssen)*cen*ctf*(qen - qcp)/qen else ct = 0. endif ! ! assemble composite params. for saltation/creep a = cen *(1. - sfssen)*qen b = (1. - sfssan)*fancan - (1. - sfssen)*cen - cbk - ci - ct c = (1. - sfssan)*fancan*1./qen ! ! solve for saltation/creep out ! collect variables s = sqrt(4.*a*c + b**2.) p = (-2.*c*qi + b)/s ! change p-values that are out of range ! math range restriction: (-1 < p < 1) edit 7-18-01 LH if (p .le. -1.) then t1 = -20 elseif (p .ge. 1.) then t1 = 20 else t1 = (s/2.0)*(-lx) + 0.5*alog((1. + p)/(1. - p)) endif ! calculate atanh qo = (s/(2.*c))*(-tanh(t1) + b/s) ! ! assemble composite params. for suspension f = sfssen*cen*qen ! added last term to intercept some ss LH edit 6-11-01 g = sfssan*fancan - sfssen*cen + cbk + cm - sfssen*ci ! ! ^^^ tmp out ! set counter ! m = (imax-1)/8 ! m = max(1,m) ! do 2 n = 1, (imax-1), m ! if (i .eq. n .and. j .eq. 2) then ! write (*,*) ! write (*,*) 'output from sbqout: i=',i,'j=',j ! write (*,*) 'a=', a, 'b=', b, 'c=',c ! write (*,*) 'g=', g, 'f=', f ! write (*,*) 'cen=', cen, 'qen=', qen ! write (*,*) 'sfssen=', sfssen, 'sfssan=',sfssan, ! & 'fancan=',fancan ! write (*,*) 'cbk=',cbk, 'ci=',ci,'ct=',ct,'cm=',cm ! endif ! 2 continue ! ^^^ end tmp out ! ! trap to prevent over range of exp if ((s*lx) .gt. 40.0) then s = 40.0/lx endif ! solve for suspension out ! collect variables t2 = alog(exp(s*lx)*(1.-p)+ p + 1.) qsso=qssi+(1./(2.*c))*((-g*s+g*b+2.*f*c)*lx+2.*g*(-alog(2.0)+t2)) ! ! assemble params. for PM-10 sf10en = sf1/(sf10 + 0.00001) sf10an = 0.67*sfcla sf10an = min(0.35,sf10an) sf10bk = 0.0015 + 0.023*(1.0 - sfsan - sfcla)**2.0 ! ! assemble composite params. for h = sf10en*f k = sf10an*sfssan*fancan - sf10en*sfssen*cen + sf10bk*cbk & & + sf10en*cm ! ! solve for PM-10 out (similar to suspension out) q10o=q10i+(1./(2.*c))*((-k*s+k*b+2.*h*c)*lx+2.*k*(-alog(2.0)+t2)) ! SURFACE UPDATE EQUATIONS !now added as a commandline argument to tsterode - LEW ! execute this section if flag is set if (flg .eq. 1) then ! calc. avg. saltation discharge in grid cell ! qav = (qi + qo)/2.0 ! ! calc change in loose surface soil for time step (+ = surface gain) ! terms: - total loss, + created by abrasion dmtlos = ((-(qo-qitmp) - (qsso-qssi))/lx & + fancan*qitmp)*time ! ! skip surface update for small dmtlos ! if (abs(dmtlos) < 0.005) go to 100 ! ! calc loose surface parameters (fsm, fdm) for later split ! of loss or gain of mobile soil on crust and other surface. ! ! fraction reservior area on crust ! asm = sfcr*(1 - svroc) ! fraction reservoir area on agg. & rock ! adm = 1 - sfcr*(1 - svroc) fsm = 1.0 fdm = 1.0 ! coef. for cover of loose mass (same as SOIL eq. s-24,2-25) sz = amax1(szrgh, 4.0*slrr) crlos = exp(-0.08*sz**0.5) ! ! execute this section if crust present if (sfcr .gt. 0.01) then ! ! update loose on crust smlos = smlos + dmtlos*fsm smlos = amax1(0.0, smlos) ! ! update cover of loose mass tmp = 3.5*smlos**1.5 if(tmp.gt.80.) then !test and prevent underflow condition sflos = crlos else sflos = (1.0 - exp(-tmp))*crlos endif ! ! update crust thickness szc = amax1(0.001, szcr) szcr = szcr - (fancr*cancr*qi/(1.2*sfcr))*time szcr = amax1(szcr, 0.00) ! ! update crust (consolidated) zone cover sfcr = sfcr*szcr/szc sfcr = amax1(0.0, sfcr) else smlos = 0. sflos = 0. endif ! ! execute this section if clods are present if (sfcr .lt. 0.99) then ! ! change in loose mass on aggregated and rock surface, (+)=increase dmlos = dmlos + dmtlos*fdm ! ! update sf84 when net loose soil gain on agg. & rock if (dmlos > 0.) then ! update cover (sf84) of loose mass when loose soil gain ! (same as SOIL eq. s-24,2-25) ! sz = amax1(szrgh, 4.0*slrr) ! crlos = exp(-0.08*sz**0.5) ! back calculate a surface mass based on initial sf84 tmp = (-alog(1-sf84ic)/3.5)**0.6667 tmp = 3.5*(dmlos+tmp)**1.5 if(tmp.gt.80.) then !test and prevent underflow condition sf84 = crlos else sf84 = (1.0 - exp(-tmp))*crlos endif ! sf84 = amax1(sf84,sf84ic) ! old code ! retain sf84old, if old > new ! if (sf84old > sf84) then ! sf84 = sf84old ! endif endif ! set limits on sf84 sf84 = amax1(0.001, sf84) sf84 = amin1(1.0, sf84) ! ! update rock cover in proportion to total of rock + ag. cover if (svrocic > 0.0) then svroc = (1-sf84)*svrocic/(1-sf84ic*(1-svrocic)-sf84*svrocic) endif ! update sf200 sf200 = (2 - sf84)*sf84 sf200 = amin1(sf200, 1.0) sf200 = amax1(0.001,sf200) ! update sf10 (crude estimate, need to update) sf10 = sf84*0.01 ! fanag = (1-sf84)*(1-sfcr)*fan endif ! ! update surface roughness ! ! rate of volume change caused by emission ! (used bulk densities of 1.2 for loose and 1.4 for crust ! to calc. mm depth per m^2 of area) szv= cen*(qen - qi)/1.2 ! ! if trapping then emission lowers roughness, ! i.e. it comes from highest areas. If (ct .gt. 0.) then szv= -szv endif ! calc. change in mean surface depth ! terms: emission, trapped, abrasion, & abrasion not emitted szv = (szv - ct*qi/1.2 - fancan*qi/1.4 & & -(1.0 - sfssan)*fancan*qi*qi/(qen*1.2))*time ! ! update ridge height if (szrgh .gt. 10.) then szrgh = szrgh + szv*2 ! set lower limit on ridge height szrgh = amax1(szrgh, 0.0) endif ! update random roughness slrr = slrr + 2.0*szv/4.0 ! set lower limit on slrr slrr = max(slrr, SLRR_MIN) ! ! !^^^^ tmp output section ! set counter ! m = (imax-1)/8 ! ! jmax = jmax ! ! do 10 n = 1, (imax-1), m ! if (i .eq. n .and. j .eq. 2) then ! write (*,19) sfssan, sfssen, cancr, canag, sfcv ! write (*,20) wus, wust, ct, fancr, fanag, fancan ! write (*,21) qo, qi, qav, qen, cen, time ! write (*,*) ! endif ! 10 continue ! 19 format(1x, 'sfssan, sfssen, cancr, canag, sfcv', 8f8.3) ! ! 20 format(1x, ' wus, wust, ct, fancr, fanag, fancan', 8f8.3) ! 21 format(1x, 'qo qi qav qen cen time', 5f8.3, f8.0) !^^^ end tmpout endif 100 continue qi = qitmp return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++