!$Author: wagner $
!$Date: 2007-09-13 02:12:56 $
!$Revision: 1.1 $
!$Source: /weru/cvs/weps/weps_ref/src/emission_parms.f95,v $
!-----------------------------------------------------------------------
!! Compute emission parameters.
!!


MODULE Emission_Parameters

  IMPLICIT NONE

    PUBLIC :: Coef_Emission

    PRIVATE :: SF84_from_lognormal_dist

    PUBLIC :: WUst_d

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



CONTAINS

!! Calculate the coefficient of emission, Cen (1/m)
REAL FUNCTION Coef_Emisson (Ceno, Renb, Renv) RESULT(Cen)

  REAL, INTENT (IN) :: Ceno    !! Coefficient of emission for a bare, loose, erodible soil (1/m)
                               !! A typical field value is about 0.06 (1/m)
  REAL, INTENT (IN) :: Renb    !! Emission ratio for flat, random vegetation
  REAL, INTENT (IN) :: Renv    !! Emission ratio for immobile cover and surface roughness

  Cen = Ceno * Renb * Renv

END FUNCTION Coef_Emission

!! Calculate the fractional emission factor to account for flat biomass cover
REAL FUNCTION EmissionRatio_flatcover (BFFcv) RESULT(Renv)

  REAL, INTENT (IN) :: BFFcv    !! biomass fraction of flat cover

  Renv = 0.075 + 0.934 * exp(BFFcv/0.149)

END FUNCTION EmissionRatio_flatcover

!! Calculate the fractional emission factor to account for immobile cover (non-vegetative) and surface roughness
REAL FUNCTION EmissionRatio_roughness (SFcv, sfa12) RESULT(Renb)

  REAL, INTENT (IN) :: SFcv    !! _non-emitting_ surface fraction, e.g., fraction not covered by clods, crust and rock 
  REAL, INTENT (IN) :: sfa12   !! fraction of area with shelter > 12 degrees

  Renb = exp(-3.0*SFcv) * (1.0 - sfa12)

END FUNCTION EmissionRatio_roughness

!! 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

!! Calculate Weibull scale factor for shelter angle (degrees) of random roughness (mm)
REAL FUNCTION WeibullSF__rr (SZrr) RESULT(SACrr)

  REAL, INTENT (IN) :: SZrr    !! Random Roughness (mm) 

  SACrr = 2.3 * sqrt(SZrr)

END FUNCTION WeibullSF__rr

!! Calculate fraction of horizontal area with shelter angles greater than 12 degrees caused by random rughness
REAL FUNCTION SFA12__rr (SACrr) RESULT(sfa12rr)

  REAL, INTENT (IN) :: SACrr    !! Weibull scale factor for shelter angle (degrees) of random roughness (mm)

  sfa12rr = exp(-12.0/SACrr)**0.77

END FUNCTION SFA12_rr

END MODULE Emission_Parameters