!
! SOI_Module
!
! Module containing the Successive Order of Interation (SOI) radiative
! transfer solution procedures used in the CRTM.
!
!
! CREATION HISTORY:
! Written by: Tom Greenwald, CIMSS/SSEC; tom.greenwald@ssec.wisc.edu
! Paul van Delst; paul.vandelst@noaa.gov
! 20-Jan-2010
MODULE SOI_Module
! ------------------
! Environment set up
! ------------------
! Module use statements
USE Type_Kinds, ONLY: fp
USE Message_Handler, ONLY: SUCCESS, FAILURE, Display_Message
USE CRTM_Parameters, ONLY: SET, ZERO, ONE, TWO, PI, &
MAX_N_LAYERS, MAX_N_ANGLES, MAX_N_LEGENDRE_TERMS, &
DEGREES_TO_RADIANS, &
SECANT_DIFFUSIVITY, &
SCATTERING_ALBEDO_THRESHOLD, &
OPTICAL_DEPTH_THRESHOLD
USE CRTM_Utility
USE RTV_Define , ONLY: RTV_type, &
DELTA_OPTICAL_DEPTH, &
MAX_N_DOUBLING, &
MAX_N_SOI_ITERATIONS
! Disable all implicit typing
IMPLICIT NONE
! --------------------
! Default visibilities
! --------------------
! Everything private by default
PRIVATE
! Procedures
PUBLIC :: CRTM_SOI
PUBLIC :: CRTM_SOI_TL
PUBLIC :: CRTM_SOI_AD
PUBLIC :: CRTM_SOI_Version
! -----------------
! Module parameters
! -----------------
! Version Id for the module
CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
'$Id: SOI_Module.f90 99117 2017-11-27 18:37:14Z tong.zhu@noaa.gov $'
CONTAINS
!################################################################################
!################################################################################
!## ##
!## ## PUBLIC MODULE ROUTINES ## ##
!## ##
!################################################################################
!################################################################################
SUBROUTINE CRTM_SOI(n_Layers, & ! Input number of atmospheric layers
w, & ! Input layer scattering albedo
T_OD, & ! Input layer optical depth
cosmic_background, & ! Input cosmic background radiance
emissivity, & ! Input surface emissivity
reflectivity, & ! Input surface reflectivity matrix
Index_Sat_Angle, & ! Input satellite angle index
RTV) ! IN/Output upward radiance and others
! ------------------------------------------------------------------------- !
! !
! FUNCTION: !
! !
! This subroutine calculates IR/MW radiance at the top of the atmosphere !
! including multiple scattering using the SOI (Successive Order of !
! Interaction) algorithm, which combines the Successive Order of !
! Scattering and doubling methods. !
! !
! Written by Tom Greenwald tomg@ssec.wisc.edu !
! ------------------------------------------------------------------------- !
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_Layers
REAL (fp), INTENT(IN), DIMENSION( : ) :: w, T_OD
REAL (fp), INTENT(IN) :: cosmic_background
REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity
REAL (fp), INTENT(IN), DIMENSION( :,: ) :: reflectivity
INTEGER, INTENT( IN ) :: Index_Sat_Angle
TYPE(RTV_type), INTENT( INOUT ) :: RTV
! -------------- internal variables --------------------------------- !
INTEGER :: i, k, niter
REAL(fp) :: rad, rad_change
REAL(fp), PARAMETER :: initial_error = 1.E10
REAL(fp), PARAMETER :: SMALL = 1.E-15
REAL(fp), PARAMETER :: SNGL_SCAT_ALB_THRESH = 0.8
REAL(fp), PARAMETER :: OPT_DEPTH_THRESH = 4.0
REAL(fp) :: radiance_thresh
REAL(fp), DIMENSION( MAX_N_ANGLES ) :: source
! Precompute layer R/T matrices and thermal sources
DO k = 1, n_Layers
! Precompute simple layer properties
RTV%e_Layer_Trans( 1 : RTV%n_Angles, k ) = EXP( -T_OD( k ) / RTV%COS_Angle( 1 : RTV%n_Angles ) )
IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD ) THEN
IF ( w( k ) < SNGL_SCAT_ALB_THRESH .AND. T_OD( k ) < OPT_DEPTH_THRESH ) THEN
CALL CRTM_Truncated_Doubling( RTV%n_Streams, RTV%n_Angles, k, w( k ), T_OD( k ), & ! Input
RTV%COS_Angle, RTV%COS_Weight, RTV%Pff( :, :, k ), & ! Input
RTV%Pbb( :, :, k ), RTV%Planck_Atmosphere( k ), & ! Input
RTV ) ! Output
ELSE
CALL CRTM_Doubling_Layer( RTV%n_Streams, RTV%n_Angles, k, w( k ), T_OD( k ), & ! Input
RTV%COS_Angle, RTV%COS_Weight, RTV%Pff( :, :, k ), & ! Input
RTV%Pbb( :, :, k ), RTV%Planck_Atmosphere( k ), & ! Input
RTV ) ! Output
END IF
! Subtract out direct transmission
DO i = 1, RTV%n_angles
RTV%s_Layer_Trans( i, i, k ) = RTV%s_Layer_Trans( i, i, k ) - RTV%e_layer_Trans( i, k )
END DO
ELSE ! Thermal sources
DO i = 1, RTV%n_Angles
RTV%s_Layer_Source_UP( i, k ) = RTV%Planck_Atmosphere( k ) * ( ONE - RTV%e_Layer_Trans( i, k ) )
RTV%s_Layer_Source_DOWN( i, k ) = RTV%s_Layer_Source_UP( i, k )
END DO
END IF
END DO
!------------------------------------
! Arrays that *must* be zeroed out
!------------------------------------
RTV%s_level_Rad_UP( Index_Sat_Angle, 0 ) = ZERO
IF ( RTV%Number_SOI_Iter > 0 ) THEN
RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, 0 : n_Layers, 1 : RTV%Number_SOI_Iter ) = ZERO
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, 0 : n_Layers, 1 : RTV%Number_SOI_Iter ) = ZERO
ELSE
RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, 0 : n_Layers, 1 ) = ZERO
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, 0 : n_Layers, 1 ) = ZERO
END IF
!---------------------------------
! Set initial/boundary conditions
!---------------------------------
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, 0, 1 ) = cosmic_background
niter = 0 ! Set that we're on the zeroth order of scatter
rad_change = initial_error ! SOI loop uses this quantity to know when to stop iteration
rad = SMALL
! ACCEL = .FALSE. ! Turn on convergence acceleration
! IF ( ACCEL ) THEN
! radiance_thresh = 5.E-7
! radsum = 0.
! q_old = 0.
! ELSE
radiance_thresh = 1.E-4
! END IF
!-----------------------------------------
! This is the Order of Interaction loop
!-----------------------------------------
soi_loop: DO WHILE ( ( ( rad_change > radiance_thresh ) .AND. ( ( niter + 1 ) <= MAX_N_SOI_ITERATIONS ) ) .OR. ( niter < 2 ) )
niter = niter + 1 ! Note: niter = 1 is the zeroth order of interaction
IF ( niter > 1 ) RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, 0, niter ) = ZERO ! No cosmic background contribution for
! higher orders of interaction
!---------------------------------------------
! Integrate down through atmospheric layers
!---------------------------------------------
layersdn_loop: DO k = 1, n_Layers
IF ( niter == 1 ) THEN
source( 1 : RTV%n_Angles ) = RTV%s_Layer_Source_DOWN( 1 : RTV%n_Angles, k )
ELSE
IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD .AND. T_OD( k ) >= OPTICAL_DEPTH_THRESHOLD ) THEN
! DO j = 1, RTV%n_Angles
! source( j ) = ZERO
! DO i = 1, RTV%n_Angles ! Integrate over angle
! source( j ) = source( j ) + RTV%s_Layer_Trans( j, i, k ) * RTV%s_Level_IterRad_DOWN( i, k - 1, niter - 1 ) &
! + RTV%s_Layer_Refl( j, i, k ) * RTV%s_Level_IterRad_UP( i, k, niter - 1 )
! END DO
! END DO
source( 1 : RTV%n_Angles ) = matmul( RTV%s_Layer_Trans( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, k - 1, niter - 1 ) ) + &
matmul( RTV%s_Layer_Refl( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, k, niter - 1 ) )
ELSE
source( 1 : RTV%n_Angles ) = ZERO
END IF
END IF
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, k, niter ) = RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, k - 1, niter ) &
* RTV%e_Layer_Trans( 1 : RTV%n_Angles, k ) &
+ source( 1 : RTV%n_Angles )
END DO layersdn_loop
!-----------------------------------------------
! Account for surface reflection and emission
!-----------------------------------------------
DO i = 1, RTV%n_Angles
RTV%s_Level_IterRad_UP( i, n_Layers, niter ) = reflectivity( i, i ) * RTV%s_Level_IterRad_DOWN( i, n_Layers, niter )
END DO
IF ( niter == 1 ) THEN
! Add surface emission only for zeroth order of interaction
RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, n_Layers, niter ) = &
RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, n_Layers, niter ) + &
emissivity( 1 : RTV%n_Angles ) * RTV%Planck_Surface
END IF
!-------------------------------------------------
! Integrate back up through atmospheric layers
!-------------------------------------------------
layersup_loop: DO k = n_Layers, 1, -1
IF ( niter == 1 ) THEN
source( 1 :RTV%n_Angles ) = RTV%s_Layer_Source_UP( 1 : RTV%n_Angles, k )
ELSE
IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD .AND. T_OD( k ) >= OPTICAL_DEPTH_THRESHOLD ) THEN
! DO j = 1, RTV%n_Angles
! source( j ) = ZERO
! DO i = 1, RTV%n_Angles ! Integrate over angle
! source( j ) = source( j ) + RTV%s_Layer_Refl( j, i, k ) * RTV%s_Level_IterRad_DOWN( i, k - 1, niter - 1 ) &
! + RTV%s_Layer_Trans( j, i, k ) * RTV%s_Level_IterRad_UP( i, k, niter - 1 )
! END DO
! END DO
source( 1 : RTV%n_Angles ) = matmul( RTV%s_Layer_Refl( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, k - 1, niter - 1 ) ) + &
matmul( RTV%s_Layer_Trans( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, k, niter - 1 ) )
ELSE
source( 1 : RTV%n_Angles ) = ZERO
END IF
END IF
RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, k - 1, niter ) = RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, k, niter ) &
* RTV%e_Layer_Trans( 1 : RTV%n_Angles, k ) &
+ source( 1 : RTV%n_Angles )
END DO layersup_loop
!-------------------------------------------------------------------
! Save solution at observation angle for this order of interaction
!-------------------------------------------------------------------
rad = RTV%s_Level_IterRad_UP( Index_Sat_Angle, 0, niter )
!---------------------------------
! Add it to cumulative solution
!---------------------------------
RTV%s_level_Rad_UP( Index_Sat_Angle, 0 ) = RTV%s_level_Rad_UP( Index_Sat_Angle, 0 ) + rad
!---------------------------------------
! Accelerate covergence of iterations?
!---------------------------------------
! IF ( ACCEL ) THEN
! radsum = radsum + rad
! IF ( niter > 1 ) THEN
! diff = rad_old - rad
! IF ( diff == 0. ) RETURN
! q = radsum + rad ** 2 / diff
! rad_change = ABS( q - q_old )
! q_old = q
! END IF
! rad_old = rad
! ELSE
! Check how much cumulative solution has changed
rad_change = rad / RTV%s_level_Rad_UP( Index_Sat_Angle, 0 )
! END IF
END DO soi_loop
RTV%Number_SOI_Iter = niter
RETURN
END SUBROUTINE CRTM_SOI
SUBROUTINE CRTM_SOI_TL(n_Layers, & ! Input number of atmospheric layers
w, & ! Input layer scattering albedo
T_OD, & ! Input layer optical depth
emissivity, & ! Input surface emissivity
reflectivity, & ! Input surface reflectivity matrix
Index_Sat_Angle, & ! Input satellite angle index
RTV, & ! Input structure containing forward part results
Planck_Atmosphere_TL, & ! Input tangent-linear atmospheric layer Planck radiance
Planck_Surface_TL, & ! Input TL surface Planck radiance
w_TL, & ! Input TL layer scattering albedo
T_OD_TL, & ! Input TL layer optical depth
emissivity_TL, & ! Input TL surface emissivity
reflectivity_TL, & ! Input TL reflectivity
Pff_TL, & ! Input TL forward phase matrix
Pbb_TL, & ! Input TL backward phase matrix
s_rad_up_TL) ! Output TL upward radiance
! ------------------------------------------------------------------------- !
! !
! FUNCTION: !
! !
! This subroutine calculates the tangent linear of the IR/MW radiance at !
! the top of the atmosphere including multiple scattering using the SOI !
! (Successive Order of Interaction) algorithm, which combines the !
! Successive Order of Scattering and doubling methods. !
! !
! Written by Tom Greenwald tomg@ssec.wisc.edu !
! ------------------------------------------------------------------------- !
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_Layers
REAL (fp), INTENT(IN), DIMENSION( : ) :: w, T_OD
REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity
REAL (fp), INTENT(IN), DIMENSION( :, : ) :: reflectivity
INTEGER, INTENT( IN ) :: Index_Sat_Angle
TYPE(RTV_type), INTENT( IN ) :: RTV
REAL (fp), INTENT(IN), DIMENSION( 0: ) :: Planck_Atmosphere_TL
REAL (fp), INTENT(IN) :: Planck_Surface_TL
REAL (fp), INTENT(IN), DIMENSION( : ) :: w_TL, T_OD_TL
REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity_TL
REAL (fp), INTENT(IN), DIMENSION( :, : ) :: reflectivity_TL
REAL (fp), INTENT(IN), DIMENSION( :, :, : ) :: Pff_TL, Pbb_TL
REAL (fp), INTENT(INOUT), DIMENSION( : ) :: s_rad_up_TL
! -------------- internal variables --------------------------------- !
INTEGER :: i, k, iter
REAL(fp), PARAMETER :: SNGL_SCAT_ALB_THRESH = 0.8
REAL(fp), PARAMETER :: OPT_DEPTH_THRESH = 4.0
REAL(fp), DIMENSION( RTV%n_Angles ) :: source_TL
REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: e_Trans_TL
REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: s_source_up_TL, s_source_down_TL
REAL(fp), DIMENSION( RTV%n_Angles, RTV%n_Angles, n_Layers ) :: s_trans_TL, s_refl_TL
REAL(fp), DIMENSION( RTV%n_Angles, 0:n_Layers, RTV%Number_SOI_Iter ) :: s_IterRad_UP_TL, s_IterRad_DOWN_TL
DO k = 1, n_Layers
e_Trans_TL( 1 : RTV%n_Angles, k ) = -T_OD_TL(k) * (EXP( -T_OD( k ) / RTV%COS_Angle( 1 : RTV%n_Angles ) ) ) / &
RTV%COS_Angle( 1 : RTV%n_Angles )
IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD ) THEN
IF ( w( k ) < SNGL_SCAT_ALB_THRESH .AND. T_OD( k ) < OPT_DEPTH_THRESH ) THEN
CALL CRTM_Truncated_Doubling_TL(RTV%n_Streams, RTV%n_Angles, k, w( k ), T_OD( k ), & !Input
RTV%COS_Angle( 1 : RTV%n_Angles ), RTV%COS_Weight( 1 : RTV%n_Angles ), & !Input
RTV%Pff( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), & !Input
RTV%Pbb( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), RTV%Planck_Atmosphere( k ), & !Input
w_TL( k ), T_OD_TL( k ), Pff_TL( :, :, k ), & !Input
Pbb_TL( :, :, k ), Planck_Atmosphere_TL( k ), RTV, & !Input
s_trans_TL( :, :, k ), s_refl_TL( :, :, k ), s_source_up_TL( :, k ), & !Output
s_source_down_TL( :, k ) ) !Output
ELSE
CALL CRTM_Doubling_layer_TL(RTV%n_Streams, RTV%n_Angles, k, w( k ), T_OD( k ), & !Input
RTV%COS_Angle( 1 : RTV%n_Angles ), RTV%COS_Weight( 1 : RTV%n_Angles ), & !Input
RTV%Pff( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), & !Input
RTV%Pbb( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), RTV%Planck_Atmosphere( k ), & !Input
w_TL( k ), T_OD_TL( k ), Pff_TL( :, :, k ), & !Input
Pbb_TL( :, :, k ), Planck_Atmosphere_TL( k ), RTV, & !Input
s_trans_TL( :, :, k), s_refl_TL( :, :, k ), s_source_up_TL( :, k ), & !Output
s_source_down_TL( :, k ) ) !Output
END IF
! Subtract out direct transmission
DO i = 1, RTV%n_angles
s_Trans_TL( i, i, k ) = s_Trans_TL( i, i, k ) - e_Trans_TL( i, k )
END DO
ELSE ! Thermal sources
DO i = 1, RTV%n_Angles
s_Source_UP_TL( i, k ) = Planck_Atmosphere_TL( k ) * ( ONE - RTV%e_Layer_Trans( i, k ) ) - &
RTV%PLanck_Atmosphere( k ) * e_Trans_TL( i, k )
s_Source_DOWN_TL( i, k ) = s_Source_UP_TL( i, k )
END DO
END IF
END DO
s_Rad_UP_TL( Index_Sat_Angle ) = ZERO
s_IterRad_UP_TL( 1 : RTV%n_Angles, 0 : n_Layers, 1 : RTV%Number_SOI_Iter ) = ZERO
s_IterRad_DOWN_TL( 1 : RTV%n_Angles, 0 : n_Layers, 1 : RTV%Number_SOI_Iter ) = ZERO
!-----------------------------------------
! This is the Order of Interaction loop
!-----------------------------------------
DO iter = 1, RTV%Number_SOI_Iter
IF ( iter > 1 ) s_IterRad_DOWN_TL( 1 : RTV%n_Angles, 0, iter ) = ZERO
!---------------------------------------------
! Integrate down through atmospheric layers
!---------------------------------------------
layersdn_loop: DO k = 1, n_Layers
IF ( iter == 1 ) THEN
source_TL( 1 : RTV%n_Angles ) = s_Source_DOWN_TL( 1 : RTV%n_Angles, k )
ELSE
IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD .AND. T_OD( k ) >= OPTICAL_DEPTH_THRESHOLD ) THEN
source_TL( 1 : RTV%n_Angles ) = matmul( RTV%s_Layer_Trans( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
s_IterRad_DOWN_TL( 1 : RTV%n_Angles, k - 1, iter - 1 ) ) + &
matmul( s_Trans_TL( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, k - 1, iter - 1 ) ) + &
matmul( RTV%s_Layer_Refl( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
s_IterRad_UP_TL( 1 : RTV%n_Angles, k, iter - 1 ) ) + &
matmul( s_Refl_TL( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
RTV%s_level_IterRad_UP( 1 : RTV%n_Angles, k, iter - 1 ) )
ELSE
source_TL( 1 : RTV%n_Angles ) = ZERO
END IF
END IF
s_IterRad_DOWN_TL( 1 : RTV%n_Angles, k, iter ) = s_IterRad_DOWN_TL( 1 : RTV%n_Angles, k - 1, iter ) &
* RTV%e_Layer_Trans( 1 : RTV%n_Angles, k ) + &
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, k - 1, iter ) &
* e_Trans_TL( 1 : RTV%n_Angles, k ) &
+ source_TL( 1 : RTV%n_Angles )
END DO layersdn_loop
!-----------------------------------------------
! Account for surface reflection and emission
!-----------------------------------------------
DO i = 1, RTV%n_Angles
s_IterRad_UP_TL( i, n_Layers, iter ) = reflectivity_TL( i, i ) * RTV%s_Level_IterRad_DOWN( i, n_Layers, iter ) + &
reflectivity( i, i ) * s_IterRad_DOWN_TL( i, n_Layers, iter )
END DO
IF ( iter == 1 ) THEN
! Add surface emission only for zeroth order of interaction
s_IterRad_UP_TL( 1 : RTV%n_Angles, n_Layers, iter ) = s_IterRad_UP_TL( 1 : RTV%n_Angles, n_Layers, iter ) + &
emissivity( 1 : RTV%n_Angles ) * Planck_Surface_TL + &
emissivity_TL( 1 : RTV%n_Angles ) * RTV%Planck_Surface
END IF
!-------------------------------------------------
! Integrate back up through atmospheric layers
!-------------------------------------------------
layersup_loop: DO k = n_Layers, 1, -1
IF ( iter == 1 ) THEN
source_TL( 1 :RTV%n_Angles ) = s_Source_UP_TL( 1 : RTV%n_Angles, k )
ELSE
IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD .AND. T_OD( k ) >= OPTICAL_DEPTH_THRESHOLD ) THEN
source_TL( 1 : RTV%n_Angles ) = matmul( RTV%s_Layer_Refl( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
s_IterRad_DOWN_TL( 1 : RTV%n_Angles, k - 1, iter - 1 ) ) + &
matmul( s_Refl_TL( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, k - 1, iter - 1 ) ) + &
matmul( RTV%s_Layer_Trans( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
s_IterRad_UP_TL( 1 : RTV%n_Angles, k, iter - 1 ) ) + &
matmul( s_Trans_TL( 1 : RTV%n_Angles, 1 : RTV%n_Angles, k ), &
RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, k, iter - 1 ) )
ELSE
source_TL( 1 : RTV%n_Angles ) = ZERO
END IF
END IF
s_IterRad_UP_TL( 1 : RTV%n_Angles, k - 1, iter ) = RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, k, iter ) &
* e_Trans_TL( 1 : RTV%n_Angles, k ) + &
s_IterRad_UP_TL( 1 : RTV%n_Angles, k, iter ) &
* RTV%e_Layer_Trans( 1 : RTV%n_Angles, k ) &
+ source_TL( 1 : RTV%n_Angles )
END DO layersup_loop
!----------------
! Add up terms
!----------------
s_Rad_UP_TL( Index_Sat_Angle ) = s_Rad_UP_TL( Index_Sat_Angle ) + s_IterRad_UP_TL( Index_Sat_Angle, 0, iter )
END DO
RETURN
END SUBROUTINE CRTM_SOI_TL
SUBROUTINE CRTM_SOI_AD(n_Layers, & ! Input number of atmospheric layers
w, & ! Input layer scattering albedo
T_OD, & ! Input layer optical depth
emissivity, & ! Input surface emissivity
reflectivity, & ! Input surface reflectivity matrix
Index_Sat_Angle, & ! Input satellite angle index
RTV, & ! Input structure containing forward results
s_rad_up_AD, & ! Input adjoint upward radiance
Planck_Atmosphere_AD, & ! Output AD atmospheric layer Planck radiance
Planck_Surface_AD, & ! Output AD surface Planck radiance
w_AD, & ! Output AD layer scattering albedo
T_OD_AD, & ! Output AD layer optical depth
emissivity_AD, & ! Output AD surface emissivity
reflectivity_AD, & ! Output AD surface reflectivity
Pff_AD, & ! Output AD forward phase matrix
Pbb_AD) ! Output AD backward phase matrix
! ------------------------------------------------------------------------- !
! FUNCTION: !
! This subroutine calculates IR/MW adjoint radiance at the top of !
! the atmosphere including atmospheric scattering. !
! !
! Tom Greenwald tomg@ssec.wisc.edu !
! ------------------------------------------------------------------------- !
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_Layers
REAL (fp), INTENT(IN), DIMENSION( : ) :: w, T_OD
REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity
REAL (fp), INTENT(IN), DIMENSION( :, : ) :: reflectivity
INTEGER, INTENT(IN) :: Index_Sat_Angle
TYPE(RTV_type), INTENT(IN) :: RTV
REAL (fp),INTENT(INOUT),DIMENSION( :, :, : ) :: Pff_AD, Pbb_AD
REAL (fp),INTENT(INOUT),DIMENSION( : ) :: w_AD,T_OD_AD
REAL (fp),INTENT(INOUT),DIMENSION( 0: ) :: Planck_Atmosphere_AD
REAL (fp),INTENT(INOUT) :: Planck_Surface_AD
REAL (fp),INTENT(INOUT),DIMENSION( : ) :: emissivity_AD
REAL (fp),INTENT(INOUT),DIMENSION( :, : ) :: reflectivity_AD
REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_AD
! Local variables
REAL(fp), PARAMETER :: SNGL_SCAT_ALB_THRESH = 0.8
REAL(fp), PARAMETER :: OPT_DEPTH_THRESH = 4.0
INTEGER :: iter, k, i, j
REAL(fp), DIMENSION( RTV%n_Angles, 0:n_Layers, RTV%Number_SOI_Iter ) :: s_IterRad_UP_AD
REAL(fp), DIMENSION( RTV%n_Angles, 0:n_Layers, RTV%Number_SOI_Iter ) :: s_IterRad_DOWN_AD
REAL(fp), DIMENSION( RTV%n_Angles ) :: source_AD
REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: s_source_UP_AD
REAL(fp), DIMENSION( RTV%n_Angles, n_layers ) :: s_source_DOWN_AD
REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: e_Trans_AD
REAL(fp), DIMENSION( RTV%n_Angles, RTV%n_Angles, n_Layers ) :: s_Refl_AD
REAL(fp), DIMENSION( RTV%n_Angles, RTV%n_Angles, n_Layers ) :: s_Trans_AD
! Zero out all local variables
s_IterRad_UP_AD = ZERO
s_IterRad_DOWN_AD = ZERO
source_AD = ZERO
s_source_UP_AD = ZERO
s_source_DOWN_AD = ZERO
e_Trans_AD = ZERO
s_Refl_AD = ZERO
s_Trans_AD = ZERO
Pff_AD = ZERO
Pbb_AD = ZERO
!--------------------------------------------------------------------
! Loop through each successive order of interaction in reverse order
!--------------------------------------------------------------------
DO iter = RTV%Number_SOI_Iter, 1, -1
s_IterRad_UP_AD( Index_Sat_Angle, 0, iter ) = s_Rad_UP_AD( Index_Sat_Angle )
!---------------------------------------
! Step down through upward integration
!---------------------------------------
DO k = 1, n_Layers
source_AD( 1 : RTV%n_Angles ) = source_AD( 1 : RTV%n_Angles ) + s_IterRad_UP_AD( 1 : RTV%n_Angles, k - 1, iter )
e_Trans_AD( 1 : RTV%n_Angles, k ) = e_Trans_AD( 1 : RTV%n_Angles, k ) + &
s_IterRad_UP_AD( 1 : RTV%n_Angles, k - 1, iter ) * &
RTV%s_Level_IterRad_UP( 1 : RTV%n_Angles, k, iter )
s_IterRad_UP_AD( 1 : RTV%n_Angles, k, iter ) = s_IterRad_UP_AD( 1 : RTV%n_Angles, k, iter ) + &
s_IterRad_UP_AD( 1 : RTV%n_Angles, k - 1, iter ) * &
RTV%e_Layer_Trans( 1 : RTV%n_Angles, k )
s_IterRad_UP_AD( 1 : RTV%n_Angles, k - 1, iter ) = ZERO
IF ( iter == 1 ) THEN
s_source_UP_AD( 1 : RTV%n_Angles, k ) = s_source_UP_AD( 1 : RTV%n_Angles, k ) + source_AD( 1 : RTV%n_Angles )
source_AD( 1 : RTV%n_Angles ) = ZERO
ELSE
IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD .AND. T_OD( k ) >= OPTICAL_DEPTH_THRESHOLD ) THEN
DO j = 1, RTV%n_Angles
DO i = 1, RTV%n_Angles ! Integrate over angle
s_Refl_AD( j, i, k ) = s_Refl_AD( j, i, k ) + source_AD( j ) * RTV%s_Level_IterRad_DOWN( i, k - 1, iter - 1 )
s_IterRad_DOWN_AD( i, k - 1, iter - 1 ) = s_IterRad_DOWN_AD( i, k - 1, iter - 1 ) + source_AD( j ) * &
RTV%s_Layer_Refl( j, i, k )
s_Trans_AD( j, i, k ) = s_Trans_AD( j, i, k ) + source_AD( j ) * RTV%s_Level_IterRad_UP( i, k, iter - 1 )
s_IterRad_UP_AD( i, k, iter - 1 ) = s_IterRad_UP_AD( i, k, iter - 1 ) + source_AD( j ) * &
RTV%s_Layer_Trans( j, i, k )
END DO
source_AD( j ) = ZERO
END DO
ELSE
source_AD( 1 : RTV%n_Angles ) = ZERO
END IF
END IF
END DO
IF ( iter == 1 ) THEN
! Add surface emission only for zeroth order of interaction
emissivity_AD( 1 : RTV%n_Angles ) = emissivity_AD( 1 : RTV%n_Angles ) + &
s_IterRad_UP_AD( 1 : RTV%n_Angles, n_Layers, iter ) * RTV%Planck_Surface
Planck_Surface_AD = Planck_Surface_AD + SUM( s_IterRad_UP_AD( :, n_Layers, iter ) * emissivity )
END IF
!-----------------------------------------------
! Account for surface reflection and emission
!-----------------------------------------------
DO i = 1, RTV%n_Angles
reflectivity_AD( i, i ) = reflectivity_AD( i, i ) + s_IterRad_UP_AD( i, n_Layers, iter ) * &
RTV%s_Level_IterRad_DOWN( i, n_Layers, iter )
s_IterRad_DOWN_AD( i, n_Layers, iter ) = s_IterRad_DOWN_AD( i, n_Layers, iter ) + &
s_IterRad_UP_AD( i, n_Layers, iter ) * reflectivity( i, i )
s_IterRad_UP_AD( i, n_Layers, iter ) = ZERO
END DO
!---------------------------------------
! Step up through downward integration
!---------------------------------------
DO k = n_Layers, 1, -1
source_AD( 1 : RTV%n_Angles ) = source_AD( 1 : RTV%n_Angles ) + s_IterRad_DOWN_AD( 1 : RTV%n_Angles, k, iter )
e_Trans_AD( 1 : RTV%n_Angles, k ) = e_Trans_AD( 1 : RTV%n_Angles, k ) + &
s_IterRad_DOWN_AD( 1 : RTV%n_Angles, k, iter ) * &
RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, k - 1, iter )
s_IterRad_DOWN_AD( 1 : RTV%n_Angles, k - 1, iter ) = s_IterRad_DOWN_AD( 1 : RTV%n_Angles, k - 1, iter ) + &
s_IterRad_DOWN_AD( 1 : RTV%n_Angles, k, iter ) * &
RTV%e_Layer_Trans( 1 : RTV%n_Angles, k )
s_IterRad_DOWN_AD( 1 : RTV%n_Angles, k, iter ) = ZERO
IF ( iter == 1 ) THEN
s_source_DOWN_AD( 1 : RTV%n_Angles, k ) = s_source_DOWN_AD( 1 : RTV%n_Angles, k ) + source_AD( 1 : RTV%n_Angles )
source_AD( 1 : RTV%n_Angles ) = ZERO
ELSE
IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD .AND. T_OD( k ) >= OPTICAL_DEPTH_THRESHOLD ) THEN
DO j = 1, RTV%n_Angles
DO i = 1, RTV%n_Angles ! Integrate over angle
s_Trans_AD( j, i, k ) = s_Trans_AD( j, i, k ) + source_AD( j ) * RTV%s_Level_IterRad_DOWN( i, k - 1, iter - 1 )
s_IterRad_DOWN_AD( i, k - 1, iter - 1 ) = s_IterRad_DOWN_AD( i, k - 1, iter - 1 ) + source_AD( j ) * &
RTV%s_Layer_Trans( j, i, k )
s_Refl_AD( j, i, k ) = s_Refl_AD( j, i, k ) + source_AD( j ) * RTV%s_Level_IterRad_UP( i, k, iter - 1 )
s_IterRad_UP_AD( i, k, iter - 1 ) = s_IterRad_UP_AD( i, k, iter - 1 ) + source_AD( j ) * &
RTV%s_Layer_Refl( j, i, k )
END DO
source_AD( j ) = ZERO
END DO
ELSE
source_AD( 1 : RTV%n_Angles ) = ZERO
END IF
END IF
END DO
IF ( iter > 1 ) s_IterRad_DOWN_AD( 1 : RTV%n_Angles, 0, iter ) = ZERO ! No cosmic background contribution for
! higher orders of interaction
END DO ! SOI iteration
s_IterRad_DOWN_AD( 1 : RTV%n_Angles, 0 : n_Layers, 1 : RTV%Number_SOI_Iter ) = ZERO
s_IterRad_UP_AD( 1 : RTV%n_Angles, 0 : n_Layers, 1 : RTV%Number_SOI_Iter ) = ZERO
s_Rad_UP_AD( Index_Sat_Angle ) = ZERO
DO k = 1, n_Layers
IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD ) THEN
DO i = 1, RTV%n_angles
e_Trans_AD( i, k ) = e_Trans_AD( i, k ) - s_Trans_AD( i, i, k )
END DO
IF ( w( k ) < SNGL_SCAT_ALB_THRESH .AND. T_OD( k ) < OPT_DEPTH_THRESH ) THEN
CALL CRTM_Truncated_Doubling_AD(RTV%n_Streams, RTV%n_Angles, k, w( k ), T_OD( k ), & !Input
RTV%COS_Angle, RTV%COS_Weight, RTV%Pff( :, :, k ), RTV%Pbb( :, :, k ), & ! Input
RTV%Planck_Atmosphere( k ), & !Input
s_trans_AD( :, :, k ), s_refl_AD( :, :, k ), s_source_up_AD( :, k ), &
s_source_down_AD( :, k ), RTV, w_AD( k ), T_OD_AD( k ), Pff_AD( :, :, k ), &
Pbb_AD( :, :, k ), Planck_Atmosphere_AD( k ) ) !Output
ELSE
CALL CRTM_Doubling_layer_AD(RTV%n_Streams, RTV%n_Angles, k, w( k ), T_OD( k ), & !Input
RTV%COS_Angle, RTV%COS_Weight, RTV%Pff( :, :, k ), RTV%Pbb( :, :, k ), & ! Input
RTV%Planck_Atmosphere( k ), & !Input
s_trans_AD( :, :, k ), s_refl_AD( :, :, k ), s_source_up_AD( :, k ), &
s_source_down_AD( :, k ), RTV, w_AD( k ), T_OD_AD( k ), Pff_AD( :, :, k ), &
Pbb_AD( :, :, k ), Planck_Atmosphere_AD( k ) ) !Output
END IF
ELSE ! Thermal sources
DO i = 1, RTV%n_Angles
s_Source_UP_AD( i, k ) = s_Source_UP_AD( i, k ) + s_Source_DOWN_AD( i, k )
s_Source_DOWN_AD( i, k ) = ZERO
Planck_Atmosphere_AD( k ) = Planck_Atmosphere_AD( k ) + s_Source_UP_AD( i, k ) * ( ONE - RTV%e_Layer_Trans( i, k ) )
e_Trans_AD( i, k ) = e_Trans_AD( i, k ) - RTV%Planck_Atmosphere( k ) * s_Source_UP_AD( i, k )
s_Source_UP_AD( i, k ) = ZERO
END DO
END IF
T_OD_AD( k ) = T_OD_AD( k ) - SUM( e_Trans_AD( :, k ) * RTV%e_Layer_Trans( 1:RTV%n_Angles, k ) / &
RTV%COS_Angle(1:RTV%n_Angles) )
e_Trans_AD( 1 : RTV%n_Angles, k ) = ZERO
END DO
END SUBROUTINE CRTM_SOI_AD
!--------------------------------------------------------------------------------
!:sdoc+:
!
! NAME:
! CRTM_SOI_Version
!
! PURPOSE:
! Subroutine to return the module version information.
!
! CALLING SEQUENCE:
! CALL CRTM_SOI_Version( Id )
!
! OUTPUT ARGUMENTS:
! Id: Character string containing the version Id information
! for the module.
! UNITS: N/A
! TYPE: CHARACTER(*)
! DIMENSION: Scalar
! ATTRIBUTES: INTENT(OUT)
!
!:sdoc-:
!--------------------------------------------------------------------------------
SUBROUTINE CRTM_SOI_Version( Id )
CHARACTER(*), INTENT(OUT) :: Id
Id = MODULE_VERSION_ID
END SUBROUTINE CRTM_SOI_Version
!################################################################################
!################################################################################
!## ##
!## ## PRIVATE MODULE ROUTINES ## ##
!## ##
!################################################################################
!################################################################################
SUBROUTINE CRTM_Truncated_Doubling( n_streams, & ! Input, number of streams
NANG, & ! Input, number of angles
KL, & ! Input, KL-th layer
single_albedo, & ! Input, single scattering albedo
optical_depth, & ! Input, layer optical depth
COS_Angle, & ! Input, COSINE of ANGLES
COS_Weight, & ! Input, GAUSSIAN Weights
ff, & ! Input, Phase matrix (forward part)
bb, & ! Input, Phase matrix (backward part)
Planck_Func, & ! Input, Planck for layer temperature
RTV) ! Output, layer transmittance, reflectance, and source
! ---------------------------------------------------------------------------------------
! FUNCTION
! Compute layer transmission, reflection matrices and source function
! at the top and bottom of the layer.
!
! Method and References
! It is a common doubling method and its theoretical basis is referred to
! Hansen, J. E., 1971: Multiple scattering of polarized light in
! planetary atmosphere. Part I. The doubling method, J. ATmos. Sci., 28, 120-125.
!
! see also ADA method.
! Quanhua Liu
! Quanhua.Liu@noaa.gov
!
! Modified from CRTM_Doubling_Layer by Tom Greenwald tomg@ssec.wisc.edu
! ----------------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_streams, NANG, KL
REAL(fp), INTENT(IN) :: single_albedo, optical_depth, Planck_Func
REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff, bb
TYPE(RTV_type), INTENT( INOUT ) :: RTV
! internal variables
REAL(fp), DIMENSION(NANG,NANG) :: term2,term3,trans,refl
REAL(fp), DIMENSION(NANG) :: C1, C2, source_up,source_down
REAL(fp) :: s, c
INTEGER :: i,j,k,L
! Check for tiny optical depth
IF ( optical_depth < OPTICAL_DEPTH_THRESHOLD ) THEN
RTV%s_Layer_Trans(1:NANG,1:NANG,KL) = ZERO
DO i = 1, NANG
RTV%s_Layer_Trans(i,i,KL) = ONE
ENDDO
RTV%s_Layer_Refl( 1 : NANG, 1 : NANG, KL ) = ZERO
RTV%s_Layer_Source_DOWN( 1 : NANG, KL ) = ZERO
RTV%s_Layer_Source_UP( 1 : NANG, KL ) = ZERO
RETURN
ENDIF
! -------------------------------------------------------------- !
! Determine number of doubling steps and construct !
! initial transmission and reflection matrix !
! --------------------------------------------------------------!
RTV%Number_Doubling( KL ) = INT( log( optical_depth / DELTA_OPTICAL_DEPTH ) / log( TWO ) ) + 1
IF ( RTV%Number_Doubling( KL ) < 1 ) RTV%Number_Doubling( KL ) = 1
IF ( RTV%Number_Doubling( KL ) <= MAX_N_DOUBLING ) THEN
RTV%Delta_Tau( KL ) = optical_depth / ( TWO**RTV%Number_Doubling( KL ) )
ELSE
RTV%Number_Doubling( KL ) = MAX_N_DOUBLING
RTV%Delta_Tau( KL ) = DELTA_OPTICAL_DEPTH
ENDIF
s = RTV%Delta_Tau( KL ) * single_albedo
DO i = 1, NANG
c = s / COS_Angle( i )
DO j = 1, NANG
RTV%Refl( i, j, 0, KL ) = c * bb( i, j ) * COS_Weight( j )
RTV%Trans( i, j, 0, KL ) = c * ff( i, j ) * COS_Weight( j )
ENDDO
RTV%Trans( i, i, 0, KL ) = RTV%Trans( i, i, 0, KL ) + ONE - RTV%Delta_Tau( KL ) / COS_Angle( i )
ENDDO
! -------------------------------------------------------------- !
! Doubling divided sub-layers !
! --------------------------------------------------------------!
DO L = 1, RTV%Number_Doubling( KL )
DO i = 1, NANG
DO j = 1, NANG
RTV%Inv_BeT( i, j, L, KL ) = ZERO
DO k = 1, NANG
! Add option to select two terms (i.e., RR+RR**2) if opd > 4 and ssa > 0.8
RTV%Inv_BeT( i, j, L, KL ) = RTV%Inv_BeT( i, j, L, KL ) + RTV%Refl( i, k, L - 1, KL ) * RTV%Refl( k, j, L - 1, KL )
ENDDO
ENDDO
RTV%Inv_BeT( i, i, L, KL ) = RTV%Inv_BeT( i, i, L, KL ) + ONE
ENDDO
term2 = matmul( RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ), RTV%Inv_BeT( 1 : NANG, 1 : NANG, L, KL ) )
term3 = matmul( term2, RTV%Refl( 1 : NANG, 1 : NANG, L - 1, KL ) )
term3 = matmul( term3, RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ) )
RTV%Refl( 1 : NANG, 1 : NANG, L, KL ) = RTV%Refl( 1 : NANG, 1 : NANG, L - 1, KL ) + term3
RTV%Trans( 1 : NANG, 1 : NANG, L, KL ) = matmul( term2, RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ) )
ENDDO
trans = RTV%Trans( 1 : NANG, 1 : NANG, RTV%Number_Doubling( KL ), KL )
refl = RTV%Refl( 1 : NANG, 1 : NANG, RTV%Number_Doubling( KL ), KL )
!
! computing source function at up and downward directions.
DO i = 1, NANG
C1( i ) = ZERO
C2( i ) = ZERO
DO j = 1, n_Streams
C1( i ) = C1( i ) + trans( i, j )
C2( i ) = C2( i ) + refl( i, j )
ENDDO
IF ( i == NANG .AND. NANG == ( n_Streams + 1 ) ) THEN
C1( i ) = C1( i ) + trans( NANG, NANG )
ENDIF
ENDDO
DO i = 1, NANG
source_up( i ) = ( ONE - C1( i ) - C2( i ) ) * Planck_Func
source_down( i ) = source_up( i )
ENDDO
RTV%C1( 1 : NANG, KL ) = C1
RTV%C2( 1 : NANG, KL ) = C2
RTV%s_Layer_Trans( 1 : NANG, 1 : NANG, KL ) = trans
RTV%s_Layer_Refl( 1 : NANG, 1 : NANG, KL ) = refl
RTV%s_Layer_Source_DOWN( 1 : NANG, KL ) = source_down
RTV%s_Layer_Source_UP( 1 : NANG, KL ) = source_up
RETURN
END SUBROUTINE CRTM_Truncated_Doubling
SUBROUTINE CRTM_Truncated_Doubling_TL(n_streams, & ! Input, number of streams
NANG, & ! Input, number of angles
KL, & ! Input, number of angles
single_albedo, & ! Input, single scattering albedo
optical_depth, & ! Input, layer optical depth
COS_Angle, & ! Input, COSINE of ANGLES
COS_Weight, & ! Input, GAUSSIAN Weights
ff, & ! Input, Phase matrix (forward part)
bb, & ! Input, Phase matrix (backward part)
Planck_Func, & ! Input, Planck for layer temperature
single_albedo_TL, & ! Input, tangent-linear single albedo
optical_depth_TL, & ! Input, TL layer optical depth
ff_TL, & ! Input, TL forward Phase matrix
bb_TL, & ! Input, TL backward Phase matrix
Planck_Func_TL, & ! Input, TL Planck for layer temperature
RTV, & ! Input, structure containing forward results
trans_TL, & ! Output, layer tangent-linear trans
refl_TL, & ! Output, layer tangent-linear refl
source_up_TL, & ! Output, layer tangent-linear source_up
source_down_TL) ! Output, layer tangent-linear source_down
! ---------------------------------------------------------------------------------------
! FUNCTION
! Compute tangent-linear layer transmission, reflection matrices and source function
! at the top and bottom of the layer.
!
! Tom Greenwald tomg@ssec.wisc.edu
! ----------------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_streams, NANG, KL
TYPE(RTV_type), INTENT(IN) :: RTV
REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff, bb
REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
REAL(fp), INTENT(IN) :: single_albedo, optical_depth, Planck_Func
! Tangent-Linear Part
REAL(fp), INTENT(OUT), DIMENSION( :, : ) :: trans_TL, refl_TL
REAL(fp), INTENT(OUT), DIMENSION( : ) :: source_up_TL, source_down_TL
REAL(fp), INTENT(IN) :: single_albedo_TL
REAL(fp), INTENT(IN) :: optical_depth_TL, Planck_Func_TL
REAL(fp), INTENT(IN), DIMENSION( : ,: ) :: ff_TL, bb_TL
! internal variables
REAL(fp), DIMENSION( NANG, NANG ) :: term2, term3, term2_TL, term3_TL, ms_term_TL
REAL(fp) :: s, c
REAL(fp) :: s_TL, c_TL, Delta_Tau_TL
REAL(fp), DIMENSION( NANG ) :: C1_TL, C2_TL
INTEGER :: i, j, k, L
! Tangent-Linear Beginning
IF( optical_depth < OPTICAL_DEPTH_THRESHOLD ) THEN
trans_TL = ZERO
refl_TL = ZERO
source_up_TL = ZERO
source_down_TL = ZERO
RETURN
ENDIF
s = RTV%Delta_Tau( KL ) * single_albedo
IF ( RTV%Number_Doubling( KL ) <= MAX_N_DOUBLING ) THEN
Delta_Tau_TL = optical_depth_TL / ( TWO**RTV%Number_Doubling( KL ) )
ELSE
Delta_Tau_TL = ZERO
ENDIF
s_TL = Delta_Tau_TL * single_albedo + RTV%Delta_Tau( KL ) * single_albedo_TL
DO i = 1, NANG
c = s/COS_Angle( i )
c_TL = s_TL/COS_Angle( i )
DO j = 1, NANG
refl_TL( i, j ) = c_TL * bb( i, j ) * COS_Weight( j ) + c * bb_TL( i, j ) * COS_Weight( j )
trans_TL( i, j ) = c_TL * ff( i, j ) * COS_Weight( j ) + c * ff_TL( i, j ) * COS_Weight( j )
ENDDO
trans_TL( i, i ) = trans_TL( i, i ) - Delta_Tau_TL / COS_Angle( i )
ENDDO
DO L = 1, RTV%Number_Doubling( KL )
DO i = 1, NANG
DO j = 1, NANG
ms_term_TL( i, j ) = ZERO
DO k = 1, NANG
ms_term_TL( i, j ) = ms_term_TL( i, j ) + RTV%Refl( i, k, L - 1, KL ) * Refl_TL( k, j ) &
+ Refl_TL( i, k ) * RTV%Refl( k, j, L - 1, KL )
END DO
END DO
END DO
term2 = matmul( RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ), RTV%Inv_BeT( 1 : NANG, 1 : NANG, L, KL ) )
term2_TL = matmul( RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ), ms_term_TL ) + &
matmul( Trans_TL, RTV%Inv_BeT( 1 : NANG, 1 : NANG, L, KL ) )
term3 = matmul( term2, RTV%Refl( 1 : NANG, 1 : NANG, L - 1, KL ) )
term3_TL = matmul( term2, Refl_TL ) + matmul( term2_TL, RTV%Refl( 1 : NANG, 1 : NANG, L - 1, KL ) )
term3 = matmul( term3, RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ) )
term3_TL = matmul( term3, Trans_TL ) + matmul( term3_TL, RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ) )
refl_TL = Refl_TL + term3_TL
trans_TL = matmul( term2, Trans_TL ) + matmul( term2_TL, RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ) )
ENDDO
! compute TL of source function at up and downward directions.
DO i = 1, NANG
C1_TL( i ) = ZERO
C2_TL( i ) = ZERO
DO j = 1, n_Streams
C1_TL( i ) = C1_TL( i ) + trans_TL( i, j )
C2_TL( i ) = C2_TL( i ) + refl_TL( i, j )
ENDDO
IF ( i == NANG .AND. NANG == ( n_Streams + 1 ) ) THEN
C1_TL( i ) = C1_TL( i ) + trans_TL( NANG, NANG )
ENDIF
ENDDO
DO i = 1, NANG
source_up_TL( i ) = - ( C1_TL( i ) + C2_TL( i ) ) * Planck_Func &
+ ( ONE - RTV%C1( i, KL ) - RTV%C2( i, KL ) ) * Planck_Func_TL
source_down_TL( i ) = source_up_TL( i )
END DO
RETURN
END SUBROUTINE CRTM_Truncated_Doubling_TL
SUBROUTINE CRTM_Truncated_Doubling_AD(n_streams, & ! Input, number of streams
NANG, & ! Input, number of angles
KL, & ! Input, number of angles
single_albedo, & ! Input, single scattering albedo
optical_depth, & ! Input, layer optical depth
COS_Angle, & ! Input, COSINE of ANGLES
COS_Weight, & ! Input, GAUSSIAN Weights
ff, & ! Input, Phase matrix (forward part)
bb, & ! Input, Phase matrix (backward part)
Planck_Func, & ! Input, Planck for layer temperature
trans_AD, & ! Input, layer tangent-linear trans
refl_AD, & ! Input, layer tangent-linear refl
source_up_AD, & ! Input, layer tangent-linear source_up
source_down_AD, & ! Input, layer tangent-linear source_down
RTV, & ! Input, structure containing forward results
single_albedo_AD, & ! Output adjoint single scattering albedo
optical_depth_AD, & ! Output AD layer optical depth
ff_AD, & ! Output AD forward Phase matrix
bb_AD, & ! Output AD backward Phase matrix
Planck_Func_AD) ! Output AD Planck for layer temperature
INTEGER, INTENT(IN) :: n_streams,NANG,KL
TYPE(RTV_type), INTENT(IN) :: RTV
REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
REAL(fp), INTENT(IN) :: single_albedo,optical_depth,Planck_Func
! Tangent-Linear Part
REAL(fp), INTENT( INOUT ), DIMENSION( :,: ) :: trans_AD,refl_AD
REAL(fp), INTENT( INOUT ), DIMENSION( : ) :: source_up_AD,source_down_AD
REAL(fp), INTENT( INOUT ) :: single_albedo_AD
REAL(fp), INTENT( INOUT ) :: optical_depth_AD,Planck_Func_AD
REAL(fp), INTENT(INOUT), DIMENSION(:,:) :: ff_AD,bb_AD
! internal variables
REAL(fp), DIMENSION( NANG, NANG ) :: term2, term3, term2_AD, term3_AD, ms_term_AD
REAL(fp) :: s, c
REAL(fp) :: s_AD, c_AD, Delta_Tau_AD
REAL(fp), DIMENSION(NANG) :: C1_AD, C2_AD
INTEGER :: i, j, L
! Adjoint Beginning
IF ( optical_depth < OPTICAL_DEPTH_THRESHOLD ) THEN
trans_AD = ZERO
refl_AD = ZERO
source_up_AD = ZERO
source_down_AD = ZERO
RETURN
ENDIF
! Remember, all output adjoint variables should be set to zero
! single_albedo_AD = 0, etc.
DO i = NANG, 1, -1
source_up_AD( i ) = source_up_AD( i ) + source_down_AD( i )
source_down_AD( i ) = ZERO
C2_AD( i ) = -source_up_AD( i ) * Planck_Func
C1_AD( i ) = -source_up_AD( i ) * Planck_Func
Planck_Func_AD = Planck_Func_AD + ( ONE - RTV%C1( i, KL ) - RTV%C2( i, KL ) ) * source_up_AD( i )
END DO
! Compute the source function in the up and downward directions.
DO i = NANG, 1, -1
IF(i == NANG .AND. NANG == ( n_Streams + 1 ) ) THEN
trans_AD( NANG, NANG ) = trans_AD( NANG, NANG ) + C1_AD( i )
ENDIF
DO j = n_Streams, 1, -1
refl_AD( i, j ) = refl_AD( i, j ) + C2_AD( i )
trans_AD( i, j ) = trans_AD( i, j ) + C1_AD( i )
END DO
END DO
term2_AD = ZERO
term3_AD = ZERO
ms_term_AD = ZERO
DO L = RTV%Number_Doubling(KL), 1, -1
! Forward calculations
term2 = matmul( RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ), RTV%Inv_BeT( 1 : NANG, 1 : NANG, L, KL ) )
term3 = matmul( term2, RTV%Refl( 1 : NANG, 1 : NANG, L - 1, KL ) )
term3 = matmul( term3, RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ) )
! Back to adjoint
term2_AD = term2_AD + matmul( trans_AD, transpose( RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ) ) )
trans_AD = matmul( transpose( term2 ), trans_AD )
term3_AD = term3_AD + refl_AD
trans_AD = trans_AD + matmul( transpose( term3 ), term3_AD )
term3_AD = matmul( term3_AD, transpose( RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ) ) )
term2_AD = term2_AD + matmul( term3_AD, transpose( RTV%Refl( 1 : NANG, 1 : NANG, L - 1, KL ) ) )
refl_AD = refl_AD + matmul( transpose( term2 ), term3_AD )
term3_AD = ZERO
ms_term_AD = ms_term_AD + matmul( transpose( RTV%Trans( 1 : NANG, 1 : NANG, L - 1, KL ) ), term2_AD )
trans_AD = trans_AD + matmul( term2_AD, transpose( RTV%Inv_BeT( 1 : NANG, 1 : NANG, L, KL ) ) )
term2_AD = ZERO
refl_AD = refl_AD + matmul( ms_term_AD, transpose( RTV%Refl( 1 : NANG, 1 : NANG, L - 1, KL ) ) ) + &
matmul( transpose( RTV%Refl( 1 : NANG, 1 : NANG, L - 1, KL ) ), ms_term_AD )
ms_term_AD = ZERO
ENDDO
s = RTV%Delta_Tau( KL ) * single_albedo
c_AD = ZERO
s_AD = ZERO
Delta_Tau_AD = ZERO
DO i = NANG, 1, -1
c = s / COS_Angle( i )
Delta_Tau_AD = Delta_Tau_AD - trans_AD( i, i ) / COS_Angle( i )
DO j = NANG, 1, -1
c_AD = c_AD + trans_AD( i, j ) * ff( i, j ) * COS_Weight( j )
ff_AD( i, j ) = ff_AD( i, j ) + trans_AD( i, j ) * c * COS_Weight( j )
c_AD = c_AD + refl_AD( i, j ) * bb( i, j ) * COS_Weight( j )
bb_AD( i, j ) = bb_AD( i, j ) + refl_AD( i, j ) * c * COS_Weight( j )
END DO
s_AD = s_AD + c_AD / COS_Angle( i )
c_AD = ZERO
ENDDO
Delta_Tau_AD = Delta_Tau_AD + s_AD * single_albedo
single_albedo_AD = single_albedo_AD + RTV%Delta_Tau( KL ) * s_AD
optical_depth_AD = optical_depth_AD + Delta_Tau_AD / ( TWO**RTV%Number_Doubling( KL ) )
END SUBROUTINE CRTM_Truncated_Doubling_AD
SUBROUTINE CRTM_Doubling_layer(n_streams, & ! Input, number of streams
NANG, & ! Input, number of angles
KL, & ! Input, KL-th layer
single_albedo, & ! Input, single scattering albedo
optical_depth, & ! Input, layer optical depth
COS_Angle, & ! Input, COSINE of ANGLES
COS_Weight, & ! Input, GAUSSIAN Weights
ff, & ! Input, Phase matrix (forward part)
bb, & ! Input, Phase matrix (backward part)
Planck_Func, & ! Input, Planck for layer temperature
RTV) ! Output, layer transmittance, reflectance, and source
! ---------------------------------------------------------------------------------------
! FUNCTION
! Compute layer transmission, reflection matrices and source function
! at the top and bottom of the layer.
!
! Method and References
! It is a common doubling method and its theoretical basis is referred to
! Hansen, J. E., 1971: Multiple scattering of polarized light in
! planetary atmosphere. Part I. The doubling method, J. ATmos. Sci., 28, 120-125.
!
! see also ADA method.
! Quanhua Liu
! Quanhua.Liu@noaa.gov
! ----------------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_streams,NANG,KL
TYPE(RTV_type), INTENT( INOUT ) :: RTV
REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
REAL(fp), INTENT(IN) :: single_albedo,optical_depth,Planck_Func
! internal variables
REAL(fp), DIMENSION(NANG,NANG) :: term2,term3,term4,trans,refl
REAL(fp), DIMENSION(NANG) :: C1, C2, source_up,source_down
REAL(fp) :: s, c
INTEGER :: i,j,k,L
INTEGER :: Error_Status
!
! Forward part beginning
IF( optical_depth < OPTICAL_DEPTH_THRESHOLD ) THEN
RTV%s_Layer_Trans(1:NANG,1:NANG,KL) = ZERO
DO i = 1, NANG
RTV%s_Layer_Trans(i,i,KL) = ONE
ENDDO
RTV%s_Layer_Refl(1:NANG,1:NANG,KL) = ZERO
RTV%s_Layer_Source_DOWN(1:NANG,KL) = ZERO
RTV%s_Layer_Source_UP(1:NANG,KL) = ZERO
RETURN
ENDIF
! -------------------------------------------------------------- !
! Determining number of doubling processes and constructing !
! initial transmission and reflection matrix
! --------------------------------------------------------------!
RTV%Number_Doubling(KL)=INT(log(optical_depth/DELTA_OPTICAL_DEPTH)/log(TWO))+1
IF( RTV%Number_Doubling(KL) < 1 ) RTV%Number_Doubling(KL) = 1
IF( RTV%Number_Doubling(KL) <= MAX_N_DOUBLING ) THEN
RTV%Delta_Tau(KL) = optical_depth/(TWO**RTV%Number_Doubling(KL))
ELSE
RTV%Number_Doubling(KL)=MAX_N_DOUBLING
RTV%Delta_Tau(KL) = DELTA_OPTICAL_DEPTH
ENDIF
s = RTV%Delta_Tau(KL) * single_albedo
DO i = 1, NANG
c = s/COS_Angle(i)
DO j = 1, NANG
RTV%Refl(i,j,0,KL) = c * bb(i,j) * COS_Weight(j)
RTV%Trans(i,j,0,KL) = c * ff(i,j) * COS_Weight(j)
ENDDO
RTV%Trans(i,i,0,KL) = RTV%Trans(i,i,0,KL) + ONE - RTV%Delta_Tau(KL)/COS_Angle(i)
ENDDO
! -------------------------------------------------------------- !
! Doubling divided sub-layers !
! --------------------------------------------------------------!
DO L = 1, RTV%Number_Doubling(KL)
DO i = 1, NANG
DO j = 1, NANG
term4(i,j) = ZERO
DO k = 1, NANG
term4(i,j) = term4(i,j) - RTV%Refl(i,k,L-1,KL)*RTV%Refl(k,j,L-1,KL)
ENDDO
ENDDO
term4(i,i) = term4(i,i) + ONE
ENDDO
RTV%Inv_BeT(1:NANG,1:NANG,L,KL) = matinv(term4, Error_Status)
IF( Error_Status /= SUCCESS ) THEN
print *,' error at matinv in CRTM_Doubling_layer '
RTV%s_Layer_Trans(1:NANG,1:NANG,KL) = ZERO
DO i = 1, NANG
RTV%s_Layer_Trans(i,i,KL) = exp(-optical_depth/COS_Angle(i))
ENDDO
RTV%s_Layer_Refl(1:NANG,1:NANG,KL) = ZERO
RTV%s_Layer_Source_DOWN(1:NANG,KL) = ZERO
RTV%s_Layer_Source_UP(1:NANG,KL) = ZERO
RETURN
ENDIF
term2 = matmul(RTV%Trans(1:NANG,1:NANG,L-1,KL), RTV%Inv_BeT(1:NANG,1:NANG,L,KL))
term3 = matmul(term2, RTV%Refl(1:NANG,1:NANG,L-1,KL))
term3 = matmul(term3, RTV%Trans(1:NANG,1:NANG,L-1,KL))
RTV%Refl(1:NANG,1:NANG,L,KL) = RTV%Refl(1:NANG,1:NANG,L-1,KL) + term3
RTV%Trans(1:NANG,1:NANG,L,KL) = matmul(term2, RTV%Trans(1:NANG,1:NANG,L-1,KL))
ENDDO
trans = RTV%Trans(1:NANG,1:NANG,RTV%Number_Doubling(KL),KL)
refl = RTV%Refl(1:NANG,1:NANG,RTV%Number_Doubling(KL),KL)
!
! computing source function at up and downward directions.
DO i = 1, NANG
C1(i) = ZERO
C2(i) = ZERO
DO j = 1, n_Streams
C1(i) = C1(i) + trans(i,j)
C2(i) = C2(i) + refl(i,j)
ENDDO
IF(i == NANG .AND. NANG == (n_Streams+1)) THEN
C1(i) = C1(i)+trans(NANG,NANG)
ENDIF
ENDDO
DO i = 1, NANG
source_up(i) = (ONE-C1(i)-C2(i))*Planck_Func
source_down(i) = source_up(i)
ENDDO
RTV%C1( 1:NANG,KL ) = C1
RTV%C2( 1:NANG,KL ) = C2
RTV%s_Layer_Trans(1:NANG,1:NANG,KL) = trans
RTV%s_Layer_Refl(1:NANG,1:NANG,KL) = refl
RTV%s_Layer_Source_DOWN(1:NANG,KL) = source_down
RTV%s_Layer_Source_UP(1:NANG,KL) = source_up
RETURN
END SUBROUTINE CRTM_Doubling_layer
!
!
SUBROUTINE CRTM_Doubling_layer_TL(n_streams, & ! Input, number of streams
NANG, & ! Input, number of angles
KL, & ! Input, number of angles
single_albedo, & ! Input, single scattering albedo
optical_depth, & ! Input, layer optical depth
COS_Angle, & ! Input, COSINE of ANGLES
COS_Weight, & ! Input, GAUSSIAN Weights
ff, & ! Input, Phase matrix (forward part)
bb, & ! Input, Phase matrix (backward part)
Planck_Func, & ! Input, Planck for layer temperature
single_albedo_TL, & ! Input, tangent-linear single albedo
optical_depth_TL, & ! Input, TL layer optical depth
ff_TL, & ! Input, TL forward Phase matrix
bb_TL, & ! Input, TL backward Phase matrix
Planck_Func_TL, & ! Input, TL Planck for layer temperature
RTV, & ! Input, structure containing forward results
trans_TL, & ! Output, layer tangent-linear trans
refl_TL, & ! Output, layer tangent-linear refl
source_up_TL, & ! Output, layer tangent-linear source_up
source_down_TL) ! Output, layer tangent-linear source_down
! ---------------------------------------------------------------------------------------
! FUNCTION
! Compute tangent-linear layer transmission, reflection matrices and source function
! at the top and bottom of the layer.
!
! Quanhua Liu
! Quanhua.Liu@noaa.gov
! ----------------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_streams,NANG,KL
TYPE(RTV_type), INTENT(IN) :: RTV
REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
REAL(fp), INTENT(IN) :: single_albedo,optical_depth,Planck_Func
! Tangent-Linear Part
REAL(fp), INTENT(OUT), DIMENSION( :,: ) :: trans_TL,refl_TL
REAL(fp), INTENT(OUT), DIMENSION( : ) :: source_up_TL,source_down_TL
REAL(fp), INTENT(IN) :: single_albedo_TL
REAL(fp), INTENT(IN) :: optical_depth_TL,Planck_Func_TL
REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff_TL,bb_TL
! internal variables
REAL(fp), DIMENSION(NANG,NANG) :: term1,term2,term3,term4,term5_TL
REAL(fp) :: s, c
REAL(fp) :: s_TL, c_TL, Delta_Tau_TL
REAL(fp), DIMENSION(NANG) :: C1_TL, C2_TL
INTEGER :: i,j,L
!
! Tangent-Linear Beginning
IF( optical_depth < OPTICAL_DEPTH_THRESHOLD ) THEN
trans_TL = ZERO
refl_TL = ZERO
source_up_TL = ZERO
source_down_TL = ZERO
RETURN
ENDIF
s = RTV%Delta_Tau(KL) * single_albedo
IF ( RTV%Number_Doubling( KL ) <= MAX_N_DOUBLING ) THEN
Delta_Tau_TL = optical_depth_TL / (TWO**RTV%Number_Doubling( KL ) )
ELSE
Delta_Tau_TL = ZERO
ENDIF
s_TL = Delta_Tau_TL * single_albedo + RTV%Delta_Tau(KL) * single_albedo_TL
DO i = 1, NANG
c = s/COS_Angle(i)
c_TL = s_TL/COS_Angle(i)
DO j = 1, NANG
refl_TL(i,j) = c_TL*bb(i,j)*COS_Weight(j)+c*bb_TL(i,j)*COS_Weight(j)
trans_TL(i,j) =c_TL*ff(i,j)*COS_Weight(j)+c*ff_TL(i,j)*COS_Weight(j)
ENDDO
trans_TL(i,i) =trans_TL(i,i) - Delta_Tau_TL/COS_Angle(i)
ENDDO
DO L = 1, RTV%Number_Doubling(KL)
term1 = matmul(RTV%Trans(1:NANG,1:NANG,L-1,KL),RTV%Inv_BeT(1:NANG,1:NANG,L,KL))
term2 = matmul(RTV%Inv_BeT(1:NANG,1:NANG,L,KL),RTV%Refl(1:NANG,1:NANG,L-1,KL))
term3 = matmul(RTV%Inv_BeT(1:NANG,1:NANG,L,KL),RTV%Trans(1:NANG,1:NANG,L-1,KL))
term4 = matmul(term2,RTV%Trans(1:NANG,1:NANG,L-1,KL))
term5_TL = matmul(refl_TL,RTV%Refl(1:NANG,1:NANG,L-1,KL)) &
+ matmul(RTV%Refl(1:NANG,1:NANG,L-1,KL),refl_TL)
refl_TL=refl_TL+matmul(matmul(term1,term5_TL),term4)+matmul(trans_TL,term4) &
+matmul(matmul(term1,refl_TL),RTV%Trans(1:NANG,1:NANG,L-1,KL))
refl_TL=refl_TL+matmul(matmul(term1,RTV%Refl(1:NANG,1:NANG,L-1,KL)),trans_TL)
trans_TL=matmul(trans_TL,term3) &
+matmul(matmul(term1,term5_TL),term3)+matmul(term1,trans_TL)
ENDDO
! computing source function at up and downward directions.
DO i = 1, NANG
C1_TL(i) = ZERO
C2_TL(i) = ZERO
DO j = 1, n_Streams
C1_TL(i) = C1_TL(i) + trans_TL(i,j)
C2_TL(i) = C2_TL(i) + refl_TL(i,j)
ENDDO
IF(i == NANG .AND. NANG == (n_Streams+1)) THEN
C1_TL(i) = C1_TL(i)+trans_TL(NANG,NANG)
ENDIF
ENDDO
DO i = 1, NANG
source_up_TL(i) = -(C1_TL(i)+C2_TL(i))*Planck_Func &
+ (ONE-RTV%C1(i,KL)-RTV%C2(i,KL))*Planck_Func_TL
source_down_TL(i) = source_up_TL(i)
END DO
RETURN
END SUBROUTINE CRTM_Doubling_layer_TL
! ---------------------------------------------------------------------------------------
! FUNCTION
! Compute layer adjoint transmission, reflection matrices and source function
! at the top and bottom of the layer.
!
! Quanhua Liu
! Quanhua.Liu@noaa.gov
! ----------------------------------------------------------------------------------------
SUBROUTINE CRTM_Doubling_layer_AD(n_streams, & ! Input, number of streams
NANG, & ! Input, number of angles
KL, & ! Input, number of angles
single_albedo, & ! Input, single scattering albedo
optical_depth, & ! Input, layer optical depth
COS_Angle, & ! Input, COSINE of ANGLES
COS_Weight, & ! Input, GAUSSIAN Weights
ff, & ! Input, Phase matrix (forward part)
bb, & ! Input, Phase matrix (backward part)
Planck_Func, & ! Input, Planck for layer temperature
trans_AD, & ! Input, layer tangent-linear trans
refl_AD, & ! Input, layer tangent-linear refl
source_up_AD, & ! Input, layer tangent-linear source_up
source_down_AD, & ! Input, layer tangent-linear source_down
RTV, & ! Input, structure containing forward results
single_albedo_AD, & ! Output adjoint single scattering albedo
optical_depth_AD, & ! Output AD layer optical depth
ff_AD, & ! Output AD forward Phase matrix
bb_AD, & ! Output AD backward Phase matrix
Planck_Func_AD) ! Output AD Planck for layer temperature
INTEGER, INTENT(IN) :: n_streams,NANG,KL
TYPE(RTV_type), INTENT(IN) :: RTV
REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
REAL(fp), INTENT(IN) :: single_albedo,optical_depth,Planck_Func
! Tangent-Linear Part
REAL(fp), INTENT( INOUT ), DIMENSION( :,: ) :: trans_AD,refl_AD
REAL(fp), INTENT( INOUT ), DIMENSION( : ) :: source_up_AD,source_down_AD
REAL(fp), INTENT( INOUT ) :: single_albedo_AD
REAL(fp), INTENT( INOUT ) :: optical_depth_AD,Planck_Func_AD
REAL(fp), INTENT(INOUT), DIMENSION(:,:) :: ff_AD,bb_AD
! internal variables
REAL(fp), DIMENSION(NANG,NANG) :: term1,term2,term3,term4,term5_AD
REAL(fp) :: s, c
REAL(fp) :: s_AD, c_AD, Delta_Tau_AD
REAL(fp), DIMENSION(NANG) :: C1_AD, C2_AD
INTEGER :: i,j,L
! Tangent-Linear Beginning
IF( optical_depth < OPTICAL_DEPTH_THRESHOLD ) THEN
trans_AD = ZERO
refl_AD = ZERO
source_up_AD = ZERO
source_down_AD = ZERO
RETURN
ENDIF
DO i = NANG, 1, -1
source_up_AD(i) = source_up_AD(i) + source_down_AD(i)
source_down_AD(i) = ZERO
C2_AD(i) = -source_up_AD(i)*Planck_Func
C1_AD(i) = -source_up_AD(i)*Planck_Func
Planck_Func_AD = Planck_Func_AD + (ONE-RTV%C1(i,KL)-RTV%C2(i,KL))*source_up_AD(i)
END DO
! Compute the source function in the up and downward directions.
DO i = NANG, 1, -1
IF(i == NANG .AND. NANG == (n_Streams+1)) THEN
trans_AD(NANG,NANG)=trans_AD(NANG,NANG)+C1_AD(i)
ENDIF
DO j = n_Streams, 1, -1
refl_AD(i,j)=refl_AD(i,j)+C2_AD(i)
trans_AD(i,j)=trans_AD(i,j)+C1_AD(i)
END DO
END DO
DO L = RTV%Number_Doubling(KL), 1, -1
term1 = MATMUL(RTV%Trans(1:NANG,1:NANG,L-1,KL),RTV%Inv_BeT(1:NANG,1:NANG,L,KL))
term2 = MATMUL(RTV%Inv_BeT(1:NANG,1:NANG,L,KL),RTV%Refl(1:NANG,1:NANG,L-1,KL))
term3 = MATMUL(RTV%Inv_BeT(1:NANG,1:NANG,L,KL),RTV%Trans(1:NANG,1:NANG,L-1,KL))
term4 = MATMUL(term2,RTV%Trans(1:NANG,1:NANG,L-1,KL))
term5_AD = MATMUL(MATMUL(TRANSPOSE(term1),trans_AD),TRANSPOSE(term3))
trans_AD = MATMUL(trans_AD,TRANSPOSE(term3))+MATMUL(TRANSPOSE(term1),trans_AD)
trans_AD=trans_AD+MATMUL(TRANSPOSE(MATMUL(term1,RTV%Refl(1:NANG,1:NANG,L-1,KL))),refl_AD)
term5_AD =term5_AD+MATMUL(MATMUL(TRANSPOSE(term1),refl_AD),TRANSPOSE(term4))
trans_AD = trans_AD+MATMUL(refl_AD,TRANSPOSE(term4))
refl_AD = refl_AD+MATMUL(MATMUL(TRANSPOSE(term1),refl_AD),TRANSPOSE(RTV%Trans(1:NANG,1:NANG,L-1,KL)))
refl_AD = refl_AD+MATMUL(term5_AD,TRANSPOSE(RTV%Refl(1:NANG,1:NANG,L-1,KL)))
refl_AD = refl_AD+MATMUL(TRANSPOSE(RTV%Refl(1:NANG,1:NANG,L-1,KL)),term5_AD)
ENDDO
s = RTV%Delta_Tau(KL) * single_albedo
c_AD = ZERO
s_AD = ZERO
Delta_Tau_AD=ZERO
DO i = NANG, 1, -1
c = s/COS_Angle(i)
Delta_Tau_AD = Delta_Tau_AD - trans_AD(i,i)/COS_Angle(i)
DO j = NANG, 1, -1
c_AD = c_AD + trans_AD(i,j)*ff(i,j)*COS_Weight(j)
ff_AD(i,j)=ff_AD(i,j)+trans_AD(i,j)*c*COS_Weight(j)
c_AD = c_AD + refl_AD(i,j)*bb(i,j)*COS_Weight(j)
bb_AD(i,j)=bb_AD(i,j) + refl_AD(i,j)*c*COS_Weight(j)
END DO
s_AD = s_AD + c_AD/COS_Angle(i)
c_AD = ZERO
ENDDO
Delta_Tau_AD = Delta_Tau_AD + s_AD* single_albedo
single_albedo_AD = single_albedo_AD+RTV%Delta_Tau(KL) * s_AD
optical_depth_AD = optical_depth_AD + Delta_Tau_AD/(TWO**RTV%Number_Doubling(KL))
END SUBROUTINE CRTM_Doubling_layer_AD
END MODULE SOI_Module