Skip to content

Add diagnostic diff reporting in regression tests + remove ifort array temporary warning #274

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
BenjaminTJohnson wants to merge 7 commits into
base: develop
Choose a base branch
from
Draft

Add diagnostic diff reporting in regression tests + remove ifort array temporary warning #274

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
d862f92
update to ctests to permit additional information on XYZ results are …
BenjaminTJohnson Dec 29, 2025
bef6476
updated symlinks to Compare_Diagnostics.inc to be relative paths
BenjaminTJohnson Dec 29, 2025
6460c23
modification to FitCoeff behavior to avoid warning messages forrtl: w…
BenjaminTJohnson Dec 29, 2025
86abb22
checkpoint
BenjaminTJohnson Dec 30, 2025
1243870
fixed RTV in Obs_4_downward usage, to RTV(nt).
BenjaminTJohnson Dec 30, 2025
efa2d72
removing check_crtm_random
BenjaminTJohnson Dec 30, 2025
8b3a25f
Zero subnormal K-matrix outputs
BenjaminTJohnson Jan 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 24 additions & 4 deletions src/AtmOptics/CRTM_AtmOptics.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,8 @@ MODULE CRTM_AtmOptics
! ---------------
! Message string length
INTEGER, PARAMETER :: ML = 256
REAL(fp), PARAMETER :: MIN_TRANSMITTANCE = TINY(ONE)
REAL(fp), PARAMETER :: MAX_OPTICAL_DEPTH = -LOG(MIN_TRANSMITTANCE)


CONTAINS
Expand Down Expand Up @@ -211,8 +213,14 @@ SUBROUTINE CRTM_Compute_Transmittance( atmoptics, transmittance )
TYPE(CRTM_AtmOptics_type), INTENT(IN) :: atmoptics
REAL(fp) , INTENT(OUT) :: transmittance
INTEGER :: k
REAL(fp) :: optical_depth
k = atmoptics%n_layers
transmittance = EXP(-ONE*SUM(atmoptics%optical_depth(1:k)))
optical_depth = SUM(atmoptics%optical_depth(1:k))
IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN
transmittance = MIN_TRANSMITTANCE
ELSE
transmittance = EXP(-optical_depth)
END IF
END SUBROUTINE CRTM_Compute_Transmittance

!--------------------------------------------------------------------------------
Expand Down Expand Up @@ -265,10 +273,16 @@ SUBROUTINE CRTM_Compute_Transmittance_TL( &
! Local variables
INTEGER :: k
REAL(fp) :: transmittance
REAL(fp) :: optical_depth

k = atmoptics%n_layers
transmittance = EXP(-ONE*SUM(atmoptics%optical_depth(1:k)))
transmittance_TL = -transmittance * SUM(atmoptics_TL%optical_depth(1:k))
optical_depth = SUM(atmoptics%optical_depth(1:k))
IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN
transmittance_TL = ZERO
ELSE
transmittance = EXP(-optical_depth)
transmittance_TL = -transmittance * SUM(atmoptics_TL%optical_depth(1:k))
END IF
END SUBROUTINE CRTM_Compute_Transmittance_TL


Expand Down Expand Up @@ -322,9 +336,15 @@ SUBROUTINE CRTM_Compute_Transmittance_AD( &
! Local variables
INTEGER :: k
REAL(fp) :: transmittance, t_delstar_t
REAL(fp) :: optical_depth

k = atmoptics%n_layers
transmittance = EXP(-ONE*SUM(atmoptics%optical_depth(1:k)))
optical_depth = SUM(atmoptics%optical_depth(1:k))
IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN
transmittance_AD = ZERO
RETURN
END IF
transmittance = EXP(-optical_depth)
t_delstar_t = transmittance*transmittance_AD
DO k = 1, atmoptics%n_layers
atmoptics_AD%optical_depth(k) = atmoptics_AD%optical_depth(k) - t_delstar_t
Expand Down
2 changes: 1 addition & 1 deletion src/CRTM_Forward_Module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -802,7 +802,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status )
WRITE( Message,'("Difference between Obs pressure level (",es22.15,&
&"hPa) and closest input profile level (",es22.15,&
&"hPa) is larger than recommended (",f4.1,"hPa) for profile #",i0)') &
Opt%Obs_4_downward_P, Atm%Level_Pressure(RTV%Obs_4_downward%idx), &
Opt%Obs_4_downward_P, Atm%Level_Pressure(RTV(nt)%Obs_4_downward%idx), &
AIRCRAFT_PRESSURE_THRESHOLD, m
CALL Display_Message( ROUTINE_NAME, Message, WARNING )
END IF
Expand Down
282 changes: 282 additions & 0 deletions src/CRTM_K_Matrix_Module.f90
View file
Open in desktop

Large diffs are not rendered by default.

4 changes: 3 additions & 1 deletion src/Coefficients/FitCoeff/FitCoeff_SetValue.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,13 @@
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'FitCoeff_SetValue'
! Local variables
CHARACTER(ML) :: msg
INTEGER :: dimensions(SIZE(self%Dimensions))
INTEGER :: i

! Allocate structure if necessary
IF ( .NOT. FitCoeff_Associated(self) ) THEN
CALL FitCoeff_Create( self, SHAPE(C) )
dimensions = SHAPE(C)
CALL FitCoeff_Create( self, dimensions )
IF ( .NOT. FitCoeff_Associated(self) ) THEN
msg = 'Allocation of FitCoeff structure failed'
CALL Display_Message( ROUTINE_NAME, msg, FAILURE ); RETURN
Expand Down
130 changes: 104 additions & 26 deletions src/RTSolution/ADA/ADA_Module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,9 +39,41 @@ MODULE ADA_Module
! -----------------
! Module parameters
! -----------------
REAL(fp), PARAMETER :: MIN_TRANSMITTANCE = TINY(ONE)
REAL(fp), PARAMETER :: MAX_OPTICAL_DEPTH = -LOG(MIN_TRANSMITTANCE)

CONTAINS

SUBROUTINE Safe_Exp_Neg(x, exp_x, is_clamped)
REAL(fp), INTENT(IN) :: x
REAL(fp), INTENT(OUT) :: exp_x
LOGICAL , INTENT(OUT) :: is_clamped

IF ( x >= MAX_OPTICAL_DEPTH ) THEN
exp_x = MIN_TRANSMITTANCE
is_clamped = .TRUE.
ELSE
exp_x = EXP(-x)
is_clamped = .FALSE.
END IF
END SUBROUTINE Safe_Exp_Neg


SUBROUTINE Safe_Exp_Neg_TL(x, x_TL, exp_x, exp_x_TL)
REAL(fp), INTENT(IN) :: x
REAL(fp), INTENT(IN) :: x_TL
REAL(fp), INTENT(OUT) :: exp_x
REAL(fp), INTENT(OUT) :: exp_x_TL
LOGICAL :: is_clamped

CALL Safe_Exp_Neg(x, exp_x, is_clamped)
IF ( is_clamped ) THEN
exp_x_TL = ZERO
ELSE
exp_x_TL = -exp_x * x_TL
END IF
END SUBROUTINE Safe_Exp_Neg_TL

!################################################################################
!################################################################################
!## ##
Expand All @@ -56,9 +88,11 @@ SUBROUTINE CRTM_SurfRef(n_Layers,total_od,dire,Index_Sat_Angle,Rsphere,Ref_0,Ref
TYPE(RTV_type), INTENT( INOUT ) :: RTV
REAL (fp), INTENT(IN) :: total_od
REAL (fp) :: Rsphere, Ref_0, Ref_1, Rfac, dire, SRef(3), S0a(3)
REAL (fp) :: optical_depth, exp_term
REAL (fp), DIMENSION(RTV%n_Angles*RTV%n_Stokes, RTV%n_Angles*RTV%n_Stokes) :: temporal_matrix,R0, T0, R2, R1
REAL (fp), DIMENSION( RTV%n_Angles*RTV%n_Stokes) :: refl_down, S0, Sun0
INTEGER :: i, j, k, L, Error_Status, nZ,Index_Sat_Angle
LOGICAL :: is_clamped
REAL(fp), PARAMETER :: rFactor(3) = (/ 0.0_fp, 0.5_fp, 1.0_fp/)
Sun0(:) = ZERO
DO i = 1, RTV%n_Angles
Expand All @@ -70,7 +104,9 @@ SUBROUTINE CRTM_SurfRef(n_Layers,total_od,dire,Index_Sat_Angle,Rsphere,Ref_0,Ref
DO 101 L = 1, 2
R0(:,:) = RTV%s_Level_Refl_UP(1:nZ,1:nZ,n_Layers)*rFactor(L)
SRef(L) = dire*rFactor(L)
S0(:)=(ONE-dire)*RTV%Planck_Surface+SRef(L)*RTV%COS_SUN*RTV%Solar_irradiance/PI*exp(-total_od/RTV%COS_SUN)*Sun0(:)
optical_depth = total_od / RTV%COS_SUN
CALL Safe_Exp_Neg(optical_depth, exp_term, is_clamped)
S0(:)=(ONE-dire)*RTV%Planck_Surface+SRef(L)*RTV%COS_SUN*RTV%Solar_irradiance/PI*exp_term*Sun0(:)
T0(:,:) = RTV%s_Layer_Trans(1:nZ,1:nZ,n_Layers)
DO k = n_Layers, 1, -1
temporal_matrix = -matmul(R0,RTV%s_Layer_Refl(1:nZ,1:nZ,k))
Expand Down Expand Up @@ -131,6 +167,8 @@ SUBROUTINE CRTM_ADA(n_Layers, & ! Input number of atmospheric layers
REAL (fp), DIMENSION(RTV%n_Angles*RTV%n_Stokes) :: temporal_vector
REAL (fp), DIMENSION(0:n_Layers) :: total_opt
INTEGER :: i, j, k, Error_Status
REAL(fp) :: optical_depth, exp_term
LOGICAL :: is_clamped
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_ADA'
CHARACTER(256) :: Message

Expand All @@ -156,8 +194,10 @@ SUBROUTINE CRTM_ADA(n_Layers, & ! Input number of atmospheric layers
END IF

IF( RTV%Solar_Flag_true ) THEN
optical_depth = total_opt(n_Layers) / RTV%COS_SUN
CALL Safe_Exp_Neg(optical_depth, exp_term, is_clamped)
RTV%s_Level_Rad_UP(1:nZ,n_Layers ) = RTV%s_Level_Rad_UP(1:nZ,n_Layers )+direct_reflectivity(1:nZ)* &
RTV%COS_SUN*RTV%Solar_irradiance/PI*exp(-total_opt(n_Layers)/RTV%COS_SUN)
RTV%COS_SUN*RTV%Solar_irradiance/PI*exp_term
END IF

! UPWARD ADDING LOOP STARTS FROM BOTTOM LAYER TO ATMOSPHERIC TOP LAYER.
Expand Down Expand Up @@ -226,7 +266,8 @@ SUBROUTINE CRTM_ADA(n_Layers, & ! Input number of atmospheric layers
matmul(RTV%Inv_GammaT(1:nZ,1:nZ,k),RTV%Refl_Trans(1:nZ,1:nZ,k))
ELSE
DO i = 1, nZ
RTV%s_Layer_Trans(i,i,k) = exp(-T_OD(k)/RTV%COS_AngleS(i))
optical_depth = T_OD(k)/RTV%COS_AngleS(i)
CALL Safe_Exp_Neg(optical_depth, RTV%s_Layer_Trans(i,i,k), is_clamped)
END DO
DO i = 1, nZ, RTV%n_Stokes
RTV%s_Layer_Source_UP(i,k) = RTV%Planck_Atmosphere(k) * (ONE - RTV%s_Layer_Trans(i,i,k) )
Expand Down Expand Up @@ -393,6 +434,7 @@ SUBROUTINE MOM(KL,nZ,optical_depth,trans,refl,RTV, Error_Status)
CHARACTER(256) :: Message
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MOM'
REAL(fp) :: xx
LOGICAL :: is_clamped
REAL(fp), DIMENSION(nZ,nZ) :: tempo

RTV%PPM(1:nZ,1:nZ,KL) = RTV%PP(1:nZ,1:nZ,KL) - RTV%PM(1:nZ,1:nZ,KL)
Expand Down Expand Up @@ -444,7 +486,7 @@ SUBROUTINE MOM(KL,nZ,optical_depth,trans,refl,RTV, Error_Status)

DO i = 1, nZ
xx = RTV%EigValue(i,KL)*optical_depth
RTV%Exp_x(i,KL) = exp(-xx)
CALL Safe_Exp_Neg(xx, RTV%Exp_x(i,KL), is_clamped)
END DO

DO i = 1, nZ
Expand Down Expand Up @@ -489,7 +531,7 @@ SUBROUTINE MOM_TL(KL,nZ,optical_depth,RTV,optical_depth_TL,PP_TL,PM_TL,trans_TL,
INTEGER :: i, j

CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MOM_TL'
REAL(fp) :: xx_TL
REAL(fp) :: xx_TL, xx
REAL(fp), DIMENSION(nZ,nZ) :: PPM_TL, i_PPM_TL, PPP_TL, HH_TL
REAL(fp), DIMENSION(nZ,nZ) :: EigVeF_TL,Gp_TL,Gm_TL,A1_TL,A2_TL,A3_TL,A4_TL,A5_TL,A6_TL
REAL(fp), DIMENSION(nZ,nZ) :: Gm_A5_TL,i_Gm_A5_TL,EigVe_TL,EigVeVa_TL,i_Gm_TL
Expand Down Expand Up @@ -529,7 +571,12 @@ SUBROUTINE MOM_TL(KL,nZ,optical_depth,RTV,optical_depth_TL,PP_TL,PM_TL,trans_TL,
i_Gm_TL = -matmul( RTV%i_Gm(1:nZ,1:nZ,KL), matmul(Gm_TL,RTV%i_Gm(1:nZ,1:nZ,KL)) )
DO i = 1, nZ
xx_TL = EigValue_TL(i)*optical_depth+RTV%EigValue(i,KL)*optical_depth_TL
Exp_x_TL(i) = -xx_TL*RTV%Exp_x(i,KL)
xx = RTV%EigValue(i,KL)*optical_depth
IF ( xx >= MAX_OPTICAL_DEPTH ) THEN
Exp_x_TL(i) = ZERO
ELSE
Exp_x_TL(i) = -xx_TL*RTV%Exp_x(i,KL)
END IF
END DO

DO i = 1, nZ
Expand Down Expand Up @@ -573,7 +620,7 @@ SUBROUTINE MOM_AD(KL,nZ,optical_depth,RTV, optical_depth_AD,PP_AD,PM_AD,trans_AD
INTEGER :: i, j

CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MOM_TL'
REAL(fp) :: xx_AD
REAL(fp) :: xx_AD, xx


REAL(fp), DIMENSION(nZ) :: Exp_x_AD,EigValue_AD,EigVa_AD
Expand Down Expand Up @@ -617,8 +664,14 @@ SUBROUTINE MOM_AD(KL,nZ,optical_depth,RTV, optical_depth_AD,PP_AD,PM_AD,trans_AD
END DO

DO i = nZ, 1, -1
xx_AD = -Exp_x_AD(i)*RTV%Exp_x(i,KL)
Exp_x_AD(i) = ZERO
xx = RTV%EigValue(i,KL)*optical_depth
IF ( xx >= MAX_OPTICAL_DEPTH ) THEN
Exp_x_AD(i) = ZERO
xx_AD = ZERO
ELSE
xx_AD = -Exp_x_AD(i)*RTV%Exp_x(i,KL)
Exp_x_AD(i) = ZERO
END IF
EigValue_AD(i) = xx_AD*optical_depth
optical_depth_AD = optical_depth_AD + RTV%EigValue(i,KL)*xx_AD
END DO
Expand Down Expand Up @@ -714,6 +767,8 @@ SUBROUTINE CRTM_AMOM_layer( n_streams, & ! Input, number of streams
INTEGER :: Error_Status
REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ)
REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ)
REAL(fp) :: optical_depth_sun, total_opt_sun
LOGICAL :: is_clamped
CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer'
CHARACTER(256) :: Message

Expand Down Expand Up @@ -807,8 +862,10 @@ SUBROUTINE CRTM_AMOM_layer( n_streams, & ! Input, number of streams
! Solar source
Sfactor = single_albedo*RTV%Solar_irradiance/PI
IF( RTV%mth_Azi == 0 ) Sfactor = Sfactor/TWO
EXPfactor = exp(-optical_depth/RTV%COS_SUN)
s_transmittance = exp(-total_opt/RTV%COS_SUN)
optical_depth_sun = optical_depth / RTV%COS_SUN
CALL Safe_Exp_Neg(optical_depth_sun, EXPfactor, is_clamped)
total_opt_sun = total_opt / RTV%COS_SUN
CALL Safe_Exp_Neg(total_opt_sun, s_transmittance, is_clamped)

DO i = 1, nZ
Solar(i) = -bb(i,nZ+1)*Sfactor
Expand Down Expand Up @@ -1290,6 +1347,7 @@ SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers
REAL (fp), DIMENSION( RTV%n_Angles*RTV%n_Stokes, RTV%n_Angles*RTV%n_Stokes ) :: s_trans_TL,s_refl_TL,Refl_Trans_TL
REAL (fp), DIMENSION( RTV%n_Angles*RTV%n_Stokes, RTV%n_Angles*RTV%n_Stokes ) :: s_refl_up_TL,Inv_Gamma_TL,Inv_GammaT_TL
REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_TL
REAL (fp) :: optical_depth, optical_depth_TL, exp_term, exp_term_TL
INTEGER :: i, j, k, nZ
!
nZ = RTV%n_Angles*RTV%n_Stokes
Expand All @@ -1309,10 +1367,12 @@ SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers
END IF

IF( RTV%Solar_Flag_true ) THEN
s_rad_up_TL = s_rad_up_TL+direct_reflectivity_TL(1:nZ)*RTV%COS_SUN*RTV%Solar_irradiance/PI &
* exp(-total_opt(n_Layers)/RTV%COS_SUN) &
- direct_reflectivity(1:nZ) * RTV%Solar_irradiance/PI &
* total_opt_TL(n_Layers) * exp(-total_opt(n_Layers)/RTV%COS_SUN)
optical_depth = total_opt(n_Layers) / RTV%COS_SUN
optical_depth_TL = total_opt_TL(n_Layers) / RTV%COS_SUN
CALL Safe_Exp_Neg_TL(optical_depth, optical_depth_TL, exp_term, exp_term_TL)
s_rad_up_TL = s_rad_up_TL + direct_reflectivity_TL(1:nZ)*RTV%COS_SUN*RTV%Solar_irradiance/PI &
* exp_term + direct_reflectivity(1:nZ) * RTV%Solar_irradiance/PI &
* RTV%COS_SUN * exp_term_TL
END IF

DO 10 k = n_Layers, 1, -1
Expand Down Expand Up @@ -1465,6 +1525,8 @@ SUBROUTINE CRTM_AMOM_layer_TL( n_streams, & ! Input, number of streams
REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ)
REAL(fp) :: EXPfactor_TL,Sfactor_TL,s_transmittance_TL,Solar_TL(2*nZ),V0_TL(2*nZ,2*nZ),Solar1_TL(2*nZ)
REAL(fp) :: Sfac2_TL
REAL(fp) :: optical_depth_sun, total_opt_sun
REAL(fp) :: optical_depth_sun_TL, total_opt_sun_TL
REAL(fp), DIMENSION( nZ ) :: thermal_up_TL,thermal_down_TL
REAL(fp), DIMENSION(nZ,nZ) :: trans, refl
INTEGER :: N2, N2_1
Expand Down Expand Up @@ -1581,11 +1643,13 @@ SUBROUTINE CRTM_AMOM_layer_TL( n_streams, & ! Input, number of streams
Sfactor = Sfactor/TWO
Sfactor_TL = Sfactor_TL/TWO
END IF
EXPfactor = exp(-optical_depth/RTV%COS_SUN)
EXPfactor_TL = -optical_depth_TL/RTV%COS_SUN*EXPfactor
optical_depth_sun = optical_depth / RTV%COS_SUN
optical_depth_sun_TL = optical_depth_TL / RTV%COS_SUN
CALL Safe_Exp_Neg_TL(optical_depth_sun, optical_depth_sun_TL, EXPfactor, EXPfactor_TL)

s_transmittance = exp(-total_opt/RTV%COS_SUN)
s_transmittance_TL = -total_opt_TL/RTV%COS_SUN*s_transmittance
total_opt_sun = total_opt / RTV%COS_SUN
total_opt_sun_TL = total_opt_TL / RTV%COS_SUN
CALL Safe_Exp_Neg_TL(total_opt_sun, total_opt_sun_TL, s_transmittance, s_transmittance_TL)

DO i = 1, nZ
Solar(i) = -bb(i,nZ+1)*Sfactor
Expand Down Expand Up @@ -1727,6 +1791,8 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers
Inv_Gamma_AD,Inv_GammaT_AD
REAL (fp) :: sum_s_AD, sums_AD, xx
REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_AD
REAL (fp) :: optical_depth, exp_term
LOGICAL :: is_clamped
INTEGER :: i, j, k,nZ
!

Expand Down Expand Up @@ -1853,9 +1919,13 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers

!
IF( RTV%Solar_Flag_true ) THEN
Copy link
Contributor

@chengdang chengdang Jan 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Ben, I'm trying to understand this set up. With such if statement, if the total AOD is too big (is_clamped), then the total_opt_AD is not updated?

xx = RTV%Solar_irradiance/PI * exp(-total_opt(n_Layers)/RTV%COS_SUN)
total_opt_AD(n_Layers) = total_opt_AD(n_Layers) &
- xx*sum(direct_reflectivity(1:nZ)*s_rad_up_AD(1:nZ))
optical_depth = total_opt(n_Layers) / RTV%COS_SUN
CALL Safe_Exp_Neg(optical_depth, exp_term, is_clamped)
xx = RTV%Solar_irradiance/PI * exp_term
IF ( .NOT. is_clamped ) THEN
total_opt_AD(n_Layers) = total_opt_AD(n_Layers) &
- xx*sum(direct_reflectivity(1:nZ)*s_rad_up_AD(1:nZ))
END IF
direct_reflectivity_AD(1:nZ) = direct_reflectivity_AD(1:nZ) &
+ s_rad_up_AD(1:nZ) * RTV%COS_SUN * xx
END IF
Expand Down Expand Up @@ -1940,6 +2010,8 @@ SUBROUTINE CRTM_AMOM_layer_AD( n_streams, & ! Input, number of streams
REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ)
REAL(fp) :: EXPfactor_AD,Sfactor_AD,s_transmittance_AD,Solar_AD(2*nZ),V0_AD(2*nZ,2*nZ),Solar1_AD(2*nZ)
REAL(fp) :: V1_AD(2*nZ,2*nZ),Sfac2_AD
REAL(fp) :: optical_depth_sun, total_opt_sun
LOGICAL :: is_clamped_exp, is_clamped_total
REAL(fp), DIMENSION(nZ,nZ) :: trans, refl
INTEGER :: N2, N2_1
REAL(fp) :: Thermal_C_AD
Expand Down Expand Up @@ -2011,8 +2083,10 @@ SUBROUTINE CRTM_AMOM_layer_AD( n_streams, & ! Input, number of streams
! Solar source
Sfactor = single_albedo*RTV%Solar_irradiance/PI
IF( RTV%mth_Azi == 0 ) Sfactor = Sfactor/TWO
EXPfactor = exp(-optical_depth/RTV%COS_SUN)
s_transmittance = exp(-total_opt/RTV%COS_SUN)
optical_depth_sun = optical_depth / RTV%COS_SUN
CALL Safe_Exp_Neg(optical_depth_sun, EXPfactor, is_clamped_exp)
total_opt_sun = total_opt / RTV%COS_SUN
CALL Safe_Exp_Neg(total_opt_sun, s_transmittance, is_clamped_total)

DO i = 1, nZ
Solar(i) = -bb(i,nZ+1)*Sfactor
Expand Down Expand Up @@ -2126,8 +2200,12 @@ SUBROUTINE CRTM_AMOM_layer_AD( n_streams, & ! Input, number of streams
bb_AD(i,nZ+1)=bb_AD(i,nZ+1) - Solar_AD(i)*Sfactor
ENDDO

total_opt_AD = total_opt_AD -s_transmittance_AD/RTV%COS_SUN*s_transmittance
optical_depth_AD = optical_depth_AD -EXPfactor_AD/RTV%COS_SUN*EXPfactor
IF ( .NOT. is_clamped_total ) THEN
total_opt_AD = total_opt_AD - s_transmittance_AD/RTV%COS_SUN*s_transmittance
END IF
IF ( .NOT. is_clamped_exp ) THEN
optical_depth_AD = optical_depth_AD -EXPfactor_AD/RTV%COS_SUN*EXPfactor
END IF

IF( RTV%mth_Azi == 0 ) THEN
Sfactor_AD = Sfactor_AD/TWO
Expand Down
Loading