CRTM_IRSSEM.f90

Go to the documentation of this file.

!
! CRTM_IRSSEM
!
! Module containing function to invoke the CRTM Infrared
! Sea Surface Emissivity Model (IRSSEM).
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 22-Jun-2005
!                       paul.vandelst@noaa.gov
!

MODULE crtm_irssem

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE type_kinds,          ONLY: fp
  USE message_handler,     ONLY: success, failure, display_message
  USE crtm_parameters,     ONLY: zero, one, degrees_to_radians
  USE crtm_interpolation,  ONLY: npts, &
                                 lpoly, &
                                 lpoly_type, &
                                 clear_lpoly, &
                                 find_index, &
                                 interp_3d, &
                                 lpoly_tl, &
                                 interp_3d_tl, &
                                 lpoly_ad, &
                                 interp_3d_ad
  USE irwatercoeff_define, ONLY: irwatercoeff_type
  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  ! Everything private by default
  PRIVATE
  ! Derived type
  PUBLIC :: ivar_type
  ! Procedures
  PUBLIC :: crtm_compute_irssem
  PUBLIC :: crtm_compute_irssem_tl
  PUBLIC :: crtm_compute_irssem_ad


  ! -----------------
  ! Module parameters
  ! -----------------
  ! Version Id for the module
  CHARACTER(*), PARAMETER :: module_version_id = &
  '$Id: CRTM_IRSSEM.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'
  ! Message string length
  INTEGER, PARAMETER :: ml = 256


  ! -------------------------------
  ! Structure definition to hold
  ! forward interpolating variables
  ! across fwd, tl and adjoint
  ! -------------------------------
  ! The interpolation routine structure
  TYPE :: einterp_type
    ! The dimensions
    INTEGER :: n_angles = 0
    INTEGER :: n_pts    = 0
    ! Allocation indicator
    LOGICAL :: is_allocated = .false.
    ! The interpolating polynomials
    TYPE(lpoly_type), ALLOCATABLE :: wlp(:)  ! Angle
    TYPE(lpoly_type)              :: xlp     ! Frequency
    TYPE(lpoly_type)              :: ylp     ! Wind Speed
    ! The LUT interpolation indices
    INTEGER, ALLOCATABLE :: i1(:), i2(:)     ! Angle
    INTEGER              :: j1    , j2       ! Frequency
    INTEGER              :: k1    , k2       ! Wind Speed
    ! The LUT interpolation boundary check
    LOGICAL, ALLOCATABLE :: a_outbound(:)    ! Angle
    LOGICAL              :: f_outbound       ! Frequency
    LOGICAL              :: v_outbound       ! Wind Speed
    ! The interpolation input
    REAL(fp), ALLOCATABLE :: a_int(:)        ! Angle
    REAL(fp)              :: f_int           ! Frequency
    REAL(fp)              :: v_int           ! Wind Speed
    ! The data to be interpolated
    REAL(fp), ALLOCATABLE :: a(:,:)          ! Angle
    REAL(fp), ALLOCATABLE :: f(:)            ! Frequency
    REAL(fp), ALLOCATABLE :: v(:)            ! Wind Speed
  END TYPE einterp_type

  ! The main internal variable structure
  TYPE :: ivar_type
    PRIVATE
    ! The interpolation data
    TYPE(einterp_type) :: ei
  END TYPE ivar_type


CONTAINS


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


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Compute_IRSSEM
!
! PURPOSE:
!       Function to compute the CRTM infrared sea surface emissivity (IRSSE)
!       for input wind speed, frequency, and angles.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Compute_IRSSEM( IRwaterCoeff, &
!                                           Wind_Speed  , &
!                                           Frequency   , &
!                                           Angle       , &
!                                           iVar        , &
!                                           Emissivity    )
!
! INPUTS:
!       IRwaterCoeff:   Infrared water emissivity model coefficient object.
!                       Load the object with the coefficients for the emissivity
!                       model to use.
!                       UNITS:      N/A
!                       TYPE:       IRwaterCoeff_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
!       Wind_Speed:     Wind speed.
!                       UNITS:      metres per second (m.s^-1)
!                       TYPE:       REAL(fp)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
!       Frequency:      Infrared frequency.
!                       UNITS:      inverse centimetres (cm^-1)
!                       TYPE:       REAL(fp)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
!       Angle:          Surface zenith angle.
!                       UNITS:      Degrees
!                       TYPE:       REAL(fp)
!                       DIMENSION:  Rank-1 (n_Angles)
!                       ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       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:       iVar_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(OUT)
!
!       Emissivity:     Sea surface emissivities for the
!                       requested wind speed, frequency, and angles.
!                       UNITS:      N/A
!                       TYPE:       REAL(fp)
!                       DIMENSION:  Same as input ANGLE argument.
!                       ATTRIBUTES: INTENT(OUT)
!
! FUNCTION RESULT:
!       Error_Status:   The return value is an integer defining the error status.
!                       The error codes are defined in the Message_Handler module.
!                       If == SUCCESS the computation was successful.
!                          == FAILURE an unrecoverable error occurred.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION crtm_compute_irssem( &
    IRwaterCoeff, &  ! Input model coefficients
    Wind_Speed  , &  ! Input
    Frequency   , &  ! Input
    Angle       , &  ! Input
    iVar        , &  ! Internal variable output
    Emissivity  ) &  ! Output
  result( err_stat )
    ! Arguments
    TYPE(irwatercoeff_type), INTENT(IN)  :: irwatercoeff
    REAL(fp)               , INTENT(IN)  :: wind_speed     ! v
    REAL(fp)               , INTENT(IN)  :: frequency      ! f
    REAL(fp)               , INTENT(IN)  :: angle(:)       ! a
    TYPE(ivar_type)        , INTENT(OUT) :: ivar
    REAL(fp)               , INTENT(OUT) :: emissivity(:)
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_Compute_IRSSEM'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: n_angles, i
    REAL(fp) :: sec_angle(size(angle))

    ! Set up
    err_stat = success
    ! ...Check dimensions
    n_angles = SIZE(angle)
    IF ( SIZE(emissivity) /= n_angles ) THEN
      err_stat = failure
      msg = 'Input Angle and output Emissivity array dimensions inconsistent.'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF
    ! ...Allocate interpolation variable structure
    CALL einterp_create( ivar%ei, npts, n_angles )
    IF ( .NOT. einterp_associated( ivar%ei ) ) THEN
      err_stat = failure
      msg = 'Error allocating interpolation variable structure.'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF
    ! ...Convert angles to secants
    sec_angle = one/cos(degrees_to_radians*angle)


    ! Compute the wind speed interpolating polynomial
    ! ...Find the LUT indices and check if input is out of bounds
    ivar%ei%v_int = wind_speed
    CALL find_index(irwatercoeff%Wind_Speed, &
                    ivar%ei%v_int, ivar%ei%k1, ivar%ei%k2, ivar%ei%v_outbound)
    ivar%ei%v = irwatercoeff%Wind_Speed(ivar%ei%k1:ivar%ei%k2)
    ! ...Compute the polynomial
    CALL lpoly( ivar%ei%v    , & ! Input
                ivar%ei%v_int, & ! Input
                ivar%ei%ylp    ) ! Output


    ! Compute the frequency interpolating polynomial
    ! ...Find the LUT indices and check if input is out of bounds
    ivar%ei%f_int = frequency
    CALL find_index(irwatercoeff%Frequency, &
                    ivar%ei%f_int, ivar%ei%j1, ivar%ei%j2, ivar%ei%f_outbound)
    ivar%ei%f = irwatercoeff%Frequency(ivar%ei%j1:ivar%ei%j2)
    ! ...Compute the polynomial
    CALL lpoly( ivar%ei%f    , & ! Input
                ivar%ei%f_int, & ! Input
                ivar%ei%xlp    ) ! Output


    ! Compute the angle interpolating polynomials
    DO i = 1, n_angles

      ! ...Find the LUT indices and check if input is out of bounds
      ivar%ei%a_int(i) = abs(angle(i))
      CALL find_index(irwatercoeff%Angle, &
                      ivar%ei%a_int(i), ivar%ei%i1(i), ivar%ei%i2(i), ivar%ei%a_outbound(i))
      ivar%ei%a(:,i) = irwatercoeff%Angle(ivar%ei%i1(i):ivar%ei%i2(i))

!! Secant interpolation test
!iVar%ei%a_int(i) = sec_angle(i)
!CALL find_index(IRwaterCoeff%Secant_Angle, &
!                iVar%ei%a_int(i), iVar%ei%i1(i), iVar%ei%i2(i), iVar%ei%a_outbound(i))
!iVar%ei%a(:,i) = IRwaterCoeff%Secant_Angle(iVar%ei%i1(i):iVar%ei%i2(i))

      ! ...Compute the polynomial
      CALL lpoly( ivar%ei%a(:,i)  , & ! Input
                  ivar%ei%a_int(i), & ! Input
                  ivar%ei%wlp(i)    ) ! Output


      ! Compute the interpolated emissivity
      CALL interp_3d( irwatercoeff%Emissivity( ivar%ei%i1(i):ivar%ei%i2(i), &
                                               ivar%ei%j1   :ivar%ei%j2   , &
                                               ivar%ei%k1   :ivar%ei%k2     ), & ! Input
                      ivar%ei%wlp(i), & ! Input
                      ivar%ei%xlp   , & ! Input
                      ivar%ei%ylp   , & ! Input
                      emissivity(i)   ) ! Output

    END DO

  END FUNCTION crtm_compute_irssem


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Compute_IRSSEM_TL
!
! PURPOSE:
!       Function to compute the tangent-linear CRTM infrared sea surface
!       emissivity (IRSSE) for input wind speed, frequency, and angles.
!
!       This function must be called *after* the forward model function,
!       CRTM_Compute_IRSSEM, has been called. The forward model function
!       populates the internal variable structure argument, iVar.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Compute_IRSSEM_TL( IRwaterCoeff , &
!                                              Wind_Speed_TL, &
!                                              iVar         , &
!                                              Emissivity_TL  )
! INPUTS:
!       IRwaterCoeff:   Infrared water emissivity model coefficient object.
!                       Load the object with the coefficients for the emissivity
!                       model to use.
!                       UNITS:      N/A
!                       TYPE:       IRwaterCoeff_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
!       Wind_Speed_TL:  The tangent-linear wind speed.
!                       UNITS:      metres per second (m.s^-1)
!                       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:       iVar_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Emissivity_TL:  Tangent-linear sea surface emissivity.
!                       UNITS:      N/A
!                       TYPE:       REAL(fp)
!                       DIMENSION:  Rank-1 (n_Angles)
!                       ATTRIBUTES: INTENT(OUT)
!
! FUNCTION RESULT:
!       Error_Status:   The return value is an integer defining the error status.
!                       The error codes are defined in the Message_Handler module.
!                       If == SUCCESS the computation was successful.
!                          == FAILURE an unrecoverable error occurred.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION crtm_compute_irssem_tl( &
    IRwaterCoeff , &  ! Input model coefficients
    Wind_Speed_TL, &  ! Input
    iVar         , &  ! Internal variable input
    Emissivity_TL) &  ! Output
  result( err_stat )
    ! Arguments
    TYPE(irwatercoeff_type), INTENT(IN)  :: irwatercoeff
    REAL(fp)               , INTENT(IN)  :: wind_speed_tl
    TYPE(ivar_type)        , INTENT(IN)  :: ivar
    REAL(fp)               , INTENT(OUT) :: emissivity_tl(:)
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_Compute_IRSSEM_TL'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER  :: i
    REAL(fp) :: v_tl(npts)
    REAL(fp) :: e_tl(npts,npts,npts)
    TYPE(lpoly_type) :: ylp_tl, xlp_tl, wlp_tl

    ! Set up
    err_stat = success
    ! ...Check internal variable allocation
    IF ( .NOT. einterp_associated( ivar%ei ) ) THEN
      err_stat = failure
      msg = 'Internal structure ei is not allocated'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF
    ! ...Check dimensions
    IF ( SIZE( emissivity_tl ) /= ivar%ei%n_Angles ) THEN
      err_stat = failure
      msg = 'Input Emissivity_TL array dimensions inconsistent with number of angles.'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF
    ! ...No TL if wind speed is out of bounds
    IF ( ivar%ei%v_outbound ) THEN
      emissivity_tl = zero
      RETURN
    END IF
    ! ...Initialise local TL variables
    v_tl = zero
    e_tl = zero
    CALL clear_lpoly(wlp_tl)
    CALL clear_lpoly(xlp_tl)


    ! Calculate the TL interpolating
    ! polynomials for wind speed
    CALL lpoly_tl( ivar%ei%v, ivar%ei%v_int, & ! FWD Input
                   ivar%ei%ylp,              & ! FWD Input
                   v_tl, wind_speed_tl,      & ! TL  Input
                   ylp_tl                    ) ! TL  Output


    ! Begin loop over angles
    DO i = 1, ivar%ei%n_Angles

      ! Perform interpolation
      CALL interp_3d_tl(irwatercoeff%Emissivity(ivar%ei%i1(i):ivar%ei%i2(i), &
                                            ivar%ei%j1   :ivar%ei%j2   , &
                                            ivar%ei%k1   :ivar%ei%k2     ), & ! FWD Emissivity input
                        ivar%ei%wlp(i), & ! FWD polynomial input
                        ivar%ei%xlp   , & ! FWD polynomial input
                        ivar%ei%ylp   , & ! FWD polynomial input
                        e_tl, wlp_tl, xlp_tl, ylp_tl, & ! TL input
                        emissivity_tl(i)              ) ! Output

    END DO

  END FUNCTION crtm_compute_irssem_tl


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Compute_IRSSEM_AD
!
! PURPOSE:
!       Function to compute the adjoint of the CRTM infrared sea surface
!       emissivity (IRSSE) for input wind speed, frequency, and angles.
!
!       This function must be called *after* the forward model function,
!       CRTM_Compute_IRSSEM, has been called. The forward model function
!       populates the internal variable structure argument, iVar.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Compute_IRSSEM_AD( IRwaterCoeff , &
!                                              Emissivity_AD, &
!                                              iVar         , &
!                                              Wind_Speed_AD  )
!
! INPUTS:
!       IRwaterCoeff:   Infrared water emissivity model coefficient object.
!                       Load the object with the coefficients for the emissivity
!                       model to use.
!                       UNITS:      N/A
!                       TYPE:       IRwaterCoeff_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
!       Emissivity_AD:  Adjoint sea surface emissivity.
!                       *** SET TO ZERO ON EXIT ***
!                       UNITS:      N/A
!                       TYPE:       REAL(fp)
!                       DIMENSION:  Rank-1 (n_Angles)
!                       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:       iVar_type
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Wind_Speed_AD:  Adjoint wind speed.
!                       *** MUST HAVE VALUE ON ENTRY ***
!                       UNITS:      per metres per second, (m.s^-1)^-1
!                       TYPE:       REAL(fp)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN OUT)
!
! FUNCTION RESULT:
!       Error_Status:   The return value is an integer defining the error status.
!                       The error codes are defined in the Message_Handler module.
!                       If == SUCCESS the computation was successful.
!                          == FAILURE an unrecoverable error occurred.
!                       UNITS:      N/A
!                       TYPE:       INTEGER
!                       DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION crtm_compute_irssem_ad( &
    IRwaterCoeff , &  ! Input model coefficients
    Emissivity_AD, &  ! Input
    iVar         , &  ! Internal Variable Input
    Wind_Speed_AD) &  ! Output
  result( err_stat )
    ! Arguments
    TYPE(irwatercoeff_type), INTENT(IN)     :: irwatercoeff
    REAL(fp)               , INTENT(IN OUT) :: emissivity_ad(:)
    TYPE(ivar_type)        , INTENT(IN)     :: ivar
    REAL(fp)               , INTENT(IN OUT) :: wind_speed_ad
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_Compute_IRSSEM_AD'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER  :: i
    REAL(fp) :: e_ad(npts,npts,npts)
    REAL(fp) :: v_ad(npts)
    TYPE(lpoly_type) :: wlp_ad, xlp_ad, ylp_ad

    ! Set Up
    err_stat = success
    e_ad = zero
    v_ad = zero
    ! ...Check internal variable allocation
    IF ( .NOT. einterp_associated( ivar%ei ) ) THEN
      err_stat = failure
      msg = 'Internal structure ei is not allocated'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF
    ! ...Check dimensions
    IF ( SIZE(emissivity_ad) /= ivar%ei%n_Angles ) THEN
      err_stat = failure
      msg = 'Input Emissivity_AD array dimensions inconsistent with number of angles.'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF
    ! ...No AD if wind speed is out of bounds
    IF ( ivar%ei%v_outbound ) RETURN
    ! ...Initialize local variables
    CALL clear_lpoly(wlp_ad)
    CALL clear_lpoly(xlp_ad)
    CALL clear_lpoly(ylp_ad)

    ! Loop over emissivity calculation angles
    DO i = 1, ivar%ei%n_Angles

      ! Get the adjoint interpoalting polynomial for wind speed
      CALL interp_3d_ad(irwatercoeff%Emissivity( ivar%ei%i1(i):ivar%ei%i2(i), &
                                                 ivar%ei%j1   :ivar%ei%j2   , &
                                                 ivar%ei%k1   :ivar%ei%k2     ), & ! FWD Input
                        ivar%ei%wlp(i)  , & ! FWD Input
                        ivar%ei%xlp     , & ! FWD Input
                        ivar%ei%ylp     , & ! FWD Input
                        emissivity_ad(i), & ! AD Input
                        e_ad, wlp_ad, xlp_ad, ylp_ad ) ! AD Output

      ! Set adjoint emissivity to zero
      emissivity_ad(i) = zero

    END DO

    ! Compute the wind speed adjoint
    CALL lpoly_ad(ivar%ei%v    , & ! FWD Input
                  ivar%ei%v_int, & ! FWD Input
                  ivar%ei%ylp  , & ! FWD Input
                  ylp_ad       , & ! AD  Input
                  v_ad         , & ! AD  Output
                  wind_speed_ad  ) ! AD  Output

  END FUNCTION crtm_compute_irssem_ad


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

  ! ----------------------------------------------
  ! Procedures to manipulate the Einterp structure
  ! ----------------------------------------------
  ELEMENTAL FUNCTION einterp_associated( ei ) RESULT( Status )
    TYPE(einterp_type), INTENT(IN) :: ei
    LOGICAL :: status
    status = ei%Is_Allocated
  END FUNCTION einterp_associated

  ELEMENTAL SUBROUTINE einterp_create( ei, n_Pts, n_Angles )
    TYPE(einterp_type), INTENT(OUT) :: ei
    INTEGER,            INTENT(IN)  :: n_pts
    INTEGER,            INTENT(IN)  :: n_angles
    INTEGER :: alloc_stat
    IF ( n_pts < 1 .OR. n_angles < 1 ) RETURN
    ALLOCATE( ei%wlp(n_angles)       , &
              ei%i1(n_angles)        , &
              ei%i2(n_angles)        , &
              ei%a_outbound(n_angles), &
              ei%a_int(n_angles)     , &
              ei%a(n_pts,n_angles)   , &
              ei%f(n_pts)            , &
              ei%v(n_pts)            , &
              stat = alloc_stat )
    IF ( alloc_stat /= 0 ) RETURN
    ei%n_Angles = n_angles
    ei%n_Pts    = n_pts
    ei%Is_Allocated = .true.
  END SUBROUTINE einterp_create

END MODULE crtm_irssem