!
! NESDIS_AMSU_SnowEM_Module
!
! Module containing the AMSU microwave snow emissivity model
!
! References:
!       Yan,B., F.Weng, and K.Okamoto, 2004, A microwave snow emissivity model,
!         8th Specialist Meeting on Microwave Radiometry and Remote Sensing Applications,
!         24-27 February, 2004, Rome, Italy.
!
!
! CREATION HISTORY:
!       Written by:     Banghua Yan, 03-Jun-2005, banghua.yan@noaa.gov
!                       Fuzhong Weng, fuzhong.weng@noaa.gov
!

MODULE NESDIS_AMSU_SnowEM_Module


  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE Type_Kinds, ONLY: fp, Double
  USE NESDIS_LandEM_Module
  USE NESDIS_SnowEM_Parameters
  ! Disable implicit typing
  IMPLICIT NONE


  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: NESDIS_AMSU_SNOWEM


  ! -----------------
  ! Module parameters
  ! -----------------
  ! Version Id for the module
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  '$Id: NESDIS_AMSU_SnowEM_Module.f90 99117 2017-11-27 18:37:14Z tong.zhu@noaa.gov $'


CONTAINS


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

!-------------------------------------------------------------------------------------------------------------
!
! NAME:
!       NESDIS_AMSU_SNOWEM
!
! PURPOSE:
!       Subroutine to simulate microwave emissivity over snow conditions from AMSU measurements at window
!       channels.
!
! REFERENCES:
!       Yan, B., F. Weng and K.Okamoto,2004: "A microwave snow emissivity model, 8th Specialist Meeting on
!       Microwave Radiometry and Remote Sension Applications,24-27 February, 2004, Rome, Italy.
!
! CATEGORY:
!       CRTM : Surface : MW SNOWEM
!
! LANGUAGE:
!       Fortran-95
!
! CALLING SEQUENCE:
!       CALL NESDIS_AMSU_SNOWEM
!
! INPUT ARGUMENTS:
!
!         Frequency                Frequency User defines
!                                  This is the "I" dimension
!                                  UNITS:      GHz
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Scalar
!
!
!         Satellite_Angle          The local zenith angle in degree for AMSU measurements.
!                                  ** NOTE: THIS IS A MANDATORY MEMBER **
!                                  **       OF THIS STRUCTURE          **
!                                  UNITS:      Degrees
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Rank-1, (I)
!
!         User_Angle               The local angle value in degree user defines.
!                                  ** NOTE: THIS IS A MANDATORY MEMBER **
!                                  **       OF THIS STRUCTURE          **
!                                  UNITS:      Degrees
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Rank-1, (I)
!
!
!         Tba                      BRIGHTNESS TEMPERATURES AT FOUR AMSU-A WINDOW CHANNELS
!                                  UNITS:      Kelvin, K
!                                  TYPE:       REAL( fp )
!                                  DIMENSION   4*1 SCALAR
!
!                        WHICH ARE
!                                  tba[1] = TB at 23.8 GHz
!                                  tba[2] = TB at: 31.4 GHz
!                                  tba[3] = TB at 50.3 GHz
!                                  tba[4] = TB at 89 GHz
!
!         Tbb                      BRIGHTNESS TEMPERATURES AT TWO AMSU-B WINDOW CHANNELS
!                                  UNITS:      Kelvin, K
!                                  TYPE:       REAL( fp )
!                                  DIMENSION   2*1 SCALAR
!
!                         WHICH ARE
!
!                                  tbb[1] = TB at 89 GHz
!                                  tbb[2] = TB at 150 GHz
!
!
!         Ts = 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
!
! **** IMPORTANT NOTES ****
!
!        When one variable among  Tba[], Tbb[] and Ts are not available, set -999.0
!
!
!
!
! OUTPUT ARGUMENTS:
!
!         Emissivity_H:            The surface emissivity at a horizontal polarization.
!                                  ** NOTE: THIS IS A MANDATORY MEMBER **
!                                  **       OF THIS STRUCTURE          **
!                                  UNITS:      N/A
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Scalar
!
!         Emissivity_V:            The surface emissivity at a vertical polarization.
!                                  ** NOTE: THIS IS A MANDATORY MEMBER **
!                                  **       OF THIS STRUCTURE          **
!                                  UNITS:      N/A
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Scalar
!
!
! IINTERNAL ARGUMENTS:
!
!       input_type (specific option index = 1 ~ 7):
!
!         input_type = 1 :  AMSU-A & B window channels of Tb and Ts are available (call AMSU_ABTs)
!
!         input_type = 2 :  AMSU-A window channels of Tb and Ts are available     (call AMSU_ATs)
!
!         input_type = 3 :  AMSU-A & B window channels of Tb are available        (call AMSU_AB)
!
!         input_type = 4 :  AMSU-A window channels of Tb are available            (call AMSU_amsua)
!
!         input_type = 5 :  AMSU-B window channels of Tb and Ts are available     (call AMSU_BTs)
!
!         input_type = 6 :  AMSU-B window channels of Tb are available            (call AMSU_amsub)
!
!         input_type = 7 :  snow depth and Ts are available                       (call ALandEM_Snow)
!
!
! OPTIONAL OUTPUT ARGUMENTS:
!
!       snow_type  -  snow type (not output here)
!                     1 : Wet Snow
!                     2 : Grass_after_Snow
!                     3 : RS_Snow (A)
!                     4 : Powder Snow
!                     5 : RS_Snow (B)
!                     6 : RS_Snow (C)
!                     7 : RS_Snow (D)
!                     8 : Thin Crust Snow
!                     9 : RS_Snow (E)
!                     10: Bottom Crust Snow (A)
!                     11: Shallow Snow
!                     12: Deep Snow
!                     13: Crust Snow
!                     14: Medium Snow
!                     15: Bottom Crust Snow (B)
!                     16: Thick Crust Snow
!                    999: AMSU measurements are not available or over non-snow conditions
!
! CALLS:
!       AMSU_ABTs          : Subroutine to calculate the microwave snow emissivity from AMSU-A/B TB and Ts
!
!       AMSU_ATs           : Subroutine to calculate the microwave snow emissivity from AMSU-A TB and Ts
!
!       AMSU_AB            : Subroutine to calculate the microwave snow emissivity from AMSU-A/B TB
!
!       AMSU_amsua         : Subroutine to calculate the microwave snow emissivity from AMSU-A TB
!
!       AMSU_BTs           : Subroutine to calculate the microwave snow emissivity from AMSU-B TB and Ts
!
!       AMSU_amsub         : Subroutine to calculate the microwave snow emissivity from AMSU-B TB
!
!       AMSU_ALandEM_Snow  : Subroutine to calculate the microwave snow emissivity from Ts and Snow Depth
!
!       em_initialization  : Subroutine to initialization snow emissivity
!
!
! SIDE EFFECTS:
!       None.
!
! RESTRICTIONS:
!       None.
!
! COMMENTS:
!       Note the INTENT on the output SensorData argument is IN OUT rather than
!       just OUT. This is necessary because the argument may be defined upon
!       input. To prevent memory leaks, the IN OUT INTENT is a must.
!
! CREATION HISTORY:
!       Written by:     Banghua Yan, QSS Group Inc., Banghua.Yan@noaa.gov (03-June-2005)
!
!
!       and             Fuzhong Weng, NOAA/NESDIS/ORA, Fuzhong.Weng@noaa.gov
!
!  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.
!
!------------------------------------------------------------------------------------------------------------

subroutine  NESDIS_AMSU_SNOWEM(Satellite_Angle,                             &  ! INPUT
                               User_Angle,                                  &  ! INPUT
                               frequency,                                   &  ! INPUT
                               Snow_Depth,                                  &  ! INPUT
                               Ts,                                          &  ! INPUT
                               tba,                                         &  ! INPUT
                               tbb,                                         &  ! INPUT
                               Emissivity_H,                                &  ! OUPUT
                               Emissivity_V)                                   ! OUTPUT

  integer,PARAMETER ::  AMSU_ABTs_ALG    = 1

  integer,PARAMETER ::  AMSU_ATs_ALG     = 2

  integer,PARAMETER ::  AMSU_AB_ALG      = 3

  integer,PARAMETER ::  AMSU_amsua_ALG   = 4

  integer,PARAMETER ::  AMSU_BTs_ALG     = 5

  integer,PARAMETER ::  AMSU_amsub_ALG   = 6

  integer,PARAMETER ::  AMSU_ALandEM_Snow_ALG  = 7

  integer, parameter  :: nch = 10, nwcha = 4, nwchb = 2, nwch = 5,nalg = 7

  integer :: snow_type,input_type,i,np,k

  real(fp)    :: Satellite_Angle,User_Angle,frequency,Ts,Snow_Depth

  real(fp)    :: em_vector(2),esh1,esv1,esh2,esv2,desh,desv,dem

  real(fp)    :: tb(nwch),tba(nwcha),tbb(nwchb)

  real(fp), intent(out) :: Emissivity_H,Emissivity_V

  logical :: INDATA(nalg)



! Initialization

  call em_initialization(frequency,em_vector)

  snow_type  = INVALID_SNOW_TYPE

  input_type = INVALID_SNOW_TYPE

  do k = 1, nalg

     INDATA(k) = .TRUE.

  end do


! Read AMSU & Ts data and set available option

! Get five AMSU-A/B window measurements

  tb(1) = tba(1); tb(2) = tba(2); tb(3) = tba(3); tb(4) = tba(4); tb(5) = tbb(2)

! Check available data

  if((Ts <= 150.0_fp) .or. (Ts >= 290.0_fp) ) then

     INDATA(1:2) = .false.;   INDATA(5)  = .false.;  INDATA(7) = .false.

  end if

  do i=1,nwcha

     if((tba(i) <= 100.0_fp) .or. (tba(i) >= 290.0_fp) ) then

        INDATA(1:4)   = .false.

        exit

     end if

  end do

  do i=1,nwchb

     if((tbb(i) <= 100.0_fp) .or. (tbb(i) >= 290.0_fp) ) then

        INDATA(1)  = .false.;  INDATA(3) = .false.;  INDATA(5:6)  = .false.

        exit

     end if

  end do

  if((Snow_Depth < 0.0_fp) .or. (Snow_Depth >= 3000.0_fp)) INDATA(7) = .false.

  if((frequency >= 80._fp) .and. (INDATA(5))) then

     INDATA(2:3) = .false.

  end if

! Check input type and call a specific Option/subroutine

  do np = 1, nalg

     if (INDATA(np)) then

        input_type = np

        exit

     end if

  end do

! GET EMISSIVITY AT SATELLITE'S ZENITH ANGLE

  GET_option: SELECT CASE (input_type)

  CASE (AMSU_ABTs_ALG)

     call AMSU_ABTs(frequency,tb,Ts,snow_type,em_vector)

  CASE (AMSU_ATs_ALG)

     call AMSU_ATs(frequency,tba,Ts,snow_type,em_vector)

  CASE (AMSU_AB_ALG)

     call AMSU_AB(frequency,tb,snow_type,em_vector)

  CASE (AMSU_amsua_ALG)

     call AMSU_amsua(frequency,tba,snow_type,em_vector)

  CASE(AMSU_BTs_ALG)

     call AMSU_BTs(frequency,tbb,Ts,snow_type,em_vector)

  CASE(AMSU_amsub_ALG)

     call AMSU_amsub(frequency,tbb,snow_type,em_vector)

  CASE(AMSU_ALandEM_Snow_ALG)

     call AMSU_ALandEM_Snow(Satellite_Angle,frequency,Snow_Depth,Ts,snow_type,em_vector)

  END SELECT GET_option

! Get the emissivity angle dependence

  call NESDIS_LandEM(Satellite_Angle,frequency,0.0_fp,0.0_fp,Ts,Ts,0.0_fp,9,13,2.0_fp,esh1,esv1)

  call NESDIS_LandEM(User_Angle,frequency,0.0_fp,0.0_fp,Ts,Ts,0.0_fp,9,13,2.0_fp,esh2,esv2)

  desh = esh1 - esh2

  desv = esv1 - esv2

  dem = ( desh + desv ) * 0.5_fp
! Emissivity at User's Angle

  Emissivity_H = em_vector(1) - dem;  Emissivity_V = em_vector(2)- dem

  if (Emissivity_H > one)         Emissivity_H = one

  if (Emissivity_V > one)         Emissivity_V = one

  if (Emissivity_H < 0.3_fp) Emissivity_H = 0.3_fp

  if (Emissivity_V < 0.3_fp) Emissivity_V = 0.3_fp


  return

end subroutine NESDIS_AMSU_SNOWEM


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


subroutine em_initialization(frequency,em_vector)

!----------------------------------------------------------------------------------------------------------!
!$$$  subprogram documentation block
!
! subprogram:   AMSU-A/B snow emissivity initialization
!
!   prgmmr:  Banghua Yan                org: nesdis              date: 2003-08-18
!
! abstract: AMSU-A/B snow emissivity initialization
!
! program history log:
!
! input argument list:
!
!      frequency   - frequency in GHz
!
! output argument list:
!
!     em_vector[1] and [2]  -  initial emissivity at two polarizations.
!
! important internal variables:
!
!      freq[1~10]  - ten frequencies for sixteen snow types of emissivity
!      em[1~16,*]  - sixteen snow emissivity spectra
!      snow_type   - snow type
!                    where it is initialized to as the type 4,i.e, Powder Snow
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------------------------------!

  integer ::  nch,ncand
  Parameter(nch = 10,ncand=16)
  real(fp)    :: frequency,em_vector(*),freq(nch)
  real(fp)    :: em(ncand,nch)
  real(fp)    :: kratio, bconst,emissivity
  integer :: ich

! Sixteen candidate snow emissivity spectra

  em(1, 1: N_FREQUENCY) = WET_SNOW_EMISS(1:N_FREQUENCY)
  em(2, 1: N_FREQUENCY) = GRASS_AFTER_SNOW_EMISS(1:N_FREQUENCY)
  em(3, 1: N_FREQUENCY) = RS_SNOW_A_EMISS(1:N_FREQUENCY)
  em(4, 1: N_FREQUENCY) = POWDER_SNOW_EMISS(1:N_FREQUENCY)
  em(5, 1: N_FREQUENCY) = RS_SNOW_B_EMISS(1:N_FREQUENCY)
  em(6, 1: N_FREQUENCY) = RS_SNOW_C_EMISS(1:N_FREQUENCY)
  em(7, 1: N_FREQUENCY) = RS_SNOW_D_EMISS(1:N_FREQUENCY)
  em(8, 1: N_FREQUENCY) = THIN_CRUST_SNOW_EMISS(1:N_FREQUENCY)
  em(9, 1: N_FREQUENCY) = RS_SNOW_E_EMISS(1:N_FREQUENCY)
  em(10, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_A_EMISS(1:N_FREQUENCY)
  em(11, 1: N_FREQUENCY) = SHALLOW_SNOW_EMISS(1:N_FREQUENCY)
  em(12, 1: N_FREQUENCY) = DEEP_SNOW_EMISS(1:N_FREQUENCY)
  em(13, 1: N_FREQUENCY) = CRUST_SNOW_EMISS(1:N_FREQUENCY)
  em(14, 1: N_FREQUENCY) = MEDIUM_SNOW_EMISS(1:N_FREQUENCY)
  em(15, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_B_EMISS(1:N_FREQUENCY)
  em(16, 1: N_FREQUENCY) = THICK_CRUST_SNOW_EMISS(1:N_FREQUENCY)


  freq = FREQUENCY_DEFAULT

! Initialization for emissivity at certain frequency
!    In case of no any inputs available for various options
!    A constant snow type & snow emissivity spectrum is assumed
!                    (e.g., powder) snow_type = 4

! Specify snow emissivity at required frequency
  do ich = 2, nch
     if(frequency <  freq(1))   exit
     if(frequency >= freq(nch)) exit
     if(frequency <  freq(ich)) then
        emissivity = em(4,ich-1) + (em(4,ich) - em(4,ich-1))     &
             *(frequency - freq(ich-1))/(freq(ich) - freq(ich-1))
        exit
     end if
  end do

! Extrapolate to lower frequencies than 4.9GHz
  if (frequency <= freq(1)) then
     kratio = (em(4,2) - em(4,1))/(freq(2) - freq(1))
     bconst = em(4,1) - kratio*freq(1)
     emissivity =  kratio*frequency + bconst
     if(emissivity >  one)         emissivity = one
     if(emissivity <= 0.8_fp) emissivity = 0.8_fp
  end if


! Assume emissivity = constant at frequencies >= 150 GHz
  if (frequency >= freq(nch)) emissivity = em(4,nch)
  em_vector(1) = emissivity
  em_vector(2) = emissivity

  return
end subroutine em_initialization



subroutine  em_interpolate(frequency,discriminator,emissivity,snow_type)

!----------------------------------------------------------------------------------------------------------!
!$$$  subprogram documentation block
!
! subprogram:  determine snow_type and calculate emissivity
!
!   prgmmr:Banghua Yan                 org: nesdis              date: 2003-08-18
!
! abstract: 1. Find one snow emissivity spectrum to mimic the emission
!              property of the realistic snow condition using a set of
!              discrminators
!           2. Interpolate/extrapolate emissivity at a required frequency
!
! program history log:
!
! input argument list:
!
!      frequency        - frequency in GHz
!      discriminators   - emissivity discriminators at five AMSU-A & B window
!                         channels
!            discriminator[1]   :  emissivity discriminator at 23.8 GHz
!            discriminator[2]   :  emissivity discriminator at 31.4 GHz
!            discriminator[3]   :  emissivity discriminator at 50.3 GHz
!            discriminator[4]   :  emissivity discriminator at 89   GHz
!            discriminator[5]   :  emissivity discriminator at 150  GHz
!
!       Note: discriminator(1) and discriminator(3) are missing value in
!            'AMSU-B & Ts','AMUS-B' and 'MODL' options., which are defined to as -999.9,
! output argument list:
!
!     emissivity  -  weighted emissivity from both V- and H- Pols.
!     snow_type             - snow type
!
! important internal variables:
!
!     freq[1 ~ 10]  -  ten frequencies for sixteen snow types of emissivity
!     em[1~16,*]    -  sixteen snow emissivity spectra
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------------------------------!

  integer,parameter:: ncand = 16,nch =10
  integer:: ich,ichmin,ichmax,i,k,snow_type
  real(fp)   :: dem,demmin0
  real(fp)   :: em(ncand,nch)
  real(fp)   :: frequency,freq(nch),emissivity,discriminator(*),emis(nch)
  real(fp)   :: cor_factor,adjust_check,kratio, bconst

! Sixteen candidate snow emissivity spectra

  em(1, 1: N_FREQUENCY) = WET_SNOW_EMISS(1:N_FREQUENCY)
  em(2, 1: N_FREQUENCY) = GRASS_AFTER_SNOW_EMISS(1:N_FREQUENCY)
  em(3, 1: N_FREQUENCY) = RS_SNOW_A_EMISS(1:N_FREQUENCY)
  em(4, 1: N_FREQUENCY) = POWDER_SNOW_EMISS(1:N_FREQUENCY)
  em(5, 1: N_FREQUENCY) = RS_SNOW_B_EMISS(1:N_FREQUENCY)
  em(6, 1: N_FREQUENCY) = RS_SNOW_C_EMISS(1:N_FREQUENCY)
  em(7, 1: N_FREQUENCY) = RS_SNOW_D_EMISS(1:N_FREQUENCY)
  em(8, 1: N_FREQUENCY) = THIN_CRUST_SNOW_EMISS(1:N_FREQUENCY)
  em(9, 1: N_FREQUENCY) = RS_SNOW_E_EMISS(1:N_FREQUENCY)
  em(10, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_A_EMISS(1:N_FREQUENCY)
  em(11, 1: N_FREQUENCY) = SHALLOW_SNOW_EMISS(1:N_FREQUENCY)
  em(12, 1: N_FREQUENCY) = DEEP_SNOW_EMISS(1:N_FREQUENCY)
  em(13, 1: N_FREQUENCY) = CRUST_SNOW_EMISS(1:N_FREQUENCY)
  em(14, 1: N_FREQUENCY) = MEDIUM_SNOW_EMISS(1:N_FREQUENCY)
  em(15, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_B_EMISS(1:N_FREQUENCY)
  em(16, 1: N_FREQUENCY) = THICK_CRUST_SNOW_EMISS(1:N_FREQUENCY)


  freq = FREQUENCY_DEFAULT



! Adjust unreasonable discriminator
  if (discriminator(4) > discriminator(2))    &
       discriminator(4) = discriminator(2) +(discriminator(5) - discriminator(2))*  &
       (150.0_fp - 89.0_fp)/(150.0_fp - 31.4_fp)
  if ( (discriminator(3) /= -999.9_fp) .and.       &
       ( ((discriminator(3)-0.01_fp) > discriminator(2)) .or.     &
       ((discriminator(3)-0.01_fp) < discriminator(4)))    )    &
       discriminator(3) = discriminator(2) +  &
       (discriminator(4) - discriminator(2))*(89.0_fp - 50.3_fp) &
       / (89.0_fp - 31.4_fp)

! Find a snow emissivity spectrum
  if(snow_type .eq. -999) then
     demmin0 = 10.0_fp
     do k = 1, ncand
        dem = zero
        ichmin = 1
        ichmax = 3
        if(discriminator(1) == -999.9_fp) then
           ichmin = 2
           ichmax = 2
        end if
        do ich = ichmin,ichmax
           dem = dem + abs(discriminator(ich) - em(k,ich+4))
        end do
        do ich = 4,5
           dem = dem + abs(discriminator(ich) - em(k,ich+5))
        end do
        if (dem < demmin0) then
           demmin0 = dem
           snow_type = k
        end if
     end do
  end if

! Shift snow emissivity according to discriminator at 31.4 GHz
  cor_factor = discriminator(2) - em(snow_type,6)
  do ich = 1, nch
     emis(ich) = em(snow_type,ich) + cor_factor
     if(emis(ich) .gt. one)         emis(ich) = one
     if(emis(ich) .lt. 0.3_fp) emis(ich) = 0.3_fp
  end do

! Emisivity data quality control
  adjust_check = zero
  do ich = 5, 9
     if (ich .le. 7) then
        if (discriminator(ich - 4) .ne. -999.9_fp) &
             adjust_check = adjust_check + abs(emis(ich) - discriminator(ich - 4))
     else
        if (discriminator(ich - 4) .ne. -999.9_fp)  &
             adjust_check = adjust_check + abs(emis(ich+1) - discriminator(ich - 4))
     end if
  end do

  if (adjust_check >= 0.04_fp) then
     if (discriminator(1) /= -999.9_fp) then
        if (discriminator(1) < emis(4)) then
           emis(5) = emis(4) + &
                (31.4_fp - 23.8_fp) * &
                (discriminator(2) - emis(4))/(31.4_fp - 18.7_fp)
        else
           emis(5) = discriminator(1)
        end if
     end if
     emis(6) = discriminator(2)
     if (discriminator(3) /= -999.9_fp) then
        emis(7) = discriminator(3)
     else
!       In case of missing the emissivity discriminator at 50.3 GHz
        emis(7) = emis(6) + (89.0_fp - 50.3_fp) * &
             (discriminator(4) - emis(6))/(89.0_fp - 31.4_fp)
     end if
     emis(8) = emis(7)
     emis(9) = discriminator(4)
     emis(10) = discriminator(5)
  end if

! Estimate snow emissivity at a required frequency
  do i = 2, nch
     if(frequency <  freq(1))   exit
     if(frequency >= freq(nch)) exit
     if(frequency <  freq(i)) then
        emissivity = emis(i-1) + (emis(i) - emis(i-1))*(frequency - freq(i-1))  &
             /(freq(i) - freq(i-1))
        exit
     end if
  end do

! Extrapolate to lower frequencies than 4.9GHz
  if (frequency <= freq(1)) then
     kratio = (emis(2) - emis(1))/(freq(2) - freq(1))
     bconst = emis(1) - kratio*freq(1)
     emissivity =  kratio*frequency + bconst
     if(emissivity > one)          emissivity = one
     if(emissivity <= 0.8_fp) emissivity = 0.8_fp
  end if

! Assume emissivity = constant at frequencies >= 150 GHz
  if (frequency >= freq(nch)) emissivity = emis(nch)

  return
end subroutine em_interpolate


subroutine AMSU_ABTs(frequency,tb,ts,snow_type,em_vector)

!----------------------------------------------------------------------------------------------------------!
!$$$  subprogram documentation block
!
! subprogram:
!
!   prgmmr:Banghua Yan                  org: nesdis              date: 2003-08-18
!
! abstract:
!         Calculate the emissivity discriminators and interpolate/extrapolate
!  emissivity at a required frequency with respect to secenery ABTs
!
! program history log:
!
! input argument list:
!
!     frequency        -  frequency in GHz
!     theta            -  local zenith angle (currently, not used here)
!     tb[1] ~ tb[5]    -  brightness temperature at five AMSU window channels:
!                              tb[1] : 23.8 GHz
!                              tb[2] : 31.4 GHz
!                              tb[3] : 50.3 GHz
!                              tb[4] : 89.0 GHz
!                              tb[5] : 150  GHz
!
! output argument list:
!
!      em_vector[1] and [2]  -  emissivity at two polarizations.
!                              set esv = esh here and will be updated
!      snow_type        -  snow type
!
! important internal variables:
!
!     nind           -  number of threshold in decision trees
!                          to identify each snow type  ( = 6)
!     em(1~16,*)     -  sixteen snow emissivity spectra
!     DI_coe         -  coefficients to generate six discriminators to describe
!                       the overall emissivity variability within a wider frequency range
!     threshold      -  thresholds in decision trees to identify snow types
!     index_in       -  six indices to discriminate snow type
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------------------------------!

  integer,parameter:: ncand = 16,nch =10,nthresh=38
  integer,parameter:: nind=6,ncoe=8,nLIcoe=6,nHIcoe=12
  integer:: i,j,k,num,npass,snow_type,md0,md1,nmodel(ncand-1)
  real(fp)   :: frequency,tb150,LI,HI,DS1,DS2,DS3
  real(fp)   :: em(ncand,nch), em_vector(*)
  real(fp)   :: tb(*),freq(nch),DTB(nind-1),DI(nind-1),       &
       DI_coe(nind-1,0:ncoe-1),threshold(nthresh,nind),       &
       index_in(nind),threshold0(nind)
  real(fp)   :: LI_coe(0:nLIcoe-1),HI_coe(0:nHIcoe-1)
  real(fp)   :: ts,emissivity
  real(fp)   :: discriminator(5)
  logical:: pick_status,tindex(nind)
  save      threshold,DI_coe,LI_coe, HI_coe,nmodel

  data  nmodel/5,10,13,16,18,24,30,31,32,33,34,35,36,37,38/

! Fitting coefficients for five discriminators
  DI_coe(1,0:ncoe-1) = (/ &
       3.285557e-002_fp,  2.677179e-005_fp,  &
       4.553101e-003_fp,  5.639352e-005_fp,  &
       -1.825188e-004_fp,  1.636145e-004_fp,  &
       1.680881e-005_fp, -1.708405e-004_fp/)
  DI_coe(2,0:ncoe-1) = (/ &
       -4.275539e-002_fp, -2.541453e-005_fp,  &
       4.154796e-004_fp,  1.703443e-004_fp,  &
       4.350142e-003_fp,  2.452873e-004_fp,  &
       -4.748506e-003_fp,  2.293836e-004_fp/)
  DI_coe(3,0:ncoe-1) = (/ &
       -1.870173e-001_fp, -1.061678e-004_fp,  &
      2.364055e-004_fp, -2.834876e-005_fp,  &
      4.899651e-003_fp, -3.418847e-004_fp,  &
      -2.312224e-004_fp,  9.498600e-004_fp/)
  DI_coe(4,0:ncoe-1) = (/ &
       -2.076519e-001_fp,  8.475901e-004_fp,  &
       -2.072679e-003_fp, -2.064717e-003_fp,  &
       2.600452e-003_fp,  2.503923e-003_fp,  &
       5.179711e-004_fp,  4.667157e-005_fp/)
  DI_coe(5,0:ncoe-1) = (/ &
       -1.442609e-001_fp, -8.075003e-005_fp,  &
       -1.790933e-004_fp, -1.986887e-004_fp,  &
       5.495115e-004_fp, -5.871732e-004_fp,  &
       4.517280e-003_fp,  7.204695e-004_fp/)

! Fitting coefficients for emissivity index at 31.4 GHz
  LI_coe = (/ &
       7.963632e-001_fp,  7.215580e-003_fp,  &
       -2.015921e-005_fp, -1.508286e-003_fp,  &
       1.731405e-005_fp, -4.105358e-003_fp/)

! Fitting coefficients for emissivity index at 150 GHz
  HI_coe = (/ &
       1.012160e+000_fp,  6.100397e-003_fp, &
       -1.774347e-005_fp, -4.028211e-003_fp, &
       1.224470e-005_fp,  2.345612e-003_fp, &
       -5.376814e-006_fp, -2.795332e-003_fp, &
       8.072756e-006_fp,  3.529615e-003_fp, &
       1.955293e-006_fp, -4.942230e-003_fp/)

! Six thresholds for sixteen candidate snow types
! Note: some snow type contains several possible
!      selections for six thresholds

!1 Wet Snow
  threshold(1,1:6) = (/0.88_fp,0.86_fp,-999.9_fp,&
       0.01_fp,0.01_fp,200._fp/)
  threshold(2,1:6) = (/0.88_fp,0.85_fp,-999.9_fp,&
       0.06_fp,0.10_fp,200._fp/)
  threshold(3,1:6) = (/0.88_fp,0.83_fp,-0.02_fp,&
       0.12_fp,0.16_fp,204._fp/)
  threshold(4,1:6) = (/0.90_fp,0.89_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)
  threshold(5,1:6) = (/0.92_fp,0.85_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!2 Grass_after_Snow
  threshold(6,1:6) = (/0.84_fp,0.83_fp,-999.9_fp,&
       0.08_fp,0.10_fp,195._fp/)
  threshold(7,1:6) = (/0.85_fp,0.85_fp,-999.9_fp,&
       0.10_fp,-999.9_fp,190._fp/)
  threshold(8,1:6) = (/0.86_fp,0.81_fp,-999.9_fp,&
       0.12_fp,-999.9_fp,200._fp/)
  threshold(9,1:6) = (/0.86_fp,0.81_fp,0.0_fp,&
       0.12_fp,-999.9_fp,189._fp/)
  threshold(10,1:6) = (/0.90_fp,0.81_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,195._fp/)

!3 RS_Snow (A)
  threshold(11,1:6) = (/0.80_fp,0.76_fp,-999.9_fp,&
       0.05_fp,-999.9_fp,185._fp/)
  threshold(12,1:6) = (/0.82_fp,0.78_fp,-999.9_fp,&
       -999.9_fp,0.25_fp,180._fp/)
  threshold(13,1:6) = (/0.90_fp,0.76_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,180._fp/)

!4 Powder  Snow
  threshold(14,1:6) = (/0.89_fp,0.73_fp,-999.9_fp,&
       0.20_fp,-999.9_fp,-999.9_fp/)
  threshold(15,1:6) = (/0.89_fp,0.75_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)
  threshold(16,1:6) = (/0.93_fp,0.72_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!5 RS_Snow (B)
  threshold(17,1:6) = (/0.82_fp,0.70_fp,-999.9_fp,&
       0.20_fp,-999.9_fp,160._fp/)
  threshold(18,1:6) = (/0.83_fp,0.70_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,160._fp/)

!6 RS_Snow (C)
  threshold(19,1:6) = (/0.75_fp,0.76_fp,-999.9_fp,&
       0.08_fp,-999.9_fp,172._fp/)
  threshold(20,1:6) = (/0.77_fp,0.72_fp,-999.9_fp,&
       0.12_fp,0.15_fp,175._fp/)
  threshold(21,1:6) = (/0.78_fp,0.74_fp,-999.9_fp,&
       -999.9_fp,0.20_fp,172._fp/)
  threshold(22,1:6) = (/0.80_fp,0.77_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,170._fp/)
  threshold(23,1:6) = (/0.82_fp,-999.9_fp,-999.9_fp,&
       0.15_fp,0.22_fp,170._fp/)
  threshold(24,1:6) = (/0.82_fp,0.73_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,170._fp/)

!7 RS_Snow (D)
  threshold(25,1:6) = (/0.75_fp,0.70_fp,-999.9_fp,&
       0.15_fp,0.25_fp,167._fp/)
  threshold(26,1:6) = (/0.77_fp,0.76_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)
  threshold(27,1:6) = (/0.80_fp,0.72_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)
  threshold(28,1:6) = (/0.77_fp,0.73_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

  threshold(29,1:6) = (/0.81_fp,0.71_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)
  threshold(30,1:6) = (/0.82_fp,0.69_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!8 Thin Crust Snow
  threshold(31,1:6) = (/0.88_fp,0.58_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!9 RS_Snow (E)
  threshold(32,1:6) = (/0.73_fp,0.67_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!10 Bottom Crust Snow (A)
  threshold(33,1:6) = (/0.83_fp,0.66_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!11 Shallow Snow
  threshold(34,1:6) = (/0.82_fp,0.60_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!12 Deep Snow
  threshold(35,1:6) = (/0.77_fp,0.60_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!13 Crust Snow
  threshold(36,1:6) = (/0.77_fp,0.7_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!14 Medium Snow
  threshold(37,1:6) = (/-999.9_fp,0.55_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!15 Bottom Crust Snow(B)
  threshold(38,1:6) = (/0.74_fp,-999.9_fp,-999.9_fp,&
       -999.9_fp,-999.9_fp,-999.9_fp/)

!16 Thick Crust Snow
! lowest priority: No constraints

! Sixteen candidate snow emissivity spectra

  em(1, 1: N_FREQUENCY) = WET_SNOW_EMISS(1:N_FREQUENCY)
  em(2, 1: N_FREQUENCY) = GRASS_AFTER_SNOW_EMISS(1:N_FREQUENCY)
  em(3, 1: N_FREQUENCY) = RS_SNOW_A_EMISS(1:N_FREQUENCY)
  em(4, 1: N_FREQUENCY) = POWDER_SNOW_EMISS(1:N_FREQUENCY)
  em(5, 1: N_FREQUENCY) = RS_SNOW_B_EMISS(1:N_FREQUENCY)
  em(6, 1: N_FREQUENCY) = RS_SNOW_C_EMISS(1:N_FREQUENCY)
  em(7, 1: N_FREQUENCY) = RS_SNOW_D_EMISS(1:N_FREQUENCY)
  em(8, 1: N_FREQUENCY) = THIN_CRUST_SNOW_EMISS(1:N_FREQUENCY)
  em(9, 1: N_FREQUENCY) = RS_SNOW_E_EMISS(1:N_FREQUENCY)
  em(10, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_A_EMISS(1:N_FREQUENCY)
  em(11, 1: N_FREQUENCY) = SHALLOW_SNOW_EMISS(1:N_FREQUENCY)
  em(12, 1: N_FREQUENCY) = DEEP_SNOW_EMISS(1:N_FREQUENCY)
  em(13, 1: N_FREQUENCY) = CRUST_SNOW_EMISS(1:N_FREQUENCY)
  em(14, 1: N_FREQUENCY) = MEDIUM_SNOW_EMISS(1:N_FREQUENCY)
  em(15, 1: N_FREQUENCY) = BOTTOM_CRUST_SNOW_B_EMISS(1:N_FREQUENCY)
  em(16, 1: N_FREQUENCY) = THICK_CRUST_SNOW_EMISS(1:N_FREQUENCY)


  freq = FREQUENCY_DEFAULT


!***  DEFINE SIX DISCRIMINATORS

  dtb(1) = tb(1) - tb(2)
  dtb(2) = tb(2) - tb(4)
  dtb(3) = tb(2) - tb(5)
  dtb(4) = tb(3) - tb(5)
  dtb(5) = tb(4) - tb(5)
  tb150  = tb(5)

  LI = LI_coe(0)
  do i=0,1
     LI = LI + LI_coe(2*i+1)*tb(i+1) + LI_coe(2*i+2)*tb(i+1)*tb(i+1)
  end do
  LI = LI + LI_coe(nLIcoe-1)*ts

  HI = HI_coe(0)
  do i=0,4
     HI = HI + HI_coe(2*i+1)*tb(i+1) + HI_coe(2*i+2)*tb(i+1)*tb(i+1)
  end do
  HI = HI + HI_coe(nHIcoe-1)*ts

  do num=1,nind-1
     DI(num) = DI_coe(num,0) + DI_coe(num,1)*tb(2)
     do i=1,5
        DI(num) = DI(num) + DI_coe(num,1+i)*DTB(i)
     end do
     DI(num) = DI(num) +  DI_coe(num,ncoe-1)*ts
  end do

!*** DEFINE FIVE INDIES
  !HI = DI(0) - DI(3)
  DS1 = DI(1) + DI(2)
  DS2 = DI(4) + DI(5)
  DS3 = DS1 + DS2 + DI(3)

  index_in(1) = LI
  index_in(2) = HI
  index_in(3) = DS1
  index_in(4) = DS2
  index_in(5) = DS3
  index_in(6) = tb150

!*** IDENTIFY SNOW TYPE


! Initialization
  md0 = 1
  snow_type = ncand
  pick_status = .false.

! Pick one snow type
! Check all possible selections for six thresholds for each snow type
  do i = 1, ncand - 1
     md1 = nmodel(i)
     do j = md0, md1
        npass = 0
        do k = 1 , nind
           threshold0(k) = threshold(j,k)
        end do
        CALL six_indices(nind,index_in,threshold0,tindex)

! Corrections
        if((i == 5)  .and. (index_in(2) >  0.75_fp)) tindex(2) = .false.
        if((i == 5)  .and. (index_in(4) >  0.20_fp)                        &
             .and. (index_in(1) >  0.88_fp)) tindex(1) = .false.
        if((i == 10) .and. (index_in(1) <= 0.83_fp)) tindex(1) = .true.
        if((i == 13) .and. (index_in(2) <  0.52_fp)) tindex(2) = .true.
        do k = 1, nind
           if(.not.tindex(k)) exit
           npass = npass + 1
        end do
        if(npass == nind) exit
     end do

     if(npass == nind) then
        pick_status = .true.
        snow_type  = i
     end if
     if(pick_status) exit
     md0 = md1 + 1
  end do

  discriminator(1) = LI + DI(1)
  discriminator(2) = LI
  discriminator(3) = DI(4) + HI
  discriminator(4) = LI - DI(2)
  discriminator(5) = HI

  call em_interpolate(frequency,discriminator,emissivity,snow_type)

  em_vector(1) = emissivity
  em_vector(2) = emissivity

  return
end subroutine AMSU_ABTs



subroutine six_indices(nind,index_in,threshold,tindex)

!----------------------------------------------------------------------------------------------------------!
!$$$  subprogram documentation block
!
! subprogram:
!
!   prgmmr: Banghua Yan                 org: nesdis              date: 2003-08-18
!
! abstract:
!
! program history log:
!
! input argument list:
!
!      nind        -  Number of threshold in decision trees
!                     to identify each snow type  ( = 6)
!      index_in    -  six indices to discriminate snow type
!      threshold   -  Thresholds in decision trees to identify snow types
!
! output argument list:
!
!      tindex      - state vaiable to show surface snow emissivity feature
!              tindex[ ] = .T.: snow emissivity feature matches the
!                                corresponding threshold for certain snow type
!              tindex[ ] = .F.: snow emissivity feature doesn't match the
!                                corresponding threshold for certain snow type
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------------------------------!

  integer ::  i,nind
  real(fp)    ::  index_in(*),threshold(*)
  logical ::  tindex(*)

  do i=1,nind
     tindex(i) = .false.
     if (threshold(i) .eq. -999.9_fp) then
        tindex(i) = .true.
     else
        if ( (i .le. 2) .or. (i .gt. (nind-1)) ) then
           if (index_in(i) .ge. threshold(i)) tindex(i) = .true.
        else
           if (index_in(i) .le. threshold(i)) tindex(i) = .true.
        end if
     end if
  end do
  return

end subroutine six_indices


subroutine AMSU_AB(frequency,tb,snow_type,em_vector)

!----------------------------------------------------------------------------------------------------------!
!$$$  subprogram documentation block
!
! subprogram:
!
!   prgmmr: Banghua Yan                 org: nesdis              date: 2003-08-18
!
! abstract:
!         Calculate the emissivity discriminators and interpolate/extrapolate
!  emissivity at required frequency with respect to option AMSUAB
!
! program history log:
!
! input argument list:
!
!      frequency    -  frequency in GHz
!      theta        -  local zenith angle (not used here)
!      tb[1]~tb[5]  -  brightness temperature at five AMSU-A & B window channels:
!                              tb[1] : 23.8 GHz
!                              tb[2] : 31.4 GHz
!                              tb[3] : 50.3 GHz
!                              tb[4] : 89   GHz
!                              tb[5] : 150  GHz
!
! output argument list:
!
!     em_vector[1] and [2] - emissivity at two polarizations.
!                            set esv = esh here and will be updated
!     snow_type       - snow type (reference [2])
!
! important internal variables:
!
!     coe    - fitting coefficients to estimate discriminator at 23.8 ~ 150  GHz
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!----------------------------------------------------------------------------------------------------------!

  integer,parameter:: nch =10,nwch = 5,ncoe = 10
  real(fp)    :: tb(*),frequency
  real(fp)    :: em_vector(*),emissivity,discriminator(nwch)
  integer :: i,snow_type,k,ich,nvalid_ch
  real(fp)  :: coe(nwch*(ncoe+1))

  save coe

! Fitting Coefficients at 23.8 GHz: Using Tb1 ~ Tb3
  coe(1:7) = (/&
       -1.326040e+000_fp,  2.475904e-002_fp, &
       -5.741361e-005_fp, -1.889650e-002_fp, &
       6.177911e-005_fp,  1.451121e-002_fp, &
       -4.925512e-005_fp/)

! Fitting Coefficients at 31.4 GHz: Using Tb1 ~ Tb3
  coe(12:18) = (/ &
       -1.250541e+000_fp,  1.911161e-002_fp, &
       -5.460238e-005_fp, -1.266388e-002_fp, &
       5.745064e-005_fp,  1.313985e-002_fp, &
       -4.574811e-005_fp/)

! Fitting Coefficients at 50.3 GHz: Using Tb1 ~ Tb3
  coe(23:29) = (/  &
       -1.246754e+000_fp,  2.368658e-002_fp, &
       -8.061774e-005_fp, -3.206323e-002_fp, &
       1.148107e-004_fp,  2.688353e-002_fp, &
       -7.358356e-005_fp/)

! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4
  coe(34:42) = (/ &
       -1.278780e+000_fp,  1.625141e-002_fp, &
       -4.764536e-005_fp, -1.475181e-002_fp, &
       5.107766e-005_fp,  1.083021e-002_fp, &
       -4.154825e-005_fp,  7.703879e-003_fp, &
       -6.351148e-006_fp/)

! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb5
  coe(45:55) = (/&
     -1.691077e+000_fp,  3.352403e-002_fp, &
     -7.310338e-005_fp, -4.396138e-002_fp, &
     1.028994e-004_fp,  2.301014e-002_fp, &
     -7.070810e-005_fp,  1.270231e-002_fp, &
     -2.139023e-005_fp, -2.257991e-003_fp, &
     1.269419e-005_fp/)

! save coe

! Calculate emissivity discriminators at five AMSU window channels
  do ich = 1, nwch
     discriminator(ich) = coe(1+(ich-1)*11)
     if (ich .le. 3) nvalid_ch = 3
     if (ich .eq. 4) nvalid_ch = 4
     if (ich .eq. 5) nvalid_ch = 5
     do i=1,nvalid_ch
        discriminator(ich) = discriminator(ich) + coe((ich-1)*11 + 2*i)*tb(i) +  &
             coe((ich-1)*11 + 2*i+1)*tb(i)*tb(i)
     end do
  end do
!  Identify one snow emissivity spectrum and interpolate/extrapolate emissivity
!  at a required frequency
  call em_interpolate(frequency,discriminator,emissivity,snow_type)

  em_vector(1) = emissivity
  em_vector(2) = emissivity

  return
end subroutine AMSU_AB


subroutine AMSU_ATs(frequency,tba,ts,snow_type,em_vector)

!-----------------------------------------------------------------------------------------------------------!
!$$$  subprogram documentation block
!
! subprogram:
!
!   prgmmr:Banghua Yan                 org: nesdis              date: 2003-08-18
!
! abstract:
!         Calculate the emissivity discriminators and interpolate/extrapolate
!  emissivity at required frequency with respect to secenery AMSUAB
!
! program history log:
!
! input argument list:
!
!      frequency        -  frequency in GHz
!      theta            -  local zenith angle (not used here)
!      ts               -  surface temperature
!      tba[1] ~ tba[4]  -  brightness temperature at five AMSU-A window channels:
!                              tba[1] : 23.8 GHz
!                              tba[2] : 31.4 GHz
!                              tba[3] : 50.3 GHz
!                              tba[4] : 89   GHz
! output argument list:
!
!     em_vector[1] and [2]  -  emissivity at two polarizations.
!                              set esv = esh here and will be updated
!     snow_type        -  snow type (reference [2])
!
! important internal variables:
!
!     coe      - fitting coefficients to estimate discriminator at 23.8 ~ 150  GHz
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!-----------------------------------------------------------------------------------------------------------!

  integer,parameter:: nch =10,nwch = 5,ncoe = 9
  real(fp)    :: tba(*)
  real(fp)    :: em_vector(*),emissivity,ts,frequency,discriminator(nwch)
  integer :: snow_type,i,k,ich,nvalid_ch
  real(fp)  :: coe(nch*(ncoe+1))

  save coe

! Fitting Coefficients at 23.8 GHz: Using Tb1, Tb2 and Ts
  coe(1:6) = (/ &
       8.210105e-001_fp,  1.216432e-002_fp,  &
       -2.113875e-005_fp, -6.416648e-003_fp, &
       1.809047e-005_fp, -4.206605e-003_fp /)

! Fitting Coefficients at 31.4 GHz: Using Tb1, Tb2 and Ts
  coe(11:16) = (/ &
       7.963632e-001_fp,  7.215580e-003_fp,  &
       -2.015921e-005_fp, -1.508286e-003_fp,  &
       1.731405e-005_fp, -4.105358e-003_fp /)

! Fitting Coefficients at 50.3 GHz: Using Tb1, Tb2, Tb3 and Ts
  coe(21:28) = (/ &
       1.724160e+000_fp,  5.556665e-003_fp, &
       -2.915872e-005_fp, -1.146713e-002_fp, &
       4.724243e-005_fp,  3.851791e-003_fp, &
       -5.581535e-008_fp, -5.413451e-003_fp /)

! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4 and Ts
  coe(31:40) = (/ &
       9.962065e-001_fp,  1.584161e-004_fp, &
       -3.988934e-006_fp,  3.427638e-003_fp, &
       -5.084836e-006_fp, -6.178904e-004_fp, &
       1.115315e-006_fp,  9.440962e-004_fp, &
       9.711384e-006_fp, -4.259102e-003_fp /)

! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb4 and Ts
  coe(41:50) = (/ &
       -5.244422e-002_fp,  2.025879e-002_fp,  &
       -3.739231e-005_fp, -2.922355e-002_fp, &
       5.810726e-005_fp,  1.376275e-002_fp, &
       -3.757061e-005_fp,  6.434187e-003_fp, &
       6.190403e-007_fp, -2.944785e-003_fp/)

! save coe

! Calculate emissivity discriminators at five AMSU window channels
  DO ich = 1, nwch
     discriminator(ich) = coe(1+(ich-1)*10)
     if (ich .le. 2) nvalid_ch = 2
     if (ich .eq. 3) nvalid_ch = 3
     if (ich .ge. 4) nvalid_ch = 4
     do i=1,nvalid_ch
        discriminator(ich) = discriminator(ich) + coe((ich-1)*10 + 2*i)*tba(i) +  &
             coe((ich-1)*10 + 2*i+1)*tba(i)*tba(i)
     end do
     discriminator(ich) = discriminator(ich) + coe( (ich-1)*10 + (nvalid_ch+1)*2 )*ts
  end do

  call em_interpolate(frequency,discriminator,emissivity,snow_type)

  em_vector(1) = emissivity
  em_vector(2) = emissivity

  return
end subroutine AMSU_ATs


subroutine AMSU_amsua(frequency,tba,snow_type,em_vector)

!------------------------------------------------------------------------------------------------------------!
!$$$  subprogram documentation block
!
! subprogram:
!
!   prgmmr: Banghua Yan                 org: nesdis              date: 2003-08-18
!
! abstract:
!         Calculate the emissivity discriminators and interpolate/extrapolate
!  emissivity at required frequency with respect to secenery AMSUA
!
! program history log:
!
! input argument list:
!
!      frequency      -  frequency in GHz
!      theta          -  local zenith angle (not used here)
!      tba[1]~tba[4]  -  brightness temperature at five AMSU-A window channels:
!                            tba[1] : 23.8 GHz
!                            tba[2] : 31.4 GHz
!                            tba[3] : 50.3 GHz
!                            tba[4] : 89   GHz
!
! output argument list:
!
!     em_vector(1) and (2)  -  emissivity at two polarizations.
!                              set esv = esh here and will be updated
!     snow_type        -  snow type
!
! important internal variables:
!
!     coe      - fitting coefficients to estimate discriminator at 23.8 ~ 150  GHz
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!--------------------------------------------------------------------------------------------------------!

  integer,parameter:: nch =10,nwch = 5,ncoe = 8
  real(fp)    :: tba(*)
  real(fp)    :: em_vector(*),emissivity,frequency,discriminator(nwch)
  integer :: snow_type,i,k,ich,nvalid_ch
  real(fp)  :: coe(50)
  save coe

! Fitting Coefficients at 23.8 GHz: Using Tb1 ~ Tb3
  coe(1:7) = (/ &
       -1.326040e+000_fp,  2.475904e-002_fp, -5.741361e-005_fp, &
       -1.889650e-002_fp,  6.177911e-005_fp,  1.451121e-002_fp, &
       -4.925512e-005_fp/)

! Fitting Coefficients at 31.4 GHz: Using Tb1 ~ Tb3
  coe(11:17) = (/ &
       -1.250541e+000_fp,  1.911161e-002_fp, -5.460238e-005_fp, &
       -1.266388e-002_fp,  5.745064e-005_fp,  1.313985e-002_fp, &
       -4.574811e-005_fp/)

! Fitting Coefficients at 50.3 GHz: Using Tb1 ~ Tb3
  coe(21:27) = (/ &
       -1.246754e+000_fp,  2.368658e-002_fp, -8.061774e-005_fp, &
       -3.206323e-002_fp,  1.148107e-004_fp,  2.688353e-002_fp, &
       -7.358356e-005_fp/)

! Fitting Coefficients at 89 GHz: Using Tb1 ~ Tb4
  coe(31:39) = (/ &
       -1.278780e+000_fp, 1.625141e-002_fp, -4.764536e-005_fp, &
       -1.475181e-002_fp, 5.107766e-005_fp,  1.083021e-002_fp, &
       -4.154825e-005_fp,  7.703879e-003_fp, -6.351148e-006_fp/)

! Fitting Coefficients at 150 GHz: Using Tb1 ~ Tb4
  coe(41:49) = (/ &
       -1.624857e+000_fp, 3.138243e-002_fp, -6.757028e-005_fp, &
       -4.178496e-002_fp, 9.691893e-005_fp,  2.165964e-002_fp, &
       -6.702349e-005_fp, 1.111658e-002_fp, -1.050708e-005_fp/)

! save coe


! Calculate emissivity discriminators at five AMSU window channels
  do ich = 1, nwch
     discriminator(ich) = coe(1+(ich-1)*10)
     if (ich .le. 2) nvalid_ch = 3
     if (ich .ge. 3) nvalid_ch = 4
     do i=1,nvalid_ch
        discriminator(ich) = discriminator(ich) + coe((ich-1)*10 + 2*i)*tba(i) +  &
             coe((ich-1)*10 + 2*i+1)*tba(i)*tba(i)
     end do
  end do

! Quality Control
  if(discriminator(4) .gt. discriminator(2))   &
       discriminator(4) = discriminator(2) + (150.0_fp - 89.0_fp)*  &
       (discriminator(5) - discriminator(2))/ &
       (150.0_fp - 31.4_fp)

! Quality control at 50.3 GHz
  if((discriminator(3) .gt. discriminator(2)) .or.  &
       (discriminator(3) .lt. discriminator(4)))      &
       discriminator(3) = discriminator(2) + (89.0_fp - 50.3_fp)*   &
       (discriminator(4) - discriminator(2))/(89.0_fp - 31.4_fp)

  call em_interpolate(frequency,discriminator,emissivity,snow_type)

  em_vector(1) = emissivity
  em_vector(2) = emissivity

  return
end subroutine AMSU_amsua


subroutine AMSU_BTs(frequency,tbb,ts,snow_type,em_vector)

!-------------------------------------------------------------------------------------------------------!
!$$$  subprogram documentation block
!
! subprogram:
!
!   prgmmr: Banghua Yan                 org: nesdis              date: 2003-08-18
!
! abstract:
!         Calculate the emissivity discriminators and interpolate/extrapolate
!  emissivity at required frequency with respect to secenery BTs
!
! program history log:
!
! input argument list:
!
!      frequency        -  frequency in GHz
!      theta            -  local zenith angle (not used here)
!      ts               -  surface temperature in degree
!      tbb[1] ~ tbb[2]  -  brightness temperature at five AMSU-B window channels:
!                              tbb[1] : 89  GHz
!                              tbb[2] : 150 GHz
!
! output argument list:
!
!     em_vector(1) and (2)  -  emissivity at two polarizations.
!                              set esv = esh here and will be updated
!     snow_type        -  snow type
!
! important internal variables:
!
!     coe      - fitting coefficients to estimate discriminator at 31.4 ~ 150  GHz
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!-------------------------------------------------------------------------------------------------------!

  integer,parameter:: nch =10,nwch = 3,ncoe = 5
  real(fp)    :: tbb(*)
  real(fp)    :: em_vector(*),emissivity,ts,frequency,ed0(nwch),discriminator(5)
  integer :: snow_type,i,k,ich,nvalid_ch
  real(fp)  :: coe(nch*(ncoe+1))
  save coe

! Fitting Coefficients at 31.4 GHz: Using Tb4, Tb5 and Ts
  coe(1:6) = (/ 3.110967e-001_fp,  1.100175e-002_fp, -1.677626e-005_fp,    &
       -4.020427e-003_fp,  9.242240e-006_fp, -2.363207e-003_fp/)
! Fitting Coefficients at 89 GHz: Using Tb4, Tb5 and Ts
  coe(11:16) = (/  1.148098e+000_fp,  1.452926e-003_fp,  1.037081e-005_fp, &
       1.340696e-003_fp, -5.185640e-006_fp, -4.546382e-003_fp /)
! Fitting Coefficients at 150 GHz: Using Tb4, Tb5 and Ts
  coe(21:26) = (/ 1.165323e+000_fp, -1.030435e-003_fp,  4.828009e-006_fp,  &
       4.851731e-003_fp, -2.588049e-006_fp, -4.990193e-003_fp/)
! save coe

! Calculate emissivity discriminators at five AMSU window channels
  do ich = 1, nwch
     ed0(ich) = coe(1+(ich-1)*10)
     nvalid_ch = 2
     do i=1,nvalid_ch
        ed0(ich) = ed0(ich) + coe((ich-1)*10 + 2*i)*tbb(i) +   &
             coe((ich-1)*10 + 2*i+1)*tbb(i)*tbb(i)
     end do
     ed0(ich) = ed0(ich) + coe( (ich-1)*10 + (nvalid_ch+1)*2 )*ts
  end do

! Quality control
  if(ed0(2) .gt. ed0(1))     &
       ed0(2) = ed0(1) + (150.0_fp - 89.0_fp)*(ed0(3) - ed0(1)) / &
       (150.0_fp - 31.4_fp)

! Match the format of the input variable
! Missing value at 23.8 GHz
  discriminator(1) = -999.9_fp;  discriminator(2) = ed0(1)
! Missing value at 50.3 GHz
  discriminator(3) = -999.9_fp; discriminator(4) = ed0(2); discriminator(5) = ed0(3)

  call em_interpolate(frequency,discriminator,emissivity,snow_type)

  em_vector(1) = emissivity
  em_vector(2) = emissivity

  return
end subroutine AMSU_BTs


subroutine AMSU_amsub(frequency,tbb,snow_type,em_vector)


!-------------------------------------------------------------------------------------------------------!
!$$$  subprogram documentation block
!
! subprogram:
!
!   prgmmr: Banghua Yan                 org: nesdis              date: 2003-08-18
!
! abstract:
!         Calculate the emissivity discriminators and interpolate/extrapolate
!  emissivity at required frequency with respect to secenery AMSUB
!
! program history log:
!
! input argument list:
!
!      frequency        -  frequency in GHz
!      theta            -  local zenith angle (not used here)
!      tbb[1] ~ tbb[2]  -  brightness temperature at five AMSU-B window channels:
!                              tbb[1] : 89  GHz
!                              tbb[2] : 150 GHz
!
! output argument list:
!     em_vector(1) and (2)  -  emissivity at two polarizations.
!                              set esv = esh here and will be updated
!     snow_type        -  snow type (reference [2])
!
! important internal variables:
!
!     coe    - fitting coefficients to estimate discriminator at 31.4 ~ 150  GHz
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!-------------------------------------------------------------------------------------------------------!

  integer,parameter:: nch =10,nwch = 3,ncoe = 4
  real(fp)    :: tbb(*)
  real(fp)    :: em_vector(*),emissivity,frequency,ed0(nwch),discriminator(5)
  integer :: snow_type,i,k,ich,nvalid_ch
  real(fp)  :: coe(50)
  save coe

! Fitting Coefficients at 31.4 GHz: Using Tb4, Tb5
  coe(1:5) = (/-4.015636e-001_fp,9.297894e-003_fp, -1.305068e-005_fp, &
       3.717131e-004_fp, -4.364877e-006_fp/)
! Fitting Coefficients at 89 GHz: Using Tb4, Tb5
  coe(11:15) = (/-2.229547e-001_fp, -1.828402e-003_fp,1.754807e-005_fp, &
       9.793681e-003_fp, -3.137189e-005_fp/)
! Fitting Coefficients at 150 GHz: Using Tb4, Tb5
  coe(21:25) = (/-3.395416e-001_fp,-4.632656e-003_fp,1.270735e-005_fp, &
       1.413038e-002_fp,-3.133239e-005_fp/)
! save coe

! Calculate emissivity discriminators at five AMSU window channels
  do ich = 1, nwch
     ed0(ich) = coe(1+(ich-1)*10)
     nvalid_ch = 2
     do i=1,nvalid_ch
        ed0(ich) = ed0(ich) + coe((ich-1)*10 + 2*i)*tbb(i) +  &
             coe((ich-1)*10 + 2*i+1)*tbb(i)*tbb(i)
     end do
  end do

! Quality Control
  if(ed0(2) .gt. ed0(1))     &
       ed0(2) = ed0(1) + (150.0_fp - 89.0_fp) * &
       (ed0(3) - ed0(1))/(150.0_fp - 31.4_fp)

! Match the format of the input variable
! Missing value at 23.8 GHz
  discriminator(1) = -999.9_fp; discriminator(2) = ed0(1)
! Missing value at 50.3 GHz
  discriminator(3) = -999.9_fp; discriminator(4) = ed0(2); discriminator(5) = ed0(3)

  call em_interpolate(frequency,discriminator,emissivity,snow_type)

  em_vector(1) = emissivity
  em_vector(2) = emissivity

  return
end subroutine AMSU_amsub


subroutine AMSU_ALandEM_Snow(theta,frequency,snow_depth,ts,snow_type,em_vector)


!------------------------------------------------------------------------------------------------------------
!$$$  subprogram documentation block
!
! subprogram:
!
!   prgmmr: Banghua Yan                 org: nesdis              date: 2003-08-18
!
! abstract:
!         Calculate the emissivity at required frequency with respect to option MODL
!   using the NESDIS_LandEM and a bias correction algorithm, where the original NESDIS_LandEM with a
!   bias correction algorithm is referred to as value-added NESDIS_LandEM or AlandEM.
!
! program history log:
!
! input argument list:
!
!      frequency        -  frequency in GHz
!      theta            -  local zenith angle in degree
!      snow_depth       -  snow depth in mm
!      ts           -  surface temperature
!
! output argument list:
!
!     em_vector(1) and (2)  -  emissivity at two polarizations.
!                              set esv = esh here and will be updated
!       snow_type        -  snow type
!
! important internal variables:
!
!    esv_3w and esh_3w   -  initial emissivity discriminator at two polarizations
!                           at three AMSU window channels computed using NESDIS_LandEM
!    esv_3w[1] and esh_3w[1] : 31.4 GHz
!    esv_3w[2] and esh_3w[2] : 89   GHz
!    esv_3w[3] and esh_3w[3] : 150  GHz
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!------------------------------------------------------------------------------------------------------------

  integer :: nw_ind
  parameter(nw_ind=3)
  real(fp) theta, frequency, freq,snow_depth, ts, em_vector(2)
  real(fp) esv,esh,esh0,esv0,theta0
  integer snow_type,ich
  real(fp)   freq_3w(nw_ind),esh_3w(nw_ind),esv_3w(nw_ind)
  complex(fp)  eair
  save freq_3w

  freq_3w = (/31.4_fp,89.0_fp,150.0_fp/)

  eair = cmplx(one,-zero,fp)

  snow_type = -999

  call NESDIS_LandEM(theta, frequency,0.0_fp,0.0_fp,ts,ts,0.0_fp,9,13,snow_depth,esh0,esv0)

  theta0 = theta
  do ich = 1, nw_ind
     freq =freq_3w(ich)
     theta = theta0
     call NESDIS_LandEM(theta, freq,0.0_fp,0.0_fp,ts,ts,0.0_fp,9,13,snow_depth,esh,esv)
     esv_3w(ich) = esv
     esh_3w(ich) = esh
  end do

  call ems_adjust(theta,frequency,snow_depth,ts,esv_3w,esh_3w,em_vector,snow_type)

  return

end subroutine AMSU_ALandEM_Snow



subroutine ems_adjust(theta,frequency,depth,ts,esv_3w,esh_3w,em_vector,snow_type)


!------------------------------------------------------------------------------------------------------------
!$$$  subprogram documentation block
!
! subprogram:
!
!   prgmmr: Banghua Yan                 org: nesdis              date: 2003-08-18
!
! abstract:
!         Calculate the emissivity discriminators and interpolate/extrapolate
!
!  emissivity at required frequency with respect to secenery MODL
!
! program history log:
!
! input argument list:
!
!      frequency   -  frequency in GHz
!      theta       -  local zenith angle in degree
!      depth       -  snow depth in mm
!      ts          -  surface temperature
!
! output argument list:
!
!     em_vector(1) and (2)  -  emissivity at two polarizations.
!                              set esv = esh here and will be updated
!     snow_type        -  snow type
!
! important internal variables:
!
!     dem_coe  -  fitting coefficients to compute discriminator correction value
!              dem_coe[1,*]   : 31.4 GHz
!              dem_coe[2,*]   : 89   GHz
!              dem_coe[3,*]   : 150  GHz
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!------------------------------------------------------------------------------------------------------------

  integer,parameter:: nch=10,nw_3=3

  integer,parameter:: ncoe=6

  real(fp),parameter  :: earthrad = 6374._fp, satheight = 833.4_fp

  integer         :: snow_type,ich,k

  real(fp)    :: theta,frequency,depth,ts,esv_3w(*),esh_3w(*)

  real(fp)    :: discriminator(5),emmod(nw_3),dem(nw_3)

  real(fp)    :: emissivity,em_vector(2)

  real(Double)  :: dem_coe(nw_3,0:ncoe-1),sinthetas,costhetas,deg2rad

  save  dem_coe

  dem_coe(1,0:ncoe-1) = (/ 2.306844e+000_Double, -7.287718e-003_Double, &

       -6.433248e-004_Double,  1.664216e-005_Double,  &

       4.766508e-007_Double, -1.754184e+000_Double/)

  dem_coe(2,0:ncoe-1) = (/ 3.152527e+000_Double, -1.823670e-002_Double, &

       -9.535361e-004_Double,  3.675516e-005_Double,  &

       9.609477e-007_Double, -1.113725e+000_Double/)

  dem_coe(3,0:ncoe-1) = (/ 3.492495e+000_Double, -2.184545e-002_Double,  &

       6.536696e-005_Double,  4.464352e-005_Double, &

       -6.305717e-008_Double, -1.221087e+000_Double/)


!

  deg2rad = 3.14159_fp*pi/180.0_fp

  sinthetas = sin(theta*deg2rad)* earthrad/(earthrad + satheight)

  sinthetas = sinthetas*sinthetas

  costhetas = one - sinthetas

  do ich = 1, nw_3

     emmod(ich) = costhetas*esv_3w(ich) + sinthetas*esh_3w(ich)

  end do

  do ich=1,nw_3

     dem(ich) = dem_coe(ich,0) + dem_coe(ich,1)*ts + dem_coe(ich,2)*depth +   &

          dem_coe(ich,3)*ts*ts + dem_coe(ich,4)*depth*depth         +   &

          dem_coe(ich,5)*emmod(ich)

  end do

  emmod(1) = emmod(1) + dem(1)

  emmod(2) = emmod(2) + dem(2)

  emmod(3) = emmod(3) + dem(3)

! Match the format of the input variable

! Missing value at 23.8 GHz

  discriminator(1) = -999.9_fp

  discriminator(2) = emmod(1)

! Missing value at 50.3 GHz

  discriminator(3) = -999.9_fp

  discriminator(4) = emmod(2)

  discriminator(5) = emmod(3)

  call em_interpolate(frequency,discriminator,emissivity,snow_type)

  em_vector(1) = emissivity

  em_vector(2) = emissivity

  return

end subroutine ems_adjust


END MODULE NESDIS_AMSU_SnowEM_Module
