! ! 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