!
! NESDIS_SEAICE_PHYEM_MODULE
!
! Module containing the NESDIS microwave seaice emissivity physical model
!
!
! CREATION HISTORY:
!       Written by:     Banghua Yan, 13-Jun-2005, banghua.yan@noaa.gov
!                       Fuzhong Weng, fuzhong.weng@noaa.gov
!

MODULE NESDIS_SEAICE_PHYEM_MODULE


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


  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: NESDIS_SIce_Phy_EM


  ! -----------------
  ! Module parameters
  ! -----------------
  ! Version Id for the module
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  '$Id: NESDIS_SEAICE_PHYEM_MODULE.f90 29405 2013-06-20 20:19:52Z paul.vandelst@noaa.gov $'


CONTAINS


  subroutine NESDIS_SIce_Phy_EM(Frequency,                                                 &  ! INPUT
                                Angle,                                                     &  ! INPUT
                                Ts_ice,                                                    &  ! INPUT
                                Salinity,                                                  &  ! INPUT
                                Emissivity_H,                                              &  ! OUTPUT
                                Emissivity_V)                                                 ! OUPTUT

!-------------------------------------------------------------------------------------------------------------
!
! NAME:
!       NESDIS_SIce_Phy_EM
!
! PURPOSE:
!       Subroutine to simulate microwave sea ice emissivity
!
! REFERENCES:
!
! Ulaby, F.T., R. K., Moore, and A. K., Fung, Microwave remote sensing: active and passive, III,
! from theory to Applications, 2020-2025, Artech House Publishers, 1990.
!
! CATEGORY:
!       CRTM : Surface : MW SEA ICE EM
!
! LANGUAGE:
!       Fortran-95
!
! CALLING SEQUENCE:
!       CALL NESDIS_SIce_Phy_EM
!
! INPUT ARGUMENTS:
!
!         Frequency                Frequency User defines
!                                  This is the "I" dimension
!                                  UNITS:      GHz
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Scalar
!
!
!         Angle                    The angle values in degree
!                                  ** NOTE: THIS IS A MANDATORY MEMBER **
!                                  **       OF THIS STRUCTURE          **
!                                  UNITS:      Degrees
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Rank-1, (I)
!
!         Ts_ice                   Sea ice temperature
!                                  UNITS:      Kelvin, K
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Scalar
!
!        Salinity                  Sea water salinity (1/thousand)
!                                  UNITS:      N/A
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Scalar
!
!
!  INTERNAL ARGUMENTS:
!
!         theta                    viewing zenith angle
!                                  UNITS:  radian
!                                  TYPE:       REAL( fp )
!                                  DIMENSION:  Scalar
!
!
! 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
! CALLS:
!
!       permitivity    : Function to calculate the dielectric constant of saline water
!
!reflection_coefficient: Function to calculate the surface reflection coefficient using Fresnel equations
!
! RESTRICTIONS:
!       None.
!
!
! CREATION HISTORY:
!       Written by:     Banghua Yan, QSS Group Inc., Banghua.Yan@noaa.gov (28-May-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.
!
!------------------------------------------------------------------------------------------------------------


     real(fp) ::  Angle,theta, Ts_ice, Salinity
     real(fp) ::  Frequency, Emissivity_H, Emissivity_V, slam, f1,epsx, epsy, pi
     integer :: ITYPE
     COMPLEX(fp) :: eair, eice
     COMPLEX(fp) :: cos_theta_i,cos_theta_t,sin_theta_i,sin_theta_t
     COMPLEX(fp) :: rv_ice, rh_ice


      pi = acos(-1.0_fp)

      theta = Angle*pi/180.0_fp

      slam = 300.0_fp/Frequency

      f1=Frequency*1.0e+9_fp

      eair = cmplx(1.0_fp, 0.0_fp, fp)

      ITYPE = 2

      call  permitivity(ITYPE,Ts_ice, Salinity, f1, epsx, epsy)

      eice = cmplx(epsx, -epsy, fp)

      sin_theta_i = sin(cmplx(theta,0.0_fp, fp))

      cos_theta_i = sqrt(1.0_fp - sin_theta_i*sin_theta_i)

      sin_theta_t = sin_theta_i*sqrt(eair)/sqrt(eice)

      cos_theta_t = sqrt(1.0_fp - sin_theta_t*sin_theta_t)

      call reflection_coefficient(eair, eice, cos_theta_i, cos_theta_t, rv_ice, rh_ice)

      Emissivity_V = 1.0_fp  - abs(rv_ice)*abs(rv_ice)

      Emissivity_H = 1.0_fp  - abs(rh_ice)*abs(rh_ice)

     return

     end subroutine NESDIS_SIce_Phy_EM


     subroutine permitivity(itype, temp, s, f, e_real, e_img)

!--------------------------------------------------------------------------------------------------
!   prgmmr:  Fuzhong Weng and Banghua Yan                org: nesdis              date: 2005-06-13
!
!   abstract: this function calculates the dielectric constant of saline water
!
!    Reference:
!
!    Microwave remote sensing by Ulaby et al (1984)  pp 2022
!
!    Concepts:
!          saline water: water containing dissolved salts. The salinity of a solution
!                        is defined as the total mass of solid salt in grams dissolved
!                        in one kilogram of solution
!          sea ice     : a heterogeneous mixture of liquid-brine inclusions and air
!                        pockets interspersed within the ice medium
!          brine inclusions: containing salt and water, exercise a strong influence
!                        on the complex dielectric constant of the mixture due to their
!                        high complex-dielectric constant when compred to that of ice
!    Input variables:
!          itype: 1 --- saline (salt) water
!                 2 --- sea ice with mixtures like brine, pure ice and air
!          default: pure water
!
!          temp   : water temperature (K)
!          s      : salinty of sea ice water
!          f      : frequency (Hz)
!
!    Output variables:
!         e_real  : real part of complex dielectric constant of sea ice/water
!         e_img   : imaginery part ..
!    Code histry:
!       Code generated: Fuzhong Weng
!       Added brine component: Banghua Yan (09/09/2003)
!------------------------------------------------------------------------------

     real(fp) ::  eo, esw, eswo, a, tswo, tsw, b, sswo, c,d, fi, ssw
     real(fp) ::  temp,t, t2, t3, s2, s3, epsp, epspp
     real(fp) ::  ESWI, EW0, TPTW, FH, X, Y, TPTWFQ
     real(fp) ::  e_real, e_img

     real(fp) ::  s,t4,t5,sigma
     real(fp) ::  ebo,eb,sb,nb,eb_real,eb_imag
     real(fp) ::   den_ice,den_sice,den_brine

     real(fp) ::   vss, cc, den_ss

     real(fp) ::   f,pi,vb,sigmao,vi

     integer :: itype

     COMPLEX(fp) :: ebrine,eice,eair,esice,xx

     pi=acos(-1.0_fp)
     ESWI=4.9_fp
     t=temp-273.15_fp
     if (t.le.-65.0_fp) t = -65.0_fp
     t2=t*t
     t3=t2*t
     t4=t3*t
     t5=t4*t
     FH=f


   GET_DielectricC: SELECT CASE (ITYPE)

   CASE (1)

! saline water

      s2 = s*s
      s3 = s**3
      eo   = 8.854e-12
      eswo = 87.134 - 1.949e-1 * t - 1.276e-2 * t2 + 2.491e-4 * t3
      a =1.0 +1.613e-5 * t*s - 3.656e-3*s + 3.210e-5*s2 - 4.232e-7*s3
      esw = eswo * a
      tswo = 1.1109e-10 - 3.824e-12 * t + 6.938e-14*t2 -5.096e-16*t3
      b = 1.0 + 2.282e-5*t*s -7.638e-4*s - 7.760e-6 *s2 + 1.105e-8*s3
      tsw = tswo*b
      epsp = ESWI + (esw - ESWI)/ (1.0 + (f*tsw)**2)
      sswo = s*(0.18252-1.4619e-3*s+2.093e-5*s2-1.282e-7*s3)
      d  = 25.0 - t
      fi = d * (2.033e-2 + 1.266e-4 * d + 2.464e-6 * d * d -   &
        s * (1.849e-5 - 2.551e-7 * d + 2.551e-8 * d * d))

      ssw = sswo * exp(-fi)
      epspp = tsw * f * (esw - ESWI) / (1.0 + (tsw * f)**2)
      epspp = epspp + ssw/(2.0 * pi * eo * f)

      e_real =   epsp
      e_img = - epspp

Case (2) ! sea ice

      eice = cmplx(3.15_fp, -0.001_fp,fp)
      eair = cmplx(1.0_fp, 0.0_fp,fp)

! brine volume fraction: vb

      if (t .le. -8.2) then
          vb = 0.001*s*(-43.795/t + 1.189)
          if (t .lt. -22.9) vb = 0.001*s*(-43.795/(-22.9) + 1.189)
      else
         if (t .le. -2.06) then
             vb = 0.001*s*(-45.917/t + 0.930)
        else
            vb = 0.001*s*(-52.56/t  - 2.28)
          if (t .gt. -0.5)  vb = 0.001*s*(-52.56/(-0.5)  - 2.28)
        endif
      endif

!salinity of liquid brine

      if (t .le. -36.8) then
         sb = 508.18 + 14.535*t + 0.2018*t2
        if (t .lt. -43.2)sb = 508.18 + 14.535*(-43.2) + 0.2018*(-43.2)*(-43.2)
      endif
      if ((t .gt. -36.8) .and. (t .le. -22.9))   &
                        sb = 242.94 + 1.5299*t + 0.0429*t2
      if ((t .gt. -22.9) .and. (t .le. 8.2))       &
                        sb = 57.041 - 9.929*t - 0.16204*t2 - 0.002396*t3
      if(t .ge. -8.2) then
         sb = 1.725 - 18.756*t - 0.3964*t2
            if(t .gt. -2.0) sb = 1.725 - 18.756*(-2.0) - 0.3964*(-2.0)*(-2.0)
      endif

! Normality of brine solution: nb

      nb  = sb*(1.707*0.01 + 1.205*1.0E-5*sb+4.058E-9*sb*sb)
      d  = 25.0 - t
      a = 1.0 - 0.255*nb + 5.15*0.01*nb*nb-6.89*1.0e-3*nb*nb*nb
      b = 1.0 + 0.146*0.01*t*nb-4.89*0.01*nb                  &
          -2.97*0.01*nb*nb +  5.64 *1.0e-3*nb*nb*nb
      c = 1.0 - 1.96*0.01*d + 8.08*1.0e-5*d*d                 &
          - nb*d*(3.02*1.0e-5 + 3.92*1.0e-5*d+nb*(1.72*1.0e-5 - 6.58*1.0e-6*d ))

      ebo = 88.045-0.4147*t+6.295E-4*t2+1.075E-5*t3

      tswo = 1.1109E-10-3.824E-12*t+6.938E-14*t2-5.096E-16*t3
      sigmao = nb*(10.39 - 2.378*nb+0.683*nb*nb-0.135*nb*nb*nb+1.01*0.01*nb**4.0)
      eb = ebo*a
      tsw = tswo*b
      sigma = sigmao*c
      eo   = 8.854e-12   ! the permittivity of free space (p. 2024)
      eb_real = ESWI + (eb - ESWI)/ (1.0 + (f*tsw)**2.)
      epspp = tsw * f * (eb - ESWI) / (1.0 + (tsw * f)**2.)
      eb_imag = epspp + sigma/(2.0 * pi * eo * f)
      ebrine=cmplx(eb_real,-eb_imag,fp)

! density of sea ice

    den_sice = 0.87
    if(t .lt. -18.0) den_sice = 0.77
    den_ice   = 0.917 - 1.403e-4*t
    den_brine = 1.0 + 0.0008*sb

!Ice volume fraction
    vi = (den_sice - vb*den_brine)/den_ice
    if ((vi+vb) .gt. 1.0) vi = 1.0 - vb
    if (vi .lt. 0.0) vi = 0.0

    xx = vb*sqrt(ebrine)+vi*sqrt(eice)+(1.0-vi-vb)*sqrt(eair)
    esice = xx*xx
    e_real = real(esice,fp)
    e_img =  -aimag(esice)
    cc = 0.387
    den_ss = 1.5
    vss =cc * den_brine*vb/den_ss

CASE Default

    EW0=88.045-0.4147*t+6.295E-4*t2+1.075E-5*t3

    TPTW=1.1109E-10-3.824E-12*t+6.938E-14*t2-5.096E-16*t3

    X=TPTW*FH

    TPTWFQ=1+X*X

    Y=(EW0-ESWI)/TPTWFQ

    e_real=ESWI+Y

    e_img=-X*Y

END SELECT GET_DielectricC

return

end subroutine permitivity


     SUBROUTINE reflection_coefficient(em1, em2, cos_theta_i, cos_theta_t, rv, rh)
!--------------------------------------------------------------------------------------------------
!   prgmmr:  Fuzhong Weng and Banghua Yan                org: nesdis              date: 2005-06-13
!
!   abstract: this function calculates the surface reflection coefficient using Fresnel equations
!
!    Input Variables
!
!      cos_theta_i        ----  cosine incident angle (degree)
!      cos_theta_t        ----  cosine_transmitted angle (degree)
!      em1                ----  dielectric constant of the medium 1
!      em2                ----  dielectric constant of the medium 2
!
!    Output Variables
!
!      rv                 ----  reflection coefficient at vertical polarization
!      rh                 ----  reflection coefficient at horizontal polarization
!
!-----------------------------------------------------------------------------------------------

    COMPLEX(fp) :: rh, rv,cos_theta_i,cos_theta_t

    COMPLEX(fp) :: em1, em2, m1, m2

    m1 = sqrt(em1)

    m2 = sqrt(em2)

    rv = (m1*cos_theta_t- m2*cos_theta_i)/(m1*cos_theta_t + m2*cos_theta_i)

    rh = (m1*cos_theta_i- m2*cos_theta_t)/(m1*cos_theta_i + m2*cos_theta_t)

    return

    end subroutine reflection_coefficient


END MODULE NESDIS_SEAICE_PHYEM_MODULE
