c********************************************************************** c subroutine sbqout ver 9/2/98 file name sbqout.for c********************************************************************** subroutine sbqout i (wus, wust, wusp, sf1, sf10, sf84, i sf200, szcr, sfcr, sflos, smlos, i szrgh, sxprg, sarrc, slrr, i seags, secr, sfcla, sfsan, i svroc, brsai, bzht, bffcv, time, i ix, qi, qssi, q10i, wustfl, o qo, qsso, q10o ) c +++ PURPOSE +++ c calculate the saltation/creep, suspension, and PM-10 discharge c from each control volume c +++ ARGUMENT DECLARATIONS +++ real*4 wus, wust, wusp, sf1, sf10, sf84 real*4 sf200, szcr, sfcr, sflos, smlos real*4 szrgh, sxprg, sarrc, slrr real*4 seags, secr, sfcla, sfsan real*4 svroc, brsai, bzht, bffcv, time real*4 ix, qi, qssi, q10i, qo, qsso, q10o C Added just to get it to compile - LEW real*4 svz integer wustfl c c +++ ARGUMENT DEFINITIONS +++ c wus = friction velocity (m/s) c wust = friction velocity threhold at emission (m/s) c wusp = friction velocity threshold at transport cap.(m/s) c sf1, sf10 soil fractions less than 0.010 and 0.10 mm c sf84, sf200 soil fractions less than 0.84 and 2.0 mm c sfcr = soil fraction area crusted c sflos = soil fraction area of loose soil (only on crust) c sarrc = soil angle r. roughness weibull c parm c szrgh = soil ridge height (mm) c sxprg = soil ridge spacing c seags = soil aggregate stability (ln(J/kg)) c secr = soil crust stability (ln(j/kg)) c sfcla = soil fraction clay by mass c sfsan = soil fraction sand by mass c svroc = soil fraction rock >2.0 mm by volume c brsai = biomass stem area index c bzht = biomass height (m) c bffcv = biomass fraction flat cover c ix = node spacing in x-direction (m) c qi, qssi, q10i=input to C.V. of saltion, creep, pm-10 (kg/m*2) c qo, qsso, q10o=output from C.V. of saltion, creep, pm-10 (kg/m*s) c +++ PARAMETERS +++ real*4 cs, cmp, ctf, cdp, c10dp parameter (cs=0.3, cmp=0.0001, ctf=0.95, cdp=0.02, & c10dp=0.001) c +++ PARAMETER DEFINITIONS +++ c cs = saltation transport coef. (kg*s^2/m^4) c cmp = mixing parameter coef. c ct = ridge trapping coef. c ctf = ridge trapping fraction c cdp = deposition coef. for saltation c c10dp = deposition coef. for pm10 C +++ LOCAL VARIABLES +++ real*4 cen, ceno, renb, renv, qen, qcp, ci real*4 sfa12, sfcv, sargc, sfsn real*4 a2, fan, fanag, fancr, canag, cancr, fancan, sfssen real*4 sfssan, sf10an, sf10en, cbk, sf10bk,a,b,c,f,g,h,k real*4 s, p, t1,t2, ct real*4 cm, tmp real*4 qav, szc, sz, crlos real*4 dmlos, dmcld, szv, sacd, sff1, sff2 integer flg c c +++ LOCAL VARIABLE DEFINITIONS +++ c cen = coef. of emission (1/m) c ceno = Coef. of emission for bare, loose, erodible sfc. c renb = reduction in emission on bare soil by cover and rough. c sargc = soil angle of ridge shelter weibull 'c' (deg.) c sfa12 = soil fraction area with shelter angle > 12 deg. c sfcv = soil fraction clod & crust cover c renv = reductin in emission on bare soil by flat biomass c qen = transport capacity using emission threshold (kg/m*s) c qcp = transport capacity using trans. cap. threshold (kg/m*s) c ci = coef. of saltation interception by biomass stems (1/m) c a2 = intermediate calc. parameter c fan = fraction saltation impacting clods & crust c fanag = fraction saltation impacting clods c fancr = fraction saltation impacting crust c canag = coef. of abrasion of clods (1/m) c cancr = coef. of abrasion of crust (1/m) c fancan = product of fan & can ie. effective abrasion coef. (1/m) c sfssen = soil fraction to suspension from emission c sfssan = soil fraction from abrasion to emission c sf10en = soil fraction of emitted suspension that is pm-10 c sf10an = soil fraction of abraded suspension that is pm-10 c a,b,c = composite parameters for saltation discharge soln eq. c f, g = composite parameters for suspension discharge soln eq. c h, k = composite parameter for pm-10 discharge eq. c s, p = intermediate calc. param. for saltation discharge soln. c t1,t2 = intermediate calc. param. for salt. and susp. soln eq. c tmp = intermdediate calc. param. for transport cap. soln eq. c c +++ END SPECIFICATIONS +++ c calc. transport cap. for emission tmp = sf84 + sflos*sfcr - sfa12 if (tmp .lt. 0.0) then tmp = 0.0 endif qen = cs * wus * wus * (wus - (wust - 0.05*tmp)) c calc. transport cap. for trapping qcp = cs * wus * wus * (wus - (wusp - 0.05*tmp)) c test for trap region with both saltation & suspension deposition if (qen .le. 0.) then qo = 0.0 qsso = qssi*exp(-cdp*ix) q10o = q10i*exp(-c10dp*ix) go to 100 c test for saltation deposition region with suspension emission elseif (qen .le. qi) then qi = qen endif C Calc. emission params for saltation/creep c emission coef. for bare, loose, erodible soil (1/m) ceno = 0.06 c fraction emission coef. reduction by flat biomass renv = 0.075 + 0.934*exp(-bffcv/0.149) c fraction emission coef. reduction by roughness and area not emitting if ((sxprg .gt. 0) .and. (szrgh .gt. 0.)) then sargc = 65.4*(szrgh/sxprg)**0.65 sfa12 =(1.- exp(-(12./sargc)**0.77))*(exp(-(12./sarrc)**0.77)) else sfa12 = exp(-(12./sarrc)**0.77) endif sfcv = ((1.-sfcr)*(1.-sf84) + sfcr - sflos)*(1.-svroc) + svroc renb = (1. - sfcv)*exp(-2.5*sfa12) cen = ceno*renv*renb c soil fraction suspension in emitted soil; lower for loose on cr sfssen = (sf10/sf200)*(1. - sfcr*0.5) c soil mixing parameter for suspension cm = cmp*sfssen c calc. trap params. for saltation/creep c calc. coef of stem interception if (bzht .eq. 0.) then ci = 0. else ci = brsai/bzht endif c Calc. abrasion params. for saltation/creep c soil fraction saltation in creep a2 = (sf84 -sf10)/(sf200 - sf10) sfsn = 1. - (1.-a2)*exp(-sfa12/20.0) c fraction of saltation abrading clods and crust fan = (1. - 4.*bffcv - 2.*svroc*(1.-bffcv))*sfsn fan = amax1(fan,0.) fanag = (1. - sf84)*(1. - sfcr)*fan fancr = (sfcr - sflos)*fan c 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)) c sum fancan = fancr*cancr + fanag*canag if (fancan .lt. 0.0001) fancan = 0.0001 c calc fraction of abraded of suspension size sfssan = (1. - sfsan*0.5) * (0.4 - 4.83*sfcla + 27.18*sfcla**2 & - 53.7*sfcla**3 + 42.25*sfcla**4 - 10.7*sfcla**5) c calc suspension breakage coef. for saltation creep c need to update cbk = 0.11*canag*(1. - sfsan) c calc total trap coef. if (qcp .lt. qen ) then ct = (1. - sfssen)*cen*ctf*(1. - qcp/qen) else ct = 0. endif c 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 c solve for saltation/creep out c collect variables s = sqrt(4.*a*c + b**2.) p = (-2.*c*qi + b)/s c change p-values that are out of range c math range restriction: (-1 < p < 1) if (p .le. -1.) then p = -0.9999 elseif (p .ge. 1.) then p = 0.9999 endif c calculate atanh (p) t1 = 0.5*alog((1. + p)/(1. - p)) qo = (s/(2.*c))*(-tanh((s/2.)*(-ix) + t1) + b/s) c assemble composite params. for suspension f = sfssen*cen*qen g = sfssan*fancan - sfssen*cen + cbk + cm c solve for suspension out c collect variables t2 = alog(exp(s*ix)*(1.-p)+ p + 1.) qsso=qssi+(1./(2.*c))*((-g*s+g*b+2.*f*c)*ix+2.*g*(-alog(2.0)+t2)) c assemble params. for PM-10 sf10en = sf1/sf10 sf10an = 0.67*sfcla sf10an = min(0.35,sf10an) sf10bk = 0.0015 + 0.023*(1.0 - sfsan - sfcla)**2.0 c assemble composite params. for h = sf10en*f k = sf10an*sfssan*fancan - sf10en*sfssen*cen + sf10bk*cbk & + sf10en*cm c solve for PM-10 out (similar to suspension out) q10o=q10i+(1./(2.*c))*((-k*s+k*b+2.*h*c)*ix+2.*k*(-alog(2.0)+t2)) 100 continue c put surface update eqns here c surface update equations wustfl = 0 flg = 0 c execute this section if flag is set if (flg .eq. 1) then c set threshold friction vel. update flag on wustfl = 1 c calc. avg. saltation discharge in grid cell qav = (qi + qo)/2.0 c execute this section if crust present if (sfcr .gt. 0.01) then c update crust thickness szc = szcr szcr = szcr - (fancr*cancr*qav/1.6)*time c update crust (consolidated) zone cover sfcr = sfcr - sfcr*(szc-szcr)/szc if (sfcr .lt. 0.) then sfcr = 0.0 endif c update loose on crust c terms: prior value, trapped, abraded but not emitted, emitted smlos = smlos +(ct*sfcr*qav+(1-sfssan)*fancr*cancr*qav**2/qen & -(1-sfssen)*cen*(qen-qav)*sflos)*time c update cover of loose mass (same as SOIL eq. s-24,2-25) sz = amax1(szrgh, 4.0*slrr) crlos = smlos**0.5/(0.24*sz) sflos = (1.0 -exp(-3.5*smlos**1.5))*crlos endif c execute this section if clods are present if (sfcr .lt. 0.99) then c change in loose mass on aggregated surface, (+)=increase c terms: trapped, abraded but not emitted, emitted c dmlos = ((1-sfcr)*ct*qav + (1-sfssan)*fanag*canag*qav**2/qen & -(1-sfssen)*cen*(qen-qav)*(1.0-sfcr)*sf84)*time c change in abraded mass from clod area dmcld = (fanag*canag*qav)*time c net volume (height) removed from loose & clod areas if(sf84 .gt. 0.99) then szv = dmlos/((1-sfcr)*sf84*1.2) elseif (sf84 .lt. 0.01) then szv = dmcld/((1-sf84*(1-sfcr))*1.6) else szv = dmlos/(sf84*1.2) + dmcld/((1-sf84)*1.6) endif c specific clod area per unit volume (crude estimate) sacd = 1.2*exp(sf84 - 1.0) c update surface loose material in zero layer sf84 = sf84 + sacd*szv endif c update surface roughness c set limits on areas in denominators sff1 = amax1((1-sf84 + sfcr),0.01) sff2 = amax1((sf84+sflos), 0.01) svz = cen*(qen - qav)/sff2 c it trapping then emission lowers roughness If (ct .gt. 0.) then svz = -svz endif c calc. change in surface height c terms: emission, trapped, abrasion, & abrasion not emitted svz = (svz - ct*qav/0.4 - fancan*qav/sff1 & -(1.0 - sfssen)*fancan*qav*qav/(qen*0.4))*time c update ridge height if (szrgh .gt. 10.) then szrgh = szrgh +szv endif c update random roughness slrr = slrr + szv/2. if (slrr .lt. 2.0) then slrr = 2.0 endif endif return end c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++