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 wus, wust, wusp, sf1, sf10, sf84 real sf200, szcr, sfcr, sflos, smlos real szrgh, sxprg, sarrc, slrr real seags, secr, sfcla, sfsan real svroc, brsai, bzht, bffcv, time real ix, qi, qssi, q10i, qo, qsso, q10o C Added just to get it to compile - LEW real 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 = weighted average 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 +++ save sfa12 data sfa12 /0/ 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++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++