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, sxrgs, sxprg, slrr,
     i seags, secr, sfcla, sfsan,sfvfs,
     i svroc, brsai, bzht, bffcv, time,
     i lx, qi, qssi, q10i, i, j, imax, jmax,
     o qo, qsso, q10o )
c
c     +++ PURPOSE +++
c     calculate the saltation/creep, suspension, and PM-10 discharge
c     from each control volume
c
c     +++ ARGUMENT DECLARATIONS +++
      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
      integer i, j, imax, jmax
c ^^^ tmp
c     integer i,j, imax, jmax, m, n

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 parallel wind direction (mm)
c     sxrg   =       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     sfvfs  =       soil fraction very fine sand by mass (0.05-0.1 mm)
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     lx     =       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
c     +++ PARAMETERS +++
      real  cs, cmp, ctf, cdp, c10dp
      parameter (cs=0.3, cmp=0.0001, ctf=1.2, 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
C     +++ LOCAL VARIABLES +++
c
      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
      real qav, szc, sz, crlos
      real dmlos, dmcld, szv, sff1, sff2, qitmp

      integer flg, kk, m,n
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     sarrc  =  soil angle of random roughness weibull 'c'(deg.)
c     sfa12  =  soil fraction area with shelter angle > 12 deg.
c     sfcv   =  soil fraction clod & crust cover
c     srrg   =  ratio of ridge spacing to ridge spacing parallel wind direction
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
c     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
c     calc. transport cap. for emission edit 5/30/01 LH
        tmp = 0
        qen = cs * wus * wus * (wus -  (wust - 0.05*tmp))
c     calc. transport cap. for trapping
        qcp = cs * wus * wus * (wus - (wusp - 0.05*tmp))
        if (qcp .lt. 0.0) then
           qcp = 0.001
        endif
c
c     store tmp qi
      qitmp = qi
c     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
c     test for saltation deposition region with suspension emission
c     (a more exact procedure would be to find new q(x) for deposition
c       and then put that solution into the suspension eqn and develop
c       another soln for the suspension eqn. For now we have set
c        qi to qen to calc. suspension and then reset qi to the real qi
c        value to calc. the deposition in sberod. This explanation added
c         12-5-2000 LH)
c
      elseif (qen .le. qi) then
        qi = qen
      endif
c
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

        sfcv = ((1.-sfcr)*(1.-sf84) + sfcr - sflos*sfcr)*(1.-svroc)
     &          + svroc
c
c         renb = exp(-4.0*sfcv)*exp(-2.5*sfa12)  edit 3-29-01 LH
        renb = exp(-3.0*sfcv)*(1.0 - sfa12)
        cen = ceno*renv*renb
c     soil fraction suspension in emitted soil; lower for loose on cr
        sfssen = (sf10/(sf200+0.001))*(1. - sfcr*0.8)
c     soil mixing parameter for suspension
        cm = cmp*sfssen
c
c     calc. trap params. for saltation/creep
c       calc. coef of stem interception
c       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
c
c     Calc. abrasion params. for saltation/creep
c       soil fraction saltation/creep in saltation
c       changed estimate  edit 5/30/01
       sfsn = sfa12**0.09
       sfsn = amax1(sfsn,0.7)
c
c       fraction of saltation abrading clods, crust and loose
      fan = (exp(- 4.0*bffcv))*(exp(-3*svroc))*sfsn
      fan = amax1(fan,0.)
c       fraction abrading aggregates
      fanag = (1. - sf84)*(1. - sfcr)*fan
c       fraction abrading crust
      fancr = (sfcr - sflos*sfcr)*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 =  (0.4 - 4.83*sfcla + 27.18*sfcla**2
     &         - 53.7*sfcla**3 + 42.25*sfcla**4 - 10.7*sfcla**5)
c
c     set upper limit on sfssan for high sand content
      sfssan = min(sfssan, (1.0 - (sfsan-sfvfs)))
c
c       calc suspension breakage coef. for saltation creep
c       edit 6-9-01 LH
      cbk = 0.11*canag*(1. - (sfsan-sfvfs))*sfsn
c
c       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
c
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
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.999999
         elseif (p .ge. 1.) then
            p = 0.999999
         endif
c        calculate atanh (p)
         t1 = 0.5*alog((1. + p)/(1. - p))
      qo = (s/(2.*c))*(-tanh((s/2.)*(-lx) + t1) + b/s)
c
c     assemble composite params. for suspension
      f = sfssen*cen*qen
c     added last term to intercept some ss LH edit  6-11-01
      g = sfssan*fancan - sfssen*cen + cbk + cm - sfssen*ci
c
c ^^^ tmp out
c     set counter
c      m = (imax-1)/8
c      m = max(1,m)
c
c      do 2 n = 1, (imax-1), m
c       if (i .eq. n  .and. j .eq. 2) then
c           write (*,*)
c           write (*,*) 'output from sbqout: i=',i,'j=',j
c           write (*,*) 'a=', a, 'b=', b, 'c=',c
c           write (*,*) 'g=', g, 'f=', f
c           write (*,*) 'cen=', cen, 'qen=', qen
c           write (*,*) 'sfssen=', sfssen, 'sfssan=',sfssan,
c     &      'fancan=',fancan
c           write (*,*) 'cbk=',cbk, 'ci=',ci,'ct=',ct,'cm=',cm
c       endif
c    2 continue

c ^^^ end tmp out
c
c     trap to prevent over range of exp
      if ((s*lx) .gt. 40.0) then
           s = 40.0/lx
      endif
c     solve for suspension out
c        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))
c
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
c     assemble composite params. for
      h = sf10en*f
      k = sf10an*sfssan*fancan - sf10en*sfssen*cen + sf10bk*cbk
     &   + sf10en*cm
c
c     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))

c     SURFACE UPDATE EQUATIONS

      flg = 1
c     execute this section if flag is set
      if (flg .eq. 1) then

c ^^^ tmp out
c     set counter
c     m = (imax-1)/8
c
c     do 2 n = 1, (imax-1), m
c       if (i .eq. n  .and. j .eq. 2) then
c           write (*,*) 'output from sbqout: i, j ', i, j
c       endif
c   2 continue
c
c ^^^ end tmp out
c
c     calc. avg. saltation discharge in grid cell
      qav = (qi + qo)/2.0

c     test for loss and abrasion masses
      szv = max (time*(qo-qi)/lx, time*fancan*qav)
c     calc no of iteration for surface update
      m = nint(min(szv/.02, 5.0))

c ^^^ tmp out
c      if (j .eq. 5) then
c        write (*,*) 'i=', i,'m=',m
c      endif
c ^^^ end tmp out

c
c     iterate to insure small step changes
      do 80 kk = 1,m
c     execute this section if crust present
      if (sfcr .gt. 0.01) then
c
c       update crust thickness
        szc = amax1(0.001, szcr)
        szcr = szcr - (fancr*cancr*qav/1.2)*time/m
        szcr = amax1(szcr, 0.00)
c
c       update crust (consolidated) zone cover
        sfcr = sfcr*szcr/szc
        sfcr = amax1(0.0, sfcr)
c
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
     &   -cen*(qen-qav)*sflos*sfcr*(1-svroc)
     &   /(1.001-sfcv))*time/m
c
         smlos = amax1(0.0, smlos)
c
c       update 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)

        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

c       update fraction saltation impacting crust
        fancr = (sfcr - sfcr*sflos)*fan
c
c ^^^ tmp out
c     set counter
c     m = (imax-1)/8
c
c     do 3 n = 1, (imax-1), m
c       if (i .eq. n  .and. j .eq. 2) then
c              write (*,200)  szc,szcr,sz,crlos
c200  format (1x, 'crust section: szc, szcr, sz, crlos', 8f8.3)
c       endif
c   3 continue
c ^^^ end tmp out
c
      endif
c


c     execute this section if clods are present
      if (sfcr .lt. 0.99) then
c
c       change in loose mass on aggregated surface, (+)=increase
c       terms: trapped, abraded but not emitted, emitted
c
        dmlos=((1-sfcr)*(1-svroc)*ct*qav
     &    + (1-sfssan)*fanag*canag*qav**2/qen
     &    -cen*(qen-qav)*(1.0-sfcr)*(1-svroc)*sf84
     &    /(1.001-sfcv))*time/m
c
c       change in abraded mass from clod area
c
        dmcld = (fanag*canag*qav)*time/m
c
c       net volume (height) removed from loose & clod areas
c         for aggregated area all erodible or all clods
c
c           szv = dmlos/((1.001-sfcr)*(1.001-svroc)*(sf84+0.001)*1.2)
c     &          + dmcld/((1.001-sf84)*(1.001-sfcr)*(1.001-svroc)*1.6)
c
c
c       specific clod area per unit volume (crude estimates)
c         1st case is all loose, but depleting to clods below
c          if (sf84 .gt. 0.99 .and. szv .lt. 0) then
c             sacd = 0.05
c          else
c             sacd = 1.2*exp(sf84 - 1.0)*(1-sf84)
c          endif
c
c         update surface loose material in zero layer
c         sf84 = sf84 + sacd*szv
c         sf84 = amax1(0.0, sf84)
c         sf84 = amin1(1.0, sf84)
c
c       revised sf84 update---
c
c      net change of loose mass at surface is then
       dmlos = dmlos + dmcld
c
c      case of depleting all loose to clods below
        if (sf84 .gt. 0.99 .and. dmlos .lt. 0.0) then
           sf84 = 0.98
        endif

c       update fraction <0.84 mm (loose)
c            based on Chepil trays at U*= 0.61 m/s
        sf84 = sf84 + ((1.0 - sf84)**2/0.4286)*dmlos
        sf84 = amin1(sf84, 1.0)
        sf84 = amax1(0.001, sf84)
c       may need to revise
        sf200 = sf84*alog(200.0)/alog(84.0)

        sf200 = amin1(sf200, 1.0)
        sf200 = amax1(0.001,sf200)

c       may need to revise
        sf10 = sf84*alog(10.0)/alog(84.0)

        fanag = (1-sf84)*(1-sfcr)*fan

      endif
c
c ^^^ tmp output
c     set counter
c     m = (imax-1)/8
c
c     do 4 n = 1, (imax-1), m
c       if (i .eq. n  .and. j .eq. 2) then
c
c      write (*,210)   dmlos, dmcld, szv, sf84
c 210 format(1x,'clod section: net_dmlos, dmcld, szv, sf84', 8f8.3)
c       endif
c   4 continue
c^^^ end tmp out
c
c
c     update surface roughness

c        set limits on areas in denominators
c         surface fraction stable that can change height
c         sff1 = amax1((sfcv-svroc), 0.01)
c         surface fraction loose
c         sff2 = amax1((1-sfcv), 0.01)
c
c     rate of volume change caused by emission
c        szv= cen*(qen - qav)/(sff2*1.2)

          szv= cen*(qen - qav)/1.0
c
c     if trapping then emission lowers roughness,
c          i.e. it comes from highest areas.
      If (ct .gt. 0.) then
        szv= -szv
      endif
c     calc. change in mean surface depth
c     terms: emission, trapped, abrasion, & abrasion not emitted
c     szv = (szv - ct*qav/0.4 - fancan*qav/sff1
c    &      -(1.0 - sfssen)*fancan*qav*qav/(qen*0.4))*time

       szv = (szv - ct*qav/1.0 - fancan*qav/1.2
     &      -(1.0 - sfssen)*fancan*qav*qav/(qen*1.0))*time/m

c
c ^^^ tmp out
c     set counter
c     m = (imax-1)/8
c
c     do 5 n = 1, (imax-1), m
c       if (i .eq. n  .and. j .eq. 2) then
c
c         write(*,220)  szv, sff1, sff2
c 220     format(1x, 'roughness section: szv, sff1, sff2', 8f8.3)
c         write(*,*)
c       endif
c   5 continue
c^^^ end tmp out
c
c     update ridge height
      if (szrgh .gt. 10.) then
          szrgh = szrgh + szv
c     set lower limit on ridge height
          szrgh = amax1(szrgh, 0.0)
      endif
c     update random roughness
         slrr = slrr + szv/4.
c    set lower limit on slrr
         slrr = amax1(slrr, 1.5)
c
c
c^^^^   tmp output section
c     set counter
c     m = (imax-1)/8
c
c     jmax = jmax
c
c     do 10 n = 1, (imax-1), m
c       if (i .eq. n  .and. j .eq. 2) then
c         write (*,19) sfssan, sfssen, cancr, canag, sfcv
c         write (*,20) wus, wust, ct, fancr, fanag, fancan
c         write (*,21) qo, qi, qav, qen,  cen,  time
c         write (*,*)
c       endif
c  10 continue
c  19 format(1x, 'sfssan, sfssen, cancr, canag, sfcv', 8f8.3)
c
c  20 format(1x, ' wus, wust, ct, fancr, fanag, fancan', 8f8.3)
c  21 format(1x, 'qo qi qav  qen  cen  time',  5f8.3, f8.0)

c^^^   end tmpout

   80 continue
      endif
  100 continue
      qi = qitmp
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++