c**********************************************************************
c     subroutine sbwus
c**********************************************************************
      subroutine sbwus
     i (anemht, awzzo, awu, wzzov,brlai, brsai,
     o  wus)
c
c     +++ PURPOSE +++
c     To calculate subregion, friction velocity, given station
c     anemometer height, surface roughness, wind speed; and subregion
c     aerodynamic roughness.
c
c     if standing biomass present, then calculate friction velocity
c     at surface below the canopy (wus).
c
c     +++ ARGUMENT DECLARATIONS +++
      real  anemht, awzzo, awu, wzzov
      real  brlai, brsai, wus
c
c     +++ ARGUMENT DEFINITIONS +++
c
c     anemht - parameter, anemometer height of input wind speed (m).
c     awzzo - parameter, surface aerodynamic roughness at input wind
c             speed location (mm).
c     awu - input wind speed driving EROSION submodel (m/s).
c     wzzov - subregion aerodynamic roughness (mm).
c     wus - subregion soil surface friction velocity (m/s)
c           i.e. below canopy, if one exists.
c
c     +++ LOCAL VARIABLES +++
      real wusst, wusv, brcd
c
c     +++ END SPECIFICATIONS +++
c     note:  in BLOCK.FOR wzoflg should be set to 1 and anemomht
c             set to correct height if anemometer is at field site
c             to obtain correct values from SBWUS or read as
c             input data in stand-alone EROSION.
c
c     Calc station (input wind speed location) friction velocity
      wusst = awu*0.4/alog(anemht*1000./awzzo)
c
c     calc subregion friction velocity
      wus = wusst * (wzzov/awzzo)**0.067
c
c     if standing biomass, calculate wus below canopy
      if (brlai .gt. 0.01 .or. brsai .gt. 0.002 ) then
         wusv = wus
c
c        calculate "effective" biomass drag coefficient
         brcd = brlai*(0.2 - 0.15*exp(-8.*brlai)) + brsai
         brcd = amax1(brcd, 0.001)
c
c        calculate friction velocity below canopy

        if( brcd.gt.2.56) then       !check to avoid underflow
            wus = wusv * 0.25*exp(-brcd/0.356)
        else
             wus = wusv*0.86*exp(-brcd/0.0298)+0.25*exp(-brcd/0.356)
        endif
c      edit 11-30-2000 LH
         wus = amin1(wus,wusv)
c ^tmp
c         write(*,*) 'awu=',awu,'wzzov=',wzzov,'wusv=',wusv,'wus=',wus
c
      endif
c
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++