CRTM_Atmosphere.f90

Go to the documentation of this file.

!
! CRTM_Atmosphere
!
! Module for adding layers to the CRTM atmosphere structure as required.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, 29-Oct-2007
!                       paul.vandelst@noaa.gov
!

MODULE crtm_atmosphere

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE type_kinds            , ONLY: fp
  USE message_handler       , ONLY: success, failure, warning, display_message
  USE crtm_parameters       , ONLY: zero, one, point_5, set, &
                                    toa_pressure           , &
                                    minimum_absorber_amount
  USE crtm_atmosphere_define, ONLY: crtm_atmosphere_type    , &
                                    OPERATOR(==), &
                                    OPERATOR(+), &
                                    crtm_atmosphere_associated, &
                                    crtm_atmosphere_create, &
                                    crtm_atmosphere_addlayercopy, &
                                    crtm_atmosphere_zero
  USE crtm_model_profiles   , ONLY: model_level_pressure, & 
                                    crtm_get_model_profile
  ! ...Internal variable definition module
  USE iatm_define,            ONLY: iatm_type      , &
                                    iatm_associated, &
                                    iatm_create    , &
                                    iatm_destroy
                                    
  ! Disable implicit typing
  IMPLICIT NONE


  ! ------------
  ! Visibilities
  ! ------------
  ! Everything private by default
  PRIVATE
  ! Module procedures
  PUBLIC :: crtm_atmosphere_addlayers
  PUBLIC :: crtm_atmosphere_addlayers_tl
  PUBLIC :: crtm_atmosphere_addlayers_ad
  ! iAtm entities
  ! ...Structure
  PUBLIC :: iatm_type      
  ! ...Procedures
  PUBLIC :: iatm_associated
  PUBLIC :: iatm_create
  PUBLIC :: iatm_destroy


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


CONTAINS


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

!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Atmosphere_AddLayers
!
! PURPOSE:
!       Function to copy an atmosphere structure and adding extra layers from
!       climatology as required to supplement the upper atmosphere profile data.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Atmosphere_AddLayers( Atm_In , &  ! Input
!                                                 Atm_Out  )  ! Output
!
! INPUTS:
!       Atm_In:          Atmosphere structure that is to be supplemented
!                        or copied.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Atm_Out:         Copy of the input atmosphere structure with extra upper
!                        atmosphere layers added as required.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        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 operation was successful
!                           == FAILURE an error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION crtm_atmosphere_addlayers( &
    Atm_In , &  ! Input
    Atm_Out) &  ! Output
  result( err_stat )
    ! Arguments
    TYPE(crtm_atmosphere_type), INTENT(IN)     :: atm_in
    TYPE(crtm_atmosphere_type), INTENT(IN OUT) :: atm_out
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_Atmosphere_AddLayers'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: i, j, k, n
    TYPE(iatm_type) :: iatm


    ! Set up
    err_stat = success


    ! If extra layers are NOT needed,
    ! then simply copy the structure
    IF ( atm_in%Level_Pressure(0) <= toa_pressure) THEN
      atm_out = atm_in
      IF ( .NOT. crtm_atmosphere_associated( atm_out ) ) THEN
        err_stat = failure
        msg = 'Error assigning Atmosphere structure with NO extra layers'
        CALL display_message( routine_name, msg, err_stat )
      END IF
      RETURN
    END IF
    
    
    ! Determine the number of extra layers required
    n = extra_layers( atm_in )
    IF ( n < 1 ) THEN
      err_stat = failure
      msg = 'Error determining extra layer count'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF
    
    
    ! Allocate the internal variable structure
    CALL iatm_create( iatm, n, atm_in%n_Absorbers )
    IF ( .NOT. iatm_associated( iatm ) ) THEN
      err_stat = failure
      msg = 'Error allocating iAtm internal structure'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF


    ! Get the extra layer profiles
    CALL crtm_get_model_profile( atm_in%Absorber_Id, &
                                 iatm%pl, iatm%tl, iatm%al, &
                                 model=atm_in%Climatology )


    ! First interpolate the extra levels to the user top pressure
    ! replacing the model data at that array index
    CALL interp_lpoly( atm_in%Level_Pressure(0), iatm%pl(n-1:n), iatm%ilpoly )
    iatm%plint_save = atm_in%Level_Pressure(0)
    iatm%pln_save   = iatm%pl(n)
    iatm%pl(n)      = iatm%plint_save
    CALL interp_linear( iatm%ilpoly, iatm%tl(n-1:n), iatm%tlint_save )
    iatm%tln_save = iatm%tl(n)
    iatm%tl(n)    = iatm%tlint_save
    DO j = 1, atm_in%n_Absorbers
      CALL interp_linear( iatm%ilpoly, iatm%al(n-1:n,j), iatm%alint_save(j) )
      iatm%aln_save(j) = iatm%al(n,j)
      iatm%al(n,j)     = iatm%alint_save(j)
    END DO
    
    ! Now compute the model profile layer averages
    DO k = 1, n
      CALL layer_p(iatm%pl(k-1:k), iatm%p(k))
      CALL layer_x(iatm%tl(k-1:k), iatm%t(k))
    END DO
    DO j = 1, atm_in%n_Absorbers
      DO k = 1, n
        CALL layer_x(iatm%al(k-1:k,j), iatm%a(k,j))
      END DO
    END DO
    
    
    ! Now, extrapolate user layer profile to get the "layer 0" value and
    ! use it to shift the model profile to avoid a discontinuity at p(n)
    CALL interp_lpoly( iatm%p(n), atm_in%Pressure(1:2), iatm%elpoly )
    CALL shift_profile( iatm%elpoly, atm_in%Temperature(1:2), iatm%t )
    DO j = 1, atm_in%n_Absorbers
      CALL shift_profile( iatm%elpoly, atm_in%Absorber(1:2,j), iatm%a(:,j) )
    END DO


    ! Make sure the absorber amounts are not negative.
    ! (Is a further, more physical, check needed here?)
    iatm%a_save = iatm%a
    WHERE (iatm%a_save < zero) iatm%a = minimum_absorber_amount
    

    ! Copy over the atmosphere structure with extra layers
    atm_out = crtm_atmosphere_addlayercopy( atm_in, n )
    IF ( .NOT. crtm_atmosphere_associated( atm_out ) ) THEN
      err_stat = failure
      msg = 'Error copying Atmosphere structure with extra layers'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF


    ! Slot the added layers into the output atmosphere structure
    ! Note: Cloud and Aerosol assignments not really needed (the
    !       zeroing is handled by the structure allocation) since
    !       at TOA, typically, there are not any clouds and/or
    !       aerosols.
    ! ...Profile
    atm_out%Level_Pressure(0:n) = iatm%pl
    atm_out%Pressure(1:n)       = iatm%p
    atm_out%Temperature(1:n)    = iatm%t
    DO j = 1, atm_out%n_Absorbers
      atm_out%Absorber(1:n,j)   = iatm%a(:,j)
    END DO
    ! ...Clouds
    IF ( atm_in%n_Clouds > 0 ) THEN
      DO i = 1, atm_in%n_Clouds
        atm_out%Cloud(i)%Effective_Radius(1:n)   = zero
        atm_out%Cloud(i)%Effective_Variance(1:n) = zero
        atm_out%Cloud(i)%Water_Content(1:n)      = zero
      END DO
    END IF
    ! ...Aerosols
    IF ( atm_in%n_Aerosols > 0 ) THEN
      DO i = 1, atm_in%n_Aerosols
        atm_out%Aerosol(i)%Effective_Radius(1:n) = zero
        atm_out%Aerosol(i)%Concentration(1:n)    = zero
      END DO
    END IF


    ! Clean up
    CALL iatm_destroy( iatm )
    
  END FUNCTION crtm_atmosphere_addlayers


!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Atmosphere_AddLayers_TL
!
! PURPOSE:
!       Function to copy a tangent-linear atmosphere structure and add extra
!       layers as required to supplement the upper atmosphere profile data.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Atmosphere_AddLayers_TL( Atm_In    , &  ! FWD Input
!                                                    Atm_In_TL , &  ! TL  Input
!                                                    Atm_Out_TL  )  ! TL  Output
!
! INPUTS:
!       Atm_In:          Forward model atmosphere structure.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Atm_In_TL:       Tangent-linear model atmosphere structure that is
!                        to be copied.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUTS:
!       Atm_Out_TL:      Copy of the input tangent-linear atmosphere structure
!                        with extra upper atmosphere layers added as required.
!                        Note that the tangent-linear values of the added layers
!                        is *always* zero.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        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 operation was successful
!                           == FAILURE an error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION crtm_atmosphere_addlayers_tl( &
    Atm_In    , &  ! FWD Input
    Atm_In_TL , &  ! TL  Input
    Atm_Out_TL) &  ! TL  Output
  result( err_stat )
    ! Arguments
    TYPE(crtm_atmosphere_type), INTENT(IN)     :: atm_in
    TYPE(crtm_atmosphere_type), INTENT(IN)     :: atm_in_tl
    TYPE(crtm_atmosphere_type), INTENT(IN OUT) :: atm_out_tl
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_Atmosphere_AddLayers_TL'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: n


    ! Set up
    err_stat = success


    ! If extra layers are NOT needed,
    ! then simply copy the structure
    IF ( atm_in%Level_Pressure(0) <= toa_pressure) THEN
      atm_out_tl = atm_in_tl
      IF ( .NOT. crtm_atmosphere_associated( atm_out_tl ) ) THEN
        err_stat = failure
        msg = 'Error assigning Atmosphere structure with NO extra layers'
        CALL display_message( routine_name, msg, err_stat )
      END IF
      RETURN
    END IF

    
    ! Determine how many extra layers are needed
    n = extra_layers( atm_in )
    IF ( n < 1 ) THEN
      err_stat = failure
      msg = 'Error determining extra layer count'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF


    ! Copy over the atmosphere structure with extra layers
    ! (which will be zero by definition)
    atm_out_tl = crtm_atmosphere_addlayercopy( atm_in_tl, n )
    IF ( .NOT. crtm_atmosphere_associated( atm_out_tl ) ) THEN
      err_stat = failure
      msg = 'Error copying Atmosphere structure with extra layers'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF
    
  END FUNCTION crtm_atmosphere_addlayers_tl
  
  
!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
!       CRTM_Atmosphere_AddLayers_AD
!
! PURPOSE:
!       Function to copy back an adjoint atmosphere structure removing added
!       extra layers as were required to supplement the upper atmosphere
!       profile data.
!
! CALLING SEQUENCE:
!       Error_Status = CRTM_Atmosphere_AddLayers_AD( Atm_In    , &  ! FWD Input
!                                                    Atm_Out_AD, &  ! AD  Input
!                                                    Atm_In_AD   )  ! AD  Output
!
! INPUTS:
!       Atm_In:          Forward model atmosphere structure.
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Atm_Out_AD:      Adjoint atmosphere structure that contains the added
!                        extra layers.
!                        ** SET TO ZERO ON EXIT ** 
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN OUT)
!
! OUTPUTS:
!       Atm_In_AD:       Adjoint atmosphere structure at the original, user
!                        specified layering.
!                        ** MUST HAVE VALUE ON ENTRY **
!                        UNITS:      N/A
!                        TYPE:       CRTM_Atmosphere_type
!                        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 operation was successful
!                           == FAILURE an error occurred
!                        UNITS:      N/A
!                        TYPE:       INTEGER
!                        DIMENSION:  Scalar
!
!:sdoc-:
!--------------------------------------------------------------------------------

  FUNCTION crtm_atmosphere_addlayers_ad( &
    Atm_In    , &  ! FWD Input
    Atm_Out_AD, &  ! AD  Input
    Atm_In_AD ) &  ! AD  Output
  result( err_stat )
    ! Arguments
    TYPE(crtm_atmosphere_type), INTENT(IN)     :: atm_in
    TYPE(crtm_atmosphere_type), INTENT(IN OUT) :: atm_out_ad
    TYPE(crtm_atmosphere_type), INTENT(IN OUT) :: atm_in_ad
    ! Function result
    INTEGER :: err_stat
    ! Local parameters
    CHARACTER(*), PARAMETER :: routine_name = 'CRTM_Atmosphere_AddLayers_AD'
    ! Local variables
    CHARACTER(ML) :: msg
    INTEGER :: i, j, n, no, nt


    ! Set up
    err_stat = success


    ! If extra layers are NOT needed, then simply perform
    ! the adjoint sum. Remember the TL form is
    !   Atm_Out_TL = Atm_In_TL
    ! so the adjoint form is
    !   Atm_In_AD  = Atm_In_AD + Atm_Out_AD
    !   Atm_Out_AD = ZERO
    IF ( atm_in%Level_Pressure(0) <= toa_pressure) THEN
      atm_in_ad = atm_in_ad + atm_out_ad
      CALL crtm_atmosphere_zero( atm_out_ad )
      RETURN
    END IF
    
    
    ! Determine how many extra layers have been used
    n = extra_layers( atm_in )
    IF ( n < 1 ) THEN
      err_stat = failure
      msg = 'Error determining extra layer count'
      CALL display_message( routine_name, msg, err_stat )
      RETURN
    END IF


    ! Perform the adjoint summations
    ! This is the adjoint equivalent of the TL Assign_Atmosphere
    no = atm_in_ad%n_Layers
    nt = n + no
    ! ...Aerosols
    IF ( atm_in_ad%n_Aerosols > 0 ) THEN
      DO i = 1, atm_in_ad%n_Aerosols
        atm_in_ad%Aerosol(i)%Concentration(1:no) = atm_in_ad%Aerosol(i)%Concentration(1:no) + &
                                                   atm_out_ad%Aerosol(i)%Concentration(n+1:nt)
        atm_in_ad%Aerosol(i)%Effective_Radius(1:no) = atm_in_ad%Aerosol(i)%Effective_Radius(1:no) + &
                                                      atm_out_ad%Aerosol(i)%Effective_Radius(n+1:nt)
        atm_in_ad%Aerosol(i)%Type = atm_out_ad%Aerosol(i)%Type
      END DO
    END IF
    ! ...Clouds    
    IF ( atm_in_ad%n_Clouds > 0 ) THEN
      DO i = 1, atm_in_ad%n_Clouds
        atm_in_ad%Cloud(i)%Water_Content(1:no) = atm_in_ad%Cloud(i)%Water_Content(1:no) + &
                                                 atm_out_ad%Cloud(i)%Water_Content(n+1:nt)
        atm_in_ad%Cloud(i)%Effective_Variance(1:no) = atm_in_ad%Cloud(i)%Effective_Variance(1:no) + &
                                                      atm_out_ad%Cloud(i)%Effective_Variance(n+1:nt)
        atm_in_ad%Cloud(i)%Effective_Radius(1:no) = atm_in_ad%Cloud(i)%Effective_Radius(1:no) + &
                                                    atm_out_ad%Cloud(i)%Effective_Radius(n+1:nt)
        atm_in_ad%Cloud(i)%Type = atm_out_ad%Cloud(i)%Type
      END DO
    END IF
    ! ...Absorber data
    DO j = 1, atm_in_ad%n_Absorbers
      atm_in_ad%Absorber(1:no,j) = atm_in_ad%Absorber(1:no,j) + atm_out_ad%Absorber(n+1:nt,j)
    END DO
    ! ...Temperature data
    atm_in_ad%Temperature(1:no) = atm_in_ad%Temperature(1:no) + atm_out_ad%Temperature(n+1:nt)
    ! ...Pressure data
    atm_in_ad%Pressure(1:no)       = atm_in_ad%Pressure(1:no) + atm_out_ad%Pressure(n+1:nt)
    atm_in_ad%Level_Pressure(0:no) = atm_in_ad%Level_Pressure(0:no) + atm_out_ad%Level_Pressure(n:nt)


    ! Zero the output atmosphere structure
    CALL crtm_atmosphere_zero( atm_out_ad )

  END FUNCTION crtm_atmosphere_addlayers_ad

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

  ! Subprogram to determine the number of extra layers required
  FUNCTION extra_layers( Atm ) RESULT( n )
    TYPE(crtm_atmosphere_type), INTENT(IN) :: atm
    INTEGER :: n
    DO n = 1, SIZE(model_level_pressure)
      IF ( model_level_pressure(n) >= atm%Level_Pressure(0) ) RETURN
    END DO
    n = 0
  END FUNCTION extra_layers


  ! Subprogram to compute the average layer pressure
  SUBROUTINE layer_p( p, p_layer )
    REAL(fp), INTENT(IN)  :: p(2)     ! Input
    REAL(fp), INTENT(OUT) :: p_layer  ! Output
    p_layer = (p(2)-p(1))/log(p(2)/p(1))
  END SUBROUTINE layer_p


  ! Subprogram to compute the average layer amount of X
  SUBROUTINE layer_x( x, x_layer )
    REAL(fp), INTENT(IN)  :: x(2)
    REAL(fp), INTENT(OUT) :: x_layer
    x_layer = point_5*(x(1)+x(2))
  END SUBROUTINE layer_x
  

  ! Subprogram to compute the interpolating polynomial linear in log(p)
  SUBROUTINE interp_lpoly( p_int, p, lpoly )
    REAL(fp), INTENT(IN)  :: p_int
    REAL(fp), INTENT(IN)  :: p(2)
    REAL(fp), INTENT(OUT) :: lpoly
    lpoly = (log(p_int)-log(p(1))) / (log(p(2))-log(p(1)))
  END SUBROUTINE interp_lpoly


  ! Subprogram to perform linear interpolation
  SUBROUTINE interp_linear( lpoly, x, x_int )
    REAL(fp), INTENT(IN)  :: lpoly
    REAL(fp), INTENT(IN)  :: x(2)
    REAL(fp), INTENT(OUT) :: x_int
    x_int = (x(2)-x(1))*lpoly + x(1)
  END SUBROUTINE interp_linear


  ! Subprogram to shifted the added profile layers
  SUBROUTINE shift_profile( lpoly, x_toa, x_shifted )
    REAL(fp), INTENT(IN)     :: lpoly
    REAL(fp), INTENT(IN)     :: x_toa(2)
    REAL(fp), INTENT(IN OUT) :: x_shifted(:)
    INTEGER :: n
    REAL(fp) :: x_int, dx
    n = SIZE(x_shifted)
    CALL interp_linear( lpoly, x_toa, x_int )
    dx = x_int - x_shifted(n)
    x_shifted = x_shifted + dx
  END SUBROUTINE shift_profile
    
END MODULE crtm_atmosphere