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 integer wustfl 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 (mm) 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 qav = average saltation/creep discharge in cell (kg/m*s) c szc = tmp variable for prior crust thickness (mm) c sz = tmp variable for roughness (mm) c crlos = factor that decreases loose cover with roughness c dmlos = change in loose mass on aggregated sfc. (kg/m^2) c dmcld = change in clod mass on aggregated sfc. (kg/m^2) c szv = change in height based on volume change (mm) c sacd = specific area of clods per unit volume (mm^2/mm^3) c sff1 = soil fraction limit on area abraded c sff2 = soil fraction limit on area emitting c +++ PARAMETERS +++ real 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 cen, ceno, renb, renv, qen, qcp, ci real sfa12, sfcv, sargc, sfsn real a2, 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 real qav, szc, sz, crlos real dmlos, dmcld, szv, sacd, sff1, sff2 integer flg 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 +++ END SPECIFICATIONS +++ c calc. fraction of area with shelter > 12 degrees 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 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 !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) cen = ceno*renv*renb c soil fraction suspension in emitted soil; lower for loose on cr sfssen = (sf10/sf200)*(1. - sfcr*0.8) 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 c set upper limit on ci ci = amin1(ci, 2.0) endif c Calc. abrasion params. for saltation/creep !soil fraction saltation/creep in saltation c sfsn = sfa12**0.25 sfsn = 1.0 !fraction of saltation abrading clods, crust and loose fan = (exp(- 4.0*bffcv))*(1 - 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 = (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)) c put surface update eqns here c surface update equations flg = 0 c execute this section if flag is set if (flg .eq. 1) then c set threshold friction vel. update flag on 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 = amax1(0.001, szcr) szcr = szcr - (fancr*cancr*qav/1.6)*time c update crust (consolidated) zone cover sfcr = sfcr - sfcr*(szc-szcr)/szc sfcr = amax1(0.0, sfcr) c update loose on crust c terms: prior value, trapped, abraded but not emitted, emitted smlos = smlos +(ct*sfcr*(1 -svroc)*qav & +(1-sfssan)*fancr*cancr*qav**2/qen & -(1-sfssen)*cen*(qen-qav)*sflos*sfcr*(1-svroc) & /(1-sfcv+0.001))*time smlos = amax1(0.0, smlos) c update cover of loose mass (same as SOIL eq. s-24,2-25) sz = amax1(szrgh, 4.0*slrr) crlos = exp(-0.003*sz) sflos = (1.0 -exp(-3.5*smlos**1.5))*crlos sflos = amax1(0.0, sflos) 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 dmlos=((1-sfcr)*(1-svroc)*ct*qav & + (1-sfssan)*fanag*canag*qav**2/qen & -(1-sfssen)*cen*(qen-qav)*(1.0-sfcr)*(1-svroc)*sf84 & /(1.001-sfcv))*time c change in abraded mass from clod area dmcld = (fanag*canag*qav)*time c net volume (height) removed from loose & clod areas c for aggregated area all erodible or all clods if((sf84 .gt. 0.99) .or. (sf84 .lt. 0.01)) then szv = 0.0 c for aggregated surface partly cloddy else szv = dmlos/((1.001-sfcr)*(1.001-svroc)*sf84*1.2) & + dmcld/((1.001-sf84)*(1.001-sfcr)*(1.001-svroc)*1.6) endif c specific clod area per unit volume (crude estimate) sacd = 1.2*exp(sf84 - 1.0)*(1-sf84) c update surface loose material in zero layer sf84 = sf84 + sacd*szv sf84 = amax1(0.0, sf84) sf84 = amin1(1.0, sf84) c ^^^ tmp output c write (*,*) 'for clods: dmlos, dmcld, szv' c write (*,*) dmlos, dmcld, szv c write (*,*) c endif c update surface roughness c set limits on areas in denominators c surface fraction stable that can change height sff1 = amax1((sfcv-svroc), 0.01) c surface fraction loose sff2 = amax1((1-sfcv), 0.01) c rate of height change caused by emission szv= cen*(qen - qav)/(sff2*1.2) c if trapping then emission lowers roughness If (ct .gt. 0.) then szv= -szv endif c calc. change in surface height c terms: emission, trapped, abrasion, & abrasion not emitted szv = (szv - 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. 1.5) then slrr = 1.5 endif endif endif 100 continue return end c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++