!$Author: wagner $
!$Date: 2007-09-13 02:12:56 $
!$Revision: 1.1 $
!$Source: /weru/cvs/weps/weps_ref/src/thresh_friction_vel.f95,v $
!-----------------------------------------------------------------------
!! Compute threshold friction velocities for specified surface conditions.
!!


MODULE Threshold_Friction_Velocity

  IMPLICIT NONE

    PUBLIC :: SFcv

    PRIVATE :: SF84_from_lognormal_dist

    PUBLIC :: WUst_d

    PUBLIC :: WUBts
    PUBLIC :: SFCcv
    PUBLIC :: WUCts
    PUBLIC :: WUCSts
    PUBLIC :: WUt



CONTAINS

!! Calculate the _non-emitting_ surface fraction, e.g., fraction not covered by clods, crust and rock
REAL FUNCTION SFcv (SFcr, SFlos, SVroc, SF84)

  REAL, INTENT (IN) :: SFcr     !! surface fraction covered by crust (but excluding fraction of rock covered area)
  REAL, INTENT (IN) :: SFlos    !! surface fraction covered with loose, erodible material on the crusted area
  REAL, INTENT (IN) :: SVroc    !! surface volume of rock > 2.0 mm in size (m^3/m^3)
  REAL, INTENT (IN) :: SF84     !! surface fraction covered with aggregates < 0.84 mm in diameter in the non-crusted area (but excluding fraction of rock covered area)

  SFcv = ((1.0-SFcr)*(1.0-SF84) + SFcr - SFlos) * (1.0 - SVroc) + SVroc

END FUNCTION SFcv

!! Calculate the fraction of aggregates less than 0.84 mm diameter in size
!! from the modified lognormal description of the aggregate size distribution
REAL FUNCTION SF84_from_lognormal_dist (SLagn, SLagx, SLagm, SOags) RESULT(SF84)

  REAL, INTENT (IN) :: SLagn    !! lower limit of size distribution (mm)
  REAL, INTENT (IN) :: SLagx    !! upper limit of size distribution (mm)
  REAL, INTENT (IN) :: SLagm    !! geometric mean of size distribution (mm)
  REAL, INTENT (IN) :: SOags    !! geometric standard deviation of size distribution

  REAL :: SLT

  SLT = ((0.84 - SLagn) * (SLagx - SLagn)) / ((SLagx - 0.84) * SLagm)

  SF84 = 0.5 * (1.0 + erf(log(SLT)/(sqrt(2.0)* log(SOags))))

END FUNCTION SF84_from_lognormal_dist

!! Calculate the _dynamic_ saltation/creep threshold friction velocity (m/s)
! in the control volume
REAL FUNCTION WUst_d (wust, sfa12)

  REAL, INTENT (IN) :: wust    ! emission threshold friction velocity (m/s)
  REAL, INTENT (IN) :: sfa12   ! fraction of area with shelter > 12 degrees


  WUst_d = wust - (0.05 * (1.0 - sfa12))  !Why was 5% value chosen?

END FUNCTION WUst_d

!! Static threshold friction velocity (m/s) for bare surface
!! Equation was derived from fitted wind tunnel data
REAL FUNCTION WUBts (SFcv, wzo)
  REAL, INTENT (IN) :: SFcv     !! _non-emitting_ surface fraction, e.g., fraction not covered by clods, crust and rock 
  REAL, INTENT (IN) :: wzo      !! aerodynamic surface roughness (mm)

  REAL :: b1, b2

  b1 = -0.179 + 0.225*((log(1.0+wzo))**0.891)
  b2 = 0.3 + 0.06*(wzo**1.2)

  WUBts = 1.7 - 1.35*exp(-b1-b2*(SFcv*SFcv))

END FUNCTION WUbts

!! Fraction change in soil surface area protected from emission
REAL FUNCTION SFCcv (SFcv, BFFcv)
  REAL, INTENT (IN) :: SFcv     !! _non-emitting_ surface fraction, e.g., fraction not covered by clods, crust and rock 
  REAL, INTENT (IN) :: BFFcv    !! biomass fraction of flat cover

  SFCcv = (1.0 - SFcv) * BFFcv

END FUNCTION SFCcv

!! Change in static threshold friction velocity caused by flat biomass cover (m/s)
REAL FUNCTION WUCts (SFC_cv)
  REAL, INTENT (IN) :: SFC_cv    !!  Fraction change in soil surface area protected from emission

  WUCts = 0.02 + SFC_cv    ! SF_cv > 0.0

END FUNCTION WUCts

!! Increase in static threshold friction velocity from surface wetness (m/s)
REAL FUNCTION WUCSts (HR0wc, HR15wc)
  REAL, INTENT (IN) :: HR0wc    !!  Surface soil water content (kg/kg)
  REAL, INTENT (IN) :: HR15wc   !!  Surface soil water content at 1.5 MPa (kg/kg)

  WUCSts = (1.0/(11.91-10.41*((HR0wc/HR15wc)**0.5)))    ! HROwc/HR15wc > 0.25

END FUNCTION WUCSts


!! Fraction change in soil surface area protected from emission
REAL FUNCTION WUt (WUts, SFA12)
  REAL, INTENT (IN) :: WUts    !! Surface static threshold friction velocity accounting for both flat biomass cover and wetness effects (m/s)
  REAL, INTENT (IN) :: SFA12   !! Fraction of area with shelter > 12 degrees

  WUt = WUts - 0.05*(1.0 - SFA12)

END FUNCTION WUt

END MODULE Threshold_Friction_Velocity