ODPS_AtmAbsorption.f90

Go to the documentation of this file.

!
! ODPS_AtmAbsorption
!
! Module containing routines to compute the optical depth profile
! due to gaseous absorption using the Optical Depth Pressure Space
! (ODPS) algorithm
!
!
! CREATION HISTORY:
!       Written by:     Yong Han & Yong Chen, JCSDA, NOAA/NESDIS 20-Jun-2008
!            TL,AD:     Tong Zhu, CIRA/CSU@NOAA/NESDIS 06-Jan-2009
!

MODULE odps_atmabsorption

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE type_kinds,               ONLY: fp
  USE crtm_parameters,          ONLY: zero, limit_exp, limit_log
  USE crtm_atmosphere_define,   ONLY: crtm_atmosphere_type
  USE crtm_geometryinfo_define, ONLY: crtm_geometryinfo_type, &
                                      crtm_geometryinfo_getvalue
  USE crtm_atmoptics_define,    ONLY: crtm_atmoptics_type
  USE odps_define,              ONLY: odps_type          , &
                                      significance_optran
  USE odps_predictor_define,    ONLY: odps_predictor_type, &
                                      pafv_associated    , &
                                      max_optran_order   , &
                                      max_optran_used_predictors
  USE odps_coordinatemapping,   ONLY: interpolate_profile      , &
                                      interpolate_profile_f1_tl, &
                                      interpolate_profile_f1_ad, &
                                      compute_interp_index

  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  ! Everything private by default
  PRIVATE
  ! Datatypes
  PUBLIC :: ivar_type
  ! Procedures
  PUBLIC :: odps_compute_atmabsorption
  PUBLIC :: odps_compute_atmabsorption_tl
  PUBLIC :: odps_compute_atmabsorption_ad


  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PRIVATE, PARAMETER :: module_version_id = &
  '$Id: $'
  ! Maximum allowed layer optical depth
  REAL(fp), PARAMETER :: max_od = 20.0_fp


  ! ------------------------------------------
  ! Structure definition to hold forward model
  ! variables across FWD, TL, and AD calls
  ! ------------------------------------------
  TYPE :: ivar_type
    PRIVATE
    INTEGER :: dummy
  END TYPE ivar_type


CONTAINS


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

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       ODPS_Compute_AtmAbsorption
!
! PURPOSE:
!       Subroutine to calculate the layer optical depths due to gaseous
!       absorption for a given sensor and channel and atmospheric profile.
!
! CALLING SEQUENCE:
!       CALL ODPS_Compute_AtmAbsorption( &
!              TC           , &
!              ChannelIndex , &
!              Predictor    , &
!              AtmAbsorption  )
!
! INPUTS:
!       TC:
!         ODPS structure holding tau coefficients
!         UNITS:      N/A
!         TYPE:       ODPS_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       ChannelIndex:
!         Channel index id. This is a unique index associated with a
!         supported sensor channel used to access the shared coefficient
!         data for a particular sensor's channel.
!         UNITS:      N/A
!         TYPE:       INTEGER
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       Predictor: 
!         Structure containing the predictor profile data.
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       AtmAbsorption:
!         Structure containing computed optical depth profile data.
!         UNITS:      N/A
!         TYPE:       CRTM_AtmOptics_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN OUT)
!
!:sdoc-:
!------------------------------------------------------------------------------

  SUBROUTINE odps_compute_atmabsorption( &
    TC           , &  ! Input
    ChannelIndex , &  ! Input
    Predictor    , &  ! Input
    AtmAbsorption  )  ! Output
    ! Arguments
    TYPE(odps_type)          , INTENT(IN)     :: tc
    INTEGER                  , INTENT(IN)     :: channelindex
    TYPE(odps_predictor_type), INTENT(IN OUT) :: predictor
    TYPE(crtm_atmoptics_type), INTENT(IN OUT) :: atmabsorption
    ! Local variables
    INTEGER  :: n_layers, n_user_layers
    INTEGER  :: i          ! coefficent index
    INTEGER  :: k          ! Layer index
    INTEGER  :: j          ! Component index
    INTEGER  :: np         ! number of predictors
    REAL(fp) :: od(predictor%n_layers)
    REAL(fp) :: od_path(0:predictor%n_layers)
    REAL(fp) :: user_od_path(0:predictor%n_user_layers)
    INTEGER  :: odps2user_idx(2, 0:predictor%n_user_layers)
    REAL(fp) :: od_tmp
    LOGICAL  :: optran
    INTEGER  :: j0, js

    ! ------
    ! Set up
    ! ------
    ! Assign the indices to a short name
    n_layers = predictor%n_Layers
    n_user_layers = predictor%n_User_Layers

    !--------------------------------------------------------
    ! Compute optical path profile
    !--------------------------------------------------------
    ! Loop over each tau component for optical depth calculation
    od = zero
    component_loop: DO j = 1, predictor%n_Components

      ! number of predictors for the current component and channel
      ! Note, TC%n_Predictors(j, ChannelIndex) <= Predictor%n_CP(j).
      ! For example, if the upper m predictors have zero coefficients,
      ! then, only coefficients with indexed 1 to (Predictor%n_CP(j) - m)
      ! are stored and used used in the OD calculations.
      np = tc%n_Predictors(j, channelindex)

      ! Check if there is any absorption for the component&channel combination.
      IF( np <= 0 ) cycle component_loop

      ! set flag for possible OPTRAN algorithm
      ! If this flag is set, this component is computed using OPTRAN algorithm.
      ! Otherwise, it is computed using ODPS algorithm.
      IF( predictor%OPTRAN .AND. j == tc%OComponent_Index)THEN
         optran = tc%OSignificance(channelindex) == significance_optran
      ELSE
         optran = .false.
      END IF

      IF(optran)THEN
        CALL add_optran_wlood(tc,          &
                             channelindex,  &
                             predictor,     &
                             od )
      ELSE

        ! ODPS algorithm
        j0 = tc%Pos_Index(j, channelindex)
        DO i = 1, np
          js = j0+(i-1)*n_layers-1
          DO k = 1, n_layers
            od(k) = od(k) + tc%C(js+k)*predictor%X(k, i, j)
          END DO
        END DO

      END IF

    END DO component_loop

    !------------------------------------------------------
    ! Compute optical path (level to space) profile
    !------------------------------------------------------
    od_path(0) = zero
    DO k = 1, n_layers
      od_tmp = od(k)
      IF(od(k) < zero)THEN
        od_tmp = zero
      ELSE IF(od(k) > max_od)THEN
        od_tmp = max_od
      END IF
      od_path(k) = od_path(k-1) + od_tmp
    END DO

    ! Save forward variables
    ! Interpolate the path profile back on the user pressure grids,
    ! Compute layer optical depths (vertical direction)
    IF ( pafv_associated(predictor%PAFV) ) THEN
      ! save forwad variables
      predictor%PAFV%OD = od
      predictor%PAFV%OD_Path = od_path
      ! If interpolation indexes are known
      user_od_path(0) = zero
      CALL interpolate_profile(predictor%PAFV%ODPS2User_Idx,    &
                               od_path,                         &
                               predictor%Ref_Level_LnPressure,  &
                               predictor%User_Level_LnPressure, &
                               user_od_path)
    ELSE
      ! interpolation indexes are not known
      CALL compute_interp_index(predictor%Ref_Level_LnPressure, &
                                predictor%User_Level_LnPressure,&
                                odps2user_idx)
      CALL interpolate_profile(odps2user_idx,                   &
                               od_path,                         &
                               predictor%Ref_Level_LnPressure,  &
                               predictor%User_Level_LnPressure, &
                               user_od_path)
    END IF

    ! Optical depth profile scaled to zenith.  Note that the scaling
    ! factor is the surface secant zenith angle.
    atmabsorption%Optical_Depth = (user_od_path(1:n_user_layers) - &
                                   user_od_path(0:n_user_layers-1)) / &
                                   predictor%Secant_Zenith_Surface

  END SUBROUTINE odps_compute_atmabsorption

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       ODPS_Compute_AtmAbsorption_TL
!
! PURPOSE:
!       Subroutine to calculate the tangent-linear layer optical depths due
!       to gaseous absorption for a given sensor and channel and atmospheric
!       profile.
!
! CALLING SEQUENCE:
!       CALL ODPS_Compute_AtmAbsorption_TL( &
!              TC           ,   &
!              ChannelIndex ,   &
!              Predictor    ,   &
!              Predictor_TL,    &
!              AtmAbsorption_TL )
!
! INPUTS:
!       TC:
!         ODPS structure holding tau coefficients
!         UNITS:      N/A
!         TYPE:       ODPS_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       ChannelIndex:
!         Channel index id. This is a unique index associated with a
!         supported sensor channel used to access the shared coefficient
!         data for a particular sensor's channel.
!         UNITS:      N/A
!         TYPE:       INTEGER
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       Predictor: 
!         Structure containing the predictor profile data.
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       Predictor_TL: 
!         Structure containing the tangent-linear predictor profile data.
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       AtmAbsorption_TL:
!         Structure containing computed tangent-linear optical depth profile
!         data.
!         UNITS:      N/A
!         TYPE:       CRTM_AtmOptics_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN OUT)
!
!:sdoc-:
!------------------------------------------------------------------------------

  SUBROUTINE odps_compute_atmabsorption_tl( &
    TC              , &  ! Input
    ChannelIndex    , &  ! Input
    Predictor       , &  ! Input
    Predictor_TL    , &  ! Input
    AtmAbsorption_TL  )  ! Output
    ! Arguments
    TYPE(odps_type)          , INTENT(IN)     :: tc
    INTEGER                  , INTENT(IN)     :: channelindex
    TYPE(odps_predictor_type), INTENT(IN)     :: predictor
    TYPE(odps_predictor_type), INTENT(IN OUT) :: predictor_tl
    TYPE(crtm_atmoptics_type), INTENT(IN OUT) :: atmabsorption_tl
    ! Local variables
    INTEGER  :: n_layers, n_user_layers
    INTEGER  :: i          ! coefficent index
    INTEGER  :: k          ! Layer index
    INTEGER  :: j          ! Component index
    INTEGER  :: np         ! number of predictors
    REAL(fp) :: od_tl(predictor%n_layers)
    REAL(fp) :: od_path_tl(0:predictor%n_layers)
    REAL(fp) :: user_od_path_tl(0:predictor%n_user_layers)
    LOGICAL  :: optran
    INTEGER  :: j0, js

    ! ------
    ! Set up
    ! ------
    ! Assign the indices to a short name
    n_layers = predictor%n_Layers
    n_user_layers = predictor%n_User_Layers

    !-----------------------------
    ! Compute optical path profile
    !-----------------------------
    ! Loop over each tau component for optical depth calculation
    od_tl = zero
    component_loop: DO j = 1, predictor%n_Components

      ! number of predictors for the current component and channel
      ! Note, TC%n_Predictors(j, ChannelIndex) <= Predictor%n_CP(j).
      ! For example, if the upper m predictors have zero coefficients,
      ! then, only coefficients with indexed 1 to (Predictor%n_CP(j) - m)
      ! are stored and used used in the OD calculations.
      np = tc%n_Predictors(j, channelindex)

      ! Check if there is any absorption for the component&channel combination.
      IF( np <= 0 ) cycle component_loop

      ! set flag for possible OPTRAN algorithm
      ! If this flag is set, this component is computed using OPTRAN algorithm.
      ! Otherwise, it is computed using ODPS algorithm.
      IF( predictor%OPTRAN .AND. j == tc%OComponent_Index)THEN
         optran = tc%OSignificance(channelindex) == significance_optran
      ELSE
         optran = .false.
      END IF

      IF(optran)THEN
        CALL add_optran_wlood_tl(tc,       &
                             channelindex,  &
                             predictor,     &
                             predictor_tl,  &
                             od_tl)
      ELSE

        ! ODPS algorithm
        j0 = tc%Pos_Index(j, channelindex)
        DO i = 1, np
          js = j0+(i-1)*n_layers-1
          DO k = 1, n_layers
            od_tl(k) = od_tl(k) + tc%C(js+k)*predictor_tl%X(k, i, j)
          END DO
        END DO

      END IF

    END DO component_loop

    !------------------------------------------------------
    ! Compute optical path (level to space) profile
    !------------------------------------------------------
    od_path_tl(0) = zero
    DO k = 1, n_layers
      IF(predictor%PAFV%OD(k) < zero)THEN
        od_tl(k) = zero
      ELSE IF(predictor%PAFV%OD(k) > max_od)THEN
        od_tl(k) = zero
      END IF
      od_path_tl(k) = od_path_tl(k-1) + od_tl(k)
    END DO

    ! Interpolate the path profile back on the user pressure grids,
    ! Compute layer optical depths (vertical direction)
    CALL interpolate_profile_f1_tl(predictor%PAFV%ODPS2User_Idx,    &
                                   predictor%Ref_Level_LnPressure,  &
                                   predictor%User_Level_LnPressure, &
                                   od_path_tl,                      &
                                   user_od_path_tl)

    atmabsorption_tl%Optical_Depth = (user_od_path_tl(1:n_user_layers) - &
                                   user_od_path_tl(0:n_user_layers-1)) / &
                                   predictor%Secant_Zenith_Surface

  END SUBROUTINE odps_compute_atmabsorption_tl

!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       ODPS_Compute_AtmAbsorption_AD
!
! PURPOSE:
!       Subroutine to calculate the layer optical depth adjoint due to gaseous
!       absorption for a given sensor and channel and atmospheric profile.
!
! CALLING SEQUENCE:
!       CALL ODPS_Compute_AtmAbsorption_AD( &
!              TC              , &
!              ChannelIndex    , &
!              Predictor       , &
!              AtmAbsorption_AD, &
!              Predictor_AD      )
!
! INPUTS:
!       TC:
!         ODPS structure holding tau coefficients
!         UNITS:      N/A
!         TYPE:       ODPS_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       ChannelIndex:
!         Channel index id. This is a unique index associated with a
!         supported sensor channel used to access the shared coefficient
!         data for a particular sensor's channel.
!         UNITS:      N/A
!         TYPE:       INTEGER
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       Predictor: 
!         Structure containing the predictor profile data.
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       AtmAbsorption_AD:
!         Structure containing optical depth adjoint profile data.
!         **NOTE: Optical depth component Set to zero upon exit **
!         UNITS:      N/A
!         TYPE:       CRTM_AtmOptics_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUTS:
!       Predictor_AD: 
!         Structure containing the adjoint predictor profile data.
!         **NOTE: Must contain data upon entry **
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN OUT)
!
!:sdoc-:
!------------------------------------------------------------------------------

  SUBROUTINE odps_compute_atmabsorption_ad( &
    TC              , &  ! Input
    ChannelIndex    , &  ! Input
    Predictor       , &  ! Input
    AtmAbsorption_AD, &  ! Input
    Predictor_AD      )  ! Output
    ! Arguments
    TYPE(odps_type)          , INTENT(IN)     :: tc
    INTEGER                  , INTENT(IN)     :: channelindex
    TYPE(odps_predictor_type), INTENT(IN)     :: predictor
    TYPE(crtm_atmoptics_type), INTENT(IN OUT) :: atmabsorption_ad
    TYPE(odps_predictor_type), INTENT(IN OUT) :: predictor_ad
    ! Local variables
    INTEGER  :: n_layers, n_user_layers
    INTEGER  :: i          ! coefficent index
    INTEGER  :: k          ! Layer index
    INTEGER  :: j          ! Component index
    INTEGER  :: np         ! number of predictors
    REAL(fp) :: od_ad(predictor%n_layers)
    REAL(fp) :: od_path_ad(0:predictor%n_layers)
    REAL(fp) :: user_od_path_ad(0:predictor%n_user_layers)
    LOGICAL  :: optran
    INTEGER  :: j0, js

    ! ------
    ! Set up
    ! ------
    ! Assign the indices to a short name
    n_layers = predictor%n_Layers
    n_user_layers = predictor%n_User_Layers

    !------- Adjoint part ---------

    ! Interpolate the path profile back on the user pressure grids,
    ! Compute layer optical depths (vertical direction)
    user_od_path_ad(n_user_layers) = zero
    DO k = n_user_layers, 1, -1
      user_od_path_ad(k) = user_od_path_ad(k) &
                           + atmabsorption_ad%Optical_Depth(k)/predictor%Secant_Zenith_Surface
      ! combined with initilization
      user_od_path_ad(k-1) = -atmabsorption_ad%Optical_Depth(k)/predictor%Secant_Zenith_Surface
    END DO
    atmabsorption_ad%Optical_Depth = zero

    od_path_ad = zero
    CALL interpolate_profile_f1_ad(predictor%PAFV%ODPS2User_Idx,       &
                                   predictor%Ref_Level_LnPressure,     &
                                   predictor%User_Level_LnPressure,    &
                                   user_od_path_ad,                    &
                                   od_path_ad )


    user_od_path_ad(0) = zero

    !-----------------------------
    ! Compute optical path profile
    !-----------------------------

    DO k = n_layers, 1, -1
      od_path_ad(k-1) = od_path_ad(k-1) + od_path_ad(k)
      ! combined with initialization
      od_ad(k) = od_path_ad(k)
      od_path_ad(k) = zero
      IF(predictor%PAFV%OD(k) < zero)THEN
        od_ad(k) = zero
      ELSE IF(predictor%PAFV%OD(k) > max_od)THEN
        od_ad(k) = zero
      END IF
    END DO
    od_path_ad(0) = zero

    ! Loop over each tau component for optical depth calculation

    component_loop_ad: DO j = 1, predictor%n_Components

      ! number of predictors for the current component and channel
      ! Note, TC%n_Predictors(j, ChannelIndex) <= Predictor%n_CP(j).
      ! For example, if the upper m predictors have zero coefficients,
      ! then, only coefficients with indexed 1 to (Predictor%n_CP(j) - m)
      ! are stored and used used in the OD calculations.
      np = tc%n_Predictors(j, channelindex)

      ! Check if there is any absorption for the component&channel combination.
      IF( np <= 0 ) cycle component_loop_ad

      IF( predictor%OPTRAN .AND. j == tc%OComponent_Index)THEN
         optran = tc%OSignificance(channelindex) == significance_optran
      ELSE
         optran = .false.
      END IF

     IF(optran)THEN

       CALL add_optran_wlood_ad(tc,          &
                                channelindex, &
                                predictor,    &
                                od_ad,        &
                                predictor_ad )
     ELSE

       ! ODPS algorithm
       j0 = tc%Pos_Index(j, channelindex)
       DO i = 1, np
         js = j0+(i-1)*n_layers-1
         DO k = n_layers, 1, -1
           predictor_ad%X(k, i, j) = predictor_ad%X(k, i, j) + tc%C(js+k)*od_ad(k)
         END DO
       END DO

     END IF

    END DO component_loop_ad

  END SUBROUTINE odps_compute_atmabsorption_ad


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

!------------------------------------------------------------------------------
!
! NAME:
!       Add_OPTRAN_wloOD
!
! PURPOSE:
!       Subroutine to calculate and add the layer optical depths due to water
!       vapor line absorption via the ODAS algorithm.
!       Note the OD argument is an in/out argument, which may hold accumulated
!       optical depths from other absorbers.
!
! CALLING SEQUENCE:
!        CALL Add_OPTRAN_wloOD( &
!               TC          , &
!               ChannelIndex, &
!               Predictor   , &
!               OD            )
!
! INPUTS:
!       TC:
!         ODPS structure holding tau coefficients
!         UNITS:      N/A
!         TYPE:       ODPS_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       ChannelIndex:
!         Channel index id. This is a unique index associated with a
!         supported sensor channel used to access the shared coefficient
!         data for a particular sensor's channel.
!         UNITS:      N/A
!         TYPE:       INTEGER
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       Predictor: 
!         Structure containing the predictor profile data.
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
! IN/OUTPUTS:
!        OD:
!          Slant path optical depth profile.
!          UNITS:      N/A
!          TYPE:       REAL(fp)
!          DIMENSION:  Rank-1 array (n_Layers)
!          ATTRIBUTES: INTENT(IN OUT)
!
!
!------------------------------------------------------------------------------

  SUBROUTINE add_optran_wlood( &
    TC          , &  ! Input
    ChannelIndex, &  ! Input
    Predictor   , &  ! Input
    OD            )  ! In/Output
    ! Arguments
    TYPE(ODPS_type),           INTENT(IN)     :: TC
    INTEGER,                   INTENT(IN)     :: ChannelIndex
    TYPE(ODPS_Predictor_type), INTENT(IN OUT) :: Predictor
    REAL(fp),                  INTENT(IN OUT) :: OD(:)
    ! Local variables
    REAL(fp) :: LN_Chi(TC%n_Layers), coeff(0:MAX_OPTRAN_ORDER)
    REAL(fp) :: b(TC%n_Layers, 0:MAX_OPTRAN_USED_PREDICTORS)
    REAL(fp) :: Chi(TC%n_Layers)
    INTEGER  :: np, n_Layers, n_orders, js, i, j, k, ii, jj

    ! -----------------------------------------
    ! Check if there is any absorption for this
    ! absorber/channel combination.
    ! -----------------------------------------
    np = tc%OP_Index(0,channelindex)  ! number of predictors
    IF ( np <= 0 ) RETURN

    n_layers = tc%n_Layers
    js = tc%OPos_Index(channelindex)
    n_orders = tc%Order(channelindex)

    DO i = 0, np
      jj = js + i*(n_orders+1)
      coeff(0:n_orders) = tc%OC(jj:jj+n_orders)
      DO k = 1, n_layers
        b(k,i) = coeff(0)
        DO j = 1, n_orders
          b(k,i) = b(k,i) + coeff(j)*predictor%Ap(k, j)
        END DO
      END DO
    END DO

    ln_chi = b(:,0)
    DO i = 1, np
      ii = tc%OP_Index(i,channelindex)
      DO k = 1, n_layers
        ln_chi(k) = ln_chi(k) + b(k, i)* predictor%OX(k, ii)
      END DO
    END DO

    DO k = 1, n_layers
      IF( ln_chi(k) > limit_exp ) THEN
        chi(k) = limit_log
      ELSE IF( ln_chi(k) < -limit_exp ) THEN
        chi(k) = zero
      ELSE
        chi(k) = exp(ln_chi(k))
      ENDIF
      od(k) = od(k) + chi(k)*predictor%dA(k)
    END DO

    IF(predictor%PAFV%OPTRAN)THEN
      predictor%PAFV%b      = b
      predictor%PAFV%LN_Chi = ln_chi
      predictor%PAFV%Chi    = chi
    END IF

  END SUBROUTINE add_optran_wlood


!------------------------------------------------------------------------------
!
! NAME:
!       Add_OPTRAN_wloOD_TL
!
! PURPOSE:
!       Subroutine to calculate and add the tangent-linear layer optical depths
!       due to water vapor line absorption via the ODAS algorithm.
!       Note the OD_TL argument is an in/out argument, which may hold accumulated
!       tangent-linear optical depths from other absorbers.
!
! CALLING SEQUENCE:
!        CALL Add_OPTRAN_wloOD_TL( &
!               TC          , &
!               ChannelIndex, &
!               Predictor   , &
!               Predictor_TL, &
!               OD_TL         )
!
! INPUTS:
!       TC:
!         ODPS structure holding tau coefficients
!         UNITS:      N/A
!         TYPE:       ODPS_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       ChannelIndex:
!         Channel index id. This is a unique index associated with a
!         supported sensor channel used to access the shared coefficient
!         data for a particular sensor's channel.
!         UNITS:      N/A
!         TYPE:       INTEGER
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       Predictor: 
!         Structure containing the predictor profile data.
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       Predictor_TL: 
!         Structure containing the tangent-linear predictor profile data.
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
! IN/OUTPUT ARGUMENTS:
!       OD_TL:
!         Slant path tangent-linear optical depth profile
!         UNITS:      N/A
!         TYPE:       REAL(fp)
!         DIMENSION:  Rank-1 array (n_Layers)
!         ATTRIBUTES: INTENT(IN OUT)
!
!------------------------------------------------------------------------------

  SUBROUTINE add_optran_wlood_tl( &
    TC          , &
    ChannelIndex, &
    Predictor   , &
    Predictor_TL, &
    OD_TL         )
    ! Arguments
    TYPE(ODPS_type),           INTENT(IN)     :: TC
    INTEGER,                   INTENT(IN)     :: ChannelIndex
    TYPE(ODPS_Predictor_type), INTENT(IN)     :: Predictor
    TYPE(ODPS_Predictor_type), INTENT(IN)     :: Predictor_TL
    REAL(fp),                  INTENT(IN OUT) :: OD_TL(:)
    ! Local variables
    REAL(fp) :: coeff(0:MAX_OPTRAN_ORDER)
    REAL(fp) :: LN_Chi_TL(TC%n_Layers), b_TL(TC%n_Layers, 0:MAX_OPTRAN_USED_PREDICTORS)
    REAL(fp) :: chi_TL(TC%n_Layers)
    INTEGER  :: np, n_Layers, n_orders, js, i, j, k, ii, jj

    ! -----------------------------------------
    ! Check if there is any absorption for this
    ! absorber/channel combination.
    ! -----------------------------------------
    np = tc%OP_Index(0,channelindex)  ! number of predictors
    IF ( np <= 0 ) RETURN

    n_layers = tc%n_Layers
    js = tc%OPos_Index(channelindex)
    n_orders = tc%Order(channelindex)

    DO i = 0, np
      jj = js + i*(n_orders+1)
      coeff(0:n_orders) = tc%OC(jj:jj+n_orders)
      DO k = 1, n_layers
        b_tl(k,i) = zero
        DO j = 1, n_orders
          b_tl(k,i) = b_tl(k,i) + coeff(j)*predictor_tl%Ap(k, j)
        END DO
      END DO
    END DO

    ln_chi_tl = b_tl(:,0)
    DO i = 1, np
      ii = tc%OP_Index(i,channelindex)
      DO k = 1, n_layers
        ln_chi_tl(k) = ln_chi_tl(k) + b_tl(k, i)* predictor%OX(k, ii) + predictor%PAFV%b(k, i)* predictor_tl%OX(k, ii)
      END DO
    END DO

    DO k = 1, n_layers
      IF( predictor%PAFV%LN_Chi(k) > limit_exp ) THEN
        chi_tl(k) = zero
      ELSE IF( predictor%PAFV%LN_Chi(k) < -limit_exp ) THEN
        chi_tl(k) = zero
      ELSE
        chi_tl(k) = predictor%PAFV%Chi(k) * ln_chi_tl(k)
      ENDIF
      od_tl(k) = od_tl(k) + chi_tl(k)*predictor%dA(k) + predictor%PAFV%Chi(k)*predictor_tl%dA(k)
    END DO

  END SUBROUTINE add_optran_wlood_tl


!------------------------------------------------------------------------------
!
! NAME:
!       Add_OPTRAN_wloOD_AD
!
! PURPOSE:
!       Subroutine to calculate and add the layer optical depth adjoints due
!       to water vapor line absorption via the ODAS algorithm.
!
! CALLING SEQUENCE:
!        CALL Add_OPTRAN_wloOD_AD( &
!               TC          , &
!               ChannelIndex, &
!               Predictor   , &
!               OD_AD       , &
!               Predictor_AD  )
!
! INPUTS:
!       TC:
!         ODPS structure holding tau coefficients
!         UNITS:      N/A
!         TYPE:       ODPS_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       ChannelIndex:
!         Channel index id. This is a unique index associated with a
!         supported sensor channel used to access the shared coefficient
!         data for a particular sensor's channel.
!         UNITS:      N/A
!         TYPE:       INTEGER
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       Predictor: 
!         Structure containing the predictor profile data.
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN)
!
!       OD_AD:
!         Slant path adjoint optical depth profile
!         UNITS:      N/A
!         TYPE:       REAL(fp)
!         DIMENSION:  Rank-1 array (n_Layers)
!         ATTRIBUTES: INTENT(IN)
!
! IN/OUTPUTS:
!       Predictor_AD: 
!         Structure containing the adjoint predictor profile data.
!         **NOTE: Must contain data upon entry **
!         UNITS:      N/A
!         TYPE:       ODPS_Predictor_type
!         DIMENSION:  Scalar
!         ATTRIBUTES: INTENT(IN OUT)
!
!------------------------------------------------------------------------------

  SUBROUTINE add_optran_wlood_ad( &
    TC          , &
    ChannelIndex, &
    Predictor   , &
    OD_AD       , &
    Predictor_AD  )
    ! Arguments
    TYPE(ODPS_type),           INTENT(IN)     :: TC
    INTEGER,                   INTENT(IN)     :: ChannelIndex
    TYPE(ODPS_Predictor_type), INTENT(IN)     :: Predictor
    REAL(fp),                  INTENT(IN OUT) :: OD_AD(:)
    TYPE(ODPS_Predictor_type), INTENT(IN OUT) :: Predictor_AD
    ! Local variables
    REAL(fp) :: coeff(0:MAX_OPTRAN_ORDER)
    REAL(fp) :: LN_Chi_AD(TC%n_Layers), b_AD(TC%n_Layers, 0:MAX_OPTRAN_USED_PREDICTORS)
    REAL(fp) :: Chi_AD(TC%n_Layers)
    INTEGER  :: np, n_Layers, n_orders, js, i, j, k, ii, jj

      !------ Forward part for LN_Chi, b -----------

      ! -----------------------------------------
      ! Check if there is any absorption for this
      ! absorber/channel combination.
      ! -----------------------------------------
      np = tc%OP_Index(0,channelindex)  ! number of predictors
      IF ( np <= 0 ) RETURN

      n_layers = tc%n_Layers
      js = tc%OPos_Index(channelindex)
      n_orders = tc%Order(channelindex)

      !------ Adjoint part ----------------------

      ! -----------------------------------------
      ! Check if there is any absorption for this
      ! absorber/channel combination.
      ! -----------------------------------------

      chi_ad = zero
      ln_chi_ad = zero
      DO k = n_layers, 1, -1

        chi_ad(k) = chi_ad(k) + od_ad(k) * predictor%dA(k)
        predictor_ad%dA(k) = predictor_ad%dA(k) + od_ad(k) * predictor%PAFV%Chi(k)
        IF( predictor%PAFV%LN_Chi(k) > limit_exp ) THEN
          chi_ad(k) = zero
        ELSE IF( predictor%PAFV%LN_Chi(k) < -limit_exp ) THEN
          chi_ad(k) = zero
        ELSE
          ln_chi_ad(k) = predictor%PAFV%Chi(k) * chi_ad(k)  ! combinded with initialization for LN_Chi_AD(k)
        ENDIF
      END DO

      DO i = 1, np
        ii = tc%OP_Index(i,channelindex)
        DO k = n_layers, 1, -1
          b_ad(k, i) = ln_chi_ad(k) * predictor%OX(k, ii)  ! Combinded with initialization for b_AD
          predictor_ad%OX(k, ii) = predictor_ad%OX(k, ii) + ln_chi_ad(k)*predictor%PAFV%b(k, i)
        END DO
      END DO
      b_ad(:,0) = ln_chi_ad

      DO i = 0, np
        jj = js + i*(n_orders+1)
        coeff(0:n_orders) = tc%OC(jj:jj+n_orders)
        DO k = n_layers, 1, -1
          DO j = 1, n_orders
            predictor_ad%Ap(k, j) = predictor_ad%Ap(k, j) + coeff(j)*b_ad(k,i)
          END DO
          b_ad(k,i) = zero
        END DO
      END DO

   END SUBROUTINE add_optran_wlood_ad

END MODULE odps_atmabsorption