Ellison.f90

Go to the documentation of this file.

!
! Ellison Ocean Permittivity module.
!
! Module containing routines to compute the complex permittivities for
! sea water based on
!
!   Ellison, W.J. et al. (2003) A comparison of ocean emissivity models
!     using the Advanced Microwave Sounding Unit, the Special Sensor
!     Microwave Imager, the TRMM Microwave Imager, and airborne radiometer
!     observations. Journal of Geophysical Research, v108, D21, Pages ACL 1,1-14
!     doi:10.1029/2002JD0032132
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 11-Apr-2007
!                       paul.vandelst@noaa.gov
!

MODULE ellison

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE type_kinds, ONLY: fp
  USE fundamental_constants, ONLY: pi, &
                                   e0 => permittivity, &  ! Permittivity of vacuum (F/m)
                                   k_to_c => standard_temperature ! Temperature units conversion
  ! Disable implicit typing
  IMPLICIT NONE


  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  ! ... Datatypes
  PUBLIC :: ivar_type
  ! ... Procedures
  PUBLIC :: ellison_ocean_permittivity
  PUBLIC :: ellison_ocean_permittivity_tl
  PUBLIC :: ellison_ocean_permittivity_ad


  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: module_version_id = &
  '$Id: Ellison.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'
  REAL(fp), PARAMETER :: zero   = 0.0_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 :: five   = 5.0_fp
  REAL(fp), PARAMETER :: twopi  = two*pi


  ! Scaling factors (used here for documenting the conversion
  ! to SI units of the double Debye model denominator)
  ! ---------------------------------------------------------
  REAL(fp), PARAMETER :: ps_to_s   = 1.0e-12_fp ! Picoseconds -> Seconds
  REAL(fp), PARAMETER :: ghz_to_hz = 1.0e+09_fp ! Gigahertz   -> Hertz
  REAL(fp), PARAMETER :: scale_factor = ps_to_s * ghz_to_hz


  ! Parameters for the Ellison et al (2003) permittivity model
  ! ----------------------------------------------------------
  ! The coefficients used to fit the Double Debye model
  REAL(fp), PARAMETER :: tau1_coeff(0:2)   = (/    17.535_fp, &
                                                 -0.61767_fp, &
                                                0.0089481_fp /)
  REAL(fp), PARAMETER :: tau2_coeff(0:3)   = (/    3.1842_fp, &
                                                 0.019189_fp, &
                                                -0.010873_fp, &
                                               0.00025818_fp /)
  REAL(fp), PARAMETER :: delta1_coeff(0:3) = (/    68.396_fp, &
                                                 -0.40643_fp, &
                                                 0.022832_fp, &
                                              -0.00053061_fp /)
  REAL(fp), PARAMETER :: delta2_coeff(0:3) = (/    4.7629_fp, &
                                                   0.1541_fp, &
                                                -0.033717_fp, &
                                               0.00084428_fp /)
  REAL(fp), PARAMETER :: einf_coeff(0:1)   = (/   5.31250_fp, &
                                               -0.0114770_fp /)
  REAL(fp), PARAMETER :: sigma_coeff(0:1) = (/      2.906_fp, &
                                                  0.09437_fp /)


  ! --------------------------------------
  ! Structure definition to hold forward
  ! variables across FWD, TL, and AD calls
  ! --------------------------------------
  TYPE :: ivar_type
    PRIVATE
    REAL(fp) :: t=zero                    ! Temperature in degC
    REAL(fp) :: f=zero, f2=zero, f0=zero  ! Frequency terms
    REAL(fp) :: tau1=zero  , tau2=zero    ! Relaxation frequencies
    REAL(fp) :: delta1=zero, delta2=zero  ! Delta terms
    REAL(fp) :: d1=zero    , d2=zero      ! Denominator terms
  END TYPE ivar_type


CONTAINS


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

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       Ellison_Ocean_Permittivity
!
! PURPOSE:
!       Subroutine to compute ocean permittivity according to the reference,
!         Ellison, W.J. et al. (2003) A comparison of ocean emissivity models
!           using the Advanced Microwave Sounding Unit, the Special Sensor
!           Microwave Imager, the TRMM Microwave Imager, and airborne radiometer
!           observations. Journal of Geophysical Research, v108, D21, ACL 1,1-14
!           doi:10.1029/2002JD0032132
!
! CALLING SEQUENCE:
!       CALL Ellison_Ocean_Permittivity( Temperature , & ! Input
!                                        Frequency   , & ! Input
!                                        Permittivity, & ! Output
!                                        iVar          ) ! Internal variable output
!
! INPUTS:
!       Temperature:   Sea surface temperature
!                      UNITS:      Kelvin (K)
!                      TYPE:       REAL(fp)
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
!       Frequency:     Frequency
!                      UNITS:      GHz
!                      TYPE:       REAL(fp)
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Permittivity:  Ocean permittivity
!                      UNITS:      N/A
!                      TYPE:       COMPLEX(fp)
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(OUT)
!
!       iVar:          Structure containing internal variables required for
!                      subsequent tangent-linear or adjoint model calls.
!                      The contents of this structure are NOT accessible
!                      outside of this module.
!                      UNITS:      N/A
!                      TYPE:       TYPE(iVar_type)
!                      DIMENSION:  Scalar
!                      ATTRIBUTES: INTENT(OUT)
!
! COMMENTS:
!       There is currently no salinity dependence.
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE ellison_ocean_permittivity( &
    Temperature , & ! Input
    Frequency   , & ! Input
    Permittivity, & ! Output
    iVar          ) ! Internal variable output
    ! Arguments
    REAL(fp),        INTENT(IN)     :: temperature
    REAL(fp),        INTENT(IN)     :: frequency
    COMPLEX(fp),     INTENT(OUT)    :: permittivity
    TYPE(ivar_type), INTENT(IN OUT) :: ivar
    ! Local variables
    REAL(fp) :: einf
    REAL(fp) :: re1, re2
    REAL(fp) :: ie1, ie2
    REAL(fp) :: sigma, iesigma
    REAL(fp) :: re, ie


    ! Compute the various polynomial components of the double Debye model
    ! -------------------------------------------------------------------
    ! Compute temperature value for polynomials
    ivar%t = temperature - k_to_c

    ! Compute the Debye model relaxation frequencies
    ! (eqn on pg ACL 1-4 of Ellison et al. 2003)
    ivar%tau1 = tau1_coeff(0) + ivar%t*(tau1_coeff(1) + &
                                  ivar%t*tau1_coeff(2))
    ivar%tau2 = tau2_coeff(0) + ivar%t*(tau2_coeff(1)   + &
                                  ivar%t*(tau2_coeff(2) + &
                                    ivar%t*tau2_coeff(3)))

    ! Compute the delta terms
    ! (eqn on pg ACL 1-4 of Ellison et al. 2003)
    ivar%delta1 = delta1_coeff(0) + ivar%t*(delta1_coeff(1) + &
                                      ivar%t*(delta1_coeff(2) + &
                                        ivar%t*delta1_coeff(3)))
    ivar%delta2 = delta2_coeff(0) + ivar%t*(delta2_coeff(1) + &
                                      ivar%t*(delta2_coeff(2) + &
                                        ivar%t*delta2_coeff(3)))

    ! Compute the "infinite" permittivity term
    ! (No coeffs provided in ref. Taken from existing code)
    einf = einf_coeff(0) + ivar%t*einf_coeff(1)


    ! Compute the permittivities using the double Debye model
    ! (eqn on pg ACL 1-3 of Ellison et al. 2003)
    ! -------------------------------------------------------
    ! The common frequency terms
    ivar%f  = twopi * frequency * scale_factor
    ivar%f2 = ivar%f**2
    ivar%f0 = twopi * frequency * ghz_to_hz * e0

    ! The denominators of the double Debye model
    ivar%d1 = one + ivar%f2*ivar%tau1**2
    ivar%d2 = one + ivar%f2*ivar%tau2**2

    ! The real parts of the "delta" terms
    re1 = ivar%delta1 / ivar%d1
    re2 = ivar%delta2 / ivar%d2

    ! The imaginary parts of the "delta" terms
    ie1 = ivar%delta1 * ivar%f * ivar%tau1 / ivar%d1
    ie2 = ivar%delta2 * ivar%f * ivar%tau2 / ivar%d2

    ! The conductivity term
    sigma   = sigma_coeff(0) + sigma_coeff(1)*ivar%t
    iesigma = sigma / ivar%f0

    ! Construct the complex permittivity, e = e' - j.e"
    re = re1 + re2 + einf
    ie = ie1 + ie2 + iesigma
    permittivity = cmplx(re, -ie, fp)

  END SUBROUTINE ellison_ocean_permittivity


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       Ellison_Ocean_Permittivity_TL
!
! PURPOSE:
!       Subroutine to compute the tangent-linear ocean permittivity according
!       to the reference,
!         Ellison, W.J. et al. (2003) A comparison of ocean emissivity models
!           using the Advanced Microwave Sounding Unit, the Special Sensor
!           Microwave Imager, the TRMM Microwave Imager, and airborne radiometer
!           observations. Journal of Geophysical Research, v108, D21, ACL 1,1-14
!           doi:10.1029/2002JD0032132
!
! CALLING SEQUENCE:
!       CALL Ellison_Ocean_Permittivity_TL( Temperature_TL , & ! Input
!                                           Permittivity_TL, & ! Output
!                                           iVar             ) ! Internal variable input
!
! INPUTS:
!       Temperature_TL:   Tangent-linear sea surface temperature
!                         UNITS:      Kelvin (K)
!                         TYPE:       REAL(fp)
!                         DIMENSION:  Scalar
!                         ATTRIBUTES: INTENT(IN)
!
!       iVar:             Structure containing internal variables required for
!                         subsequent tangent-linear or adjoint model calls.
!                         The contents of this structure are NOT accessible
!                         outside of this module.
!                         UNITS:      N/A
!                         TYPE:       TYPE(iVar_type)
!                         DIMENSION:  Scalar
!                         ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Permittivity_TL:  Tangent-linear ocean permittivity
!                         UNITS:      N/A
!                         TYPE:       COMPLEX(fp)
!                         DIMENSION:  Scalar
!                         ATTRIBUTES: INTENT(OUT)
!
! COMMENTS:
!       There is currently no salinity dependence.
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE ellison_ocean_permittivity_tl( &
    Temperature_TL , & ! Input
    Permittivity_TL, & ! Output
    iVar             ) ! Internal variable input
    ! Arguments
    REAL(fp),        INTENT(IN)  :: temperature_tl
    COMPLEX(fp),     INTENT(OUT) :: permittivity_tl
    TYPE(ivar_type), INTENT(IN)  :: ivar
    ! Local variables
    REAL(fp) :: t_tl
    REAL(fp) :: tau1_tl, tau2_tl
    REAL(fp) :: delta1_tl, delta2_tl, einf_tl
    REAL(fp) :: d1_tl, d2_tl
    REAL(fp) :: d12, d22
    REAL(fp) :: re1_tl, re2_tl
    REAL(fp) :: ie1_tl, ie2_tl
    REAL(fp) :: sigma_tl, iesigma_tl
    REAL(fp) :: re_tl, ie_tl


    ! Compute the tangent-linear for of the various
    ! polynomial components of the double Debye model
    ! -----------------------------------------------
    ! Compute temperature value for polynomials
    t_tl = temperature_tl

    ! Compute the tangent-linear Debye model relaxation frequencies
    ! (eqn on pg ACL 1-4 of Ellison et al. 2003)
    tau1_tl = (tau1_coeff(1) + ivar%t*two*tau1_coeff(2)) * t_tl
    tau2_tl = (tau2_coeff(1) + ivar%t*(two*tau2_coeff(2) + ivar%t*three*tau2_coeff(3))) * t_tl

    ! Compute the tangent-linear delta terms
    ! (eqn on pg ACL 1-4 of Ellison et al. 2003)
    delta1_tl = (delta1_coeff(1) + ivar%t*(two*delta1_coeff(2) + ivar%t*three*delta1_coeff(3))) * t_tl
    delta2_tl = (delta2_coeff(1) + ivar%t*(two*delta2_coeff(2) + ivar%t*three*delta2_coeff(3))) * t_tl

    ! Compute the tangent-liner "infinite" permittivity term
    ! (No coeffs provided in ref. Taken from existing code)
    einf_tl = einf_coeff(1) * t_tl


    ! Compute the tangent-linear permittivities
    ! using the double Debye model
    ! (eqn on pg ACL 1-3 of Ellison et al. 2003)
    ! ------------------------------------------
    ! The tangent-linear denominators of the double Debye model
    d1_tl = two*ivar%f2*ivar%tau1*tau1_tl
    d2_tl = two*ivar%f2*ivar%tau2*tau2_tl

    ! The tangent-linear real parts of the "delta" terms
    d12 = ivar%d1**2
    d22 = ivar%d2**2
    re1_tl = (ivar%d1*delta1_tl - ivar%delta1*d1_tl) / d12
    re2_tl = (ivar%d2*delta2_tl - ivar%delta2*d2_tl) / d22

    ! The tangent-linear imaginary parts of the "delta" terms
    ie1_tl = ivar%f * (delta1_tl*ivar%tau1*ivar%d1 + ivar%delta1*tau1_tl*ivar%d1 - ivar%delta1*ivar%tau1*d1_tl) / d12
    ie2_tl = ivar%f * (delta2_tl*ivar%tau2*ivar%d2 + ivar%delta2*tau2_tl*ivar%d2 - ivar%delta2*ivar%tau2*d2_tl) / d22

    ! The conductivity term
    sigma_tl   = sigma_coeff(1)*t_tl
    iesigma_tl = sigma_tl / ivar%f0

    ! Construct the complex permittivity, de = de' - j.de"
    re_tl = re1_tl + re2_tl + einf_tl
    ie_tl = ie1_tl + ie2_tl + iesigma_tl
    permittivity_tl = cmplx(re_tl, -ie_tl, fp)

  END SUBROUTINE ellison_ocean_permittivity_tl


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       Ellison_Ocean_Permittivity_AD
!
! PURPOSE:
!       Subroutine to compute the adjoint ocean permittivity according
!       to the reference,
!         Ellison, W.J. et al. (2003) A comparison of ocean emissivity models
!           using the Advanced Microwave Sounding Unit, the Special Sensor
!           Microwave Imager, the TRMM Microwave Imager, and airborne radiometer
!           observations. Journal of Geophysical Research, v108, D21, ACL 1,1-14
!           doi:10.1029/2002JD0032132
!
! CALLING SEQUENCE:
!       CALL Ellison_Ocean_Permittivity_AD( Permittivity_AD, & ! Input
!                                           Temperature_AD , & ! Output
!                                           iVar             ) ! Internal variable input
!
! INPUTS:
!       Permittivity_AD:  Adjoint ocean permittivity
!                         UNITS:      N/A
!                         TYPE:       COMPLEX(fp)
!                         DIMENSION:  Scalar
!                         ATTRIBUTES: INTENT(IN OUT)
!
!       iVar:             Structure containing internal variables required for
!                         subsequent tangent-linear or adjoint model calls.
!                         The contents of this structure are NOT accessible
!                         outside of this module.
!                         UNITS:      N/A
!                         TYPE:       TYPE(iVar_type)
!                         DIMENSION:  Scalar
!                         ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Temperature_AD:   Adjoint sea surface temperature, de/dT.
!                         UNITS:      per Kelvin (K^-1)
!                         TYPE:       REAL(fp)
!                         DIMENSION:  Scalar
!                         ATTRIBUTES: INTENT(IN OUT)
!
! SIDE EFFECTS:
!       The input adjoint variable, Permittivity_AD, is set to zero upon
!       exiting this routine.
!
! COMMENTS:
!       There is currently no salinity dependence.
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE ellison_ocean_permittivity_ad( &
    Permittivity_AD, & ! Input
    Temperature_AD , & ! Output
    iVar             ) ! Internal variable input
    ! Arguments
    COMPLEX(fp),     INTENT(IN OUT) :: permittivity_ad
    REAL(fp),        INTENT(IN OUT) :: temperature_ad
    TYPE(ivar_type), INTENT(IN)     :: ivar
    ! Local variables
    REAL(fp) :: t_ad
    REAL(fp) :: tau1_ad, tau2_ad
    REAL(fp) :: delta1_ad, delta2_ad, einf_ad
    REAL(fp) :: d1_ad, d2_ad
    REAL(fp) :: d12, d22
    REAL(fp) :: m1 , m2
    REAL(fp) :: re1_ad, re2_ad
    REAL(fp) :: ie1_ad, ie2_ad
    REAL(fp) :: sigma_ad, iesigma_ad
    REAL(fp) :: re_ad, ie_ad

    ! Compute the adjoint of the permittivities
    ! using the double Debye model
    ! (eqn on pg ACL 1-3 of Ellison et al. 2003)
    ! ------------------------------------------
    ! The complex permittivity
    ie_ad = -aimag(permittivity_ad)
    re_ad = REAL(permittivity_ad,fp)
    permittivity_ad = zero

    ! Initialise all the local adjoint variables
    iesigma_ad = ie_ad; ie2_ad = ie_ad; ie1_ad = ie_ad
    einf_ad    = re_ad; re2_ad = re_ad; re1_ad = re_ad

    ! The adjoint of the conductivity term
    sigma_ad = iesigma_ad / ivar%f0
    t_ad = sigma_coeff(1)*sigma_ad

    ! The adjoints of the imaginary parts of the "delta" terms
    d22 = ivar%d2**2; m2 = ivar%f / d22
    d12 = ivar%d1**2; m1 = ivar%f / d12

    d2_ad     = -(ivar%delta2 * ivar%tau2 * m2 * ie2_ad )
    tau2_ad   =   ivar%delta2 * ivar%d2   * m2 * ie2_ad
    delta2_ad =   ivar%tau2   * ivar%d2   * m2 * ie2_ad

    d1_ad     = -(ivar%delta1 * ivar%tau1 * m1 * ie1_ad )
    tau1_ad   =   ivar%delta1 * ivar%d1   * m1 * ie1_ad
    delta1_ad =   ivar%tau1   * ivar%d1   * m1 * ie1_ad

    ! The adjoints of the real parts of the "delta" terms
    d2_ad     = d2_ad     - (re2_ad * ivar%delta2 / d22)
    delta2_ad = delta2_ad + (re2_ad / ivar%d2)

    d1_ad     = d1_ad     - (re1_ad * ivar%delta1 / d12)
    delta1_ad = delta1_ad + (re1_ad / ivar%d1)

    ! The adjoints of the denominators of the double Debye model
    tau2_ad = tau2_ad + (two * ivar%f2 * ivar%tau2 * d2_ad)
    tau1_ad = tau1_ad + (two * ivar%f2 * ivar%tau1 * d1_ad)


    ! Compute the adjoints of the various
    ! polynomial components of the double Debye model
    ! -----------------------------------------------
    ! Compute the adjoint of the "infinite" permittivity term
    ! (No coeffs provided in ref. Taken from existing code)
    t_ad = t_ad + einf_coeff(1) * einf_ad

    ! Compute the adjoint of the delta terms
    ! (eqn on pg ACL 1-4 of Ellison et al. 2003)
    t_ad = t_ad + (delta2_coeff(1) + ivar%t*(two*delta2_coeff(2) + &
                                       ivar%t*three*delta2_coeff(3)) ) * delta2_ad
    t_ad = t_ad + (delta1_coeff(1) + ivar%t*(two*delta1_coeff(2) + &
                                       ivar%t*three*delta1_coeff(3)) ) * delta1_ad

    ! Compute the adjoint of the  Debye model relaxation frequencies
    ! (eqn on pg ACL 1-4 of Ellison et al. 2003)
    t_ad = t_ad + (tau2_coeff(1) + ivar%t*(two*tau2_coeff(2) + &
                                     ivar%t*three*tau2_coeff(3)) ) * tau2_ad
    t_ad = t_ad + (tau1_coeff(1) + ivar%t*two*tau1_coeff(2)) * tau1_ad

    ! The return value
    temperature_ad = temperature_ad + t_ad

  END SUBROUTINE ellison_ocean_permittivity_ad

END MODULE ellison