NESDIS_LandEM_Module.f90

Go to the documentation of this file.

!
! NESDIS_LandEM_Module
!
! Module containing the NESDIS microwave land emissivity model
!
!
! CREATION HISTORY:
!       Written by:     Banghua Yan, QSS Group Inc., 01-Jun-2005
!                       Banghua.Yan@noaa.gov
!                       Fuzhong Weng, NOAA/NESDIS/ORA,
!                       Fuzhong.Weng@noaa.gov
!

MODULE nesdis_landem_module
 
  ! -----------------
  ! Enviroment set up
  ! -----------------
  ! Module use
  USE type_kinds, ONLY: fp
  USE fundamental_constants, ONLY: c_2
  USE nesdis_snowem_parameters
  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  ! Procedures
  PUBLIC  :: nesdis_landem
  ! Parameters
  PUBLIC :: zero
  PUBLIC :: point1
  PUBLIC :: point5
  PUBLIC :: one
  PUBLIC :: two
  PUBLIC :: three
  PUBLIC :: four
  PUBLIC :: pi
  PUBLIC :: emissh_default
  PUBLIC :: emissv_default
  
  PUBLIC :: one_tenth
  PUBLIC :: half
  
  ! -----------------
  ! Module parameters
  ! -----------------
  REAL(fp), PARAMETER :: zero   = 0.0_fp
  REAL(fp), PARAMETER :: point1 = 0.1_fp
  REAL(fp), PARAMETER :: point5 = 0.5_fp
  REAL(fp), PARAMETER :: one    = 1.0_fp
  REAL(fp), PARAMETER :: two    = 2.0_fp
  REAL(fp), PARAMETER :: three  = 3.0_fp
  REAL(fp), PARAMETER :: four   = 4.0_fp
  REAL(fp), PARAMETER :: pi     = 3.141592653589793238462643_fp
  REAL(fp), PARAMETER :: twopi  = two*pi
  REAL(fp), PARAMETER :: emissh_default = 0.25_fp
  REAL(fp), PARAMETER :: emissv_default = 0.30_fp

  REAL(fp), PARAMETER :: one_tenth = point1
  REAL(fp), PARAMETER :: half      = point5
  
  
CONTAINS


!################################################################################
!################################################################################
!##                                                                            ##
!##                         ## PUBLIC MODULE ROUTINES ##                       ##
!##                                                                            ##
!################################################################################
!################################################################################

!--------------------------------------------------------------------------------
!
! NAME:
!       NESDIS_LandEM
!
! PURPOSE:
!       Subroutine to simulate microwave emissivity over land conditions.
!
! REFERENCES:
!       Weng, F., B. Yan, and N. Grody, 2001: "A microwave land emissivity model",
!         J. Geophys. Res., 106, 20, 115-20, 123
!
! CALLING SEQUENCE:
!       CALL NESDIS_LandEM(Angle,                 &   ! Input
!                          Frequency,             &   ! Input
!                          Soil_Moisture_Content, &   ! Input
!                          Vegetation_Fraction,   &   ! Input
!                          Soil_Temperature,      &   ! Input
!                          Land_Temperature,      &   ! Input
!                          Snow_Depth,            &   ! Input
!                          Emissivity_H,          &   ! Output
!                          Emissivity_V)              ! Output
!
! INPUT ARGUMENTS:
!         Angle:                   The angle values in degree.
!                                  UNITS:      Degrees
!                                  TYPE:       REAL(fp)
!                                  DIMENSION:  Rank-1, (I)
!
!         Frequency                Frequency User defines
!                                  This is the "I" dimension
!                                  UNITS:      GHz
!                                  TYPE:       REAL(fp)
!                                  DIMENSION:  Scalar
!
!         Soil_Moisture_Content:   The volumetric water content of the soil (0:1).
!                                  UNITS:      N/A
!                                  TYPE:       REAL(fp)
!                                  DIMENSION:  Scalar
!
!         Vegetation_Fraction:     The vegetation fraction of the surface (0:1).
!                                  UNITS:      N/A
!                                  TYPE:       REAL(fp)
!                                  DIMENSION:  Scalar
!
!         Soil_Temperature:        The soil temperature.
!                                  UNITS:      Kelvin, K
!                                  TYPE:       REAL(fp)
!                                  DIMENSION:  Scalar
!
!         Land_Temperature:        The land surface temperature.
!                                  UNITS:      Kelvin, K
!                                  TYPE:       REAL(fp)
!                                  DIMENSION:  Scalar
!
!         Snow_Depth:              The snow depth.
!                                  UNITS:      mm
!                                  TYPE:       REAL(fp)
!                                  DIMENSION:  Scalar
!
! OUTPUT ARGUMENTS:
!         Emissivity_H:            The surface emissivity at a horizontal
!                                  polarization.
!                                  UNITS:      N/A
!                                  TYPE:       REAL(fp)
!                                  DIMENSION:  Scalar
!
!         Emissivity_V:            The surface emissivity at a vertical polarization.
!                                  UNITS:      N/A
!                                  TYPE:       REAL(fp)
!                                  DIMENSION:  Scalar
!
!
! INTERNAL ARGUMENTS:
!       theta       -  local zenith angle in radian
!       rhob        -  bulk volume density of the soil (1.18-1.12)
!       rhos        -  density of the solids (2.65 g.cm^3 for solid soil material)
!       sand        -  sand fraction (sand + clay = 1.0)
!       clay        -  clay fraction
!       lai         -  leaf area index (eg. lai = 4.0 for corn leaves)
!       sigma       -  surface roughness formed between medium 1 and 2,
!                      expressed as the standard deviation of roughtness height (mm)
!       leaf_thick  --  leaf thickness (mm)
!       rad         -  radius of dense medium scatterers (mm)
!       va          -  fraction volume of dense medium scatterers(0.0 - 1.0)
!       ep          -  dielectric constant of ice or sand particles, complex value
!                               (e.g, 3.0+i0.0)
!
! CREATION HISTORY:
!       Written by:     Banghua Yan, QSS Group Inc., 16-May-2005
!                       Banghua.Yan@noaa.gov
!                       Fuzhong Weng, NOAA/NESDIS/ORA,
!                       Fuzhong.Weng@noaa.gov
!
!------------------------------------------------------------------------------------------------------------

  SUBROUTINE nesdis_landem(Angle,                 &   ! Input
                           Frequency,             &   ! Input
                           Soil_Moisture_Content, &   ! Input
                           Vegetation_Fraction,   &   ! Input
                           Soil_Temperature,      &   ! Input
                           t_skin,                &   ! Input
                           Lai,                   &   ! Input
                           Soil_Type,             &   ! Input
                           Vegetation_Type,       &   ! Input
                           Snow_Depth,            &   ! Input
                           Emissivity_H,          &   ! Output
                           Emissivity_V)              ! Output
    ! Arguments
    REAL(fp), intent(in) :: angle
    REAL(fp), intent(in) :: frequency
    REAL(fp), intent(in) :: soil_moisture_content
    REAL(fp), intent(in) :: vegetation_fraction
    REAL(fp), intent(in) :: soil_temperature
    REAL(fp), intent(in) :: t_skin
    REAL(fp), intent(in) :: lai
    INTEGER,  intent(in) :: soil_type
    INTEGER,  intent(in) :: vegetation_type
    REAL(fp), intent(in) :: snow_depth
    REAL(fp), intent(out):: emissivity_v,emissivity_h
    ! Local parameters
    REAL(fp), PARAMETER :: snow_depth_c     = 10.0_fp
    REAL(fp), PARAMETER :: tsoilc_undersnow = 280.0_fp
    REAL(fp), PARAMETER :: rhos = 2.65_fp
    REAL(fp)            :: sand, clay, rhob
    REAL(fp), PARAMETER, dimension(0:9) :: frac_sand = (/ 0.80_fp,     &
                          0.92_fp, 0.10_fp, 0.20_fp, 0.51_fp, 0.50_fp, &
                          0.35_fp, 0.60_fp, 0.42_fp,  0.92_fp  /)
    REAL(fp), PARAMETER, dimension(0:9) :: frac_clay = (/ 0.20_fp,     &
                          0.06_fp, 0.34_fp, 0.63_fp, 0.14_fp, 0.43_fp, &
                          0.34_fp, 0.28_fp, 0.085_fp, 0.06_fp /)
    REAL(fp), PARAMETER, dimension(0:9) :: rhob_soil = (/ 1.48_fp,     &
                          1.68_fp, 1.27_fp, 1.21_fp, 1.48_fp, 1.31_fp, &
                          1.32_fp, 1.40_fp, 1.54_fp, 1.68_fp /)
! Specific Density
    REAL(fp), PARAMETER, dimension(0:13) :: veg_rho  = (/ 0.33_fp,     &
                          0.40_fp, 0.40_fp, 0.40_fp, 0.40_fp, 0.40_fp, &
                          0.25_fp, 0.25_fp, 0.40_fp, 0.40_fp, 0.40_fp, &
                          0.40_fp, 0.33_fp, 0.33_fp            /)
! MGE
    REAL(fp), PARAMETER, dimension(0:13) :: veg_mge  = (/ 0.50_fp,     &
                          0.45_fp, 0.45_fp, 0.45_fp, 0.40_fp, 0.40_fp, &
                          0.30_fp, 0.35_fp, 0.30_fp, 0.30_fp, 0.40_fp, &
                          0.30_fp, 0.50_fp, 0.40_fp            /)
! LAI
    REAL(fp), PARAMETER, dimension(0:13) :: lai_min  = (/ 0.52_fp,     &
                          3.08_fp, 1.85_fp, 2.80_fp, 5.00_fp, 1.00_fp, &
                          0.50_fp, 0.52_fp, 0.60_fp, 0.50_fp, 0.60_fp, &
                          0.10_fp, 1.56_fp, 0.01_fp            /)
    REAL(fp), PARAMETER, dimension(0:13) :: lai_max  = (/ 2.90_fp,     &
                          6.48_fp, 3.31_fp, 5.50_fp, 6.40_fp, 5.16_fp, &
                          3.66_fp, 2.90_fp, 2.60_fp, 3.66_fp, 2.60_fp, &
                          0.75_fp, 5.68_fp, 0.01_fp            /)
! Leaf_thickness
    REAL(fp), PARAMETER, dimension(0:13) :: leaf_th  = (/ 0.07_fp,     &
                          0.18_fp, 0.18_fp, 0.18_fp, 0.18_fp, 0.18_fp, &
                          0.12_fp, 0.12_fp, 0.12_fp, 0.12_fp, 0.12_fp, &
                          0.12_fp, 0.15_fp, 0.12_fp            /)
    ! Local variables
    REAL(fp) :: mv,veg_frac,theta,theta_i,theta_t,mu,r21_h,r21_v,r23_h,r23_v,  &
                t21_v,t21_h,gv,gh,ssalb_h,ssalb_v,tau_h,tau_v,mge, &
                    leaf_thick,rad,sigma,va,ep_real,ep_imag
    REAL(fp) :: t_soil
    REAL(fp) :: rhoveg, vlai
    REAL(fp) :: local_snow_depth
    COMPLEX(fp) :: esoil, eveg, esnow, eair
    LOGICAL :: snowem_physical_model

    eair = cmplx(one,-zero,fp)
    theta = angle*pi/180.0_fp

    ! By default use the 
    ! Assign local variable
    mv               = soil_moisture_content
    veg_frac         = vegetation_fraction
    t_soil           = soil_temperature
    sand = frac_sand(soil_type)
    clay = frac_clay(soil_type )
    rhob = rhob_soil(soil_type )
    local_snow_depth = snow_depth

    ! Check soil/skin temperature
    if ( (t_soil <= 100.0_fp .OR.  t_soil >= 350.0_fp) .AND. &
         (t_skin >= 100.0_fp .AND. t_skin <= 350.0_fp) ) t_soil = t_skin

    ! Check soil moisture content range
    mv = max(min(mv,one),zero)

    ! Surface type based on snow depth
    IF (local_snow_depth > point1) THEN

      ! O.k.; we're going to compute snow emissivities....
      
      ! By default use the physical model for snow
      snowem_physical_model = .true.
      if (local_snow_depth > snow_depth_c) snowem_physical_model = .false.

      ! Compute the snow emissivity
      IF ( snowem_physical_model ) THEN

        ep_real = 3.2_fp
        ep_imag = -0.0005_fp
        sigma = one

        ! For deep snow, the performance of the model is poor
        local_snow_depth = min(local_snow_depth,1000.0_fp)

        ! The fraction volume of dense medium
        ! scatterers must be in the range (0-1)
        va = 0.4_fp + 0.0004_fp*local_snow_depth
        va = max(min(va,one),zero)

        ! Limit for snow grain size
        rad = min((point5 + 0.005_fp*local_snow_depth),one)

        ! Limit for soil temperature
        t_soil = min(t_soil,tsoilc_undersnow)

        CALL snow_diel(frequency, ep_real, ep_imag, rad, va, esnow)
        CALL soil_diel(frequency, t_soil, mv, rhob, rhos, sand, clay, esoil)

        theta_i = asin(REAL(SIN(theta)*SQRT(eair)/SQRT(esnow),fp))

        CALL reflectance(esnow, eair, theta_i,  theta, r21_v, r21_h)
        CALL transmittance(esnow, eair, theta_i, theta, t21_v, t21_h)

        mu      = cos(theta_i)
        theta_t = asin(REAL(SIN(theta_i)*SQRT(esnow)/SQRT(esoil),fp))

        CALL reflectance(esnow, esoil, theta_i, theta_t, r23_v, r23_h)
        CALL roughness_reflectance(frequency, sigma, r23_v, r23_h)
        CALL snow_optic(frequency,rad,local_snow_depth,va,ep_real, ep_imag,gv,gh,&
                        ssalb_v,ssalb_h,tau_v,tau_h)
        CALL two_stream_solution(mu,gv,gh,ssalb_h,ssalb_v,tau_h,tau_v, &
                                 r21_h,r21_v,r23_h,r23_v,t21_v,t21_h,emissivity_v,emissivity_h, &
                               frequency, t_soil, t_skin)
      ELSE
        ! Use the empirical method 
        CALL snowem_default(frequency,t_skin, local_snow_depth,emissivity_v,emissivity_h)
      END IF

    ELSE

      ! No snow, so we're going to compute canopy emissivities....
      
      ! Limit for vegetation fraction
      veg_frac = max(min(veg_frac,one),zero)

!     lai = THREE*veg_frac + POINT5
!     mge = POINT5*veg_frac
!     leaf_thick = 0.07_fp
      mu  = cos(theta)
      sigma = point5
    
      vlai = lai*veg_frac
      mge = veg_mge(vegetation_type)
      rhoveg = veg_rho(vegetation_type)
      leaf_thick = leaf_th(vegetation_type)

      r21_h    = zero
      r21_v    = zero
      t21_h    = one
      t21_v    = one

      CALL soil_diel(frequency, t_soil, mv, rhob, rhos, sand, clay, esoil)
      theta_t = asin(REAL(SIN(theta)*SQRT(eair)/SQRT(esoil),fp))
      CALL reflectance(eair, esoil, theta, theta_t, r23_v, r23_h)
      CALL roughness_reflectance(frequency, sigma, r23_v, r23_h)
      CALL canopy_diel(frequency, mge, eveg, rhoveg)
      CALL canopy_optic(vlai,frequency,theta,eveg,leaf_thick,gv,gh,ssalb_v,ssalb_h,tau_v,tau_h)
      CALL two_stream_solution(mu,gv,gh,ssalb_h,ssalb_v,tau_h,tau_v, &
                               r21_h,r21_v,r23_h,r23_v,t21_v,t21_h,emissivity_v,emissivity_h, &
                               frequency, t_soil, t_skin)
    END IF

  END SUBROUTINE nesdis_landem



!##################################################################################
!##################################################################################
!##                                                                              ##
!##                          ## PRIVATE MODULE ROUTINES ##                       ##
!##                                                                              ##
!##################################################################################
!##################################################################################



subroutine snowem_default(frequency,ts, depth,Emissivity_V,Emissivity_H)

!----------------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .
!   prgmmr:  Banghua Yan and Fuzhong Weng               org: nesdis              date: 2005-12-01
!
! abstract: preliminary estimate of snow emissivity using  surface temperature and snow depth
!
! input argument list:
!
!      ts         -  surface temperature
!      frequency   -  frequency (ghz)
!
! output argument list:
!
!      Emissivity_V         -  snow emissivty at V-POL
!      Emissivity_H         -  snow emissivty at H-POL
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!------------------------------------------------------------------------------------------------------------

  ! Arguments
  REAL(fp) :: frequency,ts, depth,Emissivity_V,Emissivity_H
  ! Local parameters  
  INTEGER , PARAMETER :: new = 7
  INTEGER , PARAMETER :: NFRESH_SHALLOW_SNOW = 1
  INTEGER , PARAMETER :: NPOWDER_SNOW        = 2
  INTEGER , PARAMETER :: NWET_SNOW           = 3
  INTEGER , PARAMETER :: NDEEP_SNOW          = 4
  REAL(fp), PARAMETER :: twet    = 270.0_fp
  REAL(fp), PARAMETER :: tcrust  = 235.0_fp
  REAL(fp), PARAMETER :: depth_s =  50.0_fp
  REAL(fp), PARAMETER :: depth_c = 100.0_fp
  ! Local variables
  INTEGER :: ich,basic_snow_type
  REAL(fp), DIMENSION(new) :: ev, eh, freq
  REAL(fp) :: df, df0


  freq = frequency_amsre(1:new)
  
  ! Determine the snow type based on temperatures
  basic_snow_type = nfresh_shallow_snow
  if (ts >= twet .and. depth <= depth_s) then
    basic_snow_type = nwet_snow
  else
    if (depth <= depth_s) then
      basic_snow_type = nfresh_shallow_snow
    else
      basic_snow_type = npowder_snow
    endif
  endif
  if (ts <= tcrust .and. depth >= depth_c) basic_snow_type = ndeep_snow

  ! Assign the emissivity spectrum
  SELECT CASE (basic_snow_type)
    CASE (nfresh_shallow_snow)
      ev = grass_after_snow_ev_amsre(1:new)
      eh = grass_after_snow_eh_amsre(1:new)
    CASE (npowder_snow)
      ev = powder_snow_ev_amsre(1:new)
      eh = powder_snow_eh_amsre(1:new)
    CASE (nwet_snow)
      ev = wet_snow_ev_amsre(1:new)
      eh = wet_snow_eh_amsre(1:new)
    CASE (ndeep_snow)
      ev = deep_snow_ev_amsre(1:new)
      eh = deep_snow_eh_amsre(1:new)
  END SELECT

  ! Handle possible extrapolation
  if (frequency <= freq(1)) then
    emissivity_h = eh(1)
    emissivity_v = ev(1)
    return
  end if
  if (frequency >= freq(new)) then
    emissivity_h = eh(new)
    emissivity_v = ev(new)
    return
  end if

  ! Interpolate emissivity at a certain frequency
  channel_loop: do ich=2,new
    if (frequency <= freq(ich)) then
      df  = frequency-freq(ich-1)
      df0 = freq(ich)-freq(ich-1)
      emissivity_h = eh(ich-1) + (df*(eh(ich)-eh(ich-1))/df0)
      emissivity_v = ev(ich-1) + (df*(ev(ich)-ev(ich-1))/df0)
      exit channel_loop
    end if
  end do channel_loop

end subroutine snowem_default


subroutine canopy_optic(vlai,frequency,theta,esv,d,gv,gh,&
                        ssalb_v,ssalb_h,tau_v, tau_h)

!----------------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    canopy_optic compute optic parameters for canopy
!
!   prgmmr:  Fuzhong Weng and Banghua Yan                org: nesdis              date: 2000-11-28
!
! abstract: compute optic parameters for canopy
!
! program history log:
!
! input argument list:
!
!      lai         -  leaf area index
!      frequency   -  frequency (ghz)
!      theta       -  incident angle
!      esv         -  leaf dielectric constant
!      d           -  leaf thickness (mm)
!
! output argument list:
!
!      gv           -  asymmetry factor for v pol
!      gh           -  asymmetry factor for h pol
!      ssalb_v      -  single scattering albedo at v. polarization
!      ssalb_h      -  single scattering albedo at h. polarization
!      tau_v        -  optical depth at v. polarization
!      tau_h        -  optical depth at h. polarization
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!------------------------------------------------------------------------------------------------------------

  REAL(fp) :: frequency,theta,d,vlai,ssalb_v,ssalb_h,tau_v,tau_h,gv, gh, mu
  COMPLEX(fp) :: ix,k0,kz0,kz1,rhc,rvc,esv,expval1,factt,factrvc,factrhc
  REAL(fp) :: rh,rv,th,tv
  REAL(fp), PARAMETER :: threshold = 0.999_fp

  mu = cos(theta)
  ix = cmplx(zero, one, fp)

  k0  = cmplx(twopi*frequency/300.0_fp, zero, fp)   ! 1/mm
  kz0 = k0*mu
  kz1 = k0*sqrt(esv - sin(theta)**2)

  rhc = (kz0 - kz1)/(kz0 + kz1)
  rvc = (esv*kz0 - kz1)/(esv*kz0 + kz1)

  expval1 = exp(-two*ix*kz1*d)
  factrvc = one-rvc**2*expval1
  factrhc = one-rhc**2*expval1
  factt   = four*kz0*kz1*exp(ix*(kz0-kz1)*d)

  rv = abs(rvc*(one - expval1)/factrvc)**2
  rh = abs(rhc*(one - expval1)/factrhc)**2

  th = abs(factt/((kz1+kz0)**2*factrhc))**2
  tv = abs(esv*factt/((kz1+esv*kz0)**2*factrvc))**2

  gv = point5
  gh = point5

  tau_v = point5*vlai*(two-tv-th)
  tau_h = tau_v

  ssalb_v = min((rv+rh)/(two-tv-th),threshold)
  ssalb_h = ssalb_v

end subroutine canopy_optic


subroutine snow_optic(frequency,a,h,f,ep_real,ep_imag,gv,gh, ssalb_v,ssalb_h,tau_v,tau_h)

!-------------------------------------------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    landem      comput optic parameters for snow
!
!   prgmmr: Fuzhong Weng and Banghua Yan                 org: nesdis              date: 2000-11-28
!
! abstract: compute optic parameters for snow
!
! program history log:
!
! input argument list:
!
!      theta        -  local zenith angle (degree)
!      frequency    -  frequency (ghz)
!      ep_real      -  real part of dielectric constant of particles
!      ep_imag      -  imaginary part of dielectric constant of particles
!      a            -  particle radiu (mm)
!      h            -  snow depth(mm)
!      f            -  fraction volume of snow (0.0 - 1.0)
!
! output argument list:
!
!       ssalb       -  single scattering albedo
!       tau         -  optical depth
!       g           -  asymmetry factor
!
!   important internal variables:
!
!       ks          -  scattering coeffcient (/mm)
!       ka          -  absorption coeffient (/mm)
!       kp          -  eigenvalue of two-stream approximation
!       y           -  = yr+iyi
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------

  REAL(fp) :: yr,yi,ep_real,ep_imag
  REAL(fp) :: frequency,a,h,f,ssalb_v,ssalb_h,tau_v,tau_h,gv,gh,k
  REAL(fp) :: ks1,ks2,ks3,ks,kr1,kr2,kr3,kr,ki1,ki2,ki3,ki
  REAL(fp) :: fact1,fact2,fact3,fact4,fact5

  k = twopi/(300._fp/frequency)

  yr = (ep_real - one)/(ep_real + two)
  yi = -ep_imag/(ep_real + two)

  fact1 = (one+two*f)**2
  fact2 = one-f*yr
  fact3 = (one-f)**4
  fact4 = f*(k*a)**3
  fact5 = one+two*f*yr

  ks1 = k*sqrt(fact2/fact5)
  ks2 = fact4*fact3/fact1
  ks3 = (yr/fact2)**2
  ks = ks1*ks2*ks3

  kr1 = fact5/fact2
  kr2 = two*ks2
  kr3 = two*yi*yr/(fact2**3)
  kr = k*sqrt(kr1+kr2*kr3)

  ki1 = three*f*yi/fact2**2
  ki2 = kr2
  ki3 = ks3
  ki  = k**2/(two*kr)*(ki1+ki2*ki3)

  gv = point5
  gh = point5

  ssalb_v = min(ks/ki, 0.999_fp)
  ssalb_h = ssalb_v
  tau_v = two*ki*h
  tau_h = tau_v

end subroutine snow_optic


subroutine soil_diel(freq,t_soil,vmc,rhob,rhos,sand,clay,esm)

!----------------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    Soil_Diel   calculate the dielectric properties of soil
!
!   prgmmr: Fuzhong Weng and Banghua Yan                 org: nesdis              date: 2000-11-28
!
! abstract: compute the dilectric constant of the bare soil
!
! program history log:
!
! input argument list:
!
!      theta        -  local zenith angle (degree)
!      frequency    -  frequency (ghz)
!      t_soil       -  soil temperature
!      vmc          -  volumetric moisture content (demensionless)
!      rhob         -  bulk volume density of the soil (1.18-1.12)
!      rhos         -  density of the solids (2.65 g.cm^3 for
!                       solid soil material)
!      sand         -  sand fraction (sand + clay = 1.0)
!      clay         -  clay fraction
!
! output argument list:
!
!      esm          -  dielectric constant for bare soil
!
! important internal variables:
!
!      esof         -  the permittivity of free space
!      eswo         -  static dieletric constant
!      tauw         -  relaxation time of water
!      s            -  salinity
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------

  REAL(fp) :: f,tauw,freq,t_soil,vmc,rhob,rhos,sand,clay
  REAL(fp) :: alpha,beta,ess,rhoef,t,eswi,eswo
  REAL(fp) :: esof
  COMPLEX(fp) :: esm,esw,es1,es2

  alpha = 0.65_fp
  beta  = 1.09_fp - 0.11_fp*sand + 0.18_fp*clay
  ess = (1.01_fp + 0.44_fp*rhos)**2 - 0.062_fp
  rhoef = -1.645_fp + 1.939_fp*rhob - 0.020213_fp*sand + 0.01594_fp*clay
  t = t_soil - 273.0_fp
  f = freq*1.0e9_fp

  ! the permittivity at the high frequency limit
  eswi = 5.5_fp

  ! the permittivity of free space (esof)
  esof = 8.854e-12_fp

  ! static dieletric constant (eswo)
  eswo = 87.134_fp+(-1.949e-1_fp+(-1.276e-2_fp+2.491e-4_fp*t)*t)*t
  tauw = 1.1109e-10_fp+(-3.824e-12_fp+(6.938e-14_fp-5.096e-16_fp*t)*t)*t

  if (vmc > zero) then
     es1 = cmplx(eswi, -rhoef*(rhos-rhob)/(twopi*f*esof*rhos*vmc), fp)
  else
     es1 = cmplx(eswi, zero, fp)
  endif

  es2 = cmplx(eswo-eswi, zero, fp)/cmplx(one, f*tauw, fp)
  esw = es1 + es2
  esm = one + (ess**alpha - one)*rhob/rhos + vmc**beta*esw**alpha - vmc
  esm = esm**(one/alpha)

  if(aimag(esm) >= zero) esm = cmplx(REAL(esm,fp),-0.0001_fp, fp)

end subroutine soil_diel


subroutine snow_diel(frequency,ep_real,ep_imag,rad,frac,ep_eff)

!----------------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    Snow_Diel   compute dielectric constant of snow
!
!   prgmmr: Fuzhong Weng and Banghua Yan                 org: nesdis              date: 2000-11-28
!
! abstract: compute dielectric constant of snow
!
!
! program history log:
!
! input argument list:
!
!       frequency   -  frequency (ghz)
!       ep_real     -  real part of dielectric constant of particle
!       ep_imag     -  imaginary part of dielectric constant of particle
!       rad         -  particle radiu (mm)
!       frac        -  fraction volume of snow (0.0 - 1.0)
!
! output argument list:
!
!       ep_eff      -  dielectric constant of the dense medium
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!  Copyright (C) 2005 Fuzhong Weng and Banghua Yan
!
!  This program is free software; you can redistribute it and/or modify it under the terms of the GNU
!  General Public License as published by the Free Software Foundation; either version 2 of the License,
!  or (at your option) any later version.
!
!  This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
!  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
!  License for more details.
!
!  You should have received a copy of the GNU General Public License along with this program; if not, write
!  to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!----------------------------------------------------------------------------------

  REAL(fp) :: ep_imag,ep_real
  REAL(fp) :: frequency,rad,frac,k0,yr,yi
  COMPLEX(fp) :: y,ep_r,ep_i,ep_eff,fracy

  k0 = twopi/(300.0_fp/frequency)

  yr = (ep_real - one)/(ep_real + two)
  yi = ep_imag/(ep_real + two)

  y = cmplx(yr, yi, fp)
  fracy=frac*y

  ep_r = (one + two*fracy)/(one - fracy)
  ep_i = two*fracy*y*(k0*rad)**3*(one-frac)**4/((one-fracy)**2*(one+two*frac)**2)
  ep_eff = ep_r - cmplx(zero,one,fp)*ep_i

  if (aimag(ep_eff) >= zero) ep_eff = cmplx(REAL(ep_eff), -0.0001_fp, fp)

end subroutine snow_diel


subroutine canopy_diel(frequency,mg,esv,rhoveg)

!----------------------------------------------------------------------------------
!
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:   canopy_diel compute the dielectric constant of the vegetation canopy
!
!   prgmmr:  Fuzhong Weng and Banghua Yan                org: nesdis              date: 2000-11-28
!
! abstract: compute the dielectric constant of the vegetation canopy geomatrical optics approximation
!
!           for vegetation canopy work for horizontal leaves
!
! program history log:
!
! input argument list:
!
!      frequency    -  frequency (ghz)
!      mg           -  gravimetric water content
!
! output argument list:
!
!      esv          -  dielectric constant of leaves
!
! remarks:
!
! references:
!
!     ulaby and el-rayer, 1987: microwave dielectric spectrum of vegetation part ii,
!           dual-dispersion model, ieee trans geosci. remote sensing, 25, 550-557
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------

  REAL(fp) :: frequency,  mg, en, vf, vb
  REAL(fp) :: rhoveg, vmv
  COMPLEX(fp) :: esv, xx

  vmv = mg*rhoveg/( one - mg*(one-rhoveg) )
  en = 1.7_fp + (3.2_fp + 6.5_fp*vmv)*vmv

  vf = vmv*(0.82_fp*vmv + 0.166_fp)
  vb = 31.4_fp*vmv*vmv/( one + 59.5_fp*vmv*vmv)

  xx = cmplx(zero,one,fp)

  esv = en + vf*(4.9_fp + 75.0_fp/(one + xx*frequency/18.0_fp)-xx*(18.0_fp/frequency)) + &
        vb*(2.9_fp + 55.0_fp/(one + sqrt(xx*frequency/0.18_fp)))

  if (aimag(esv) >= zero) esv = cmplx(REAL(esv), -0.0001_fp, fp)

end subroutine canopy_diel


subroutine reflectance(em1, em2, theta_i, theta_t, rv, rh)

!----------------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    Reflectance compute the surface reflectivity
!
!   prgmmr:                  org: nesdis              date: 2000-11-28
!
! abstract: compute the surface reflectivety using fresnel equations
!    for a rough surface having a standard deviation of height of sigma
!
! program history log:
!
! input argument list:
!      theta_i      -  incident angle (degree)
!      theta_t      -  transmitted angle (degree)
!      em1          -  dielectric constant of the medium 1
!      em2          -  dielectric constant of the medium 2
!
! output argument list:
!
!      rv           -  reflectivity at vertical polarization
!      rh           -  reflectivity at horizontal polarization
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------

  REAL(fp) :: theta_i, theta_t
  REAL(fp) :: rh, rv,cos_i,cos_t
  COMPLEX(fp) :: em1, em2, m1, m2, angle_i, angle_t

  ! compute the refractive index ratio between medium 2 and 1
  ! using dielectric constant (n = SQRT(e))
  cos_i = cos(theta_i)
  cos_t = cos(theta_t)

  angle_i = cmplx(cos_i, zero, fp)
  angle_t = cmplx(cos_t, zero, fp)

  m1 = sqrt(em1)
  m2 = sqrt(em2)

  rv = (abs((m1*angle_t-m2*angle_i)/(m1*angle_t+m2*angle_i)))**2
  rh = (abs((m1*angle_i-m2*angle_t)/(m1*angle_i+m2*angle_t)))**2

end subroutine reflectance


subroutine transmittance(em1,em2,theta_i,theta_t,tv,th)

!----------------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    Transmittance    calculate Transmittance
!
!   prgmmr:  Banghua Yan and Fuzhong Weng               org: nesdis              date: 2000-11-28
!
! abstract: compute Transmittance
!
! program history log:
!
! input argument list:
!
!      theta        -  local zenith angle (degree)
!      theta_i      -  incident angle (degree)
!      theta_t      -  transmitted angle (degree)
!      em1          -  dielectric constant of the medium 1
!      em2          -  dielectric constant of the medium 2
!
! output argument list:
!
!      tv           -  transmisivity at vertical polarization
!      th           -  transmisivity at horizontal polarization
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------

  REAL(fp) :: theta_i, theta_t
  REAL(fp) :: th, tv, rr, cos_i,cos_t
  COMPLEX(fp) :: em1, em2, m1, m2, angle_i, angle_t

  ! compute the refractive index ratio between medium 2 and 1
  ! using dielectric constant (n = SQRT(e))
  cos_i = cos(theta_i)
  cos_t = cos(theta_t)

  angle_i = cmplx(cos_i, zero, fp)
  angle_t = cmplx(cos_t, zero, fp)

  m1 = sqrt(em1)
  m2 = sqrt(em2)

  rr = abs(m2/m1)*cos_t/cos_i
  tv = rr*(abs(two*m1*angle_i/(m1*angle_t + m2*angle_i)))**2
  th = rr*(abs(two*m1*angle_i/(m1*angle_i + m2*angle_t)))**2

end subroutine transmittance


subroutine roughness_reflectance(frequency,sigma,rv,rh)

!-------------------------------------------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    rought_reflectance    calculate surface relectivity
!
!   prgmmr: Banghua Yan and Fuzhong Weng                 org: nesdis              date: 2000-11-28
!
! abstract: compute the surface reflectivety for a rough surface having a standard devoation of height of sigma
!
!
! program history log:
!
! input argument list:
!
!      frequency    -  frequency (ghz)
!
!      theta        -  local zenith angle (degree) (currently, not used here)
!
!      sigma        -  standard deviation of rough surface height
!
!                      smooth surface:0.38, medium: 1.10, rough:2.15 cm
!
!    internal variables
!
!
! output argument list:
!
!      rv            -  reflectivity at vertical polarization
!      rh            -  reflectivity at horizontal polarization
!
!
!   important internal variables:
!
!      k0           -  a propagation constant or wavenumber in a free space
!
! remarks:
!
! references:
!
!   wang, j. and b. j. choudhury, 1992: passive microwave radiation from soil: examples...
!    passive microwave remote sensing of .. ed. b. j. choudhury, etal vsp.
!    also wang and choudhury (1982)
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!-------------------------------------------------------------------------------------------------------------

  REAL(fp) :: frequency
  REAL(fp) :: q, rh, rv, rh_s, rv_s, sigma

  rh_s = 0.3_fp*rh
  rv_s = 0.3_fp*rv
  q = 0.35_fp*(one - exp(-0.60_fp*frequency*sigma**two))
  rh = rh_s + q*(rv_s-rh_s)
  rv = rv_s + q*(rh_s-rv_s)

end subroutine roughness_reflectance


subroutine two_stream_solution(mu,gv,gh,ssalb_h,ssalb_v,tau_h,tau_v, &
      r21_h,r21_v,r23_h,r23_v,t21_v,t21_h,esv,esh,frequency,t_soil,t_skin)

!-------------------------------------------------------------------------------------------------------------
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    two_stream_solution
!
!   prgmmr: Banghua Yan and Fuzhong Weng                 org: nesdis              date: 2000-11-28
!
! abstract: two stream solution
!         Updated with the more accurate formula of total upwelling radiance emanating from the surface.
!
! REFERENCES:
!       Weng, F., B. Yan, and N. Grody, 2001: "A microwave land emissivity model", J. Geophys. Res., 106,
!                                             20, 115-20, 123
!   version: beta
!
! program history log:
!
! input argument list:
!
!      b            -  scattering layer temperature (k)         (gdas)   (not used here)
!      mu           -  cos(theta)
!      gv           -  asymmetry factor for v pol
!      gh           -  asymmetry factor for h pol
!      ssalb_v      -  single scattering albedo at v. polarization
!      ssalb_h      -  single scattering albedo at h. polarization
!      tau_v        -  optical depth at v. polarization
!      tau_h        -  optical depth at h. polarization
!      r12_v        -  reflectivity at vertical polarization   (not used here)
!      r12_h        -  reflectivity at horizontal polarization (not used here)
!      r21_v        -  reflectivity at vertical polarization
!      r21_h        -  reflectivity at horizontal polarization
!      r23_v        -  reflectivity at vertical polarization
!      r23_h        -  reflectivity at horizontal polarization
!      t21_v        -  transmisivity at vertical polarization
!      t21_h        -  transmisivity at horizontal polarization
!      t12_v        -  transmisivity at vertical polarization   (not used here)
!      t12_h        -  transmisivity at horizontal polarization (not used here)
!      Frequency    -  frequency
!      t_soil       -  soil temperature
!      t_skin       -  land surface temperature
!
! output argument list:
!
!       esh         -  emissivity for horizontal polarization
!       esv         -  emissivity for vertical polarization
!
! Local variables:
!       gsect0, gsect1_h, gsect1_v, gsect2_h, gsect2_v
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!-------------------------------------------------------------------------------------------------------------

  REAL(fp) :: mu, gv, gh, ssalb_h, ssalb_v, tau_h,tau_v,                 &
              r21_h, r21_v, r23_h, r23_v, t21_v, t21_h, esv, esh
  REAL(fp) :: alfa_v, alfa_h, kk_h, kk_v, gamma_h, gamma_v, beta_v, beta_h
  REAL(fp) :: fact1,fact2
  REAL(fp) :: frequency, t_soil, t_skin
  REAL(fp) :: gsect0, gsect1_h, gsect1_v, gsect2_h, gsect2_v

  alfa_h  = sqrt((one - ssalb_h)/(one - gh*ssalb_h))
  kk_h    = sqrt((one - ssalb_h)*(one -  gh*ssalb_h))/mu
  beta_h  = (one - alfa_h)/(one + alfa_h)
  gamma_h = (beta_h -r23_h)/(one-beta_h*r23_h)

  alfa_v  = sqrt((one-ssalb_v)/(one - gv*ssalb_v))
  kk_v    = sqrt((one-ssalb_v)*(one - gv*ssalb_v))/mu
  beta_v  = (one - alfa_v)/(one + alfa_v)
  gamma_v = (beta_v -r23_v)/(one-beta_v*r23_v)

  fact1=gamma_h*exp(-two*kk_h*tau_h)
  fact2=gamma_v*exp(-two*kk_v*tau_v)

  gsect0  =(exp(c_2*frequency/t_skin) -one)/(exp(c_2*frequency/t_soil) -one)

  gsect1_h=(one-r23_h)*(gsect0-one)
  gsect2_h=((one-beta_h*beta_h)/(one-beta_h*r23_h))*exp(-kk_h*tau_h)

  gsect1_v=(one-r23_v)*(gsect0-one)
  gsect2_v=((one-beta_v*beta_v)/(one-beta_v*r23_v))*exp(-kk_h*tau_v)

  esh  = t21_h*((one - beta_h)*(one + fact1)+gsect1_h*gsect2_h) /(one-beta_h*r21_h-(beta_h-r21_h)*fact1)
  esv  = t21_v*((one - beta_v)*(one + fact2)+gsect1_v*gsect2_v) /(one-beta_v*r21_v-(beta_v-r21_v)*fact2)

  if (esh < emissh_default) esh = emissh_default
  if (esv < emissv_default) esv = emissv_default

  if (esh > one) esh = one
  if (esv > one) esv = one

end subroutine two_stream_solution

END MODULE nesdis_landem_module