Liu.f90

Go to the documentation of this file.

!
! Liu Ocean Permittivity module.
!
! Module containing routines to compute the complex permittivities for
! sea water based on
!
!   Liu, Q. et al. (2010) An improved fast microwave water emissivity model.
!      IEEE Trans. Geosci. Remote Sensing, accepted June 25, 2010
!
!
! CREATION HISTORY:
!       Written by:     Quanhua (Mark) Liu, JCSDA 30-Jul-2009
!                       Quanhua.Liu@noaa.gov

MODULE liu

  ! -----------------
  ! 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 :: liu_ocean_permittivity
  PUBLIC :: liu_ocean_permittivity_tl
  PUBLIC :: liu_ocean_permittivity_ad


  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: module_version_id = &
  '$Id: Liu.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'

  ! Literal constants
  ! -----------------
  REAL(fp), PARAMETER :: zero   = 0.0_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

  ! Scaling factors
  ! ---------------
  REAL(fp), PARAMETER :: ghz_to_hz = 1.0e+09_fp ! Gigahertz   -> Hertz

  ! Fixed value for ionic conductivity denominator term, scaled
  ! for frequency. See the last term of eqn.(3) in reference.
  ! -----------------------------------------------------------
  REAL(fp), PARAMETER :: tau0 = two*pi*e0*ghz_to_hz
  
  ! Parameters for the Liu et al (2010) permittivity model
  ! ------------------------------------------------------
  ! The coefficients for the high-frequency permittivity temperature
  ! polynomial. Eqn.(4a) in reference.
  REAL(fp), PARAMETER :: einf_coeff(0:1) = (/ 3.8_fp, &
                                              2.48033e-02_fp /)
  
  ! The coefficients for the static permittivity temperature
  ! and salinity polynomials. Eqn.(4b) in reference.
  REAL(fp), PARAMETER :: es_t_coeff(0:3) = (/ 87.9181727_fp      , &
                                              -4.031592248e-01_fp, &
                                               9.493088010e-04_fp, &
                                              -1.930858348e-06_fp /)
  REAL(fp), PARAMETER :: es_s_coeff(0:2) = (/-2.697e-03_fp, &
                                             -7.3e-06_fp  , &
                                             -8.9e-06_fp   /)

  ! The coefficients for the intermediate frequency permittivity
  ! temperature and salinity polynomials. Eqn.(4c) in reference.
  REAL(fp), PARAMETER :: e1_t_coeff(0:2) = (/ 5.723_fp     , &
                                              2.2379e-02_fp, &
                                             -7.1237e-04_fp /)
  REAL(fp), PARAMETER :: e1_s_coeff(0:2) = (/-6.28908e-03_fp, &
                                              1.76032e-04_fp, &
                                             -9.22144e-05_fp /)

  ! The coefficients for the relaxation time temperature and
  ! salinity polynomials. Eqns.(4e) and (4f) in reference.
  REAL(fp), PARAMETER :: tau1_t_coeff(0:3) = (/ 1.124465e-01_fp , &
                                               -3.9815727e-03_fp, &
                                                8.113381e-05_fp , &
                                               -7.1824242e-07_fp /)
  REAL(fp), PARAMETER :: tau1_s_coeff(0:2) = (/-2.39357e-03_fp, &  
                                                3.1353e-05_fp , &  
                                               -2.52477e-07_fp /)  

  REAL(fp), PARAMETER :: tau2_t_coeff(0:3) = (/ 3.049979018e-03_fp, &
                                               -3.010041629e-05_fp, &
                                                4.811910733e-06_fp, &
                                               -4.259775841e-08_fp /)
  REAL(fp), PARAMETER :: tau2_s_coeff(0:2) = (/ 1.49e-01_fp, &
                                               -8.8e-04_fp , &
                                               -1.05e-04_fp /)

  ! The coefficients for the ionic conductivity exponential.
  ! Eqn.(4g) in reference.
  REAL(fp), PARAMETER :: alpha_coeff = -4.259775841e-08_fp

  ! The coefficients for the ionic conductivity exponent term
  ! polynomial. Eqn.(4i) in reference.
  REAL(fp), PARAMETER :: beta_coeff(0:5) = (/ 2.033e-02_fp, &
                                              1.266e-04_fp, &
                                              2.464e-06_fp, &
                                             -1.849e-05_fp, &
                                              2.551e-07_fp, &
                                             -2.551e-08_fp /)
     
  ! The coefficients for the ionic conductivity at 25C polynomial.
  ! Eqn.(4j) in reference.
  REAL(fp), PARAMETER :: alpha25_coeff(0:3) = (/ 1.82521e-01_fp, &
                                                -1.46192e-03_fp, &
                                                 2.09324e-05_fp, &
                                                -1.28205e-07_fp /)


  ! --------------------------------------
  ! Structure definition to hold forward
  ! variables across FWD, TL, and AD calls
  ! --------------------------------------
  TYPE :: ivar_type
    PRIVATE
    REAL(fp) :: t=zero                             ! Temperature in deg.C
    REAL(fp) :: s=zero                             ! Salinity
    REAL(fp) :: delta=zero, beta=zero              ! Ionic conductivity components
    REAL(fp) :: alpha25=zero, alpha=zero           ! Ionic conductivity terms
    REAL(fp) :: es_t=zero, es_s=zero               ! The temperature and salinity es terms
    REAL(fp) :: e1_t=zero, e1_s=zero               ! The temperature and salinity e1 terms
    REAL(fp) :: tau1_t=zero, tau1_s=zero           ! The temperature and salinity tau1 terms
    REAL(fp) :: tau2_t=zero, tau2_s=zero           ! The temperature and salinity tau2 terms
    REAL(fp) :: f1=zero, f2=zero                   ! The relaxation compound terms, f.tau
    REAL(fp) :: del1=zero, del2=zero               ! The permittivity differences
  END TYPE ivar_type


CONTAINS


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       Liu_Ocean_Permittivity
!
! PURPOSE:
!       Subroutine to compute ocean permittivity according to the reference,
!         Liu, Q. et al. (2010) An improved fast microwave water emissivity model.
!            IEEE Trans. Geosci. Remote Sensing, accepted June 25, 2010
!
! CALLING SEQUENCE:
!       CALL Liu_Ocean_Permittivity( Temperature , & ! Input
!                                    Salinity    , & ! 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)
!
!       Salinity:      Water salinity
!                      UNITS:      ppt (parts per thousand)
!                      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)
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE liu_ocean_permittivity( &
    Temperature , & ! Input
    Salinity    , & ! Input
    Frequency   , & ! Input
    Permittivity, & ! Output
    iVar          ) ! Internal variable output                   
    ! Arguments
    REAL(fp),        INTENT(IN)     :: temperature
    REAL(fp),        INTENT(IN)     :: salinity
    REAL(fp),        INTENT(IN)     :: frequency
    COMPLEX(fp),     INTENT(OUT)    :: permittivity
    TYPE(ivar_type), INTENT(IN OUT) :: ivar
    ! Local variables
    REAL(fp) :: einf
    REAL(fp) :: tau1, tau2, es, e1
    REAL(fp) :: re, ie

    ! Setup
    ! -----
    ! ...Initialise imaginary component of result
    ie = zero
    ! ...Save the inputs
    ivar%t = temperature - k_to_c
    ivar%s = salinity


    ! Compute the TEMPERATURE polynomial parameterisations
    ! ----------------------------------------------------
    ! ...The high-frequency permittivity temperature polynomial (eqn.4a)
    einf = einf_coeff(0) + ivar%t*einf_coeff(1)
    ! ...The static permittivity temperature polynomial (eqn.4b)
    es = es_t_coeff(0) + ivar%t*(es_t_coeff(1) + &
                           ivar%t*(es_t_coeff(2) + &
                             ivar%t*es_t_coeff(3)))
    ivar%es_t = es  ! Save it
    ! ...The intermediate frequency permittivity temperature polynomial (eqn.4c)
    e1 = e1_t_coeff(0) + ivar%t*(e1_t_coeff(1) + &
                           ivar%t*e1_t_coeff(2))
    ivar%e1_t = e1  ! Save it
    ! ...The Debye relaxation time constants temperature polynomials (eqns.4e & 4f)
    ! ...Units of tau: nanoseconds (for use with GHz frequencies)
    tau1 = tau1_t_coeff(0) + ivar%t*(tau1_t_coeff(1) + &
                               ivar%t*(tau1_t_coeff(2) + &
                                 ivar%t*tau1_t_coeff(3)))
    ivar%tau1_t = tau1  ! Save it
    tau2 = tau2_t_coeff(0) + ivar%t*(tau2_t_coeff(1) + &
                               ivar%t*(tau2_t_coeff(2) + &
                                 ivar%t*tau2_t_coeff(3)))
    ivar%tau2_t = tau2  ! Save it 


    ! Compute the SALINITY polynomial parameterisations
    ! -------------------------------------------------
    IF ( ivar%s > zero ) THEN
      ! ...The temperature difference from 25C (eqn.4h) used to compute ionic conductivity.
      ivar%delta = 25.0_fp - ivar%t
      ! ...The beta term (eqn.4i) used to compute ionic conductivity
      ivar%beta = beta_coeff(0) + ivar%delta*(beta_coeff(1) + &
                                    ivar%delta*beta_coeff(2)) + &
                  (beta_coeff(3) + ivar%delta*(beta_coeff(4) + &
                                     ivar%delta*beta_coeff(5)))*ivar%S
      ! ...The ionic conductivity at 25C (eqn.4j)
      ivar%alpha25 = ivar%s*(alpha25_coeff(0) + &
                       ivar%s*(alpha25_coeff(1) + &
                         ivar%s*(alpha25_coeff(2) + &
                           ivar%s*alpha25_coeff(3))))
      ! ...The ionic conductivity (eqn.4g)
      ivar%alpha = ivar%alpha25*exp(-ivar%delta*ivar%beta)
      !  ...The imaginary component dependent on ionic conductivity (eqn.3)
      ie = -ivar%alpha/(frequency*tau0)


      ! ...The static permittivity salinity polynomial (eqn.4b)
      ivar%es_s = one + ivar%s*(es_s_coeff(0) + ivar%s*es_s_coeff(1) + ivar%t*es_s_coeff(2))
      es = es * ivar%es_s 


      ! ...The intermediate frequency permittivity salinity polynomial (eqn.4c)
      ivar%e1_s = one + ivar%s*(e1_s_coeff(0) + ivar%s*e1_s_coeff(1) + ivar%t*e1_s_coeff(2))
      e1 = e1 * ivar%e1_s 


      ! ...The Debye relaxation time constants salinity polynomials (eqns.4e & 4f)
      ! ...Units of tau: nanoseconds (for use with GHz frequencies)
      ivar%tau1_s = one + ivar%s*(tau1_s_coeff(0) + ivar%t*(tau1_s_coeff(1) + &
                                                      ivar%t*tau1_s_coeff(2)))
      tau1 = tau1 * ivar%tau1_s
      ivar%tau2_s = one + ivar%s*(tau2_s_coeff(0) + ivar%t*tau2_s_coeff(1) + &
                                                    (ivar%s**2)*tau2_s_coeff(2))
      tau2 = tau2 * ivar%tau2_s

    END IF
    
    
    ! Compute the complex permittivity
    ! --------------------------------
    ! ...The compound terms
    ivar%f1 = frequency*tau1  ! Note there is no GHz->Hz conversion.
    ivar%f2 = frequency*tau2  ! That is embedded in the tau values.
    ivar%del1 = es - e1
    ivar%del2 = e1 - einf
    ! ...The real part
    re = einf + ivar%del1/(one + ivar%f1**2) + &
                ivar%del2/(one + ivar%f2**2)
    ! ...The imaginary part
    ie = -ie + ivar%del1*ivar%f1/(one + ivar%f1**2) + &
               ivar%del2*ivar%f2/(one + ivar%f2**2)
    ! ...Combine them (e = re - j.ie)
    permittivity = cmplx(re,-ie,fp)                                  

  END SUBROUTINE liu_ocean_permittivity


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       Liu_Ocean_Permittivity_TL
!
! PURPOSE:
!       Subroutine to compute the tangent-linear ocean permittivity according
!       to the model in the reference,
!         Liu, Q. et al. (2010) An improved fast microwave water emissivity model.
!            IEEE Trans. Geosci. Remote Sensing, accepted June 25, 2010
!
! CALLING SEQUENCE:
!       CALL Liu_Ocean_Permittivity_TL( Temperature_TL , &  ! TL Input
!                                       Salinity_TL    , &  ! TL Input
!                                       Frequency      , &  ! Invariant Input
!                                       Permittivity_TL, &  ! 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)
!
!       Salinity_TL:     Tangent-linear water salinity.
!                        UNITS:      ppt (parts per thousand)
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Frequency:       Frequency
!                        UNITS:      GHz
!                        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 permittivity
!                        UNITS:      N/A
!                        TYPE:       COMPLEX(fp)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(OUT)
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE liu_ocean_permittivity_tl( &
    Temperature_TL , &  ! TL Input
    Salinity_TL    , &  ! TL Input                             
    Frequency      , &  ! Invariant input
    Permittivity_TL, &  ! TL Output
    iVar             )  ! Internal variable input
    ! Arguments
    REAL(fp),        INTENT(IN)  :: temperature_tl
    REAL(fp),        INTENT(IN)  :: salinity_tl
    REAL(fp),        INTENT(IN)  :: frequency
    COMPLEX(fp),     INTENT(OUT) :: permittivity_tl
    TYPE(ivar_type), INTENT(IN)  :: ivar
 
    ! Local variables
    REAL(fp) :: einf_tl
    REAL(fp) :: es_s_tl, es_tl
    REAL(fp) :: e1_s_tl, e1_tl
    REAL(fp) :: tau1_s_tl, tau1_tl
    REAL(fp) :: tau2_s_tl, tau2_tl
    REAL(fp) :: delta_tl, beta_tl
    REAL(fp) :: alpha25_tl, alpha_tl
    REAL(fp) :: f1_tl, f2_tl, del1_tl, del2_tl
    REAL(fp) :: x1, x2
    REAL(fp) :: re_tl, ie_tl


    ! Setup
    ! -----
    ! ...Initialise imaginary component of result
    ie_tl = zero
    

    ! Compute the tangent-linear TEMPERATURE polynomial parameterisations
    ! -------------------------------------------------------------------
    ! ...The high-frequency permittivity temperature polynomial (eqn.4a)
    einf_tl = einf_coeff(1)*temperature_tl
    ! ...The static permittivity temperature polynomial (eqn.4b)
    es_tl = (es_t_coeff(1) + ivar%t*(two*es_t_coeff(2) + &
                               ivar%t*three*es_t_coeff(3))) * temperature_tl
    ! ...The intermediate frequency permittivity temperature polynomial (eqn.4c)
    e1_tl = (e1_t_coeff(1) + ivar%t*two*e1_t_coeff(2)) * temperature_tl
    ! ...The Debye relaxation time constants temperature polynomials (eqns.4e & 4f)
    ! ...Units of tau: nanoseconds (for use with GHz frequencies)
    tau1_tl = (tau1_t_coeff(1) + ivar%t*(two*tau1_t_coeff(2) + &
                                   ivar%t*three*tau1_t_coeff(3))) * temperature_tl
    tau2_tl = (tau2_t_coeff(1) + ivar%t*(two*tau2_t_coeff(2) + &
                                   ivar%t*three*tau2_t_coeff(3))) * temperature_tl


    ! Compute the SALINITY polynomial parameterisations
    ! -------------------------------------------------
    IF ( ivar%s > zero ) THEN
      ! ...The temperature difference from 25C (eqn.4h) used to compute ionic conductivity.
      delta_tl = -temperature_tl
      ! ...The beta term (eqn.4i) used to compute ionic conductivity
      beta_tl = (beta_coeff(1) + ivar%delta*two*beta_coeff(2) + &
                 ivar%s*(beta_coeff(4) + ivar%delta*two*beta_coeff(5)))*delta_tl + &
                (beta_coeff(3) + ivar%delta*(beta_coeff(4) + &
                                   ivar%delta*beta_coeff(5)))*salinity_tl
      ! ...The ionic conductivity at 25C (eqn.4j)
      alpha25_tl = (alpha25_coeff(0) + &
                     ivar%s*(two*alpha25_coeff(1) + &
                       ivar%s*(three*alpha25_coeff(2) + &
                         ivar%s*four*alpha25_coeff(3))))*salinity_tl
      ! ...The ionic conductivity (eqn.4g)
      alpha_tl = (alpha25_tl/ivar%alpha25 - ivar%delta*beta_tl - ivar%beta*delta_tl)*ivar%alpha
      !  ...The imaginary component dependent on ionic conductivity (eqn.3)
      ie_tl = -alpha_tl/(frequency*tau0)


      ! ...The static permittivity salinity polynomial (eqn.4b)
      es_s_tl = (es_s_coeff(0) + ivar%s*two*es_s_coeff(1) + ivar%t*es_s_coeff(2))*salinity_tl + &
                ivar%s*es_s_coeff(2)*temperature_tl
      es_tl = ivar%es_t*es_s_tl + ivar%es_s*es_tl


      ! ...The intermediate frequency permittivity salinity polynomial (eqn.4c)
      e1_s_tl = (e1_s_coeff(0) + ivar%s*two*e1_s_coeff(1) + ivar%t*e1_s_coeff(2))*salinity_tl + &
                ivar%s*e1_s_coeff(2)*temperature_tl
      e1_tl = ivar%e1_t*e1_s_tl + ivar%e1_s*e1_tl


      ! ...The Debye relaxation time constants salinity polynomials (eqns.4e & 4f)
      ! ...Units of tau: nanoseconds (for use with GHz frequencies)
      tau1_s_tl = (tau1_s_coeff(0) + ivar%t*(tau1_s_coeff(1) + &
                                       ivar%t*tau1_s_coeff(2)))*salinity_tl + &
                  ivar%s*(tau1_s_coeff(1) + ivar%t*two*tau1_s_coeff(2))*temperature_tl
      tau1_tl = ivar%tau1_t*tau1_s_tl + ivar%tau1_s*tau1_tl

      tau2_s_tl = (tau2_s_coeff(0) + ivar%t*tau2_s_coeff(1) + &
                                     (ivar%s**2)*three*tau2_s_coeff(2))*salinity_tl + &
                  ivar%s*tau2_s_coeff(1)*temperature_tl
      tau2_tl = ivar%tau2_t*tau2_s_tl + ivar%tau2_s*tau2_tl

    END IF


    ! Compute the complex permittivity
    ! --------------------------------
    ! ...The compound terms
    f1_tl = frequency*tau1_tl  ! Note there is no GHz->Hz conversion.
    f2_tl = frequency*tau2_tl  ! That is embedded in the tau values.
    del1_tl = es_tl - e1_tl
    del2_tl = e1_tl - einf_tl
    ! ...Calculate the denominator terms
    x1 = one + ivar%f1**2
    x2 = one + ivar%f2**2
    ! ...The real part
    re_tl = einf_tl + &
            (x1*del1_tl - ivar%del1*ivar%f1*two*f1_tl)/x1**2 + &
            (x2*del2_tl - ivar%del2*ivar%f2*two*f2_tl)/x2**2
    ! ...The imaginary part
    ie_tl = ivar%f1*del1_tl/x1 + &
            ivar%del1*(one - ivar%f1**2)*f1_tl/x1**2 + &
            ivar%f2*del2_tl/x2 + &
            ivar%del2*(one - ivar%f2**2)*f2_tl/x2**2 - &
            ie_tl
    ! ...Combine them (e = re - j.ie)
    permittivity_tl = cmplx(re_tl,-ie_tl,fp)                                     
 
  END SUBROUTINE liu_ocean_permittivity_tl


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       Liu_Ocean_Permittivity_AD
!
! PURPOSE:
!       Subroutine to compute the ocean permittivity adjoint according to the
!       model in the reference,
!         Liu, Q. et al. (2010) An improved fast microwave water emissivity model.
!            IEEE Trans. Geosci. Remote Sensing, accepted June 25, 2010
!
! CALLING SEQUENCE:
!       CALL Liu_Ocean_Permittivity_AD( Permittivity_AD, &  ! AD Input
!                                       Frequency      , &  ! Invariant Input
!                                       Temperature_AD , &  ! AD Output
!                                       Salinity_AD    , &  ! AD Output
!                                       iVar             )  ! Internal variable input
!
! INPUTS:
!       Permittivity_AD:  Adjoint ocean permittivity
!                         *** SET TO ZERO UPON EXIT ***
!                         UNITS:      N/A
!                         TYPE:       COMPLEX(fp)
!                         DIMENSION:  Scalar
!                         ATTRIBUTES: INTENT(IN OUT)
!
!       Frequency:        Frequency
!                         UNITS:      GHz
!                         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:
!       Temperature_AD:   Adjoint sea surface temperature, de/dT.
!                         *** MUST HAVE VALUE UPON ENTRY ***
!                         UNITS:      per Kelvin (K^-1)
!                         TYPE:       REAL(fp)
!                         DIMENSION:  Scalar
!                         ATTRIBUTES: INTENT(IN OUT)
!
!       Salinity_AD:      Adjoint water salinity, de/dS
!                         *** MUST HAVE VALUE UPON ENTRY ***
!                         UNITS:      per ppt (parts-per-thousand^-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.
!
!:sdoc-:
!--------------------------------------------------------------------------------

  SUBROUTINE liu_ocean_permittivity_ad( &
    Permittivity_AD, & ! AD Input
    Frequency      , & ! Invariant Input
    Temperature_AD , & ! AD Output
    Salinity_AD    , & ! AD Output 
    iVar             ) ! Internal variable input  
    IMPLICIT NONE                                                           
    ! Arguments
    COMPLEX(fp),     INTENT(IN OUT) :: permittivity_ad
    REAL(fp),        INTENT(IN)     :: frequency
    REAL(fp),        INTENT(IN OUT) :: temperature_ad
    REAL(fp),        INTENT(IN OUT) :: salinity_ad
    TYPE(ivar_type), INTENT(IN)     :: ivar
    ! Local variables
    REAL(fp) :: re_ad, ie_ad
    REAL(fp) :: x1, x2
    REAL(fp) :: f1_ad, f2_ad, del1_ad, del2_ad
    REAL(fp) :: delta_ad, beta_ad
    REAL(fp) :: alpha25_ad, alpha_ad
    REAL(fp) :: tau1_s_ad, tau1_ad
    REAL(fp) :: tau2_s_ad, tau2_ad
    REAL(fp) :: einf_ad
    REAL(fp) :: es_s_ad, es_ad
    REAL(fp) :: e1_s_ad, e1_ad


    ! Complex permittivity
    ! --------------------
    ! ...Separate the real and imaginary parts
    ie_ad = -aimag(permittivity_ad)
    re_ad = REAL(permittivity_ad, fp)
    permittivity_ad = zero
    ! ...Calculate the denominator terms
    x1 = one + ivar%f1**2
    x2 = one + ivar%f2**2
    ! ...The imaginary part
    f2_ad = ivar%del2*(one - ivar%f2**2)*ie_ad/x2**2
    del2_ad = ivar%f2*ie_ad/x2
    f1_ad = ivar%del1*(one - ivar%f1**2)*ie_ad/x1**2
    del1_ad = ivar%f1*ie_ad/x1
    ie_ad = -ie_ad
    ! ...The real part
    f2_ad = f2_ad - (ivar%del2*ivar%f2*two*re_ad)/x2**2
    del2_ad = del2_ad + re_ad/x2
    f1_ad = f1_ad - (ivar%del1*ivar%f1*two*re_ad)/x1**2
    del1_ad = del1_ad + re_ad/x1
    einf_ad = re_ad
    ! ...The compound terms
    einf_ad = einf_ad - del2_ad
    e1_ad = del2_ad
    e1_ad = e1_ad - del1_ad
    es_ad = del1_ad
    tau2_ad = frequency*f2_ad
    tau1_ad = frequency*f1_ad


    ! Compute the SALINITY polynomial parameterisations
    ! -------------------------------------------------
    IF ( ivar%s > zero ) THEN
      ! ...The Debye relaxation time constants salinity polynomials (eqns.4e & 4f)
      ! ...Units of tau: nanoseconds (for use with GHz frequencies)
      tau2_s_ad = ivar%tau2_t*tau2_ad
      tau2_ad   = ivar%tau2_s*tau2_ad
      salinity_ad = salinity_ad + &
                    (tau2_s_coeff(0) + ivar%t*tau2_s_coeff(1) + &
                                       (ivar%s**2)*three*tau2_s_coeff(2))*tau2_s_ad
      temperature_ad = temperature_ad + &
                       ivar%s*tau2_s_coeff(1)*tau2_s_ad
      
      tau1_s_ad = ivar%tau1_t*tau1_ad
      tau1_ad   = ivar%tau1_s*tau1_ad
      salinity_ad = salinity_ad + &
                    (tau1_s_coeff(0) + ivar%t*(tau1_s_coeff(1) + &
                                       ivar%t*tau1_s_coeff(2)))*tau1_s_ad
      temperature_ad = temperature_ad + &
                       ivar%s*(tau1_s_coeff(1) + ivar%t*two*tau1_s_coeff(2))*tau1_s_ad


      ! ...The intermediate frequency permittivity salinity polynomial (eqn.4c)
      e1_s_ad = ivar%e1_t*e1_ad
      e1_ad   = ivar%e1_s*e1_ad
      salinity_ad = salinity_ad + &
                    (e1_s_coeff(0) + ivar%s*two*e1_s_coeff(1) + ivar%t*e1_s_coeff(2))*e1_s_ad
      temperature_ad = temperature_ad + &
                       ivar%s*e1_s_coeff(2)*e1_s_ad


      ! ...The static permittivity salinity polynomial (eqn.4b)
      es_s_ad = ivar%es_t*es_ad
      es_ad   = ivar%es_s*es_ad
      salinity_ad = salinity_ad + &
                    (es_s_coeff(0) + ivar%s*two*es_s_coeff(1) + ivar%t*es_s_coeff(2))*es_s_ad
      temperature_ad = temperature_ad + &
                       ivar%s*es_s_coeff(2)*es_s_ad
            

      !  ...The imaginary component dependent on ionic conductivity (eqn.3)
      alpha_ad = -ie_ad/(frequency*tau0)
      ! ...The ionic conductivity (eqn.4g)
      alpha25_ad = ivar%alpha*alpha_ad/ivar%alpha25
      beta_ad  = -ivar%delta*ivar%alpha*alpha_ad
      delta_ad = -ivar%beta*ivar%alpha*alpha_ad
      ! ...The ionic conductivity at 25C (eqn.4j)
      salinity_ad = salinity_ad + &
                    (alpha25_coeff(0) + &
                     ivar%s*(two*alpha25_coeff(1) + &
                       ivar%s*(three*alpha25_coeff(2) + &
                         ivar%s*four*alpha25_coeff(3))))*alpha25_ad
      ! ...The beta term (eqn.4i) used to compute ionic conductivity
      salinity_ad = salinity_ad + &
                    (beta_coeff(3) + ivar%delta*(beta_coeff(4) + &
                                       ivar%delta*beta_coeff(5)))*beta_ad
      delta_ad = delta_ad + &
                 (beta_coeff(1) + ivar%delta*two*beta_coeff(2) + &
                 (beta_coeff(4) + ivar%delta*two*beta_coeff(5))*ivar%s)*beta_ad
      ! ...The temperature difference from 25C (eqn.4h) used to compute ionic conductivity.
      temperature_ad = temperature_ad - delta_ad
    END IF


    ! ...The Debye relaxation time constants temperature polynomials (eqns.4e & 4f)
    ! ...Units of tau: nanoseconds (for use with GHz frequencies)
    temperature_ad = temperature_ad + &
                     (tau2_t_coeff(1) + ivar%t*(two*tau2_t_coeff(2) + &
                                          ivar%t*three*tau2_t_coeff(3)))*tau2_ad
    temperature_ad = temperature_ad + &
                     (tau1_t_coeff(1) + ivar%t*(two*tau1_t_coeff(2) + &
                                          ivar%t*three*tau1_t_coeff(3)))*tau1_ad
    ! ...The intermediate frequency permittivity temperature polynomial (eqn.4c)
    temperature_ad = temperature_ad + &
                     (e1_t_coeff(1) + ivar%t*two*e1_t_coeff(2))*e1_ad
    ! ...The static permittivity salinity polynomial (eqn.4b)
    temperature_ad = temperature_ad + &
                     (es_t_coeff(1) + ivar%t*(two*es_t_coeff(2) + &
                                        ivar%t*three*es_t_coeff(3)))*es_ad
    ! ...The high-frequency permittivity temperature polynomial (eqn.4a)
    temperature_ad = temperature_ad + einf_coeff(1)*einf_ad

  END SUBROUTINE liu_ocean_permittivity_ad

END MODULE liu