From be26002d160536179196b8042aff5ee864535be9 Mon Sep 17 00:00:00 2001 From: Cheng Dang Date: Wed, 2 Jul 2025 23:46:28 -0600 Subject: [PATCH 1/4] Add backscatter profile --- src/AtmOptics/CRTM_AtmOptics_Define.f90 | 5 +- .../AerosolScatter/ASvar_Define.f90 | 2 + src/AtmScatter/CRTM_AOD_Module.f90 | 11 ++- src/AtmScatter/CRTM_AerosolScatter.f90 | 86 ++++++++++++++++--- src/AtmScatter/CRTM_AtmScatter_Define.f90 | 14 ++- .../AerosolCoeff/AerosolCoeff_netCDF_IO.f90 | 21 +++++ src/RTSolution/CRTM_RTSolution_Define.f90 | 80 +++++++++++++++-- .../regression/forward/test_AOD/test_AOD.f90 | 3 + .../regression/k_matrix/test_AOD/test_AOD.f90 | 3 + 9 files changed, 199 insertions(+), 26 deletions(-) diff --git a/src/AtmOptics/CRTM_AtmOptics_Define.f90 b/src/AtmOptics/CRTM_AtmOptics_Define.f90 index b1f3cec1..8ba21616 100644 --- a/src/AtmOptics/CRTM_AtmOptics_Define.f90 +++ b/src/AtmOptics/CRTM_AtmOptics_Define.f90 @@ -289,6 +289,9 @@ ELEMENTAL SUBROUTINE CRTM_AtmOptics_Create( & CALL AtmOptics_Allocate( self, alloc_stat ) IF ( alloc_stat /= 0 ) RETURN END IF + + ! Set the default value + CALL CRTM_AtmOptics_Zero( self ) ! Initialise dimensions (but not arrays!) self%n_Layers = n_Layers @@ -400,7 +403,7 @@ SUBROUTINE Scalar_Inspect(self) WRITE(*,'(5(1x,es22.15,:))') self%Single_Scatter_Albedo(1:self%n_Layers) WRITE(*,'(3x,"Asymmetry_Factor :")') WRITE(*,'(5(1x,es22.15,:))') self%Asymmetry_Factor(1:self%n_Layers) - WRITE(*,'(3x,"Backscat_Coefficient Backscattering coefficient :")') + WRITE(*,'(3x,"Backscat_Coefficient :")') WRITE(*,'(5(1x,es22.15,:))') self%Backscat_Coefficient(1:self%n_Layers) WRITE(*,'(3x,"Delta_Truncation :")') WRITE(*,'(5(1x,es22.15,:))') self%Delta_Truncation(1:self%n_Layers) diff --git a/src/AtmScatter/AerosolScatter/ASvar_Define.f90 b/src/AtmScatter/AerosolScatter/ASvar_Define.f90 index 586f2690..cd8c2901 100644 --- a/src/AtmScatter/AerosolScatter/ASvar_Define.f90 +++ b/src/AtmScatter/AerosolScatter/ASvar_Define.f90 @@ -131,6 +131,7 @@ MODULE ASvar_Define REAL(fp), ALLOCATABLE :: ke(:,:) ! I3 x I4 Mass extinction coefficient REAL(fp), ALLOCATABLE :: w(:,:) ! I3 x I4 Single Scatter Albedo REAL(fp), ALLOCATABLE :: g(:,:) ! I3 x I4 Asymmetry factor + REAL(fp), ALLOCATABLE :: kb(:,:) ! I3 x I4 Mass backscatter coefficient REAL(fp), ALLOCATABLE :: pcoeff(:,:,:,:) ! 0:I1 x I2 x I3 x I4 Phase coefficients ! The accumulated scattering coefficient REAL(fp), ALLOCATABLE :: total_bs(:) ! I3 Volume scattering coefficient @@ -191,6 +192,7 @@ ELEMENTAL SUBROUTINE ASvar_Create( & self%ke(n_Layers, n_Aerosols), & self%w(n_Layers, n_Aerosols), & self%g(n_Layers, n_Aerosols), & + self%kb(n_Layers, n_Aerosols), & self%pcoeff(0:n_Legendre_Terms,n_Phase_Elements,n_Layers, n_Aerosols), & self%total_bs(n_Layers), & STAT = alloc_stat ) diff --git a/src/AtmScatter/CRTM_AOD_Module.f90 b/src/AtmScatter/CRTM_AOD_Module.f90 index 084e917a..77e70e87 100644 --- a/src/AtmScatter/CRTM_AOD_Module.f90 +++ b/src/AtmScatter/CRTM_AOD_Module.f90 @@ -19,7 +19,7 @@ MODULE CRTM_AOD_Module ! Module usage ! ------------ USE Message_Handler, ONLY: SUCCESS, FAILURE, Display_Message - USE CRTM_Parameters, ONLY: ZERO, & + USE CRTM_Parameters, ONLY: ZERO, TWO, PI, & MAX_N_PHASE_ELEMENTS, & MAX_N_LEGENDRE_TERMS USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type, & @@ -344,6 +344,10 @@ FUNCTION CRTM_AOD( & ! Save the nadir optical depth RTSolution(ln,m)%Layer_Optical_Depth(1:Atmosphere(m)%n_Layers) = AtmOptics%Optical_Depth + ! Save the Backscat_Coefficient, only non-zero for GOCART-GEOS5 + ! CD: double check if should be divided by 2PI + RTSolution(ln,m)%Backscat_Coefficient(1:Atmosphere(m)%n_Layers) = AtmOptics%Backscat_Coefficient/(TWO*PI) + END DO Channel_Loop END DO Sensor_Loop @@ -669,6 +673,11 @@ FUNCTION CRTM_AOD_TL( & RTSolution(ln,m)%Layer_Optical_Depth(1:Atmosphere(m)%n_Layers) = AtmOptics%Optical_Depth RTSolution_TL(ln,m)%Layer_Optical_Depth(1:Atmosphere(m)%n_Layers) = AtmOptics_TL%Optical_Depth + ! Save the backscatter coefficient, only non-zero for GOCART-GEOS5 + ! CD: double check if should be divided by 2PI + RTSolution(ln,m)%Backscat_Coefficient(1:Atmosphere(m)%n_Layers) = AtmOptics%Backscat_Coefficient/(TWO*PI) + RTSolution_TL(ln,m)%Backscat_Coefficient(1:Atmosphere(m)%n_Layers) = AtmOptics_TL%Backscat_Coefficient/(TWO*PI) + END DO Channel_Loop END DO Sensor_Loop diff --git a/src/AtmScatter/CRTM_AerosolScatter.f90 b/src/AtmScatter/CRTM_AerosolScatter.f90 index 0c46f915..bfad500b 100644 --- a/src/AtmScatter/CRTM_AerosolScatter.f90 +++ b/src/AtmScatter/CRTM_AerosolScatter.f90 @@ -251,10 +251,12 @@ FUNCTION CRTM_Compute_AerosolScatter( & Atm%Relative_Humidity(ka) , & ! Input ASV%ke(ka,n) , & ! Output ASV%w(ka,n) , & ! Output + ASV%kb(ka,n) , & ! Output ASV%pcoeff(:,:,ka,n) , & ! Output ASV%asi(ka,n) ) ! Interpolation ! interpolation quality control + ! CD: kb quality control? IF( ASV%ke(ka,n) <= ZERO ) THEN ASV%ke(ka,n) = ZERO ASV%w(ka,n) = ZERO @@ -302,6 +304,9 @@ FUNCTION CRTM_Compute_AerosolScatter( & bs = Atm%Aerosol(n)%Concentration(ka) * ASV%ke(ka,n) * ASV%w(ka,n) ASV%Total_bs(ka) = ASV%Total_bs(ka) + bs AScat%Single_Scatter_Albedo(ka) = AScat%Single_Scatter_Albedo(ka) + bs + ! Back scatter coefficient, only non zero for GOCART-GEOS5 + AScat%Backscat_Coefficient(ka) = AScat%Backscat_Coefficient(ka) + & + (ASV%kb(ka,n) * Atm%Aerosol(n)%Concentration(ka)) DO m = 1, AScat%n_Phase_Elements DO l = 0, AScat%n_Legendre_Terms @@ -432,7 +437,7 @@ FUNCTION CRTM_Compute_AerosolScatter_TL( & LOGICAL :: Layer_Mask(Atm%n_Layers) INTEGER :: Layer_Index(Atm%n_Layers) INTEGER :: nAerosol_Layers - REAL(fp) :: ke_TL, w_TL + REAL(fp) :: ke_TL, w_TL, kb_TL REAL(fp) :: pcoeff_TL(0:AScat%n_Legendre_Terms, AScat%n_Phase_Elements) REAL(fp) :: bs, bs_TL @@ -473,11 +478,13 @@ FUNCTION CRTM_Compute_AerosolScatter_TL( & Atm%Aerosol(n)%Type , & ! Input ASV%ke(ka,n) , & ! Input ASV%w(ka,n) , & ! Input + ASV%kb(ka,n) , & ! Input Atm_TL%Aerosol(n)%Effective_Radius(ka) , & ! TL Input Atm_TL%Aerosol(n)%Effective_Variance(ka), & ! TL Input Atm_TL%Relative_Humidity(ka) , & ! TL Input ke_TL , & ! TL Output w_TL , & ! TL Output + kb_TL , & ! TL Output pcoeff_TL , & ! TL Output ASV%asi(ka,n) ) ! Interpolation @@ -498,6 +505,10 @@ FUNCTION CRTM_Compute_AerosolScatter_TL( & AScat_TL%Optical_Depth(ka) = AScat_TL%Optical_Depth(ka) + & (ke_TL * Atm%Aerosol(n)%Concentration(ka)) + & (ASV%ke(ka,n) * Atm_TL%Aerosol(n)%Concentration(ka)) + ! Compute the back scatter coefficient + AScat_TL%Backscat_Coefficient(ka) = AScat_TL%Backscat_Coefficient(ka) + & + (kb_TL * Atm%Aerosol(n)%Concentration(ka)) + & + (ASV%kb(ka,n) * Atm_TL%Aerosol(n)%Concentration(ka)) ! Compute the phase matrix coefficients IF( n_Phase_Elements > 0 .and. AScat%Include_Scattering ) THEN @@ -639,7 +650,7 @@ FUNCTION CRTM_Compute_AerosolScatter_AD( & LOGICAL :: Layer_Mask(Atm%n_Layers) INTEGER :: Layer_Index(Atm%n_Layers) INTEGER :: nAerosol_Layers - REAL(fp) :: ke_AD, w_AD + REAL(fp) :: ke_AD, w_AD, kb_AD REAL(fp) :: pcoeff_AD(0:AScat%n_Legendre_Terms, AScat%n_Phase_Elements) REAL(fp) :: bs, bs_AD @@ -683,10 +694,11 @@ FUNCTION CRTM_Compute_AerosolScatter_AD( & ! Initialize the individual ! Aerosol adjoint variables - bs_AD = ZERO + bs_AD = ZERO pcoeff_AD = ZERO ke_AD = ZERO w_AD = ZERO + kb_AD = ZERO ! Compute the adjoint of the ! phase matrix coefficients @@ -715,11 +727,13 @@ FUNCTION CRTM_Compute_AerosolScatter_AD( & ! Compute the adjoint of the volume ! scattering coefficient. - ke_AD = ke_AD + (Atm%Aerosol(n)%Concentration(ka) * bs_AD * ASV%w(ka,n) ) Atm_AD%Aerosol(n)%Concentration(ka) = Atm_AD%Aerosol(n)%Concentration(ka) + & ( bs_AD * ASV%ke(ka,n) * ASV%w(ka,n) ) + ! ! Compute the adjoint of the backscattering coefficient + ! CD: placeholder for now + ! interpolation quality control IF( ASV%w(ka,n) >= ONE ) THEN w_AD = ZERO @@ -737,8 +751,10 @@ FUNCTION CRTM_Compute_AerosolScatter_AD( & Atm%Aerosol(n)%Type , & ! Input ASV%ke(ka,n) , & ! input ASV%w(ka,n) , & ! input + ASV%kb(ka,n) , & ! input ke_AD , & ! AD Input w_AD , & ! AD Input + kb_AD , & ! AD Input pcoeff_AD , & ! AD Input Atm_AD%Aerosol(n)%Effective_Radius(ka) , & ! AD Output Atm_AD%Aerosol(n)%Effective_Variance(ka), & ! AD Output @@ -777,6 +793,7 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure RH , & ! Input relative humidity, unitless ke , & ! Output extinction coefficient (=~ optical depth) w , & ! Output single scattering albedo + kb , & ! Output backscattering coefficient pcoeff , & ! Output spherical Legendre coefficients asi ) ! Output interpolation data ! Arguments @@ -787,6 +804,7 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure REAL(fp) , INTENT(IN) :: RH REAL(fp) , INTENT(OUT) :: ke REAL(fp) , INTENT(OUT) :: w + REAL(fp) , INTENT(OUT) :: kb REAL(fp) , INTENT(IN OUT) :: pcoeff(0:,:) TYPE(ASinterp_type) , INTENT(IN OUT) :: asi ! Local variables @@ -801,6 +819,7 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure IF ( Aerosol_Type == AerosolCoeff_BYPASS_AEROSOL ) THEN ke = ZERO w = ZERO + kb = ZERO pcoeff = ZERO RETURN END IF @@ -823,6 +842,8 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure ! Fixed indices fix_rh = 1 fix_sig = 1 + ! kb is not used in CRTM scheme + kb = ZERO ! Find effective radius indices for Interpolation asi%r_int = MAX(MIN(AeroC%Reff(AeroC%n_Radii,k),Reff),AeroC%Reff(1,k)) @@ -851,6 +872,8 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure ELSE IF ( AeroC%Scheme == 'CMAQ' ) THEN ! Fixed indices fix_rh = 1 + ! kb is not used in CMAQ scheme + kb = ZERO ! Find effective radius indices for Interpolation asi%r_int = MAX(MIN(AeroC%Reff(AeroC%n_Radii,k),Reff),AeroC%Reff(1,k)) @@ -884,7 +907,7 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure ! Absorption coefficient ke = ke * (ONE- w) END IF - + ELSE IF ( AeroC%Scheme == 'GOCART-GEOS5' .OR. AeroC%Scheme == 'NAAPS' ) THEN ! Fixed indices fix_sig = 1 @@ -901,6 +924,13 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure ! Perform Interpolation CALL interp_2D( AeroC%ke( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k), asi%wlp, asi%hlp, ke ) CALL interp_2D( AeroC%w( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k), asi%wlp, asi%hlp, w ) + ! CD: This is for GOCART-GEOS5 only, for NAAPS, kb is zero + IF (AeroC%Scheme == 'GOCART-GEOS5') THEN + CALL interp_2D( AeroC%kb( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k), asi%wlp, asi%hlp, kb ) + ELSE + kb = ZERO + END IF + CALL interp_2D( AeroC%kb( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k), asi%wlp, asi%hlp, kb ) IF (AerosolScatter%n_Phase_Elements > 0 .and. AerosolScatter%Include_Scattering ) THEN pcoeff(0,1) = POINT_5 DO m = 1, AerosolScatter%n_Phase_Elements @@ -914,7 +944,6 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure ke = ke * (ONE- w) END IF - END IF !IF ( AeroC%Scheme == 'CRTM') THEN END SUBROUTINE Get_Aerosol_Opt @@ -933,20 +962,22 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st Aerosol_Type , & ! Input see CRTM_Aerosol_Define.f90 ke , & ! Input w , & ! Input + kb , & ! Input Reff_TL , & ! Input TL effective radius (mm) Rsig_TL , & ! Input TL effective radius (mm) RH_TL , & ! Input TL effective radius (mm) ke_TL , & ! Output TL extinction coefficient (=~ optical depth) w_TL , & ! Output TL single scattering albedo + kb_TL , & ! Output TL backscattering coefficient pcoeff_TL , & ! TL Output spherical Legendre coefficients asi ) ! Input interpolation data ! Arguments TYPE(CRTM_AtmOptics_type), INTENT(IN) :: AerosolScatter_TL INTEGER , INTENT(IN) :: Aerosol_Type - REAL(fp), INTENT(IN) :: ke, w + REAL(fp), INTENT(IN) :: ke, w, kb REAL(fp), INTENT(IN) :: Reff_TL, Rsig_TL, RH_TL - REAL(fp), INTENT(OUT) :: ke_TL, w_TL + REAL(fp), INTENT(OUT) :: ke_TL, w_TL, kb_TL REAL(fp), INTENT(IN OUT) :: pcoeff_TL(0:,:) TYPE(ASinterp_type), INTENT(IN) :: asi ! Local variables @@ -969,6 +1000,7 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st IF ( Aerosol_Type == AerosolCoeff_BYPASS_AEROSOL ) THEN ke_TL = ZERO w_TL = ZERO + kb_TL = ZERO pcoeff_TL = ZERO RETURN END IF @@ -1017,6 +1049,8 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st ! Fixed indices fix_rh = 1 fix_sig = 1 + ! kb is not used in CRTM scheme + kb_TL = ZERO ! Effective radius term CALL LPoly_TL( asi%r, asi%r_int, & ! FWD Input @@ -1058,6 +1092,8 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st ELSE IF ( AeroC%Scheme == 'CMAQ' ) THEN ! Fixed indices fix_rh = 1 + ! kb is not used in CMAQ scheme + kb_TL = ZERO ! Effective radius term CALL LPoly_TL( asi%r, asi%r_int, & ! FWD Input @@ -1123,6 +1159,15 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st CALL interp_2D_TL( zg , asi%wlp, asi%hlp, & ! FWD Input zg_TL, wlp_TL , hlp_TL , & ! TL Input w_TL ) ! TL Output + ! Backscattering coefficient + IF (AeroC%Scheme == 'GOCART-GEOS5') THEN + zg => AeroC%kb( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k) + CALL interp_2D_TL( zg , asi%wlp, asi%hlp, & ! FWD Input + zg_TL, wlp_TL , hlp_TL , & ! TL Input + kb_TL ) ! TL Output + ELSE + kb_TL = ZERO + END IF ! Phase matrix coefficients IF (AerosolScatter_TL%n_Phase_Elements > 0 .and. AerosolScatter_TL%Include_Scattering ) THEN pcoeff_TL(0,1) = ZERO @@ -1161,8 +1206,10 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st Aerosol_Type , & ! Input see CRTM_Aerosol_Define.f90 ke , & ! Input w , & ! Input + kb , & ! Input ke_AD , & ! AD Input extinction cross section w_AD , & ! AD Input single scatter albedo + kb_AD , & ! AD Input backscattering coefficient pcoeff_AD , & ! AD Input spherical Legendre coefficients Reff_AD , & ! AD Outputmode radius Rsig_AD , & ! AD Output radius deviation @@ -1171,9 +1218,10 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st ! Arguments TYPE(CRTM_AtmOptics_type), INTENT(IN) :: AerosolScatter_AD INTEGER , INTENT(IN) :: Aerosol_Type - REAL(fp), INTENT(IN) :: ke, w + REAL(fp), INTENT(IN) :: ke, w, kb REAL(fp), INTENT(IN OUT) :: ke_AD ! AD Input REAL(fp), INTENT(IN OUT) :: w_AD ! AD Input + REAL(fp), INTENT(IN OUT) :: kb_AD ! AD Input REAL(fp), INTENT(IN OUT) :: pcoeff_AD(0:,:) ! AD Input REAL(fp), INTENT(IN OUT) :: Reff_AD ! AD Output REAL(fp), INTENT(IN OUT) :: Rsig_AD ! AD Output @@ -1200,6 +1248,7 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st Reff_AD = ZERO ke_AD = ZERO w_AD = ZERO + kb_AD = ZERO pcoeff_AD = ZERO Rsig_AD = ZERO RH_AD = ZERO @@ -1214,6 +1263,7 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st Reff_AD = ZERO ke_AD = ZERO w_AD = ZERO + kb_AD = ZERO pcoeff_AD = ZERO Rsig_AD = ZERO RH_AD = ZERO @@ -1246,6 +1296,8 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st ! Fixed indices fix_rh = 1 fix_sig = 1 + ! kb is not used in CRTM scheme + kb_AD = ZERO ! Phase matrix coefficients IF (AerosolScatter_AD%n_Phase_Elements > 0 .and. AerosolScatter_AD%Include_Scattering ) THEN @@ -1300,6 +1352,8 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st ELSE IF ( AeroC%Scheme == 'CMAQ' ) THEN ! Fixed indices fix_rh = 1 + ! kb is not used in CMAQ scheme + kb_AD = ZERO ! Phase matrix coefficients IF (AerosolScatter_AD%n_Phase_Elements > 0 .and. AerosolScatter_AD%Include_Scattering ) THEN @@ -1333,6 +1387,7 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st CALL interp_3D_AD( zc , asi%wlp, asi%xlp, asi%vlp, & ! FWD Input ke_AD , & ! AD Input zc_AD, wlp_AD, xlp_AD, vlp_AD ) ! AD Output + NULLIFY(zc) ! Compute the AD of the interpolating polynomials @@ -1385,16 +1440,25 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st END IF ! Single scatter albedo - zg => AeroC%w( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k) + zg => AeroC%w( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k ) CALL interp_2D_AD( zg , asi%wlp, asi%hlp, & ! FWD Input w_AD , & ! AD Input zg_AD, wlp_AD , hlp_AD ) ! AD Output ! Extinction coefficient - zg => AeroC%ke( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k) + zg => AeroC%ke( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k ) CALL interp_2D_AD( zg , asi%wlp, asi%hlp, & ! FWD Input ke_AD , & ! AD Input zg_AD, wlp_AD, hlp_AD ) ! AD Output + ! Backscattering coefficient + IF (AeroC%Scheme == 'GOCART-GEOS5') THEN + zg => AeroC%kb( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k) + CALL interp_2D_AD( zg , asi%wlp, asi%hlp, & ! FWD Input + kb_AD , & ! AD Input + zg_AD, wlp_AD, hlp_AD ) ! AD Output + ELSE + kb_AD = ZERO + END IF NULLIFY(zg) ! Compute the AD of the interpolating polynomials diff --git a/src/AtmScatter/CRTM_AtmScatter_Define.f90 b/src/AtmScatter/CRTM_AtmScatter_Define.f90 index 42bc080e..11c002ef 100644 --- a/src/AtmScatter/CRTM_AtmScatter_Define.f90 +++ b/src/AtmScatter/CRTM_AtmScatter_Define.f90 @@ -43,8 +43,6 @@ MODULE CRTM_AtmScatter_Define ! ----------------- ! Module parameters ! ----------------- - ! RCS Id for the module - CHARACTER(*), PARAMETER :: MODULE_RCS_ID = & ! Message string length INTEGER, PARAMETER :: ML = 256 @@ -66,6 +64,7 @@ MODULE CRTM_AtmScatter_Define REAL(fp), POINTER :: Optical_Depth(:) => NULL() ! K REAL(fp), POINTER :: Single_Scatter_Albedo(:) => NULL() ! K REAL(fp), POINTER :: Asymmetry_Factor(:) => NULL() ! K + REAL(fp), POINTER :: Backscat_Coefficient(:) => NULL() ! K REAL(fp), POINTER :: Delta_Truncation(:) => NULL() ! K REAL(fp), POINTER :: Phase_Coefficient(:,:,:) => NULL() ! 0:Ic x Ip x K END TYPE CRTM_AtmScatter_type @@ -162,6 +161,7 @@ FUNCTION CRTM_Associated_AtmScatter( AtmScatter, & ! Input IF ( ALL_Test ) THEN IF ( ASSOCIATED(AtmScatter%Optical_Depth ) .AND. & ASSOCIATED(AtmScatter%Single_Scatter_Albedo) .AND. & + ASSOCIATED(AtmScatter%Backscat_Coefficient ) .AND. & ASSOCIATED(AtmScatter%Asymmetry_Factor ) .AND. & ASSOCIATED(AtmScatter%Delta_Truncation ) .AND. & ASSOCIATED(AtmScatter%Phase_Coefficient ) ) THEN @@ -171,6 +171,7 @@ FUNCTION CRTM_Associated_AtmScatter( AtmScatter, & ! Input IF ( ASSOCIATED(AtmScatter%Optical_Depth ) .OR. & ASSOCIATED(AtmScatter%Single_Scatter_Albedo) .OR. & ASSOCIATED(AtmScatter%Asymmetry_Factor ) .OR. & + ASSOCIATED(AtmScatter%Backscat_Coefficient ) .AND. & ASSOCIATED(AtmScatter%Delta_Truncation ) .OR. & ASSOCIATED(AtmScatter%Phase_Coefficient ) ) THEN Association_Status = .TRUE. @@ -276,6 +277,7 @@ FUNCTION CRTM_Destroy_AtmScatter( AtmScatter, & ! Output DEALLOCATE( AtmScatter%Optical_Depth , & AtmScatter%Single_Scatter_Albedo, & AtmScatter%Asymmetry_Factor , & + AtmScatter%Backscat_Coefficient , & AtmScatter%Delta_Truncation , & AtmScatter%Phase_Coefficient , & STAT=Allocate_Status ) @@ -439,6 +441,7 @@ FUNCTION CRTM_Allocate_AtmScatter( n_Layers , & ! Input ALLOCATE( AtmScatter%Optical_Depth( n_Layers ), & AtmScatter%Single_Scatter_Albedo( n_Layers ), & AtmScatter%Asymmetry_Factor( n_Layers ), & + AtmScatter%Backscat_Coefficient( n_Layers ), & AtmScatter%Delta_Truncation( n_Layers ), & AtmScatter%Phase_Coefficient( 0:n_Legendre_Terms, & n_Phase_Elements, & @@ -467,6 +470,7 @@ FUNCTION CRTM_Allocate_AtmScatter( n_Layers , & ! Input AtmScatter%Optical_Depth = ZERO AtmScatter%Single_Scatter_Albedo = ZERO AtmScatter%Asymmetry_Factor = ZERO + AtmScatter%Backscat_Coefficient = ZERO AtmScatter%Delta_Truncation = ZERO AtmScatter%Phase_Coefficient = ZERO @@ -598,7 +602,8 @@ FUNCTION CRTM_Assign_AtmScatter( AtmScatter_in , & ! Input ! --------------- AtmScatter_out%Optical_Depth = AtmScatter_in%Optical_Depth AtmScatter_out%Single_Scatter_Albedo = AtmScatter_in%Single_Scatter_Albedo - AtmScatter_out%Asymmetry_Factor = AtmScatter_in%Asymmetry_Factor + AtmScatter_out%Asymmetry_Factor = AtmScatter_in%Asymmetry_Factor + AtmScatter_out%Backscat_Coefficient = AtmScatter_in%Backscat_Coefficient AtmScatter_out%Delta_Truncation = AtmScatter_in%Delta_Truncation AtmScatter_out%Phase_Coefficient = AtmScatter_in%Phase_Coefficient @@ -644,7 +649,8 @@ SUBROUTINE CRTM_Zero_AtmScatter( AtmScatter ) ! Reset the array components AtmScatter%Optical_Depth = ZERO AtmScatter%Single_Scatter_Albedo = ZERO - AtmScatter%Asymmetry_Factor = ZERO + AtmScatter%Asymmetry_Factor = ZERO + AtmScatter%Backscat_Coefficient = ZERO AtmScatter%Delta_Truncation = ZERO AtmScatter%Phase_Coefficient = ZERO END SUBROUTINE CRTM_Zero_AtmScatter diff --git a/src/Coefficients/AerosolCoeff/AerosolCoeff_netCDF_IO.f90 b/src/Coefficients/AerosolCoeff/AerosolCoeff_netCDF_IO.f90 index 69e59030..1bcb7460 100644 --- a/src/Coefficients/AerosolCoeff/AerosolCoeff_netCDF_IO.f90 +++ b/src/Coefficients/AerosolCoeff/AerosolCoeff_netCDF_IO.f90 @@ -90,6 +90,7 @@ MODULE AerosolCoeff_netCDF_IO CHARACTER(*), PARAMETER :: RSIG_VARNAME = 'Rsig' CHARACTER(*), PARAMETER :: RH_VARNAME = 'RH' CHARACTER(*), PARAMETER :: KE_VARNAME = 'ke' + CHARACTER(*), PARAMETER :: KB_VARNAME = 'kb' CHARACTER(*), PARAMETER :: W_VARNAME = 'w' CHARACTER(*), PARAMETER :: G_VARNAME = 'g' CHARACTER(*), PARAMETER :: PCOEFF_VARNAME = 'pcoeff' @@ -104,6 +105,7 @@ MODULE AerosolCoeff_netCDF_IO CHARACTER(*), PARAMETER :: RSIG_LONGNAME = 'Mode radius deviation' CHARACTER(*), PARAMETER :: RH_LONGNAME = 'Relative humidity' CHARACTER(*), PARAMETER :: KE_LONGNAME = 'ke' + CHARACTER(*), PARAMETER :: KB_LONGNAME = 'kb' CHARACTER(*), PARAMETER :: W_LONGNAME = 'w' CHARACTER(*), PARAMETER :: G_LONGNAME = 'g' CHARACTER(*), PARAMETER :: PCOEFF_LONGNAME = 'pcoeff' @@ -117,6 +119,7 @@ MODULE AerosolCoeff_netCDF_IO CHARACTER(*), PARAMETER :: RSIG_DESCRIPTION = 'Mode radius standard deviation LUT dimension vector' CHARACTER(*), PARAMETER :: RH_DESCRIPTION = 'Relative humidity LUT dimension vector' CHARACTER(*), PARAMETER :: KE_DESCRIPTION = 'Mass extinction coefficient for aerosol scatterers' + CHARACTER(*), PARAMETER :: KB_DESCRIPTION = 'Mass backscatter coefficient for aerosol scatterers' CHARACTER(*), PARAMETER :: W_DESCRIPTION = 'Single scatter albedo for aerosol scatterers' CHARACTER(*), PARAMETER :: G_DESCRIPTION = 'Asymmetry parameter for aerosol scatterers' CHARACTER(*), PARAMETER :: PCOEFF_DESCRIPTION = 'Phase coefficients for aerosol scatterers' @@ -130,6 +133,7 @@ MODULE AerosolCoeff_netCDF_IO CHARACTER(*), PARAMETER :: RSIG_UNITS = 'N/A' CHARACTER(*), PARAMETER :: RH_UNITS = 'fraction (0->1)' CHARACTER(*), PARAMETER :: KE_UNITS = 'Metres squared per kilogram (m^2.kg^-1)' + CHARACTER(*), PARAMETER :: KB_UNITS = 'Metres squared per kilogram (m^2.kg^-1)' CHARACTER(*), PARAMETER :: W_UNITS = 'N/A' CHARACTER(*), PARAMETER :: G_UNITS = 'N/A' CHARACTER(*), PARAMETER :: PCOEFF_UNITS = 'N/A' @@ -144,6 +148,7 @@ MODULE AerosolCoeff_netCDF_IO REAL(Double) , PARAMETER :: RSIG_FILLVALUE = FILL_FLOAT REAL(Double) , PARAMETER :: RH_FILLVALUE = FILL_FLOAT REAL(Double) , PARAMETER :: KE_FILLVALUE = FILL_FLOAT + REAL(Double) , PARAMETER :: KB_FILLVALUE = FILL_FLOAT REAL(Double) , PARAMETER :: W_FILLVALUE = FILL_FLOAT REAL(Double) , PARAMETER :: G_FILLVALUE = FILL_FLOAT REAL(Double) , PARAMETER :: PCOEFF_FILLVALUE = FILL_FLOAT @@ -157,6 +162,7 @@ MODULE AerosolCoeff_netCDF_IO INTEGER, PARAMETER :: RSIG_TYPE = NF90_DOUBLE INTEGER, PARAMETER :: RH_TYPE = NF90_DOUBLE INTEGER, PARAMETER :: KE_TYPE = NF90_DOUBLE + INTEGER, PARAMETER :: KB_TYPE = NF90_DOUBLE INTEGER, PARAMETER :: W_TYPE = NF90_DOUBLE INTEGER, PARAMETER :: G_TYPE = NF90_DOUBLE INTEGER, PARAMETER :: PCOEFF_TYPE = NF90_DOUBLE @@ -1143,6 +1149,21 @@ FUNCTION AerosolCoeff_netCDF_ReadFile( & ' - '//TRIM(NF90_STRERROR( NF90_Status )) CALL Read_Cleanup(); RETURN END IF + ! ...kb variable, only for GOCART-GEOS5 scheme + IF ( Aerosol_Model == 'GOCART-GEOS5' ) THEN + NF90_Status = NF90_INQ_VARID( FileId,KB_VARNAME,VarId ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error inquiring '//TRIM(Filename)//' for '//KB_VARNAME//& + ' variable ID - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Read_Cleanup(); RETURN + END IF + NF90_Status = NF90_GET_VAR( FileId,VarID,AerosolCoeff%kb ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error reading '//KB_VARNAME//' from '//TRIM(Filename)//& + ' - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Read_Cleanup(); RETURN + END IF + END IF ! ...ke variable NF90_Status = NF90_INQ_VARID( FileId,KE_VARNAME,VarId ) IF ( NF90_Status /= NF90_NOERR ) THEN diff --git a/src/RTSolution/CRTM_RTSolution_Define.f90 b/src/RTSolution/CRTM_RTSolution_Define.f90 index 7a205a2d..d6c72944 100644 --- a/src/RTSolution/CRTM_RTSolution_Define.f90 +++ b/src/RTSolution/CRTM_RTSolution_Define.f90 @@ -155,6 +155,7 @@ MODULE CRTM_RTSolution_Define CHARACTER(*), PARAMETER :: SSA_VARNAME = 'Single_Scatter_Albedo' CHARACTER(*), PARAMETER :: ACREFL_VARNAME = 'Reflectivity' ! Active sensor CHARACTER(*), PARAMETER :: ACRATT_VARNAME = 'Reflectivity_Attenuated' ! Active sensor + CHARACTER(*), PARAMETER :: BS_VARNAME = 'Backscat_Coefficient' ! Variable units attribute. CHARACTER(*), PARAMETER :: UNITS_ATTNAME = 'units' @@ -181,7 +182,9 @@ MODULE CRTM_RTSolution_Define ! ...Brightness_Temperature, Tb_clear CHARACTER(*), PARAMETER :: BT_UNITS = 'Kelvin' ! ...Visible (or UV) reflectance - CHARACTER(*), PARAMETER :: RF_UNITS = 'fraction (0->1)' + CHARACTER(*), PARAMETER :: RF_UNITS = 'fraction (0->1)' + ! ...Backscat_Coefficient + CHARACTER(*), PARAMETER :: BS_UNITS = 'per Rad (m^-1 rad^-1)' ! Variable _FillValue attribute. CHARACTER(*), PARAMETER :: FILLVALUE_ATTNAME = '_FillValue' @@ -234,6 +237,7 @@ MODULE CRTM_RTSolution_Define REAL(fp), ALLOCATABLE :: Upwelling_Radiance(:) ! K REAL(fp), ALLOCATABLE :: Layer_Optical_Depth(:) ! K REAL(fp), ALLOCATABLE :: Single_Scatter_Albedo(:) ! K + REAL(fp), ALLOCATABLE :: Backscat_Coefficient(:) ! K ! Radiative transfer results for a single channel REAL(fp) :: Radiance = ZERO REAL(fp) :: Brightness_Temperature = ZERO @@ -373,6 +377,7 @@ ELEMENTAL SUBROUTINE CRTM_RTSolution_Create( RTSolution, n_Layers ) RTSolution%Single_Scatter_Albedo(n_Layers), & RTSolution%Reflectivity(n_Layers), & RTSolution%Reflectivity_Attenuated(n_Layers), & + RTSolution%Backscat_Coefficient(n_Layers), & STAT = alloc_stat ) IF ( alloc_stat /= 0 ) RETURN @@ -386,6 +391,7 @@ ELEMENTAL SUBROUTINE CRTM_RTSolution_Create( RTSolution, n_Layers ) RTSolution%Single_Scatter_Albedo = ZERO RTSolution%Reflectivity = ZERO RTSolution%Reflectivity_Attenuated = ZERO + RTSolution%Backscat_Coefficient = ZERO ! Set allocation indicator RTSolution%Is_Allocated = .TRUE. @@ -451,6 +457,7 @@ ELEMENTAL SUBROUTINE CRTM_RTSolution_Zero( RTSolution ) RTSolution%Single_Scatter_Albedo = ZERO RTSolution%Reflectivity = ZERO RTSolution%Reflectivity_Attenuated = ZERO + RTSolution%Backscat_Coefficient = ZERO END IF END SUBROUTINE CRTM_RTSolution_Zero @@ -540,6 +547,8 @@ SUBROUTINE Scalar_Inspect( RTSolution, Unit ) WRITE(fid,'(5(1x,es22.15,:))') RTSolution%Reflectivity WRITE(fid,'(3x,"Reflectivity_Attenuated :")') WRITE(fid,'(5(1x,es22.15,:))') RTSolution%Reflectivity_Attenuated + WRITE(fid,'(3x,"Backscat_Coefficient :")') + WRITE(fid,'(5(1x,es22.15,:))') RTSolution%Backscat_Coefficient END IF FLUSH(fid) END SUBROUTINE Scalar_Inspect @@ -690,7 +699,8 @@ ELEMENTAL FUNCTION CRTM_RTSolution_Compare( & (.NOT. ALL(Compares_Within_Tolerance(x%Upwelling_Radiance , y%Upwelling_Radiance , n))) .OR. & (.NOT. ALL(Compares_Within_Tolerance(x%Layer_Optical_Depth , y%Layer_Optical_Depth , n))) .OR. & (.NOT. ALL(Compares_Within_Tolerance(x%Reflectivity , y%Reflectivity , n))) .OR. & - (.NOT. ALL(Compares_Within_Tolerance(x%Reflectivity_Attenuated , y%Reflectivity_Attenuated , n)))) RETURN + (.NOT. ALL(Compares_Within_Tolerance(x%Reflectivity_Attenuated , y%Reflectivity_Attenuated , n))) .OR. & + (.NOT. ALL(Compares_Within_Tolerance(x%Backscat_Coefficient , y%Backscat_Coefficient , n)))) RETURN END IF ! If we get here, the structures are comparable @@ -1624,6 +1634,7 @@ FUNCTION CRTM_RTSolution_ReadFile_NetCDF( & REAL(fp), ALLOCATABLE :: Single_Scatter_Albedo(:,:,:) REAL(fp), ALLOCATABLE :: Reflectivity(:,:,:) REAL(fp), ALLOCATABLE :: Reflectivity_Attenuated(:,:,:) + REAL(fp), ALLOCATABLE :: Backscat_Coefficient(:,:,:) ! Set up @@ -1671,6 +1682,7 @@ FUNCTION CRTM_RTSolution_ReadFile_NetCDF( & Single_Scatter_Albedo( n_Channels, n_Layers, n_Profiles ), & Reflectivity( n_Channels, n_Layers, n_Profiles ), & Reflectivity_Attenuated( n_Channels, n_Layers, n_Profiles ), & + Backscat_Coefficient( n_Channels, n_Layers, n_Profiles ), & STAT = alloc_stat ) IF ( alloc_stat /= 0 ) THEN msg = 'Error allocating RTSolution output arrays' @@ -2011,6 +2023,19 @@ FUNCTION CRTM_RTSolution_ReadFile_NetCDF( & ' - '//TRIM(NF90_STRERROR( NF90_Status )) CALL Read_Cleanup(); RETURN END IF + ! ... Backscat_Coefficient variable + NF90_Status = NF90_INQ_VARID( FileId,BS_VARNAME,VarId ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error inquiring '//TRIM(Filename)//' for '//BS_VARNAME//& + ' variable ID - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Read_Cleanup(); RETURN + END IF + NF90_Status = NF90_GET_VAR( FileId,VarID,Backscat_Coefficient) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error reading '//BS_VARNAME//' from '//TRIM(Filename)//& + ' - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Read_Cleanup(); RETURN + END IF ! Close the file NF90_Status = NF90_CLOSE( FileId ); Close_File = .FALSE. @@ -2073,6 +2098,7 @@ FUNCTION CRTM_RTSolution_ReadFile_NetCDF( & RTSolution(l,m)%Single_Scatter_Albedo(c) = Single_Scatter_Albedo(l,c,m) RTSolution(l,m)%Reflectivity(c) = Reflectivity(l,c,m) RTSolution(l,m)%Reflectivity_Attenuated(c) = Reflectivity_Attenuated(l,c,m) + RTSolution(l,m)%Backscat_Coefficient(c) = Backscat_Coefficient(l,c,m) END DO END DO Channel_Loop END DO Profile_Loop @@ -2465,7 +2491,8 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & REAL(fp), ALLOCATABLE :: Single_Scatter_Albedo(:,:,:) REAL(fp), ALLOCATABLE :: Reflectivity(:,:,:) REAL(fp), ALLOCATABLE :: Reflectivity_Attenuated(:,:,:) - + REAL(fp), ALLOCATABLE :: Backscat_Coefficient(:,:,:) + ! Set up err_stat = SUCCESS @@ -2509,6 +2536,7 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & Single_Scatter_Albedo( n_Channels, n_Layers, n_Profiles ), & Reflectivity( n_Channels, n_Layers, n_Profiles ), & Reflectivity_Attenuated( n_Channels, n_Layers, n_Profiles ), & + Backscat_Coefficient( n_Channels, n_Layers, n_Profiles ), & STAT = alloc_stat ) IF ( alloc_stat /= 0 ) THEN msg = 'Error allocating RTSolution output arrays' @@ -2548,11 +2576,11 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & Single_Scatter_Albedo(l,c,m) = RTSolution(l,m)%Single_Scatter_Albedo(c) Reflectivity(l,c,m) = RTSolution(l,m)%Reflectivity(c) Reflectivity_Attenuated(l,c,m) = RTSolution(l,m)%Reflectivity_Attenuated(c) + Backscat_Coefficient(l,c,m) = RTSolution(l,m)%Backscat_Coefficient(c) END DO END DO Channel_Loop END DO Profile_Loop - ! Write to output file ! ...Create output netCDF file err_stat = CreateFile_netCDF( & @@ -2902,8 +2930,19 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & ' - '//TRIM(NF90_STRERROR( NF90_Status )) CALL Write_Cleanup(); RETURN END IF - - + ! ... Backscat_Coefficient variable + NF90_Status = NF90_INQ_VARID( FileId,BS_VARNAME,VarId ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error inquiring '//TRIM(Filename)//' for '//BS_VARNAME//& + ' variable ID - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Write_Cleanup(); RETURN + END IF + NF90_Status = NF90_PUT_VAR( FileId,VarID, Backscat_Coefficient) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error writing '//BS_VARNAME//' to '//TRIM(Filename)//& + ' - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Write_Cleanup(); RETURN + END IF ! Close the file NF90_Status = NF90_CLOSE( FileId ) @@ -2946,6 +2985,7 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & Single_Scatter_Albedo, & Reflectivity, & Reflectivity_Attenuated, & + Backscat_Coefficient, & STAT = alloc_stat ) IF ( alloc_stat /= 0 ) THEN msg = 'Error deallocating RTSolution output arrays' @@ -3052,8 +3092,10 @@ ELEMENTAL FUNCTION CRTM_RTSolution_Equal( x, y ) RESULT( is_equal ) ALL(x%Upwelling_Overcast_Radiance .EqualTo. y%Upwelling_Overcast_Radiance ) .AND. & ALL(x%Upwelling_Radiance .EqualTo. y%Upwelling_Radiance ) .AND. & ALL(x%Layer_Optical_Depth .EqualTo. y%Layer_Optical_Depth ) .AND. & + ALL(x%Single_Scatter_Albedo .EqualTo. y%Single_Scatter_Albedo ) .AND. & ALL(x%Reflectivity .EqualTo. y%Reflectivity ) .AND. & - ALL(x%Reflectivity_Attenuated .EqualTo. y%Reflectivity_Attenuated ) + ALL(x%Reflectivity_Attenuated .EqualTo. y%Reflectivity_Attenuated ) .AND. & + ALL(x%Backscat_Coefficient .EqualTo. y%Backscat_Coefficient ) END IF END FUNCTION CRTM_RTSolution_Equal @@ -3563,7 +3605,8 @@ FUNCTION Read_Record( & rts%Upwelling_Radiance, & rts%Layer_Optical_Depth, & rts%Reflectivity, & - rts%Reflectivity_Attenuated + rts%Reflectivity_Attenuated, & + rts%Backscat_Coefficient IF ( io_stat /= 0 ) THEN msg = 'Error reading array intermediate results - '//TRIM(io_msg) CALL Read_Record_Cleanup(); RETURN @@ -3677,7 +3720,8 @@ FUNCTION Write_Record( & rts%Upwelling_Radiance, & rts%Layer_Optical_Depth, & rts%Reflectivity, & - rts%Reflectivity_Attenuated + rts%Reflectivity_Attenuated, & + rts%Backscat_Coefficient IF ( io_stat /= 0 ) THEN msg = 'Error writing array intermediate results - '//TRIM(io_msg) CALL Write_Record_Cleanup(); RETURN @@ -4268,6 +4312,24 @@ FUNCTION CreateFile_netCDF( & CALL Create_Cleanup(); RETURN END IF + ! ... Backscat_Coefficient variable + NF90_Status = NF90_DEF_VAR( FileID, & + BS_VARNAME, & + FLOAT_TYPE, & + dimIDs=(/n_Channels_DimID, n_Layers_DimID, n_Profiles_DimID/), & + varID=VarID ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error defining '//BS_VARNAME//' variable in '//& + TRIM(Filename)//' - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Create_Cleanup(); RETURN + END IF + Put_Status(1) = NF90_PUT_ATT( FileID,VarID,UNITS_ATTNAME ,BS_UNITS ) + Put_Status(2) = NF90_PUT_ATT( FileID,VarID,FILLVALUE_ATTNAME ,FILL_FLOAT ) + IF ( ANY(Put_Status /= NF90_NOERR) ) THEN + msg = 'Error writing '//BS_VARNAME//' variable attributes to '//TRIM(Filename) + CALL Create_Cleanup(); RETURN + END IF + ! Take netCDF file out of define mode NF90_Status = NF90_ENDDEF( FileId ) IF ( NF90_Status /= NF90_NOERR ) THEN diff --git a/test/mains/regression/forward/test_AOD/test_AOD.f90 b/test/mains/regression/forward/test_AOD/test_AOD.f90 index 950a3e22..0e169835 100644 --- a/test/mains/regression/forward/test_AOD/test_AOD.f90 +++ b/test/mains/regression/forward/test_AOD/test_AOD.f90 @@ -94,6 +94,9 @@ PROGRAM test_AOD WRITE( *,'(/5x,"Initializing the CRTM...")' ) Error_Status = CRTM_Init( (/Sensor_Id/), & ChannelInfo, & + Aerosol_Model = 'GOCART-GEOS5', & + AerosolCoeff_Format = 'netCDF', & + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc', & File_Path=COEFFICIENTS_PATH) IF ( Error_Status /= SUCCESS ) THEN Message = 'Error initializing CRTM' diff --git a/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 b/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 index 0de8e78c..ddb4b9a6 100644 --- a/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 +++ b/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 @@ -97,6 +97,9 @@ PROGRAM test_AOD WRITE( *,'(/5x,"Initializing the CRTM...")' ) Error_Status = CRTM_Init( (/Sensor_Id/), & ChannelInfo, & + Aerosol_Model = 'GOCART-GEOS5', & + AerosolCoeff_Format = 'netCDF', & + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc', & File_Path=COEFFICIENTS_PATH) IF ( Error_Status /= SUCCESS ) THEN Message = 'Error initializing CRTM' From e30abbb55406dce4fb9b4c267deda38100caebf0 Mon Sep 17 00:00:00 2001 From: Cheng Dang Date: Fri, 11 Jul 2025 16:09:01 -0600 Subject: [PATCH 2/4] Quality control on kb --- src/AtmScatter/CRTM_AerosolScatter.f90 | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/AtmScatter/CRTM_AerosolScatter.f90 b/src/AtmScatter/CRTM_AerosolScatter.f90 index bfad500b..c3239845 100644 --- a/src/AtmScatter/CRTM_AerosolScatter.f90 +++ b/src/AtmScatter/CRTM_AerosolScatter.f90 @@ -256,10 +256,10 @@ FUNCTION CRTM_Compute_AerosolScatter( & ASV%asi(ka,n) ) ! Interpolation ! interpolation quality control - ! CD: kb quality control? IF( ASV%ke(ka,n) <= ZERO ) THEN ASV%ke(ka,n) = ZERO ASV%w(ka,n) = ZERO + ASV%kb(ka,n) = ZERO END IF IF( ASV%w(ka,n) <= ZERO ) THEN ASV%w(ka,n) = ZERO @@ -492,6 +492,7 @@ FUNCTION CRTM_Compute_AerosolScatter_TL( & IF( ASV%ke(ka,n) <= ZERO ) THEN ke_TL = ZERO w_TL = ZERO + kb_TL = ZERO END IF IF( ASV%w(ka,n) <= ZERO ) THEN w_TL = ZERO @@ -732,20 +733,22 @@ FUNCTION CRTM_Compute_AerosolScatter_AD( & ( bs_AD * ASV%ke(ka,n) * ASV%w(ka,n) ) ! ! Compute the adjoint of the backscattering coefficient - ! CD: placeholder for now + ! CD: to be implemented + ! interpolation quality control IF( ASV%w(ka,n) >= ONE ) THEN w_AD = ZERO END IF - IF( ASV%ke(ka,n) <= ZERO ) THEN - ke_AD = ZERO - w_AD = ZERO - END IF IF( ASV%w(ka,n) <= ZERO ) THEN w_AD = ZERO pcoeff_AD = ZERO END IF + IF( ASV%ke(ka,n) <= ZERO ) THEN + ke_AD = ZERO + w_AD = ZERO + kb_AD = ZERO + END IF CALL Get_Aerosol_Opt_AD(AScat_AD , & ! Input Atm%Aerosol(n)%Type , & ! Input @@ -924,7 +927,6 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure ! Perform Interpolation CALL interp_2D( AeroC%ke( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k), asi%wlp, asi%hlp, ke ) CALL interp_2D( AeroC%w( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k), asi%wlp, asi%hlp, w ) - ! CD: This is for GOCART-GEOS5 only, for NAAPS, kb is zero IF (AeroC%Scheme == 'GOCART-GEOS5') THEN CALL interp_2D( AeroC%kb( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k), asi%wlp, asi%hlp, kb ) ELSE From 910dc3f44b94b8435d231da7f30193ead1beb813 Mon Sep 17 00:00:00 2001 From: Cheng Dang Date: Mon, 14 Jul 2025 20:51:54 -0600 Subject: [PATCH 3/4] Remove 2pi --- src/AtmScatter/CRTM_AOD_Module.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/AtmScatter/CRTM_AOD_Module.f90 b/src/AtmScatter/CRTM_AOD_Module.f90 index 77e70e87..53a485a3 100644 --- a/src/AtmScatter/CRTM_AOD_Module.f90 +++ b/src/AtmScatter/CRTM_AOD_Module.f90 @@ -19,7 +19,7 @@ MODULE CRTM_AOD_Module ! Module usage ! ------------ USE Message_Handler, ONLY: SUCCESS, FAILURE, Display_Message - USE CRTM_Parameters, ONLY: ZERO, TWO, PI, & + USE CRTM_Parameters, ONLY: ZERO, & MAX_N_PHASE_ELEMENTS, & MAX_N_LEGENDRE_TERMS USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type, & @@ -346,7 +346,7 @@ FUNCTION CRTM_AOD( & ! Save the Backscat_Coefficient, only non-zero for GOCART-GEOS5 ! CD: double check if should be divided by 2PI - RTSolution(ln,m)%Backscat_Coefficient(1:Atmosphere(m)%n_Layers) = AtmOptics%Backscat_Coefficient/(TWO*PI) + RTSolution(ln,m)%Backscat_Coefficient(1:Atmosphere(m)%n_Layers) = AtmOptics%Backscat_Coefficient END DO Channel_Loop @@ -675,8 +675,8 @@ FUNCTION CRTM_AOD_TL( & ! Save the backscatter coefficient, only non-zero for GOCART-GEOS5 ! CD: double check if should be divided by 2PI - RTSolution(ln,m)%Backscat_Coefficient(1:Atmosphere(m)%n_Layers) = AtmOptics%Backscat_Coefficient/(TWO*PI) - RTSolution_TL(ln,m)%Backscat_Coefficient(1:Atmosphere(m)%n_Layers) = AtmOptics_TL%Backscat_Coefficient/(TWO*PI) + RTSolution(ln,m)%Backscat_Coefficient(1:Atmosphere(m)%n_Layers) = AtmOptics%Backscat_Coefficient + RTSolution_TL(ln,m)%Backscat_Coefficient(1:Atmosphere(m)%n_Layers) = AtmOptics_TL%Backscat_Coefficient END DO Channel_Loop From a445d401a35327bf16d2e5f54823cf581be59716 Mon Sep 17 00:00:00 2001 From: Cheng Dang Date: Wed, 14 Jan 2026 21:37:23 -0700 Subject: [PATCH 4/4] change nc filename to nc4 for now --- test/mains/regression/forward/test_AOD/test_AOD.f90 | 2 +- test/mains/regression/k_matrix/test_AOD/test_AOD.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/mains/regression/forward/test_AOD/test_AOD.f90 b/test/mains/regression/forward/test_AOD/test_AOD.f90 index 0e169835..7944b373 100644 --- a/test/mains/regression/forward/test_AOD/test_AOD.f90 +++ b/test/mains/regression/forward/test_AOD/test_AOD.f90 @@ -96,7 +96,7 @@ PROGRAM test_AOD ChannelInfo, & Aerosol_Model = 'GOCART-GEOS5', & AerosolCoeff_Format = 'netCDF', & - AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc', & + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc4', & File_Path=COEFFICIENTS_PATH) IF ( Error_Status /= SUCCESS ) THEN Message = 'Error initializing CRTM' diff --git a/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 b/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 index ddb4b9a6..e1f437a5 100644 --- a/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 +++ b/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 @@ -99,7 +99,7 @@ PROGRAM test_AOD ChannelInfo, & Aerosol_Model = 'GOCART-GEOS5', & AerosolCoeff_Format = 'netCDF', & - AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc', & + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc4', & File_Path=COEFFICIENTS_PATH) IF ( Error_Status /= SUCCESS ) THEN Message = 'Error initializing CRTM'