Azimuth_Emissivity_F6_Module.f90

Go to the documentation of this file.

!
! Helper module containing the azimuth emissivity routines for the
! CRTM implementation of FASTEM6
!
!
! CREATION HISTORY:
!       Written by:     Original FASTEM1-5 authors, and Masahiro Kazumori
!                       for the new azimuthal emissivity model.
!
!       Refactored by:  Paul van Delst, January 2015
!                       paul.vandelst@noaa.gov
!

MODULE azimuth_emissivity_f6_module

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE type_kinds     , ONLY: fp
  USE fitcoeff_define, ONLY: fitcoeff_3d_type
  USE search_utility , ONLY: bisection_search
  ! Disable implicit typing
  IMPLICIT NONE


  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  ! Data types
  PUBLIC :: ivar_type
  ! Science routines
  PUBLIC :: azimuth_emissivity_f6
  PUBLIC :: azimuth_emissivity_f6_tl
  PUBLIC :: azimuth_emissivity_f6_ad


  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: module_version_id = &
  '$Id: Azimuth_Emissivity_F6_Module.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 :: pi = 3.141592653589793238462643383279_fp
  REAL(fp), PARAMETER :: degrees_to_radians = pi / 180.0_fp

  ! Dimensions
  ! ...Number of Stokes parameters handled
  INTEGER, PARAMETER :: n_stokes = 2
  INTEGER, PARAMETER :: ivpol = 1  ! Index for vertical polarisation
  INTEGER, PARAMETER :: ihpol = 2  ! Index for horizontal polarisation
  ! ...The number of fitting frequencies
  INTEGER, PARAMETER :: n_frequencies = 6

  ! Parameters used in the model
  ! ...The fitting frequencies
  REAL(fp), PARAMETER :: fit_frequency(n_frequencies) = &
  [ 6.925_fp, 10.65_fp, 18.7_fp, 23.8_fp, 36.5_fp, 89.0_fp ]
  ! ...Wind speed limits
  REAL(fp), PARAMETER :: wind_speed_max18 = 18.0_fp
  REAL(fp), PARAMETER :: wind_speed_max15 = 15.0_fp
  ! ...Frequency limits
  REAL(fp), PARAMETER :: frequency_max37 = 37.0_fp
  ! ...Exponents for the AxSy_theta terms
  REAL(fp), PARAMETER :: xs11 = two
  REAL(fp), PARAMETER :: xs12 = two
  REAL(fp), PARAMETER :: xs21 = one
  REAL(fp), PARAMETER :: xs22 = four
  ! ...Reference zenith angle
  REAL(fp), PARAMETER :: theta_ref = 55.2_fp


  ! --------------------------------------
  ! Structure definition to hold internal
  ! variables across FWD, TL, and AD calls
  ! --------------------------------------
  TYPE :: ivar_type
    PRIVATE
    ! Direct inputs
    REAL(fp) :: wind_speed   = zero
    REAL(fp) :: frequency    = zero
    REAL(fp) :: zenith_angle = zero
    ! Derived inputs
    REAL(fp) :: phi  = zero    ! Azimuth angle in radians
    LOGICAL  :: lw18 = .false. ! Logical to flag wind speed > 18m/s
    REAL(fp) :: w18  = zero    ! Wind speed with 18m/s maximum
    LOGICAL  :: lw15 = .false. ! Logical to flag wind speed > 15m/s
    REAL(fp) :: w15  = zero    ! Wind speed with 15m/s maximum
    REAL(fp) :: f37  = zero    ! Frequency with 37GHz maximum
    ! Intermediate variables
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1v , a2v , a1h , a2h
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1s1, a1s2, a2s1, a2s2
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a2s2_theta0
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1s1_theta, a1s2_theta, a2s1_theta, a2s2_theta
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1v_theta , a1h_theta , a2v_theta , a2h_theta
    REAL(fp) :: azimuth_component(n_frequencies, n_stokes)
    ! Interpolation variables
    INTEGER :: i1 = 0
    INTEGER :: i2 = 0
    REAL(fp) :: lpoly = zero
   END TYPE ivar_type


CONTAINS


  ! ===========================================================
  ! Compute emissivity as a function of relative azimuth angle.
  ! ===========================================================

  ! Forward model
  SUBROUTINE azimuth_emissivity_f6( &
    AZCoeff      , &  ! Input
    Wind_Speed   , &  ! Input
    Azimuth_Angle, &  ! Input
    Frequency    , &  ! Input
    Zenith_Angle , &  ! Input
    e_Azimuth    , &  ! Output
    iVar           )  ! Internal variable output
    ! Arguments
    TYPE(fitcoeff_3d_type), INTENT(IN)     :: azcoeff
    REAL(fp)              , INTENT(IN)     :: wind_speed
    REAL(fp)              , INTENT(IN)     :: azimuth_angle
    REAL(fp)              , INTENT(IN)     :: frequency
    REAL(fp)              , INTENT(IN)     :: zenith_angle
    REAL(fp)              , INTENT(OUT)    :: e_azimuth(:)
    TYPE(ivar_type)       , INTENT(IN OUT) :: ivar
    ! Local variables
    INTEGER :: j

    ! Initialise output
    e_azimuth = zero

    ! Save inputs for TL and AD calls
    ivar%wind_speed   = wind_speed
    ivar%frequency    = frequency
    ivar%zenith_angle = zenith_angle

    ! Convert inputs
    ivar%phi  = azimuth_angle * degrees_to_radians
    ivar%lw18 = wind_speed > wind_speed_max18
    ivar%w18  = min(wind_speed, wind_speed_max18)
    ivar%lw15 = wind_speed > wind_speed_max15
    ivar%w15  = min(wind_speed, wind_speed_max15)
    ivar%f37  = min(frequency, frequency_max37)


    ! Loop over frequencies to compute the intermediate terms
    frequency_loop: DO j = 1, n_frequencies

      ivar%A1v(j) = azcoeff%C(1,j,ivpol) * ( exp(-azcoeff%C(5,j,ivpol) * ivar%w18**2 ) - one ) * &
                    ( azcoeff%C(2,j,ivpol) * ivar%w18    + &
                      azcoeff%C(3,j,ivpol) * ivar%w18**2 + &
                      azcoeff%C(4,j,ivpol) * ivar%w18**3   )
      ivar%A2v(j) = azcoeff%C(6,j,ivpol) * ivar%w18

      ivar%A1h(j) = azcoeff%C(1,j,ihpol) * ivar%w18
      ivar%A2h(j) = azcoeff%C(2,j,ihpol) * ( exp(-azcoeff%C(6,j,ihpol) * ivar%w18**2 ) - one ) * &
                    ( azcoeff%C(3,j,ihpol) * ivar%w18    + &
                      azcoeff%C(4,j,ihpol) * ivar%w18**2 + &
                      azcoeff%C(5,j,ihpol) * ivar%w18**3   )

      ivar%A1s1(j) = (ivar%A1v(j) + ivar%A1h(j))/two
      ivar%A1s2(j) =  ivar%A1v(j) - ivar%A1h(j)
      ivar%A2s1(j) = (ivar%A2v(j) + ivar%A2h(j))/two
      ivar%A2s2(j) =  ivar%A2v(j) - ivar%A2h(j)

      ivar%A2s2_theta0(j) = (ivar%w15**2 - (ivar%w15**3)/22.5_fp)/55.5556_fp * &
                            (two/290.0_fp) * &
                            (one - log10(30.0_fp/ivar%f37) )

      ivar%A1s1_theta(j) = ivar%A1s1(j)*((ivar%zenith_angle/theta_ref)**xs11)
      ivar%A2s1_theta(j) = ivar%A2s1(j)*((ivar%zenith_angle/theta_ref)**xs12)
      ivar%A1s2_theta(j) = ivar%A1s2(j)*((ivar%zenith_angle/theta_ref)**xs21)
      ivar%A2s2_theta(j) = ivar%A2s2_theta0(j) + &
                           (ivar%A2s2(j) - ivar%A2s2_theta0(j))*((ivar%zenith_angle/theta_ref)**xs22)

      ivar%A1v_theta(j) = point5*(two*ivar%A1s1_theta(j) + ivar%A1s2_theta(j))
      ivar%A1h_theta(j) = point5*(two*ivar%A1s1_theta(j) - ivar%A1s2_theta(j))
      ivar%A2v_theta(j) = point5*(two*ivar%A2s1_theta(j) + ivar%A2s2_theta(j))
      ivar%A2h_theta(j) = point5*(two*ivar%A2s1_theta(j) - ivar%A2s2_theta(j))

      ivar%azimuth_component(j,ivpol) = (ivar%A1v_theta(j) * cos(ivar%phi)) + (ivar%A2v_theta(j) * cos(two*ivar%phi))
      ivar%azimuth_component(j,ihpol) = (ivar%A1h_theta(j) * cos(ivar%phi)) + (ivar%A2h_theta(j) * cos(two*ivar%phi))

    END DO frequency_loop


    ! Interpolate to input frequency for result. Only V and H polarisation.
    ! ...Check for lower out of bounds frequency
    IF ( frequency < fit_frequency(1) ) THEN
      e_azimuth(ivpol) = ivar%azimuth_component(1,ivpol)
      e_azimuth(ihpol) = ivar%azimuth_component(1,ihpol)
      RETURN
    END IF
    ! ...Check for upper out of bounds frequency
    IF ( frequency > fit_frequency(n_frequencies) ) THEN
      e_azimuth(ivpol) = ivar%azimuth_component(n_frequencies,ivpol)
      e_azimuth(ihpol) = ivar%azimuth_component(n_frequencies,ihpol)
      RETURN
    END IF
    ! ...In bounds, so interpolate
    ivar%i1 = bisection_search(fit_frequency, frequency)
    ivar%i2 = ivar%i1 + 1
    ivar%lpoly = (frequency - fit_frequency(ivar%i1))/(fit_frequency(ivar%i2)-fit_frequency(ivar%i1))

    e_azimuth(ivpol) = (   ivar%lpoly      * ivar%azimuth_component(ivar%i2,ivpol)) + &
                       ((one - ivar%lpoly) * ivar%azimuth_component(ivar%i1,ivpol))
    e_azimuth(ihpol) = (   ivar%lpoly      * ivar%azimuth_component(ivar%i2,ihpol)) + &
                       ((one - ivar%lpoly) * ivar%azimuth_component(ivar%i1,ihpol))

  END SUBROUTINE azimuth_emissivity_f6



  ! Tangent-linear model
  SUBROUTINE azimuth_emissivity_f6_tl( &
    AZCoeff         , &  ! Input
    Wind_Speed_TL   , &  ! Input
    Azimuth_Angle_TL, &  ! Input
    e_Azimuth_TL    , &  ! Output
    iVar              )  ! Internal variable input
    ! Arguments
    TYPE(fitcoeff_3d_type), INTENT(IN)  :: azcoeff
    REAL(fp)              , INTENT(IN)  :: wind_speed_tl
    REAL(fp)              , INTENT(IN)  :: azimuth_angle_tl
    REAL(fp)              , INTENT(OUT) :: e_azimuth_tl(:)
    TYPE(ivar_type)       , INTENT(IN)  :: ivar
    ! Local variables
    INTEGER  :: j
    REAL(fp) :: phi_tl
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1v_tl , a2v_tl , a1h_tl , a2h_tl
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1s1_tl, a1s2_tl, a2s1_tl, a2s2_tl
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a2s2_theta0_tl
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1s1_theta_tl, a1s2_theta_tl, a2s1_theta_tl, a2s2_theta_tl
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1v_theta_tl , a1h_theta_tl , a2v_theta_tl , a2h_theta_tl
    REAL(fp) :: azimuth_component_tl(n_frequencies, n_stokes)

    ! Initialise output
    e_azimuth_tl = zero

    ! Convert inputs
    phi_tl = azimuth_angle_tl * degrees_to_radians

    ! Loop over frequencies to compute the intermediate terms
    frequency_loop: DO j = 1, n_frequencies

      ! Only compute TL values if wind speed is not maxed out
      IF ( ivar%lw18 ) THEN
        a1v_tl(j) = zero
        a2v_tl(j) = zero
        a1h_tl(j) = zero
        a2h_tl(j) = zero
      ELSE
        a1v_tl(j) = ( azcoeff%C(1,j,ivpol) * ( exp(-azcoeff%C(5,j,ivpol) * ivar%w18**2 ) - one ) * &
                      (         azcoeff%C(2,j,ivpol)               + &
                        two   * azcoeff%C(3,j,ivpol) * ivar%w18    + &
                        three * azcoeff%C(4,j,ivpol) * ivar%w18**2   ) - &
                      two * azcoeff%C(1,j,ivpol) * azcoeff%C(5,j,ivpol) * ivar%w18 * &
                      exp(-azcoeff%C(5,j,ivpol) * ivar%w18**2 ) * &
                      ( azcoeff%C(2,j,ivpol) * ivar%w18    + &
                        azcoeff%C(3,j,ivpol) * ivar%w18**2 + &
                        azcoeff%C(4,j,ivpol) * ivar%w18**3   ) ) * wind_speed_tl
        a2v_tl(j) = azcoeff%C(6,j,ivpol) * wind_speed_tl

        a1h_tl(j) = azcoeff%C(1,j,ihpol) * wind_speed_tl

        a2h_tl(j) = ( azcoeff%C(2,j,ihpol) * ( exp(-azcoeff%C(6,j,ihpol) * ivar%w18**2 ) - one ) * &
                      (         azcoeff%C(3,j,ihpol)               + &
                        two   * azcoeff%C(4,j,ihpol) * ivar%w18    + &
                        three * azcoeff%C(5,j,ihpol) * ivar%w18**2   ) - &
                      two * azcoeff%C(2,j,ihpol) * azcoeff%C(6,j,ihpol) * ivar%w18 * &
                      exp(-azcoeff%C(6,j,ihpol) * ivar%w18**2 ) * &
                      ( azcoeff%C(3,j,ihpol) * ivar%w18    + &
                        azcoeff%C(4,j,ihpol) * ivar%w18**2 + &
                        azcoeff%C(5,j,ihpol) * ivar%w18**3   ) ) * wind_speed_tl

      END IF

      a1s1_tl(j) = (a1v_tl(j) + a1h_tl(j))/two
      a1s2_tl(j) =  a1v_tl(j) - a1h_tl(j)
      a2s1_tl(j) = (a2v_tl(j) + a2h_tl(j))/two
      a2s2_tl(j) =  a2v_tl(j) - a2h_tl(j)


      ! Only compute TL value if wind speed is not maxed out
      IF ( ivar%lw15 ) THEN
        a2s2_theta0_tl(j) = zero
      ELSE
        a2s2_theta0_tl(j) = (two*ivar%w15 - (three*ivar%w15**2)/22.5_fp)/55.5556_fp * &
                            (two/290.0_fp) * &
                            (one - log10(30.0_fp/ivar%f37) ) * wind_speed_tl
      END IF


      a1s1_theta_tl(j) = a1s1_tl(j)*((ivar%zenith_angle/theta_ref)**xs11)
      a2s1_theta_tl(j) = a2s1_tl(j)*((ivar%zenith_angle/theta_ref)**xs12)
      a1s2_theta_tl(j) = a1s2_tl(j)*((ivar%zenith_angle/theta_ref)**xs21)
      a2s2_theta_tl(j) = a2s2_theta0_tl(j) + &
                         (a2s2_tl(j) - a2s2_theta0_tl(j))*((ivar%zenith_angle/theta_ref)**xs22)


      a1v_theta_tl(j) = point5*(two*a1s1_theta_tl(j) + a1s2_theta_tl(j))
      a1h_theta_tl(j) = point5*(two*a1s1_theta_tl(j) - a1s2_theta_tl(j))
      a2v_theta_tl(j) = point5*(two*a2s1_theta_tl(j) + a2s2_theta_tl(j))
      a2h_theta_tl(j) = point5*(two*a2s1_theta_tl(j) - a2s2_theta_tl(j))


      azimuth_component_tl(j,ivpol) = (cos(  ivar%phi  ) * a1v_theta_tl(j)) + &
                                      (cos(two*ivar%phi) * a2v_theta_tl(j)) - &
                                      ( (      ivar%A1v_theta(j) * sin(  ivar%phi  )) + &
                                        (two * ivar%A2v_theta(j) * sin(two*ivar%phi))   ) * phi_tl
      azimuth_component_tl(j,ihpol) = (cos(  ivar%phi  ) * a1h_theta_tl(j)) + &
                                      (cos(two*ivar%phi) * a2h_theta_tl(j)) - &
                                      ( (      ivar%A1h_theta(j) * sin(  ivar%phi  )) + &
                                        (two * ivar%A2h_theta(j) * sin(two*ivar%phi))   ) * phi_tl

    END DO frequency_loop


    ! Interpolate to input frequency for result. Only V and H polarisation.
    ! ...Check for lower out of bounds frequency
    IF ( ivar%frequency < fit_frequency(1) ) THEN
      e_azimuth_tl(ivpol) = azimuth_component_tl(1,ivpol)
      e_azimuth_tl(ihpol) = azimuth_component_tl(1,ihpol)
      RETURN
    END IF
    ! ...Check for upper out of bounds frequency
    IF ( ivar%frequency > fit_frequency(n_frequencies) ) THEN
      e_azimuth_tl(ivpol) = azimuth_component_tl(n_frequencies,ivpol)
      e_azimuth_tl(ihpol) = azimuth_component_tl(n_frequencies,ihpol)
      RETURN
    END IF
    ! ...In bounds, so interpolate
    e_azimuth_tl(ivpol) = (   ivar%lpoly      * azimuth_component_tl(ivar%i2,ivpol)) + &
                          ((one - ivar%lpoly) * azimuth_component_tl(ivar%i1,ivpol))
    e_azimuth_tl(ihpol) = (   ivar%lpoly      * azimuth_component_tl(ivar%i2,ihpol)) + &
                          ((one - ivar%lpoly) * azimuth_component_tl(ivar%i1,ihpol))

  END SUBROUTINE azimuth_emissivity_f6_tl


  ! Adjoint model
  SUBROUTINE azimuth_emissivity_f6_ad( &
    AZCoeff         , &  ! Input
    e_Azimuth_AD    , &  ! AD Input
    Wind_Speed_AD   , &  ! AD Output
    Azimuth_Angle_AD, &  ! AD Output
    iVar              )  ! Internal variable input
    ! Arguments
    TYPE(fitcoeff_3d_type), INTENT(IN)     :: azcoeff
    REAL(fp)              , INTENT(IN OUT) :: e_azimuth_ad(:)
    REAL(fp)              , INTENT(IN OUT) :: wind_speed_ad
    REAL(fp)              , INTENT(IN OUT) :: azimuth_angle_ad
    TYPE(ivar_type)       , INTENT(IN)     :: ivar
    ! Local variables
    INTEGER  :: j
    REAL(fp) :: phi_ad
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1v_ad , a2v_ad , a1h_ad , a2h_ad
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1s1_ad, a1s2_ad, a2s1_ad, a2s2_ad
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a2s2_theta0_ad
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1s1_theta_ad, a1s2_theta_ad, a2s1_theta_ad, a2s2_theta_ad
    REAL(fp), DIMENSION(N_FREQUENCIES) :: a1v_theta_ad , a1h_theta_ad , a2v_theta_ad , a2h_theta_ad
    REAL(fp) :: azimuth_component_ad(n_frequencies, n_stokes)

    ! Initialise local adjoint variables
    phi_ad = zero
    a1v_ad  = zero;  a2v_ad  = zero;  a1h_ad  = zero;  a2h_ad  = zero
    a1s1_ad = zero;  a1s2_ad = zero;  a2s1_ad = zero;  a2s2_ad = zero
    a2s2_theta0_ad = zero
    a1s1_theta_ad = zero; a1s2_theta_ad = zero; a2s1_theta_ad = zero; a2s2_theta_ad = zero
    a1v_theta_ad  = zero; a1h_theta_ad  = zero; a2v_theta_ad  = zero; a2h_theta_ad  = zero
    azimuth_component_ad = zero


    ! Adjoint of frequency interpolation. Only V and H polarisation.
    IF ( ivar%frequency < fit_frequency(1) ) THEN
      ! ...Lower out of bounds frequency
      azimuth_component_ad(1,ihpol) = azimuth_component_ad(1,ihpol) + e_azimuth_ad(ihpol)
      azimuth_component_ad(1,ivpol) = azimuth_component_ad(1,ivpol) + e_azimuth_ad(ivpol)
    ELSE IF ( ivar%frequency > fit_frequency(n_frequencies) ) THEN
      ! ...Upper out of bounds frequency
      azimuth_component_ad(n_frequencies,ihpol) = azimuth_component_ad(n_frequencies,ihpol) + e_azimuth_ad(ihpol)
      azimuth_component_ad(n_frequencies,ivpol) = azimuth_component_ad(n_frequencies,ivpol) + e_azimuth_ad(ivpol)
    ELSE
      ! ...In bounds, so interpolate
      azimuth_component_ad(ivar%i1,ihpol) = ((one - ivar%lpoly) * e_azimuth_ad(ihpol)) + azimuth_component_ad(ivar%i1,ihpol)
      azimuth_component_ad(ivar%i2,ihpol) = (   ivar%lpoly      * e_azimuth_ad(ihpol)) + azimuth_component_ad(ivar%i2,ihpol)
      azimuth_component_ad(ivar%i1,ivpol) = ((one - ivar%lpoly) * e_azimuth_ad(ivpol)) + azimuth_component_ad(ivar%i1,ivpol)
      azimuth_component_ad(ivar%i2,ivpol) = (   ivar%lpoly      * e_azimuth_ad(ivpol)) + azimuth_component_ad(ivar%i2,ivpol)
    END IF
    e_azimuth_ad(ihpol) = zero
    e_azimuth_ad(ivpol) = zero


    ! Loop over frequencies to compute the intermediate term adjoints
    frequency_loop: DO j = 1, n_frequencies

      phi_ad = phi_ad - &
               ( (      ivar%A1h_theta(j) * sin(  ivar%phi  )) + &
                 (two * ivar%A2h_theta(j) * sin(two*ivar%phi))   ) * azimuth_component_ad(j,ihpol)
      a2h_theta_ad(j) = (cos(two*ivar%phi) * azimuth_component_ad(j,ihpol)) + a2h_theta_ad(j)
      a1h_theta_ad(j) = (cos(  ivar%phi  ) * azimuth_component_ad(j,ihpol)) + a1h_theta_ad(j)
      azimuth_component_ad(j,ihpol) = zero

      phi_ad = phi_ad - &
               ( (      ivar%A1v_theta(j) * sin(  ivar%phi  )) + &
                 (two * ivar%A2v_theta(j) * sin(two*ivar%phi))   ) * azimuth_component_ad(j,ivpol)
      a2v_theta_ad(j) = (cos(two*ivar%phi) * azimuth_component_ad(j,ivpol)) + a2v_theta_ad(j)
      a1v_theta_ad(j) = (cos(  ivar%phi  ) * azimuth_component_ad(j,ivpol)) + a1v_theta_ad(j)
      azimuth_component_ad(j,ivpol) = zero


      a2s2_theta_ad(j) = a2s2_theta_ad(j) - point5*a2h_theta_ad(j)
      a2s1_theta_ad(j) = a2s1_theta_ad(j) +        a2h_theta_ad(j)
      a2h_theta_ad(j)  = zero

      a2s2_theta_ad(j) = a2s2_theta_ad(j) + point5*a2v_theta_ad(j)
      a2s1_theta_ad(j) = a2s1_theta_ad(j) +        a2v_theta_ad(j)
      a2v_theta_ad(j)  = zero

      a1s2_theta_ad(j) = a1s2_theta_ad(j) - point5*a1h_theta_ad(j)
      a1s1_theta_ad(j) = a1s1_theta_ad(j) +        a1h_theta_ad(j)
      a1h_theta_ad(j)  = zero

      a1s2_theta_ad(j) = a1s2_theta_ad(j) + point5*a1v_theta_ad(j)
      a1s1_theta_ad(j) = a1s1_theta_ad(j) +        a1v_theta_ad(j)
      a1v_theta_ad(j)  = zero


      a2s2_ad(j)        = a2s2_ad(j)        + a2s2_theta_ad(j)*((ivar%zenith_angle/theta_ref)**xs22)
      a2s2_theta0_ad(j) = a2s2_theta0_ad(j) + a2s2_theta_ad(j)*(one - (ivar%zenith_angle/theta_ref)**xs22)
      a2s2_theta_ad(j)  = zero

      a1s2_ad(j)       = a1s2_ad(j) + a1s2_theta_ad(j)*((ivar%zenith_angle/theta_ref)**xs21)
      a1s2_theta_ad(j) = zero

      a2s1_ad(j)       = a2s1_ad(j) + a2s1_theta_ad(j)*((ivar%zenith_angle/theta_ref)**xs12)
      a2s1_theta_ad(j) = zero

      a1s1_ad(j)       = a1s1_ad(j) + a1s1_theta_ad(j)*((ivar%zenith_angle/theta_ref)**xs11)
      a1s1_theta_ad(j) = zero


      ! Only compute AD value if wind speed is not maxed out
      IF ( ivar%lw15 ) THEN
        a2s2_theta0_ad(j) = zero
      ELSE
        wind_speed_ad = wind_speed_ad + &
                        ((two*ivar%w15 - (three*ivar%w15**2)/22.5_fp)/55.5556_fp * &
                         (two/290.0_fp) * &
                         (one - log10(30.0_fp/ivar%f37)))*a2s2_theta0_ad(j)
        a2s2_theta0_ad(j) = zero
      END IF


      a2h_ad(j)  = a2h_ad(j) - a2s2_ad(j)
      a2v_ad(j)  = a2v_ad(j) + a2s2_ad(j)
      a2s2_ad(j) = zero

      a2h_ad(j)  = a2h_ad(j) + point5*a2s1_ad(j)
      a2v_ad(j)  = a2v_ad(j) + point5*a2s1_ad(j)
      a2s1_ad(j) = zero

      a1h_ad(j)  = a1h_ad(j) - a1s2_ad(j)
      a1v_ad(j)  = a1v_ad(j) + a1s2_ad(j)
      a1s2_ad(j) = zero

      a1h_ad(j)  = a1h_ad(j) + point5*a1s1_ad(j)
      a1v_ad(j)  = a1v_ad(j) + point5*a1s1_ad(j)
      a1s1_ad(j) = zero
      
      
      ! Only compute AD values if wind speed is not maxed out
      IF ( ivar%lw18 ) THEN
        a1v_ad(j) = zero
        a2v_ad(j) = zero
        a1h_ad(j) = zero
        a2h_ad(j) = zero
      ELSE
        wind_speed_ad = wind_speed_ad + &
                        ( azcoeff%C(2,j,ihpol) * ( exp(-azcoeff%C(6,j,ihpol) * ivar%w18**2 ) - one ) * &
                          (         azcoeff%C(3,j,ihpol)               + &
                            two   * azcoeff%C(4,j,ihpol) * ivar%w18    + &
                            three * azcoeff%C(5,j,ihpol) * ivar%w18**2   ) - &
                          two * azcoeff%C(2,j,ihpol) * azcoeff%C(6,j,ihpol) * ivar%w18 * &
                          exp(-azcoeff%C(6,j,ihpol) * ivar%w18**2 ) * &
                          ( azcoeff%C(3,j,ihpol) * ivar%w18    + &
                            azcoeff%C(4,j,ihpol) * ivar%w18**2 + &
                            azcoeff%C(5,j,ihpol) * ivar%w18**3   ) ) * a2h_ad(j)
        a2h_ad(j) = zero
        
        wind_speed_ad = wind_speed_ad + azcoeff%C(1,j,ihpol)*a1h_ad(j) 
        a1h_ad(j) = zero
      
        wind_speed_ad = wind_speed_ad + azcoeff%C(6,j,ivpol)*a2v_ad(j) 
        a2v_ad(j) = zero
      
        wind_speed_ad = wind_speed_ad + &
                        ( azcoeff%C(1,j,ivpol) * ( exp(-azcoeff%C(5,j,ivpol) * ivar%w18**2 ) - one ) * &
                          (         azcoeff%C(2,j,ivpol)               + &
                            two   * azcoeff%C(3,j,ivpol) * ivar%w18    + &
                            three * azcoeff%C(4,j,ivpol) * ivar%w18**2   ) - &
                          two * azcoeff%C(1,j,ivpol) * azcoeff%C(5,j,ivpol) * ivar%w18 * &
                          exp(-azcoeff%C(5,j,ivpol) * ivar%w18**2 ) * &
                          ( azcoeff%C(2,j,ivpol) * ivar%w18    + &
                            azcoeff%C(3,j,ivpol) * ivar%w18**2 + &
                            azcoeff%C(4,j,ivpol) * ivar%w18**3   ) ) * a1v_ad(j)
        a1v_ad(j) = zero
      END IF

    END DO frequency_loop

    
    ! Adjoint of the angle perturbation in radians
    azimuth_angle_ad = azimuth_angle_ad + degrees_to_radians*phi_ad
    
  END SUBROUTINE azimuth_emissivity_f6_ad


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


END MODULE azimuth_emissivity_f6_module