!
! CRTM_AerosolScatter
!
!
! Module to compute the aerosol absorption and scattering properties
! required for radiative transfer in an atmosphere with aerosols.
!
!
! CREATION HISTORY:
! Written by: Paul van Delst, 15-Feb-2005
! paul.vandelst@noaa.gov
! Modified by Quanhua Liu, 03-Oct-2006
! Quanhua.Liu@noaa.gov
!
MODULE CRTM_AerosolScatter
! -----------------
! Environment setup
! -----------------
! Module use
USE Type_Kinds, ONLY: fp
USE Message_Handler, ONLY: SUCCESS, FAILURE, Display_Message
USE CRTM_Parameters, ONLY: ZERO, ONE, POINT_5, ONEpointFIVE, &
MAX_N_LAYERS, &
MAX_N_AEROSOLS, &
AEROSOL_CONTENT_THRESHOLD, &
BS_THRESHOLD, &
MAX_N_LEGENDRE_TERMS, &
MAX_N_PHASE_ELEMENTS, &
HGPHASE ! <<< NEED TO REMOVE THIS IN FUTURE
USE CRTM_SpcCoeff, ONLY: SC, &
SpcCoeff_IsMicrowaveSensor , &
SpcCoeff_IsInfraredSensor , &
SpcCoeff_IsVisibleSensor , &
SpcCoeff_IsUltravioletSensor
USE CRTM_AerosolCoeff, ONLY: AeroC
USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type , &
DUST_AEROSOL , &
SEASALT_SSAM_AEROSOL , &
SEASALT_SSCM1_AEROSOL , &
SEASALT_SSCM2_AEROSOL , &
SEASALT_SSCM3_AEROSOL , &
ORGANIC_CARBON_AEROSOL , &
BLACK_CARBON_AEROSOL , &
SULFATE_AEROSOL
USE CRTM_GeometryInfo_Define, ONLY: CRTM_GeometryInfo_type
USE CRTM_Interpolation, ONLY: NPTS , &
LPoly_type , &
find_index , &
interp_1D , &
interp_2D , &
interp_3D , &
interp_2D_TL, &
interp_3D_TL, &
interp_2D_AD, &
interp_3D_AD, &
Clear_LPoly , &
LPoly , &
LPoly_TL , &
LPoly_AD
USE CRTM_AtmOptics_Define, ONLY: CRTM_AtmOptics_type
! Internal variable definition module
USE ASvar_Define, ONLY: ASvar_type, &
ASinterp_type, &
ASvar_Associated, &
ASvar_Destroy , &
ASvar_Create
! Disable implicit typing
IMPLICIT NONE
! ------------
! Visibilities
! ------------
! Everything private by default
PRIVATE
! Procedures
PUBLIC :: CRTM_Compute_AerosolScatter
PUBLIC :: CRTM_Compute_AerosolScatter_TL
PUBLIC :: CRTM_Compute_AerosolScatter_AD
! -----------------
! Module parameters
! -----------------
! Version Id for the module
CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
'$Id: CRTM_AerosolScatter.f90 99117 2017-11-27 18:37:14Z tong.zhu@noaa.gov $'
! Message string length
INTEGER, PARAMETER :: ML = 256
! Number of stream angle definitions
INTEGER, PARAMETER :: TWO_STREAMS = 2
INTEGER, PARAMETER :: FOUR_STREAMS = 4
INTEGER, PARAMETER :: SIX_STREAMS = 6
INTEGER, PARAMETER :: EIGHT_STREAMS = 8
INTEGER, PARAMETER :: SIXTEEN_STREAMS = 16
INTEGER, PARAMETER :: THIRTYTWO_STREAMS = 32
CONTAINS
!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
! CRTM_Compute_AerosolScatter
!
! PURPOSE:
! Function to compute the aerosol absorption and scattering properties
! and populate the output AerosolScatter structure for a single channel.
!
! CALLING SEQUENCE:
! Error_Status = CRTM_Compute_AerosolScatter( Atmosphere , &
! SensorIndex , &
! ChannelIndex , &
! AerosolScatter, &
! ASvar )
!
! INPUT ARGUMENTS:
! Atmosphere: CRTM_Atmosphere structure containing the atmospheric
! profile data.
! UNITS: N/A
! TYPE: CRTM_Atmosphere_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.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! SensorIndex: Sensor index id. This ia a unique index associated
! with a sensor.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
! AerosolScatter: CRTM_AtmOptics structure containing the aerosol
! absorption and scattering properties required by
! the radiative transfer.
! UNITS: N/A
! TYPE: CRTM_AtmOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN OUT)
!
! ASvar: Structure containing internal variables required for
! subsequent tangent-linear or adjoint model calls.
! UNITS: N/A
! TYPE: ASvar_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(OUT)
!
! FUNCTION RESULT:
! Error_Status: The return value is an integer defining the error status.
! The error codes are defined in the ERROR_HANDLER module.
! If == SUCCESS the computation was sucessful
! == FAILURE an unrecoverable error occurred
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------
FUNCTION CRTM_Compute_AerosolScatter( &
Atm , & ! Input
SensorIndex , & ! Input
ChannelIndex, & ! Input
AScat , & ! Output
ASV ) & ! Internal variable output
RESULT( Error_Status )
!Arguments
TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm
INTEGER , INTENT(IN) :: ChannelIndex
INTEGER , INTENT(IN) :: SensorIndex
TYPE(CRTM_AtmOptics_type) , INTENT(IN OUT) :: AScat
TYPE(ASvar_type) , INTENT(IN OUT) :: ASV
! Function result
INTEGER :: Error_Status
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Compute_AerosolScatter'
! Local Variables
CHARACTER(ML) :: Message
INTEGER :: k, ka, l, m, n
REAL(fp) :: Frequency
LOGICAL :: Layer_Mask(Atm%n_Layers)
INTEGER :: Layer_Index(Atm%n_Layers)
INTEGER :: nAerosol_Layers
REAL(fp) :: bs
! ------
! Set up
! ------
Error_Status = SUCCESS
IF ( SpcCoeff_IsMicrowaveSensor(SC(SensorIndex)) ) RETURN
IF ( Atm%n_Aerosols == 0 ) RETURN
ASV%Total_bs = ZERO
! Frequency in inverse centimetres
Frequency = SC(SensorIndex)%Wavenumber(ChannelIndex)
! Determine offset for Legendre coefficients in the
! AeroC lookup table corresponding to the
! number of streams
SELECT CASE(AScat%n_Legendre_Terms)
CASE (TWO_STREAMS) ; AScat%lOffset = 0
CASE (FOUR_STREAMS) ; AScat%lOffset = 0
CASE (SIX_STREAMS) ; AScat%lOffset = 5
CASE (EIGHT_STREAMS) ; AScat%lOffset = 12
CASE (SIXTEEN_STREAMS); AScat%lOffset = 21
CASE DEFAULT
AScat%lOffset = 0
! Use two-stream model or HG and RAYLEIGH Phase function
IF( HGPHASE ) THEN
AScat%n_Legendre_Terms = 0
ELSE
Error_Status = FAILURE
WRITE(Message,'("The n_Legendre_Terms in AerosolScatter, ",i0,", do not fit model")') &
AScat%n_Legendre_Terms
CALL Display_Message(ROUTINE_NAME, Message, Error_Status )
RETURN
END IF
END SELECT
! -----------------------------------------------
! Loop over the different Aerosols in the profile
! -----------------------------------------------
Aerosol_loop: DO n = 1, Atm%n_Aerosols
! Only process aerosols with more
! than the threshold Aerosol amount
Layer_Mask = Atm%Aerosol(n)%Concentration > AEROSOL_CONTENT_THRESHOLD
nAerosol_Layers = COUNT(Layer_Mask)
IF ( nAerosol_Layers == 0 ) CYCLE Aerosol_loop
! --------------------------------------
! Loop over the current Aerosol's Layers
! --------------------------------------
Layer_Index(1:nAerosol_Layers) = PACK((/(k,k=1,Atm%Aerosol(n)%n_Layers)/), Layer_Mask)
Aerosol_Layer_loop: DO k = 1, nAerosol_Layers
ka = Layer_Index(k)
! Get the aerosol optical properties
CALL Get_Aerosol_Opt( AScat , & ! Input
Frequency , & ! Input
Atm%Aerosol(n)%Type , & ! Input
Atm%Aerosol(n)%Effective_Radius(ka), & ! Input
ASV%ke(ka,n) , & ! Output
ASV%w(ka,n) , & ! Output
ASV%pcoeff(:,:,ka,n) , & ! Output
ASV%asi(ka,n) ) ! Interpolation
! interpolation quality control
IF( ASV%ke(ka,n) <= ZERO ) THEN
ASV%ke(ka,n) = ZERO
ASV%w(ka,n) = ZERO
END IF
IF( ASV%w(ka,n) <= ZERO ) THEN
ASV%w(ka,n) = ZERO
ASV%pcoeff(:,:,ka,n) = ZERO
END IF
IF( ASV%w(ka,n) >= ONE ) THEN
ASV%w(ka,n) = ONE
END IF
! Compute the optical depth (absorption + scattering)
! tau = rho.ke
! where
! rho = Integrated Aerosol Concentration for a layer(kg/m^2) [M.L^-2]
! ke = mass extintion coefficient (m^2/kg) [L^2.M^-1]
! Note that since all these computations are done for a given
! layer, the optical depth is the same as the volume extinction
! coefficient, be. Usually,
! tau = be.d(z)
! but we are working with height/thickness independent quantities
! so that
! tau = be
! This is why the optical depth is used in the denominator to
! compute the single scatter albedo in the Layer_loop below.
AScat%Optical_Depth(ka) = AScat%Optical_Depth(ka) + &
(ASV%ke(ka,n)*Atm%Aerosol(n)%Concentration(ka))
! Compute the phase matrix coefficients
! p = p + p(LUT)*bs
! where
! p(LUT) = the phase coefficient from the LUT
IF( AScat%n_Phase_Elements > 0 .and. AScat%Include_Scattering ) THEN
! Compute the volume scattering coefficient for the current
! aerosol layer and accumulate it for the layer total for the
! profile (i.e. all aerosols)
! bs = rho.w.ke
! where
! bs = volume scattering coefficient for a layer [dimensionless]
! rho = integrated aerosol concentration for a layer (kg/m^2) [M.L^-2]
! w = single scatter albedo [dimensionless]
! ke = mass extintion coefficient (m^2/kg) [L^2.M^-1]
! qliu
bs = Atm%Aerosol(n)%Concentration(ka) * ASV%ke(ka,n) * ASV%w(ka,n)
ASV%Total_bs(ka) = ASV%Total_bs(ka) + bs
AScat%Single_Scatter_Albedo(ka) = AScat%Single_Scatter_Albedo(ka) + bs
DO m = 1, AScat%n_Phase_Elements
DO l = 0, AScat%n_Legendre_Terms
AScat%Phase_Coefficient(l,m,ka) = AScat%Phase_Coefficient(l,m,ka) + &
(ASV%pcoeff(l,m,ka,n) * bs)
END DO
END DO
END IF
END DO Aerosol_Layer_loop
END DO Aerosol_loop
END FUNCTION CRTM_Compute_AerosolScatter
!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
! CRTM_Compute_AerosolScatter_TL
!
! PURPOSE:
! Function to compute the tangent-linear aerosol absorption and
! scattering properties and populate the output AerosolScatter_TL
! structure for a single channel.
!
! CALLING SEQUENCE:
! Error_Status = CRTM_Compute_AerosolScatter_TL( Atmosphere , &
! AerosolScatter , &
! Atmosphere_TL , &
! SensorIndex , &
! ChannelIndex , &
! AerosolScatter_TL, &
! ASvar )
!
! INPUT ARGUMENTS:
! Atmosphere: CRTM_Atmosphere structure containing the atmospheric
! profile data.
! UNITS: N/A
! TYPE: CRTM_Atmosphere_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! AerosolScatter: CRTM_AtmOptics structure containing the forward model
! aerosol absorption and scattering properties required
! for radiative transfer.
! UNITS: N/A
! TYPE: CRTM_AtmOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! Atmosphere_TL: CRTM Atmosphere structure containing the tangent-linear
! atmospheric state data.
! UNITS: N/A
! TYPE: CRTM_Atmosphere_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! SensorIndex: Sensor index id. This ia a unique index associated
! with a sensor.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
!
! Channel_Index: Channel index id. This is a unique index associated
! with a (supported) sensor channel used to access the
! shared coefficient data.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! ASvar: Structure containing internal variables required for
! subsequent tangent-linear or adjoint model calls.
! UNITS: N/A
! TYPE: ASvar_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
! AerosolScatter_TL: CRTM_AtmOptics structure containing the tangent-linear
! aerosol absorption and scattering properties required
! for radiative transfer.
! UNITS: N/A
! TYPE: CRTM_AtmOptics_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 ERROR_HANDLER module.
! If == SUCCESS the computation was sucessful
! == FAILURE an unrecoverable error occurred
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------
FUNCTION CRTM_Compute_AerosolScatter_TL( &
Atm , & ! FWD Input
AScat , & ! FWD Input
Atm_TL , & ! TL Input
SensorIndex , & ! Input
ChannelIndex, & ! Input
AScat_TL , & ! TL Input
ASV ) & ! Internal variable input
RESULT( Error_Status )
! Arguments
TYPE(CRTM_Atmosphere_type) , INTENT(IN) :: Atm
TYPE(CRTM_AtmOptics_type) , INTENT(IN) :: AScat
TYPE(CRTM_Atmosphere_type) , INTENT(IN) :: Atm_TL
INTEGER , INTENT(IN) :: SensorIndex
INTEGER , INTENT(IN) :: ChannelIndex
TYPE(CRTM_AtmOptics_type) , INTENT(IN OUT) :: AScat_TL
TYPE(ASvar_type) , INTENT(IN) :: ASV
! Function result
INTEGER :: Error_Status
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Compute_AerosolScatter_TL'
! Local variables
INTEGER :: k, ka, l, m, n
INTEGER :: n_Legendre_Terms, n_Phase_Elements
REAL(fp) :: Frequency
LOGICAL :: Layer_Mask(Atm%n_Layers)
INTEGER :: Layer_Index(Atm%n_Layers)
INTEGER :: nAerosol_Layers
REAL(fp) :: ke_TL, w_TL
REAL(fp) :: pcoeff_TL(0:AScat%n_Legendre_Terms, AScat%n_Phase_Elements)
REAL(fp) :: bs, bs_TL
! ------
! Set up
! ------
Error_Status = SUCCESS
IF ( SpcCoeff_IsMicrowaveSensor(SC(SensorIndex)) ) RETURN
IF ( Atm%n_Aerosols == 0 ) RETURN
! Frequency
Frequency = SC(SensorIndex)%Wavenumber(ChannelIndex)
! Phase matrix dimensions
n_Legendre_Terms = AScat_TL%n_Legendre_Terms
n_Phase_Elements = AScat_TL%n_Phase_Elements
AScat_TL%lOffset = AScat%lOffset
! -----------------------------------------------
! Loop over the different Aerosols in the profile
! -----------------------------------------------
Aerosol_loop: DO n = 1, Atm%n_Aerosols
! Only process aerosols with more
! than the threshold aerosol amount
Layer_Mask = Atm%Aerosol(n)%Concentration > AEROSOL_CONTENT_THRESHOLD
nAerosol_Layers = COUNT(Layer_Mask)
IF (nAerosol_Layers == 0) CYCLE Aerosol_loop
! --------------------------------------
! Loop over the current aerosol's layers
! --------------------------------------
Layer_Index(1:nAerosol_Layers) = PACK((/(k,k=1,Atm%Aerosol(n)%n_Layers)/), Layer_Mask)
Aerosol_Layer_loop: DO k = 1, nAerosol_Layers
ka = Layer_Index(k)
! Obtain bulk aerosol optical properties
CALL Get_Aerosol_Opt_TL(AScat_TL , & ! Input
Atm%Aerosol(n)%Type , & ! Input
ASV%ke(ka,n) , & ! Input
ASV%w(ka,n) , & ! Input
Atm_TL%Aerosol(n)%Effective_Radius(ka), & ! TL Input
ke_TL , & ! TL Output
w_TL , & ! TL Output
pcoeff_TL , & ! TL Output
ASV%asi(ka,n) ) ! Interpolation
! interpolation quality control
IF( ASV%ke(ka,n) <= ZERO ) THEN
ke_TL = ZERO
w_TL = ZERO
END IF
IF( ASV%w(ka,n) <= ZERO ) THEN
w_TL = ZERO
pcoeff_TL = ZERO
END IF
IF( ASV%w(ka,n) >= ONE ) THEN
w_TL = ZERO
END IF
! Compute the optical depth (absorption + scattering)
AScat_TL%Optical_Depth(ka) = AScat_TL%Optical_Depth(ka) + &
(ke_TL * Atm%Aerosol(n)%Concentration(ka)) + &
(ASV%ke(ka,n) * Atm_TL%Aerosol(n)%Concentration(ka))
! Compute the phase matrix coefficients
IF( n_Phase_Elements > 0 .and. AScat%Include_Scattering ) THEN
! Compute the volume scattering coefficient
bs = Atm%Aerosol(n)%Concentration(ka) * ASV%ke(ka,n) * ASV%w(ka,n)
bs_TL = (Atm_TL%Aerosol(n)%Concentration(ka) * ASV%ke(ka,n) * ASV%w(ka,n) ) + &
(Atm%Aerosol(n)%Concentration(ka) * ke_TL * ASV%w(ka,n) ) + &
(Atm%Aerosol(n)%Concentration(ka) * ASV%ke(ka,n) * w_TL )
AScat_TL%Single_Scatter_Albedo(ka) = AScat_TL%Single_Scatter_Albedo(ka) + bs_TL
DO m = 1, n_Phase_Elements
DO l = 0, n_Legendre_Terms
AScat_TL%Phase_Coefficient(l,m,ka) = AScat_TL%Phase_Coefficient(l,m,ka) + &
(pcoeff_TL(l,m) * bs ) + &
(ASV%pcoeff(l,m,ka,n) * bs_TL)
END DO
END DO
END IF
END DO Aerosol_Layer_loop
END DO Aerosol_loop
END FUNCTION CRTM_Compute_AerosolScatter_TL
!------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
! CRTM_Compute_AerosolScatter_AD
!
! PURPOSE:
! Function to compute the adjoint aerosol absorption and scattering
! properties for a single channel.
!
! CALLING SEQUENCE:
! Error_Status = CRTM_Compute_AerosolScatter_AD( Atmosphere , &
! AerosolScatter , &
! AerosolScatter_AD, &
! SensorIndex , &
! ChannelIndex , &
! Atmosphere_AD , &
! ASVariables )
!
! INPUT ARGUMENTS:
! Atmosphere: CRTM_Atmosphere structure containing the atmospheric
! profile data.
! UNITS: N/A
! TYPE: CRTM_Atmosphere_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! AerosolScatter: CRTM_AtmOptics structure containing the forward model
! aerosol absorption and scattering properties required
! for radiative transfer.
! UNITS: N/A
! TYPE: CRTM_AtmOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! AerosolScatter_AD: CRTM_AtmOptics structure containing the adjoint
! aerosol absorption and scattering properties.
! **NOTE: On EXIT from this function, the contents of
! this structure may be modified (e.g. set to
! zero.)
! UNITS: N/A
! TYPE: CRTM_AtmOptics_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN OUT)
!
! SensorIndex: Sensor index id. This ia a unique index associated
! with a sensor.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! Channel_Index: Channel index id. This is a unique index associated
! with a (supported) sensor channel used to access the
! shared coefficient data.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! ASvar: Structure containing internal variables required for
! subsequent tangent-linear or adjoint model calls.
! UNITS: N/A
! TYPE: ASvar_type
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
! Atmosphere_AD: CRTM Atmosphere structure containing the adjoint
! atmospheric state data.
! 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
! ERROR_HANDLER module.
! If == SUCCESS the computation was sucessful
! == FAILURE an unrecoverable error occurred
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
!
!:sdoc-:
!------------------------------------------------------------------------------
FUNCTION CRTM_Compute_AerosolScatter_AD( &
Atm , & ! FWD Input
AScat , & ! FWD Input
AScat_AD , & ! AD Input
SensorIndex , & ! Input
ChannelIndex, & ! Input
Atm_AD , & ! AD Output
ASV ) & ! Internal Variable input
RESULT( Error_Status )
! Arguments
TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm
TYPE(CRTM_AtmOptics_type) , INTENT(IN) :: AScat
TYPE(CRTM_AtmOptics_type) , INTENT(IN OUT) :: AScat_AD
INTEGER , INTENT(IN) :: SensorIndex
INTEGER , INTENT(IN) :: ChannelIndex
TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_AD
TYPE(ASvar_type) , INTENT(IN) :: ASV
! Function result
INTEGER :: Error_Status
! Local parameters
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Compute_AerosolScatter_AD'
! Local variables
INTEGER :: k, ka, l, m, n
INTEGER :: n_Legendre_Terms, n_Phase_Elements
REAL(fp) :: Frequency
LOGICAL :: Layer_Mask(Atm%n_Layers)
INTEGER :: Layer_Index(Atm%n_Layers)
INTEGER :: nAerosol_Layers
REAL(fp) :: ke_AD, w_AD
REAL(fp) :: pcoeff_AD(0:AScat%n_Legendre_Terms, AScat%n_Phase_Elements)
REAL(fp) :: bs, bs_AD
! ------
! Set up
! ------
Error_Status = SUCCESS
IF ( SpcCoeff_IsMicrowaveSensor(SC(SensorIndex)) ) RETURN
IF ( Atm%n_Aerosols == 0 ) RETURN
! Frequency
Frequency = SC(SensorIndex)%Wavenumber(ChannelIndex)
! Phase matrix dimensions
n_Legendre_Terms = AScat_AD%n_Legendre_Terms
n_Phase_Elements = AScat_AD%n_Phase_Elements
AScat_AD%lOffset = AScat%lOffset
! ----------------------------------------------------------
! Adjoint of accumulated optical properties for all aerosols
! ----------------------------------------------------------
! Shorten variable name
l = n_Legendre_Terms
! ------------------------------------------------
! Loop over different types of aerosols in profile
! ------------------------------------------------
Aerosol_loop: DO n = 1, Atm%n_Aerosols
! Only process aerosols with more than
! the threshold aerosol concentration
Layer_Mask = Atm%Aerosol(n)%Concentration > AEROSOL_CONTENT_THRESHOLD
nAerosol_Layers = COUNT(Layer_Mask)
IF ( nAerosol_Layers == 0 ) CYCLE Aerosol_loop
! --------------------------------------
! Loop over the current aerosol's layers
! --------------------------------------
Layer_Index(1:nAerosol_Layers) = PACK((/(k,k=1,Atm%Aerosol(n)%n_Layers)/), Layer_Mask)
Aerosol_Layer_loop: DO k = 1, nAerosol_Layers
ka = Layer_Index(k)
! Initialize the individual
! Aerosol adjoint variables
bs_AD = ZERO
pcoeff_AD = ZERO
ke_AD = ZERO
w_AD = ZERO
! Compute the adjoint of the
! phase matrix coefficients
IF( n_Phase_Elements > 0 .and. AScat%Include_Scattering ) THEN
! Recompute the forward model volume scattering
! coefficient for the current aerosol type ONLY
bs = Atm%Aerosol(n)%Concentration(ka) * ASV%ke(ka,n) * ASV%w(ka,n)
DO m = 1, n_Phase_Elements
DO l = 0, n_Legendre_Terms
bs_AD = bs_AD + (ASV%pcoeff(l,m,ka,n) * AScat_AD%Phase_Coefficient(l,m,ka))
pcoeff_AD(l,m) = pcoeff_AD(l,m) + (bs * AScat_AD%Phase_Coefficient(l,m,ka))
END DO
END DO
! NOTE: bs_AD is not reinitialized after this
! point since it is reinitialized at the
! start of the Aerosol_Layer_loop
bs_AD = bs_AD + AScat_AD%Single_Scatter_Albedo(ka)
w_AD = w_AD + (Atm%Aerosol(n)%Concentration(ka) * ASV%ke(ka,n)* bs_AD )
END IF
! Compute the adjoint of the optical
! depth (absorption + scattering)
Atm_AD%Aerosol(n)%Concentration(ka) = Atm_AD%Aerosol(n)%Concentration(ka) + &
(ASV%ke(ka,n) * AScat_AD%Optical_Depth(ka))
ke_AD = ke_AD + (Atm%Aerosol(n)%Concentration(ka) * AScat_AD%Optical_Depth(ka))
! Compute the adjoint of the volume
! scattering coefficient.
ke_AD = ke_AD + (Atm%Aerosol(n)%Concentration(ka) * bs_AD * ASV%w(ka,n) )
Atm_AD%Aerosol(n)%Concentration(ka) = Atm_AD%Aerosol(n)%Concentration(ka) + &
( bs_AD * ASV%ke(ka,n) * ASV%w(ka,n) )
! interpolation quality control
IF( ASV%w(ka,n) >= ONE ) THEN
w_AD = ZERO
END IF
IF( ASV%ke(ka,n) <= ZERO ) THEN
ke_AD = ZERO
w_AD = ZERO
END IF
IF( ASV%w(ka,n) <= ZERO ) THEN
w_AD = ZERO
pcoeff_AD = ZERO
END IF
! interpolation quality control
IF( ASV%w(ka,n) >= ONE ) THEN
w_AD = ZERO
END IF
IF( ASV%ke(ka,n) <= ZERO ) THEN
ke_AD = ZERO
w_AD = ZERO
END IF
IF( ASV%w(ka,n) <= ZERO ) THEN
w_AD = ZERO
pcoeff_AD = ZERO
END IF
! Adjoint AScat interpolation routine
CALL Get_Aerosol_Opt_AD(AScat_AD , & ! Input
Atm%Aerosol(n)%Type , & ! Input
ASV%ke(ka,n) , & ! Output
ASV%w(ka,n) , & ! Output
ke_AD , & ! AD Input
w_AD , & ! AD Input
pcoeff_AD , & ! AD Input
Atm_AD%Aerosol(n)%Effective_Radius(ka), & ! AD Output
ASV%asi(ka,n) ) ! Interpolation
END DO Aerosol_Layer_loop
END DO Aerosol_loop
END FUNCTION CRTM_Compute_AerosolScatter_AD
!################################################################################
!################################################################################
!## ##
!## ## PRIVATE MODULE ROUTINES ## ##
!## ##
!################################################################################
!################################################################################
! --------------------------------------------
! Determine the aerosol type index, k, for the
! aerosols based on AeroC LUT organisation
! --------------------------------------------
FUNCTION Aerosol_Type_Index( Aerosol_Type ) RESULT( k )
INTEGER, INTENT(IN) :: Aerosol_Type
INTEGER :: k
SELECT CASE (Aerosol_Type)
CASE(DUST_AEROSOL) ; k=1
CASE(SEASALT_SSAM_AEROSOL) ; k=2
CASE(SEASALT_SSCM1_AEROSOL) ; k=3
CASE(SEASALT_SSCM2_AEROSOL) ; k=4
CASE(SEASALT_SSCM3_AEROSOL) ; k=5
CASE(ORGANIC_CARBON_AEROSOL) ; k=6
CASE(BLACK_CARBON_AEROSOL) ; k=7
CASE(SULFATE_AEROSOL) ; k=8
END SELECT
END FUNCTION Aerosol_Type_Index
! ----------------------------------------
! Subroutine to obtain the bulk
! optical properties of aerosol
! extinction coefficient (ke),
! scattering coefficient (w)
! asymmetry factor (g), and
! spherical Legendre coefficients (pcoeff)
! ----------------------------------------
SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure
Frequency , & ! Input in cm^-1
Aerosol_Type , & ! Input see CRTM_Aerosol_Define.f90
Reff , & ! Input effective radius (mm)
ke , & ! Output extinction coefficient (=~ optical depth)
w , & ! Output single scattering albedo
pcoeff , & ! Output spherical Legendre coefficients
asi ) ! Output interpolation data
! Arguments
TYPE(CRTM_AtmOptics_type), INTENT(IN) :: AerosolScatter
REAL(fp) , INTENT(IN) :: Frequency
INTEGER , INTENT(IN) :: Aerosol_Type
REAL(fp) , INTENT(IN) :: Reff
REAL(fp) , INTENT(OUT) :: ke
REAL(fp) , INTENT(OUT) :: w
REAL(fp) , INTENT(IN OUT) :: pcoeff(0:,:)
TYPE(ASinterp_type) , INTENT(IN OUT) :: asi
! Local variables
INTEGER :: k, l, m
! Get the aerosol type LUT index
! ------------------------------
k = Aerosol_Type_Index( Aerosol_Type )
! Find the Frequency and effective
! radius indices for interpolation
! --------------------------------
asi%f_int = MAX(MIN(AeroC%Frequency(AeroC%n_Wavelengths),Frequency),AeroC%Frequency(1))
CALL find_index(AeroC%Frequency(:),asi%f_int,asi%i1,asi%i2, asi%f_outbound)
asi%f = AeroC%Frequency(asi%i1:asi%i2)
asi%r_int = MAX(MIN(AeroC%Reff(AeroC%n_Radii,k),Reff),AeroC%Reff(1,k))
CALL find_index(AeroC%Reff(:,k), asi%r_int, asi%j1,asi%j2, asi%r_outbound)
asi%r = AeroC%Reff(asi%j1:asi%j2,k)
! Calculate the interpolating polynomials
! ---------------------------------------
! Frequency term
CALL LPoly( asi%f, asi%f_int, & ! Input
asi%wlp ) ! Output
! Effective radius term
CALL LPoly( asi%r, asi%r_int, & ! Input
asi%xlp ) ! Output
! Perform Interpolation
! ---------------------
CALL interp_2D( AeroC%ke(asi%i1:asi%i2,asi%j1:asi%j2,k), asi%wlp, asi%xlp, ke )
CALL interp_2D( AeroC%w(asi%i1:asi%i2,asi%j1:asi%j2,k) , asi%wlp, asi%xlp, w )
IF (AerosolScatter%n_Phase_Elements > 0 .and. AerosolScatter%Include_Scattering ) THEN
pcoeff(0,1) = POINT_5
DO m = 1, AerosolScatter%n_Phase_Elements
DO l = 1, AerosolScatter%n_Legendre_Terms
CALL interp_2D( AeroC%pcoeff(asi%i1:asi%i2,asi%j1:asi%j2,k,l+AerosolScatter%lOffset,m), &
asi%wlp, asi%xlp, pcoeff(l,m) )
END DO
END DO
ELSE
! Absorption coefficient
ke = ke * (ONE- w)
END IF
END SUBROUTINE Get_Aerosol_Opt
! ---------------------------------------------
! Subroutine to obtain the tangent-linear
! bulk optical properties of a aerosol:
! extinction coefficient (ke_TL),
! scattereing coefficient (w_TL)
! asymmetry factor (g_TL), and
! spherical Legendre coefficients (pcoeff_TL)
! ---------------------------------------------
SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL structure
Aerosol_Type , & ! Input see CRTM_Aerosol_Define.f90
ke , & ! Input
w , & ! Input
Reff_TL , & ! Input TL effective radius (mm)
ke_TL , & ! Output TL extinction coefficient (=~ optical depth)
w_TL , & ! Output TL single scattering albedo
pcoeff_TL , & ! TL Output spherical Legendre coefficients
asi ) ! Input interpolation data
! Arguments
TYPE(CRTM_AtmOptics_type), INTENT(IN) :: AerosolScatter_TL
INTEGER , INTENT(IN) :: Aerosol_Type
REAL(fp), INTENT(IN) :: w, ke, Reff_TL
REAL(fp), INTENT(OUT) :: ke_TL
REAL(fp), INTENT(OUT) :: w_TL
REAL(fp), INTENT(IN OUT) :: pcoeff_TL(0:,:)
TYPE(ASinterp_type), INTENT(IN) :: asi
! Local variables
INTEGER :: k, l, m
REAL(fp) :: f_int_TL, r_int_TL
REAL(fp) :: f_TL(NPTS), r_TL(NPTS)
REAL(fp) :: z_TL(NPTS,NPTS)
TYPE(LPoly_type) :: wlp_TL, xlp_TL
REAL(fp), POINTER :: z(:,:) => NULL()
! Setup
! -----
! No TL output when all dimensions
! are outside LUT bounds
IF ( asi%f_outbound .AND. asi%r_outbound ) THEN
ke_TL = ZERO
w_TL = ZERO
pcoeff_TL = ZERO
RETURN
END IF
! The TL inputs
f_int_TL = ZERO
f_TL = ZERO
r_int_TL = Reff_TL
r_TL = ZERO
z_TL = ZERO
! Calculate the TL interpolating polynomials
! ------------------------------------------
! Frequency term (always zero. This is a placeholder for testing)
CALL LPoly_TL( asi%f, asi%f_int, & ! FWD Input
asi%wlp, & ! FWD Input
f_TL, f_int_TL, & ! TL Input
wlp_TL ) ! TL Output
! Effective radius term
CALL LPoly_TL( asi%r, asi%r_int, & ! FWD Input
asi%xlp, & ! FWD Input
r_TL, r_int_TL, & ! TL Input
xlp_TL ) ! TL Output
! Get the aerosol type LUT index
! ------------------------------
k = Aerosol_Type_Index( Aerosol_Type )
! Perform Interpolation
! ---------------------
! Extinction coefficient
z => AeroC%ke(asi%i1:asi%i2,asi%j1:asi%j2,k)
CALL interp_2D_TL( z , asi%wlp, asi%xlp, & ! FWD Input
z_TL, wlp_TL , xlp_TL , & ! TL Input
ke_TL ) ! TL Output
! Single scatter albedo
z => AeroC%w(asi%i1:asi%i2,asi%j1:asi%j2,k)
CALL interp_2D_TL( z , asi%wlp, asi%xlp, & ! FWD Input
z_TL, wlp_TL , xlp_TL , & ! TL Input
w_TL ) ! TL Output
! Phase matrix coefficients
IF (AerosolScatter_TL%n_Phase_Elements > 0 .and. AerosolScatter_TL%Include_Scattering ) THEN
pcoeff_TL(0,1) = ZERO
DO m = 1, AerosolScatter_TL%n_Phase_Elements
DO l = 1, AerosolScatter_TL%n_Legendre_Terms
z => AeroC%pcoeff(asi%i1:asi%i2,asi%j1:asi%j2,k,l+AerosolScatter_TL%lOffset,m)
CALL interp_2D_TL( z , asi%wlp, asi%xlp, & ! FWD Input
z_TL, wlp_TL , xlp_TL , & ! TL Input
pcoeff_TL(l,m) ) ! TL Output
END DO
END DO
ELSE
! Absorption coefficient
IF( w < ONE ) THEN
ke_TL = ke_TL * (ONE - w) - ke/(ONE -w) * w_TL
ELSE
ke_TL = ZERO
END IF
END IF
NULLIFY(z)
END SUBROUTINE Get_Aerosol_Opt_TL
! ---------------------------------------------
! Subroutine to obtain the adjoint
! bulk optical properties of a aerosol:
! extinction coefficient (ke_AD),
! scattereing coefficient (w_AD)
! asymmetry factor (g_AD), and
! spherical Legendre coefficients (pcoeff_AD)
! ---------------------------------------------
SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD structure
Aerosol_Type , & ! Input see CRTM_Aerosol_Define.f90
ke , & ! Input
w , & ! Input
ke_AD , & ! AD Input extinction cross section
w_AD , & ! AD Input single scatter albedo
pcoeff_AD , & ! AD Input spherical Legendre coefficients
Reff_AD , & ! AD Output effective radius
asi ) ! Input interpolation data
! Arguments
TYPE(CRTM_AtmOptics_type), INTENT(IN) :: AerosolScatter_AD
INTEGER , INTENT(IN) :: Aerosol_Type
REAL(fp), INTENT(IN) :: ke, w
REAL(fp), INTENT(IN OUT) :: ke_AD ! AD Input
REAL(fp), INTENT(IN OUT) :: w_AD ! AD Input
REAL(fp), INTENT(IN OUT) :: pcoeff_AD(0:,:) ! AD Input
REAL(fp), INTENT(IN OUT) :: Reff_AD ! AD Output
TYPE(ASinterp_type), INTENT(IN) :: asi
! Local variables
INTEGER :: k, l, m
REAL(fp) :: f_int_AD, r_int_AD
REAL(fp) :: f_AD(NPTS), r_AD(NPTS)
REAL(fp) :: z_AD(NPTS,NPTS)
TYPE(LPoly_type) :: wlp_AD, xlp_AD
REAL(fp), POINTER :: z(:,:) => NULL()
! Setup
! -----
! No AD output all dimensions
! are outside LUT bounds
IF ( asi%f_outbound .AND. asi%r_outbound ) THEN
Reff_AD = ZERO
ke_AD = ZERO
w_AD = ZERO
pcoeff_AD = ZERO
RETURN
END IF
! Initialise local adjoint variables
f_int_AD = ZERO
r_int_AD = ZERO
f_AD = ZERO
r_AD = ZERO
z_AD = ZERO
CALL Clear_LPoly(wlp_AD)
CALL Clear_LPoly(xlp_AD)
! Get the aerosol type LUT index
! ------------------------------
k = Aerosol_Type_Index( Aerosol_Type )
! Perform interpolation
! ---------------------
! Phase matrix coefficients
IF (AerosolScatter_AD%n_Phase_Elements > 0 .and. AerosolScatter_AD%Include_Scattering ) THEN
DO m = 1, AerosolScatter_AD%n_Phase_Elements
DO l = 1, AerosolScatter_AD%n_Legendre_Terms
z => AeroC%pcoeff(asi%i1:asi%i2,asi%j1:asi%j2,k,l+AerosolScatter_AD%lOffset,m)
CALL interp_2D_AD( z , asi%wlp, asi%xlp, & ! FWD Input
pcoeff_AD(l,m) , & ! AD Input
z_AD, wlp_AD, xlp_AD ) ! AD Output
END DO
END DO
pcoeff_AD(0,1) = ZERO
ELSE
! Absorption coefficient
IF( w < ONE ) THEN
w_AD = w_AD - ke/(ONE -w) * ke_AD
ke_AD = ke_AD * (ONE - w)
ELSE
ke_AD = ZERO
END IF
END IF
! Single scatter albedo
z => AeroC%w(asi%i1:asi%i2,asi%j1:asi%j2,k)
CALL interp_2D_AD( z , asi%wlp, asi%xlp, & ! FWD Input
w_AD , & ! AD Input
z_AD, wlp_AD , xlp_AD ) ! AD Output
! Extinction coefficient
z => AeroC%ke(asi%i1:asi%i2,asi%j1:asi%j2,k)
CALL interp_2D_AD( z , asi%wlp, asi%xlp, & ! FWD Input
ke_AD , & ! AD Input
z_AD, wlp_AD , xlp_AD ) ! AD Output
NULLIFY(z)
! Compute the AD of the interpolating polynomials
! -----------------------------------------------
! Efective radius term
CALL LPoly_AD( asi%r, asi%r_int, & ! FWD Input
asi%xlp, & ! FWD Input
xlp_AD, & ! AD Input
r_AD, r_int_AD ) ! AD Output
! Frequency term (always zero. This is a placeholder for testing)
CALL LPoly_AD( asi%f, asi%f_int, & ! FWD Input
asi%wlp, & ! FWD Input
wlp_AD, & ! AD Input
f_AD, f_int_AD ) ! AD Output
! The AD outputs
! --------------
Reff_AD = Reff_AD + r_int_AD
END SUBROUTINE Get_Aerosol_Opt_AD
END MODULE CRTM_AerosolScatter