Azimuth_Emissivity_Module.f90

Go to the documentation of this file.

!
! Helper module containing the azimuth emissivity routines for the
! CRTM implementation of FASTEM4 and FASTEM5
!
!
! CREATION HISTORY:
!       Written by:     Original FASTEM1/2/3 authors
!
!       Modified by:    Quanhua Liu, Quanhua.Liu@noaa.gov
!                       Stephen English, Stephen.English@metoffice.gov.uk
!                       July, 2009
!
!       Refactored by:  Paul van Delst, December 2011
!                       paul.vandelst@noaa.gov
!

MODULE azimuth_emissivity_module

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


  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  ! Data types
  PUBLIC :: ivar_type
  ! Science routines
  PUBLIC :: azimuth_emissivity
  PUBLIC :: azimuth_emissivity_tl
  PUBLIC :: azimuth_emissivity_ad


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

  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 :: pi = 3.141592653589793238462643383279_fp
  REAL(fp), PARAMETER :: degrees_to_radians = pi / 180.0_fp

  ! Dimensions
  ! ...Number of component predictors for harmonic coefficients
  INTEGER, PARAMETER :: n_predictors = 10
  ! ...Number of Stokes parameters
  INTEGER, PARAMETER :: n_stokes = 4
  ! ...The number of harmonics considered in the trignometric parameterisation
  INTEGER, PARAMETER :: n_harmonics = 3
  
  
  ! --------------------------------------
  ! 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
    ! Derived inputs
    REAL(fp) :: sec_z = zero
    ! Cosine and sine of the harmonic angle, m*phi
    REAL(fp) :: cos_angle(n_harmonics) = zero
    REAL(fp) :: sin_angle(n_harmonics) = zero
    ! Intermediate variables
    REAL(fp) :: trig_coeff(n_stokes, n_harmonics) = zero
   END TYPE ivar_type

  
CONTAINS


  ! ===========================================================
  ! Compute emissivity as a function of relative azimuth angle.
  ! ===========================================================
  
  ! Forward model
  SUBROUTINE azimuth_emissivity( &
    AZCoeff      , &  ! Input
    Wind_Speed   , &  ! Input
    Azimuth_Angle, &  ! Input
    Frequency    , &  ! Input
    cos_z        , &  ! 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)     :: cos_z        
    REAL(fp)              , INTENT(OUT)    :: e_azimuth(:)
    TYPE(ivar_type)       , INTENT(IN OUT) :: ivar
    ! Local variables
    INTEGER :: i, m
    REAL(fp) :: phi, angle
    REAL(fp) :: predictor(n_predictors)

    ! Initialise output
    e_azimuth = zero

    ! Save inputs for TL and AD calls
    ivar%wind_speed = wind_speed
    ivar%frequency  = frequency
    ivar%sec_z      = one/cos_z
    
    ! Convert angle
    phi = azimuth_angle * degrees_to_radians

    ! Compute the azimuth emissivity component predictors
    CALL compute_predictors( wind_speed, frequency, ivar%sec_z, predictor )

    ! Compute the azimuth emissivity vector
    harmonic_loop: DO m = 1, n_harmonics

      ! Compute the angles
      angle = REAL(m,fp) * phi
      ivar%cos_angle(m) = cos(angle)
      ivar%sin_angle(m) = sin(angle)

      ! Compute the coefficients
      DO i = 1, n_stokes
        CALL compute_coefficient( &
               azcoeff%C(:,i,m), &
               predictor, &
               ivar%trig_coeff(i,m) )
      END DO

      ! Compute the emissivities
      e_azimuth(1) = e_azimuth(1) + ivar%trig_coeff(1,m)*ivar%cos_angle(m) ! Vertical
      e_azimuth(2) = e_azimuth(2) + ivar%trig_coeff(2,m)*ivar%cos_angle(m) ! Horizontal
      e_azimuth(3) = e_azimuth(3) + ivar%trig_coeff(3,m)*ivar%sin_angle(m) ! +/- 45deg.
      e_azimuth(4) = e_azimuth(4) + ivar%trig_coeff(4,m)*ivar%sin_angle(m) ! Circular

    END DO harmonic_loop

    ! Apply frequency correction 
    e_azimuth = e_azimuth * azimuth_freq_correction(frequency)

  END SUBROUTINE azimuth_emissivity


  ! Tangent-linear model
  SUBROUTINE azimuth_emissivity_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 :: i, m
    REAL(fp) :: phi_tl, angle_tl
    REAL(fp) :: predictor_tl(n_predictors)
    REAL(fp) :: trig_coeff_tl(n_stokes)
    
    ! Initialise output
    e_azimuth_tl = zero

    ! Compute angle perturbation in radians    
    phi_tl = azimuth_angle_tl * degrees_to_radians
    
    ! Compute the azimuth emissivity component tangent-linear predictors
    CALL compute_predictors_tl( ivar%wind_speed, ivar%frequency, ivar%sec_z, wind_speed_tl, predictor_tl )

    ! Compute the tangent-linear azimuth emissivity vector
    harmonic_loop: DO m = 1, n_harmonics
      
      ! Compute the angles
      angle_tl = REAL(m,fp) * phi_tl

      ! Compute the coefficients
      DO i = 1, n_stokes
        CALL compute_coefficient_tl( &
               azcoeff%C(:,i,m), &
               predictor_tl, &
               trig_coeff_tl(i) )
      END DO
      
      ! Compute the tangent-linear emissivities
      ! ...Vertical polarisation
      e_azimuth_tl(1) = e_azimuth_tl(1) + ivar%cos_angle(m)*trig_coeff_tl(1) - &
                                          ivar%trig_coeff(1,m)*ivar%sin_angle(m)*angle_tl
      ! ...Horizontal polarisation
      e_azimuth_tl(2) = e_azimuth_tl(2) + ivar%cos_angle(m)*trig_coeff_tl(2) - &
                                          ivar%trig_coeff(2,m)*ivar%sin_angle(m)*angle_tl
      ! ...+/- 45deg. polarisation
      e_azimuth_tl(3) = e_azimuth_tl(3) + ivar%sin_angle(m)*trig_coeff_tl(3) + &
                                          ivar%trig_coeff(3,m)*ivar%cos_angle(m)*angle_tl 
      ! ...Circular polarisation
      e_azimuth_tl(4) = e_azimuth_tl(4) + ivar%sin_angle(m)*trig_coeff_tl(4) + &
                                          ivar%trig_coeff(4,m)*ivar%cos_angle(m)*angle_tl
    END DO harmonic_loop

    ! Apply tangent-linear frequency correction
    e_azimuth_tl = e_azimuth_tl * azimuth_freq_correction(ivar%frequency)

  END SUBROUTINE azimuth_emissivity_tl


  ! Adjoint model
  SUBROUTINE azimuth_emissivity_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 :: i, m
    REAL(fp) :: phi_ad, angle_ad
    REAL(fp) :: predictor_ad(n_predictors)
    REAL(fp) :: trig_coeff_ad(n_stokes)

    ! Initialise local adjoints
    phi_ad       = zero
    predictor_ad = zero
    
    ! Adjoint of frequency correction
    e_azimuth_ad = e_azimuth_ad * azimuth_freq_correction(ivar%frequency)
    
    ! Compute the azimuth emissivity vector adjoint
    harmonic_loop: DO m = 1, n_harmonics
      
      ! Compute the adjoint of the emissivities
      ! ...Circular polarisation
      angle_ad         = ivar%trig_coeff(4,m)*ivar%cos_angle(m)*e_azimuth_ad(4)
      trig_coeff_ad(4) = ivar%sin_angle(m)*e_azimuth_ad(4)
      ! ...+/- 45deg. polarisation
      angle_ad         = angle_ad + ivar%trig_coeff(3,m)*ivar%cos_angle(m)*e_azimuth_ad(3)
      trig_coeff_ad(3) = ivar%sin_angle(m)*e_azimuth_ad(3)
      ! ...Horizontal polarisation
      angle_ad         = angle_ad - ivar%trig_coeff(2,m)*ivar%sin_angle(m)*e_azimuth_ad(2)
      trig_coeff_ad(2) = ivar%cos_angle(m)*e_azimuth_ad(2)
      ! ...Vertical polarisation
      angle_ad         = angle_ad - ivar%trig_coeff(1,m)*ivar%sin_angle(m)*e_azimuth_ad(1)
      trig_coeff_ad(1) = ivar%cos_angle(m)*e_azimuth_ad(1)

      ! Compute the adjoint of the coefficients
      DO i = n_stokes, 1, -1
        CALL compute_coefficient_ad( &
               azcoeff%C(:,i,m), &
               trig_coeff_ad(i), &
               predictor_ad      )
      END DO

      ! Compute the angle adjoint
      phi_ad = phi_ad + REAL(m,fp)*angle_ad
      
    END DO harmonic_loop
    
    ! Compute the azimuth emissivity component predictor adjoint
    CALL compute_predictors_ad( ivar%wind_speed, ivar%frequency, ivar%sec_z, predictor_ad, wind_speed_ad )

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


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

  ! =============================================
  ! Compute predictors for the azimuth components
  ! =============================================
  
  ! Forward model
  SUBROUTINE compute_predictors( &
    Wind_Speed, &  ! Input
    Frequency , &  ! Input
    sec_z     , &  ! Input
    Predictor   )  ! Output
    ! Arguments
    REAL(fp), INTENT(IN)  :: Wind_Speed
    REAL(fp), INTENT(IN)  :: Frequency 
    REAL(fp), INTENT(IN)  :: sec_z     
    REAL(fp), INTENT(OUT) :: Predictor(N_PREDICTORS) 
    ! Compute the predictors.
    predictor( 1) = one
    predictor( 2) = frequency
    predictor( 3) = sec_z
    predictor( 4) = sec_z * frequency
    predictor( 5) = wind_speed
    predictor( 6) = wind_speed * frequency
    predictor( 7) = wind_speed**2
    predictor( 8) = frequency * wind_speed**2
    predictor( 9) = wind_speed * sec_z
    predictor(10) = wind_speed * sec_z * frequency
  END SUBROUTINE compute_predictors
    
  
  ! Tangent-linear model  
  SUBROUTINE compute_predictors_tl(  &
    Wind_Speed   , &  ! FWD Input
    Frequency    , &  ! FWD Input
    sec_z        , &  ! FWD Input
    Wind_Speed_TL, &  ! TL  Input
    Predictor_TL   )  ! TL  Output
    ! Arguments
    REAL(fp), INTENT(IN)  :: Wind_Speed
    REAL(fp), INTENT(IN)  :: Frequency 
    REAL(fp), INTENT(IN)  :: sec_z     
    REAL(fp), INTENT(IN)  :: Wind_Speed_TL
    REAL(fp), INTENT(OUT) :: Predictor_TL(N_PREDICTORS)
    ! Compute the tangent-linear predictors.
    predictor_tl( 1) = zero
    predictor_tl( 2) = zero
    predictor_tl( 3) = zero
    predictor_tl( 4) = zero
    predictor_tl( 5) = wind_speed_tl
    predictor_tl( 6) = frequency * wind_speed_tl
    predictor_tl( 7) = two * wind_speed * wind_speed_tl
    predictor_tl( 8) = two * frequency * wind_speed * wind_speed_tl
    predictor_tl( 9) = sec_z * wind_speed_tl
    predictor_tl(10) = sec_z * frequency * wind_speed_tl
  END SUBROUTINE compute_predictors_tl
    
  
  ! Adjoint model
  SUBROUTINE compute_predictors_ad(  &
    Wind_Speed   , &  ! FWD Input
    Frequency    , &  ! FWD Input
    sec_z        , &  ! FWD Input
    Predictor_AD , &  ! AD  Input
    Wind_Speed_AD  )  ! AD  Output
    ! Arguments
    REAL(fp), INTENT(IN)     :: Wind_Speed
    REAL(fp), INTENT(IN)     :: Frequency 
    REAL(fp), INTENT(IN)     :: sec_z     
    REAL(fp), INTENT(IN OUT) :: Predictor_AD(N_PREDICTORS)
    REAL(fp), INTENT(IN OUT) :: Wind_Speed_AD
    ! Compute the predictor adjoints
    wind_speed_ad = wind_speed_ad + &
                    sec_z * frequency            * predictor_ad(10) + &
                    sec_z                        * predictor_ad( 9) + &
                    two * frequency * wind_speed * predictor_ad( 8) + &
                    two             * wind_speed * predictor_ad( 7) + &
                          frequency              * predictor_ad( 6) + &
                                                   predictor_ad( 5)
    predictor_ad = zero
  END SUBROUTINE compute_predictors_ad
    
    
  ! ==============================================================
  ! Compute the component coefficient from the regression equation
  ! ==============================================================
  
  ! Forward model
  SUBROUTINE compute_coefficient( &
    c           , &  ! Input
    X           , &  ! Input
    Coefficient   )  ! Output
    ! Arguments
    REAL(fp), INTENT(IN)  :: c(:)  ! regression coefficient
    REAL(fp), INTENT(IN)  :: X(:)  ! predictor
    REAL(fp), INTENT(OUT) :: Coefficient
    ! Local variables
    INTEGER :: i
    ! Compute component coefficient
    coefficient = zero
    DO i = 1, n_predictors
      coefficient = coefficient + c(i)*x(i)
    END DO
  END SUBROUTINE compute_coefficient


  ! Tangent-linear model
  SUBROUTINE compute_coefficient_tl( &
    c             , &  ! Input
    X_TL          , &  ! Input
    Coefficient_TL  )  ! Output
    ! Arguments
    REAL(fp), INTENT(IN)  :: c(:)     ! regression coefficient
    REAL(fp), INTENT(IN)  :: X_TL(:)  ! predictor
    REAL(fp), INTENT(OUT) :: Coefficient_TL
    ! Local variables
    INTEGER :: i
    ! Compute tangent-linear component coefficient
    coefficient_tl = zero
    DO i = 1, n_predictors
      coefficient_tl = coefficient_tl + c(i)*x_tl(i)
    END DO
  END SUBROUTINE compute_coefficient_tl


  ! Adjoint model
  SUBROUTINE compute_coefficient_ad( &
    c             , &  ! Input
    Coefficient_AD, &  ! Input
    X_AD            )  ! Output
    ! Arguments
    REAL(fp), INTENT(IN)     :: c(:)     ! regression coefficient
    REAL(fp), INTENT(IN OUT) :: Coefficient_AD
    REAL(fp), INTENT(IN OUT) :: X_AD(:)  ! predictor
    ! Local variables
    INTEGER :: i
    ! Compute adjoint of the component coefficient
    DO i = 1, n_predictors
      x_ad(i) = x_ad(i) + c(i)*coefficient_ad
    END DO
    coefficient_ad = zero
  END SUBROUTINE compute_coefficient_ad


  PURE FUNCTION  azimuth_freq_correction( Frequency ) RESULT( Fre_C )
    IMPLICIT NONE
    REAL( fp ), INTENT(IN) :: frequency
    REAL( fp ) :: fre_c
    INTEGER :: i
      ! Data for the frequency correction
    REAL(fp), PARAMETER :: x(9) = (/ 0.0_fp, 1.4_fp, 6.8_fp, 10.7_fp, 19.35_fp, &
                                   37._fp, 89._fp, 150._fp, 200._fp/)
    REAL(fp), PARAMETER :: y(9) = (/ 0.0_fp, 0.1_fp, 0.6_fp, 0.9_fp, 1._fp, &
                                   1.0_fp, 0.4_fp, 0.2_fp, 0.0_fp/)
    ! 
    IF( frequency <= zero .or. frequency >= 200.0_fp ) THEN
      fre_c = zero
      RETURN
    ELSE
      DO i = 1, 8
        IF( frequency >= x(i) .and. frequency <= x(i+1) ) THEN
          fre_c = y(i) + (y(i+1)-y(i))/(x(i+1)-x(i))*(frequency-x(i))
          RETURN
        END IF
      END DO
    END IF
    fre_c = zero
    
  END FUNCTION  azimuth_freq_correction

END MODULE azimuth_emissivity_module