diff --git a/.gitignore b/.gitignore index 62f2735a..d2322c0c 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,4 @@ build.bash rebuild.bash .gitignore *~ +conductor/ diff --git a/src/AtmScatter/AerosolScatter/ASvar_Define.f90 b/src/AtmScatter/AerosolScatter/ASvar_Define.f90 index cd8c2901..345566b3 100644 --- a/src/AtmScatter/AerosolScatter/ASvar_Define.f90 +++ b/src/AtmScatter/AerosolScatter/ASvar_Define.f90 @@ -241,6 +241,11 @@ SUBROUTINE ASvar_Inspect( self) WRITE(*,'(5x,"g Aerosol index #",i0)') i4 WRITE(*,'(5(1x,es22.15,:))') self%g(:,i4) END DO + WRITE(*,'(3x,"Backscatter coefficient (kb) :")') + DO i4 = 1, self%n_Aerosols + WRITE(*,'(5x,"kb Aerosol index #",i0)') i4 + WRITE(*,'(5(1x,es22.15,:))') self%kb(:,i4) + END DO WRITE(*,'(3x,"Phase coefficients (pcoeff) :")') DO i4 = 1, self%n_Aerosols WRITE(*,'(5x,"pcoeff Aerosol index #",i0)') i4 @@ -588,6 +593,13 @@ FUNCTION ASvar_ReadFile( & msg = 'Error reading asymmetry factor - '//TRIM(io_msg) CALL Read_Cleanup(); RETURN END IF + ! ...Backscatter coefficient + READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) & + ASvar%kb + IF ( io_stat /= 0 ) THEN + msg = 'Error reading backscatter coefficient - '//TRIM(io_msg) + CALL Read_Cleanup(); RETURN + END IF ! ...Phase coefficients READ( fid, IOSTAT=io_stat, IOMSG=io_msg ) & ASvar%pcoeff @@ -763,6 +775,13 @@ FUNCTION ASvar_WriteFile( & msg = 'Error writing asymmetry factor - '//TRIM(io_msg) CALL Write_Cleanup(); RETURN END IF + ! ...Backscatter coefficient + WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) & + ASvar%kb + IF ( io_stat /= 0 ) THEN + msg = 'Error writing backscatter coefficient - '//TRIM(io_msg) + CALL Write_Cleanup(); RETURN + END IF ! ...Phase coefficients WRITE( fid, IOSTAT=io_stat, IOMSG=io_msg ) & ASvar%pcoeff @@ -842,6 +861,7 @@ ELEMENTAL FUNCTION ASvar_Equal( x, y ) RESULT( is_equal ) IF ( ALL(x%ke .EqualTo. y%ke ) .AND. & ALL(x%w .EqualTo. y%w ) .AND. & ALL(x%g .EqualTo. y%g ) .AND. & + ALL(x%kb .EqualTo. y%kb ) .AND. & ALL(x%pcoeff .EqualTo. y%pcoeff ) .AND. & ALL(x%total_bs .EqualTo. y%total_bs ) ) & is_equal = .TRUE. diff --git a/src/AtmScatter/CRTM_AerosolScatter.f90 b/src/AtmScatter/CRTM_AerosolScatter.f90 index c3239845..26353925 100644 --- a/src/AtmScatter/CRTM_AerosolScatter.f90 +++ b/src/AtmScatter/CRTM_AerosolScatter.f90 @@ -251,6 +251,7 @@ FUNCTION CRTM_Compute_AerosolScatter( & Atm%Relative_Humidity(ka) , & ! Input ASV%ke(ka,n) , & ! Output ASV%w(ka,n) , & ! Output + ASV%g(ka,n) , & ! Output ASV%kb(ka,n) , & ! Output ASV%pcoeff(:,:,ka,n) , & ! Output ASV%asi(ka,n) ) ! Interpolation @@ -304,9 +305,10 @@ 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%Asymmetry_Factor(ka) = AScat%Asymmetry_Factor(ka) + (ASV%g(ka,n) * bs) AScat%Backscat_Coefficient(ka) = AScat%Backscat_Coefficient(ka) + & - (ASV%kb(ka,n) * Atm%Aerosol(n)%Concentration(ka)) + (ASV%kb(ka,n) * Atm%Aerosol(n)%Concentration(ka)) + DO m = 1, AScat%n_Phase_Elements DO l = 0, AScat%n_Legendre_Terms @@ -437,7 +439,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, kb_TL + REAL(fp) :: ke_TL, w_TL, g_TL, kb_TL REAL(fp) :: pcoeff_TL(0:AScat%n_Legendre_Terms, AScat%n_Phase_Elements) REAL(fp) :: bs, bs_TL @@ -484,6 +486,7 @@ FUNCTION CRTM_Compute_AerosolScatter_TL( & Atm_TL%Relative_Humidity(ka) , & ! TL Input ke_TL , & ! TL Output w_TL , & ! TL Output + g_TL , & ! TL Output kb_TL , & ! TL Output pcoeff_TL , & ! TL Output ASV%asi(ka,n) ) ! Interpolation @@ -519,6 +522,11 @@ FUNCTION CRTM_Compute_AerosolScatter_TL( & (Atm%Aerosol(n)%Concentration(ka) * ke_TL * ASV%w(ka,n) ) + & (Atm%Aerosol(n)%Concentration(ka) * ASV%ke(ka,n) * w_TL ) AScat_TL%Single_Scatter_Albedo(ka) = AScat_TL%Single_Scatter_Albedo(ka) + bs_TL + AScat_TL%Asymmetry_Factor(ka) = AScat_TL%Asymmetry_Factor(ka) + & + (g_TL * bs) + (ASV%g(ka,n) * bs_TL) + 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)) DO m = 1, n_Phase_Elements DO l = 0, n_Legendre_Terms @@ -651,7 +659,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, kb_AD + REAL(fp) :: ke_AD, w_AD, g_AD, kb_AD REAL(fp) :: pcoeff_AD(0:AScat%n_Legendre_Terms, AScat%n_Phase_Elements) REAL(fp) :: bs, bs_AD @@ -699,6 +707,7 @@ FUNCTION CRTM_Compute_AerosolScatter_AD( & pcoeff_AD = ZERO ke_AD = ZERO w_AD = ZERO + g_AD = ZERO kb_AD = ZERO ! Compute the adjoint of the @@ -713,6 +722,10 @@ FUNCTION CRTM_Compute_AerosolScatter_AD( & pcoeff_AD(l,m) = pcoeff_AD(l,m) + (bs * AScat_AD%Phase_Coefficient(l,m,ka)) END DO END DO + ! Adjoint of asymmetry factor + g_AD = g_AD + (bs * AScat_AD%Asymmetry_Factor(ka)) + bs_AD = bs_AD + (ASV%g(ka,n) * AScat_AD%Asymmetry_Factor(ka)) + ! NOTE: bs_AD is not reinitialized after this ! point since it is reinitialized at the ! start of the Aerosol_Layer_loop @@ -726,6 +739,11 @@ FUNCTION CRTM_Compute_AerosolScatter_AD( & (ASV%ke(ka,n) * AScat_AD%Optical_Depth(ka)) ke_AD = ke_AD + (Atm%Aerosol(n)%Concentration(ka) * AScat_AD%Optical_Depth(ka)) + ! Compute the adjoint of the backscatter coefficient + Atm_AD%Aerosol(n)%Concentration(ka) = Atm_AD%Aerosol(n)%Concentration(ka) + & + (ASV%kb(ka,n) * AScat_AD%Backscat_Coefficient(ka)) + kb_AD = kb_AD + (Atm%Aerosol(n)%Concentration(ka) * AScat_AD%Backscat_Coefficient(ka)) + ! Compute the adjoint of the volume ! scattering coefficient. ke_AD = ke_AD + (Atm%Aerosol(n)%Concentration(ka) * bs_AD * ASV%w(ka,n) ) @@ -757,6 +775,7 @@ FUNCTION CRTM_Compute_AerosolScatter_AD( & ASV%kb(ka,n) , & ! input ke_AD , & ! AD Input w_AD , & ! AD Input + g_AD , & ! AD Input kb_AD , & ! AD Input pcoeff_AD , & ! AD Input Atm_AD%Aerosol(n)%Effective_Radius(ka) , & ! AD Output @@ -796,7 +815,8 @@ 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 + g , & ! Output asymmetry factor + kb , & ! Output backscatter coefficient pcoeff , & ! Output spherical Legendre coefficients asi ) ! Output interpolation data ! Arguments @@ -807,6 +827,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) :: g REAL(fp) , INTENT(OUT) :: kb REAL(fp) , INTENT(IN OUT) :: pcoeff(0:,:) TYPE(ASinterp_type) , INTENT(IN OUT) :: asi @@ -822,6 +843,7 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure IF ( Aerosol_Type == AerosolCoeff_BYPASS_AEROSOL ) THEN ke = ZERO w = ZERO + g = ZERO kb = ZERO pcoeff = ZERO RETURN @@ -859,6 +881,8 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure ! Perform Interpolation CALL interp_2D( AeroC%ke( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, fix_sig, k), asi%wlp, asi%xlp, ke ) CALL interp_2D( AeroC%w( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, fix_sig, k), asi%wlp, asi%xlp, w ) + CALL interp_2D( AeroC%g( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, fix_sig, k), asi%wlp, asi%xlp, g ) + CALL interp_2D( AeroC%kb( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, fix_sig, k), asi%wlp, asi%xlp, kb ) IF (AerosolScatter%n_Phase_Elements > 0 .and. AerosolScatter%Include_Scattering ) THEN pcoeff(0,1) = POINT_5 DO m = 1, AerosolScatter%n_Phase_Elements @@ -898,6 +922,8 @@ SUBROUTINE Get_Aerosol_Opt( AerosolScatter, & ! Input AerosolScatter structure ! Perform Interpolation CALL interp_3D( AeroC%ke( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, asi%k1:asi%k2, k), asi%wlp, asi%xlp, asi%vlp, ke ) CALL interp_3D( AeroC%w( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, asi%k1:asi%k2, k), asi%wlp, asi%xlp, asi%vlp, w ) + CALL interp_3D( AeroC%g( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, asi%k1:asi%k2, k), asi%wlp, asi%xlp, asi%vlp, g ) + CALL interp_3D( AeroC%kb( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, asi%k1:asi%k2, k), asi%wlp, asi%xlp, asi%vlp, kb ) IF (AerosolScatter%n_Phase_Elements > 0 .and. AerosolScatter%Include_Scattering ) THEN pcoeff(0,1) = POINT_5 DO m = 1, AerosolScatter%n_Phase_Elements @@ -927,11 +953,7 @@ 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 ) - 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%g( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k), asi%wlp, asi%hlp, g ) 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 @@ -970,7 +992,8 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st 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 + g_TL , & ! Output TL asymmetry factor + kb_TL , & ! Output TL backscatter coefficient pcoeff_TL , & ! TL Output spherical Legendre coefficients asi ) ! Input interpolation data @@ -979,7 +1002,7 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st INTEGER , INTENT(IN) :: Aerosol_Type 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, kb_TL + REAL(fp), INTENT(OUT) :: ke_TL, w_TL, g_TL, kb_TL REAL(fp), INTENT(IN OUT) :: pcoeff_TL(0:,:) TYPE(ASinterp_type), INTENT(IN) :: asi ! Local variables @@ -1002,6 +1025,7 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st IF ( Aerosol_Type == AerosolCoeff_BYPASS_AEROSOL ) THEN ke_TL = ZERO w_TL = ZERO + g_TL = ZERO kb_TL = ZERO pcoeff_TL = ZERO RETURN @@ -1014,6 +1038,8 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st IF ( asi%f_outbound .AND. asi%r_outbound .AND. asi%v_outbound .AND. asi%h_outbound) THEN ke_TL = ZERO w_TL = ZERO + g_TL = ZERO + kb_TL = ZERO pcoeff_TL = ZERO RETURN END IF @@ -1070,6 +1096,16 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st CALL interp_2D_TL( zg , asi%wlp, asi%xlp, & ! FWD Input zg_TL, wlp_TL , xlp_TL , & ! TL Input w_TL ) ! TL Output + ! Asymmetry factor + zg => AeroC%g( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, fix_sig, k) + CALL interp_2D_TL( zg , asi%wlp, asi%xlp, & ! FWD Input + zg_TL, wlp_TL , xlp_TL , & ! TL Input + g_TL ) ! TL Output + ! Backscatter coefficient + zg => AeroC%kb( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, fix_sig, k) + CALL interp_2D_TL( zg , asi%wlp, asi%xlp, & ! FWD Input + zg_TL, wlp_TL , xlp_TL , & ! TL Input + kb_TL ) ! TL Output ! Phase matrix coefficients IF (AerosolScatter_TL%n_Phase_Elements > 0 .and. AerosolScatter_TL%Include_Scattering ) THEN pcoeff_TL(0,1) = ZERO @@ -1119,6 +1155,16 @@ SUBROUTINE Get_Aerosol_Opt_TL(AerosolScatter_TL, & ! Input AerosolScatterTL st CALL interp_3D_TL( zc , asi%wlp, asi%xlp, asi%vlp, & ! FWD Input zc_TL, wlp_TL , xlp_TL , vlp_TL , & ! TL Input w_TL ) ! TL Output + ! Asymmetry factor + zc => AeroC%g( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, asi%k1:asi%k2, k) + CALL interp_3D_TL( zc , asi%wlp, asi%xlp, asi%vlp, & ! FWD Input + zc_TL, wlp_TL , xlp_TL , vlp_TL , & ! TL Input + g_TL ) ! TL Output + ! Backscatter coefficient + zc => AeroC%kb( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, asi%k1:asi%k2, k) + CALL interp_3D_TL( zc , asi%wlp, asi%xlp, asi%vlp, & ! FWD Input + zc_TL, wlp_TL , xlp_TL , vlp_TL , & ! TL Input + kb_TL ) ! TL Output ! Phase matrix coefficients IF (AerosolScatter_TL%n_Phase_Elements > 0 .and. AerosolScatter_TL%Include_Scattering ) THEN pcoeff_TL(0,1) = ZERO @@ -1161,15 +1207,21 @@ 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 + ! Asymmetry factor + zg => AeroC%g( 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 + g_TL ) ! TL Output + ! Backscatter coefficient + 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 ! 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 + 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 ! Phase matrix coefficients IF (AerosolScatter_TL%n_Phase_Elements > 0 .and. AerosolScatter_TL%Include_Scattering ) THEN pcoeff_TL(0,1) = ZERO @@ -1211,7 +1263,8 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st kb , & ! Input ke_AD , & ! AD Input extinction cross section w_AD , & ! AD Input single scatter albedo - kb_AD , & ! AD Input backscattering coefficient + g_AD , & ! AD Input asymmetry factor + kb_AD , & ! AD Input backscatter coefficient pcoeff_AD , & ! AD Input spherical Legendre coefficients Reff_AD , & ! AD Outputmode radius Rsig_AD , & ! AD Output radius deviation @@ -1223,6 +1276,7 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st 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) :: g_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 @@ -1250,6 +1304,7 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st Reff_AD = ZERO ke_AD = ZERO w_AD = ZERO + g_AD = ZERO kb_AD = ZERO pcoeff_AD = ZERO Rsig_AD = ZERO @@ -1265,6 +1320,7 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st Reff_AD = ZERO ke_AD = ZERO w_AD = ZERO + g_AD = ZERO kb_AD = ZERO pcoeff_AD = ZERO Rsig_AD = ZERO @@ -1322,6 +1378,16 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st END IF END IF + ! Backscatter coefficient + zg => AeroC%kb( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, fix_sig, k) + CALL interp_2D_AD( zg , asi%wlp, asi%xlp, & ! FWD Input + kb_AD , & ! AD Input + zg_AD, wlp_AD , xlp_AD ) ! AD Output + ! Asymmetry factor + zg => AeroC%g( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, fix_sig, k) + CALL interp_2D_AD( zg , asi%wlp, asi%xlp, & ! FWD Input + g_AD , & ! AD Input + zg_AD, wlp_AD , xlp_AD ) ! AD Output ! Single scatter albedo zg => AeroC%w( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, fix_sig, k) CALL interp_2D_AD( zg , asi%wlp, asi%xlp, & ! FWD Input @@ -1378,6 +1444,16 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st END IF END IF + ! Backscatter coefficient + zc => AeroC%kb( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, asi%k1:asi%k2, k) + CALL interp_3D_AD( zc , asi%wlp, asi%xlp, asi%vlp, & ! FWD Input + kb_AD , & ! AD Input + zc_AD, wlp_AD, xlp_AD, vlp_AD ) ! AD Output + ! Asymmetry factor + zc => AeroC%g( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, asi%k1:asi%k2, k) + CALL interp_3D_AD( zc , asi%wlp, asi%xlp, asi%vlp, & ! FWD Input + g_AD , & ! AD Input + zc_AD, wlp_AD, xlp_AD, vlp_AD ) ! AD Output ! Single scatter albedo zc => AeroC%w( asi%i1:asi%i2, fix_rh, asi%j1:asi%j2, asi%k1:asi%k2, k) CALL interp_3D_AD( zc , asi%wlp, asi%xlp, asi%vlp, & ! FWD Input @@ -1441,6 +1517,16 @@ SUBROUTINE Get_Aerosol_Opt_AD( AerosolScatter_AD, & ! Input AerosolScatter AD st END IF END IF + ! Backscatter coefficient + 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 + ! Asymmetry factor + zg => AeroC%g( asi%i1:asi%i2, asi%h1:asi%h2, fix_r, fix_sig, k) + CALL interp_2D_AD( zg , asi%wlp, asi%hlp, & ! FWD Input + g_AD , & ! AD Input + zg_AD, wlp_AD , hlp_AD ) ! AD Output ! Single scatter albedo 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 diff --git a/src/Atmosphere/CRTM_Atmosphere.f90 b/src/Atmosphere/CRTM_Atmosphere.f90 index 4951161c..44cace4d 100644 --- a/src/Atmosphere/CRTM_Atmosphere.f90 +++ b/src/Atmosphere/CRTM_Atmosphere.f90 @@ -389,6 +389,7 @@ FUNCTION CRTM_Atmosphere_AddLayers( & IF ( Atm_In%n_Aerosols > 0 ) THEN DO i = 1, Atm_In%n_Aerosols Atm_Out%Aerosol(i)%Effective_Radius(1:n) = ZERO + Atm_Out%Aerosol(i)%Effective_Variance(1:n) = ZERO Atm_Out%Aerosol(i)%Concentration(1:n) = ZERO END DO END IF @@ -622,6 +623,8 @@ FUNCTION CRTM_Atmosphere_AddLayers_AD( & Atm_Out_AD%Aerosol(i)%Concentration(n+1:nt) Atm_In_AD%Aerosol(i)%Effective_Radius(1:no) = Atm_In_AD%Aerosol(i)%Effective_Radius(1:no) + & Atm_Out_AD%Aerosol(i)%Effective_Radius(n+1:nt) + Atm_In_AD%Aerosol(i)%Effective_Variance(1:no) = Atm_In_AD%Aerosol(i)%Effective_Variance(1:no) + & + Atm_Out_AD%Aerosol(i)%Effective_Variance(n+1:nt) Atm_In_AD%Aerosol(i)%Type = Atm_Out_AD%Aerosol(i)%Type END DO END IF @@ -645,6 +648,9 @@ FUNCTION CRTM_Atmosphere_AddLayers_AD( & END DO ! ...Temperature data Atm_In_AD%Temperature(1:no) = Atm_In_AD%Temperature(1:no) + Atm_Out_AD%Temperature(n+1:nt) + ! ...Relative humidity data + Atm_In_AD%Relative_Humidity(1:no) = Atm_In_AD%Relative_Humidity(1:no) + & + Atm_Out_AD%Relative_Humidity(n+1:nt) ! ...Pressure data Atm_In_AD%Pressure(1:no) = Atm_In_AD%Pressure(1:no) + Atm_Out_AD%Pressure(n+1:nt) Atm_In_AD%Level_Pressure(0:no) = Atm_In_AD%Level_Pressure(0:no) + Atm_Out_AD%Level_Pressure(n:nt) @@ -748,6 +754,7 @@ FUNCTION CRTM_Atmosphere_ClearSkyCopy( atm, atm_clear ) RESULT( err_stat ) atm_clear%Level_Pressure = atm%Level_Pressure(0:k) atm_clear%Pressure = atm%Pressure(1:k) atm_clear%Temperature = atm%Temperature(1:k) + atm_clear%Relative_Humidity = atm%Relative_Humidity(1:k) atm_clear%Absorber = atm%Absorber(1:k,:) atm_clear%Cloud_Fraction = atm%Cloud_Fraction(1:k) ! ...Aerosol components @@ -884,6 +891,7 @@ FUNCTION CRTM_Atmosphere_ClearSkyCopy_TL( & atm_clear_TL%Level_Pressure = atm_TL%Level_Pressure(0:k) atm_clear_TL%Pressure = atm_TL%Pressure(1:k) atm_clear_TL%Temperature = atm_TL%Temperature(1:k) + atm_clear_TL%Relative_Humidity = atm_TL%Relative_Humidity(1:k) atm_clear_TL%Absorber = atm_TL%Absorber(1:k,:) atm_clear_TL%Cloud_Fraction = atm_TL%Cloud_Fraction(1:k) ! ...Aerosol components @@ -1019,6 +1027,7 @@ FUNCTION CRTM_Atmosphere_ClearSkyCopy_AD( & atm_AD%Level_Pressure(0:k) = atm_AD%Level_Pressure(0:k) + atm_clear_AD%Level_Pressure(0:k) atm_AD%Pressure(1:k) = atm_AD%Pressure(1:k) + atm_clear_AD%Pressure(1:k) atm_AD%Temperature(1:k) = atm_AD%Temperature(1:k) + atm_clear_AD%Temperature(1:k) + atm_AD%Relative_Humidity(1:k) = atm_AD%Relative_Humidity(1:k) + atm_clear_AD%Relative_Humidity(1:k) atm_AD%Absorber(1:k,:) = atm_AD%Absorber(1:k,:) + atm_clear_AD%Absorber(1:k,:) atm_AD%Cloud_Fraction(1:k) = atm_AD%Cloud_Fraction(1:k) + atm_clear_AD%Cloud_Fraction(1:k) diff --git a/src/Coefficients/AerosolCoeff/AerosolCoeff_Inspect/AerosolCoeff_Inspect.f90 b/src/Coefficients/AerosolCoeff/AerosolCoeff_Inspect/AerosolCoeff_Inspect.f90 index 4867675b..0b864587 100644 --- a/src/Coefficients/AerosolCoeff/AerosolCoeff_Inspect/AerosolCoeff_Inspect.f90 +++ b/src/Coefficients/AerosolCoeff/AerosolCoeff_Inspect/AerosolCoeff_Inspect.f90 @@ -27,20 +27,26 @@ PROGRAM AerosolCoeff_Inspect ! Parameters ! ---------- CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'AerosolCoeff_Inspect' - CHARACTER(*), PARAMETER :: PROGRAM_RCS_ID = & + CHARACTER(*), PARAMETER :: PROGRAM_RCS_ID = '$Id$' ! --------- ! Variables ! --------- INTEGER :: Error_Status CHARACTER(256) :: Filename + CHARACTER(256) :: Aerosol_Model TYPE(AerosolCoeff_type) :: A ! Output program header CALL Program_Message( PROGRAM_NAME, & 'Program to display the contents of a CRTM Binary format '//& 'AerosolCoeff file to stdout.', & - '$Revision$' ) + PROGRAM_RCS_ID ) + + ! Get the aerosol model + WRITE( *,FMT='(/5x,"Enter the Aerosol Model (CRTM, CMAQ, GOCART-GEOS5, or NAAPS): ")',ADVANCE='NO' ) + READ( *,'(a)' ) Aerosol_Model + Aerosol_Model = ADJUSTL( Aerosol_Model ) ! Get the filename WRITE( *,FMT='(/5x,"Enter the Binary AerosolCoeff filename: ")',ADVANCE='NO' ) @@ -54,7 +60,7 @@ PROGRAM AerosolCoeff_Inspect END IF ! Read the binary data file - Error_Status = AerosolCoeff_Binary_ReadFile( Filename, A ) + Error_Status = AerosolCoeff_Binary_ReadFile( Aerosol_Model, Filename, A ) IF ( Error_Status /= SUCCESS ) THEN CALL Display_Message( PROGRAM_NAME, & 'Error reading Binary AerosolCoeff file '//& diff --git a/src/Coefficients/AerosolCoeff/AerosolCoeff_netCDF_IO.f90 b/src/Coefficients/AerosolCoeff/AerosolCoeff_netCDF_IO.f90 index 1bcb7460..806f3619 100644 --- a/src/Coefficients/AerosolCoeff/AerosolCoeff_netCDF_IO.f90 +++ b/src/Coefficients/AerosolCoeff/AerosolCoeff_netCDF_IO.f90 @@ -1152,16 +1152,21 @@ FUNCTION AerosolCoeff_netCDF_ReadFile( & ! ...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 + ! Older GOCART-GEOS5 files may not include kb, so treat it as optional. + IF ( NF90_Status .EQ. NF90_ENOTVAR ) THEN + AerosolCoeff%kb = 0.0_fp + ELSE + 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 END IF ! ...ke variable diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b49c8fd1..40859bc6 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -308,6 +308,24 @@ add_test(NAME test_Unit_Aerosol_Bypass COMMAND $) set_tests_properties(test_Unit_Aerosol_Bypass PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") +add_executable(Unit_AerosolScatter_TL mains/unit/Unit_Test/test_AerosolScatter_TL.f90) +target_link_libraries(Unit_AerosolScatter_TL PRIVATE crtm) +add_test(NAME test_Unit_AerosolScatter_TL + COMMAND $) +set_tests_properties(test_Unit_AerosolScatter_TL PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") + +add_executable(Unit_AerosolScatter_AD mains/unit/Unit_Test/test_AerosolScatter_AD.f90) +target_link_libraries(Unit_AerosolScatter_AD PRIVATE crtm) +add_test(NAME test_Unit_AerosolScatter_AD + COMMAND $) +set_tests_properties(test_Unit_AerosolScatter_AD PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") + +add_executable(Unit_AerosolScatter_K mains/unit/Unit_Test/test_AerosolScatter_K.f90) +target_link_libraries(Unit_AerosolScatter_K PRIVATE crtm) +add_test(NAME test_Unit_AerosolScatter_K + COMMAND $) +set_tests_properties(test_Unit_AerosolScatter_K PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") + add_executable(Unit_Aerosol_Bypass_TL mains/unit/Unit_Test/test_Aerosol_Bypass_TL.f90) target_link_libraries(Unit_Aerosol_Bypass_TL PRIVATE crtm) add_test(NAME test_Unit_Aerosol_Bypass_TL @@ -378,6 +396,15 @@ foreach(testtype IN LISTS CCoef_Utils) target_link_libraries(${testtype} PRIVATE crtm) endforeach() +#AerosolCoeff utilities +list (APPEND ACoeff_Utils + AerosolCoeff_Inspect +) +foreach(testtype IN LISTS ACoeff_Utils) + add_executable(${testtype} ${CRTM_SOURCE_DIR}/src/Coefficients/AerosolCoeff/${testtype}/${testtype}.f90) + target_link_libraries(${testtype} PRIVATE crtm) +endforeach() + list (APPEND SEcategoryCoef_Utils SEcategory_NC2BIN SEcategory_BIN2NC @@ -546,6 +573,10 @@ set_tests_properties(test_TL_convergence_active_sensor PROPERTIES ENVIRONMENT OM list( APPEND crtm_test_input AerosolCoeff/Little_Endian/AerosolCoeff.bin +AerosolCoeff/Little_Endian/AerosolCoeff.GOCART-GEOS5.bin +AerosolCoeff/Little_Endian/AerosolCoeff.CMAQ.bin +AerosolCoeff/netCDF/AerosolCoeff.GOCART-GEOS5.nc4 +AerosolCoeff/netCDF/AerosolCoeff.GOCART-GEOS5.BRC.kb.v2.nc AerosolCoeff/netCDF/AerosolCoeff.nc4 CloudCoeff/Little_Endian/CloudCoeff.bin CloudCoeff/netCDF/CloudCoeff.nc4 diff --git a/test/mains/regression/forward/test_AOD/test_AOD.f90 b/test/mains/regression/forward/test_AOD/test_AOD.f90 index 7944b373..5e00d632 100644 --- a/test/mains/regression/forward/test_AOD/test_AOD.f90 +++ b/test/mains/regression/forward/test_AOD/test_AOD.f90 @@ -12,6 +12,7 @@ PROGRAM test_AOD ! ! Module usage USE CRTM_Module + USE File_Utility, ONLY: File_Exists ! Disable all implicit typing IMPLICIT NONE ! ============================================================================ @@ -45,10 +46,13 @@ PROGRAM test_AOD CHARACTER(256) :: Message CHARACTER(256) :: Version CHARACTER(256) :: Sensor_Id + CHARACTER(256) :: AerosolCoeff_File INTEGER :: Error_Status INTEGER :: Allocate_Status INTEGER :: n_Channels INTEGER :: l, m + LOGICAL :: has_new_coeff + LOGICAL :: has_old_coeff ! Declarations for RTSolution comparison INTEGER :: n_l, n_m CHARACTER(256) :: rts_File @@ -92,12 +96,28 @@ PROGRAM test_AOD ! 2a. Initialise for the requested sensor ! --------------------------------------- WRITE( *,'(/5x,"Initializing the CRTM...")' ) + has_new_coeff = File_Exists(COEFFICIENTS_PATH//'AerosolCoeff.GOCART-GEOS5.BRC.kb.v2.nc') + has_old_coeff = File_Exists(COEFFICIENTS_PATH//'AerosolCoeff.GOCART-GEOS5.nc4') + IF ( has_new_coeff ) THEN + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.BRC.kb.v2.nc' + ELSE + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc4' + END IF Error_Status = CRTM_Init( (/Sensor_Id/), & ChannelInfo, & Aerosol_Model = 'GOCART-GEOS5', & AerosolCoeff_Format = 'netCDF', & - AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc4', & + AerosolCoeff_File = TRIM(AerosolCoeff_File), & File_Path=COEFFICIENTS_PATH) + IF ( Error_Status /= SUCCESS .AND. has_old_coeff .AND. TRIM(AerosolCoeff_File) /= 'AerosolCoeff.GOCART-GEOS5.nc4' ) THEN + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc4' + Error_Status = CRTM_Init( (/Sensor_Id/), & + ChannelInfo, & + Aerosol_Model = 'GOCART-GEOS5', & + AerosolCoeff_Format = 'netCDF', & + AerosolCoeff_File = TRIM(AerosolCoeff_File), & + File_Path=COEFFICIENTS_PATH) + END IF IF ( Error_Status /= SUCCESS ) THEN Message = 'Error initializing CRTM' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) 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 e1f437a5..d4a82428 100644 --- a/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 +++ b/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 @@ -12,6 +12,7 @@ PROGRAM test_AOD ! ! Module usage USE CRTM_Module + USE File_Utility, ONLY: File_Exists ! Disable all implicit typing IMPLICIT NONE ! ============================================================================ @@ -43,10 +44,13 @@ PROGRAM test_AOD CHARACTER(256) :: Message CHARACTER(256) :: Version CHARACTER(256) :: Sensor_Id + CHARACTER(256) :: AerosolCoeff_File INTEGER :: Error_Status INTEGER :: Allocate_Status INTEGER :: n_Channels INTEGER :: l, m + LOGICAL :: has_new_coeff + LOGICAL :: has_old_coeff ! Declarations for Jacobian comparisons INTEGER :: n_l, n_m CHARACTER(256) :: atmk_File @@ -95,12 +99,28 @@ PROGRAM test_AOD ! 2a. Initialise for the requested sensor ! --------------------------------------- WRITE( *,'(/5x,"Initializing the CRTM...")' ) + has_new_coeff = File_Exists(COEFFICIENTS_PATH//'AerosolCoeff.GOCART-GEOS5.BRC.kb.v2.nc') + has_old_coeff = File_Exists(COEFFICIENTS_PATH//'AerosolCoeff.GOCART-GEOS5.nc4') + IF ( has_new_coeff ) THEN + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.BRC.kb.v2.nc' + ELSE + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc4' + END IF Error_Status = CRTM_Init( (/Sensor_Id/), & ChannelInfo, & Aerosol_Model = 'GOCART-GEOS5', & AerosolCoeff_Format = 'netCDF', & - AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc4', & + AerosolCoeff_File = TRIM(AerosolCoeff_File), & File_Path=COEFFICIENTS_PATH) + IF ( Error_Status /= SUCCESS .AND. has_old_coeff .AND. TRIM(AerosolCoeff_File) /= 'AerosolCoeff.GOCART-GEOS5.nc4' ) THEN + AerosolCoeff_File = 'AerosolCoeff.GOCART-GEOS5.nc4' + Error_Status = CRTM_Init( (/Sensor_Id/), & + ChannelInfo, & + Aerosol_Model = 'GOCART-GEOS5', & + AerosolCoeff_Format = 'netCDF', & + AerosolCoeff_File = TRIM(AerosolCoeff_File), & + File_Path=COEFFICIENTS_PATH) + END IF IF ( Error_Status /= SUCCESS ) THEN Message = 'Error initializing CRTM' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) diff --git a/test/mains/unit/Unit_Test/test_AerosolScatter_AD.f90 b/test/mains/unit/Unit_Test/test_AerosolScatter_AD.f90 new file mode 100644 index 00000000..8b32c5c7 --- /dev/null +++ b/test/mains/unit/Unit_Test/test_AerosolScatter_AD.f90 @@ -0,0 +1,147 @@ + +! +! test_AerosolScatter_AD +! +! Unit test for CRTM_Compute_AerosolScatter_AD +! Verifies Adjoint model consistency with Tangent Linear model. +! + +PROGRAM test_AerosolScatter_AD + + USE Type_Kinds, ONLY: fp + USE Message_Handler, ONLY: SUCCESS, FAILURE, Display_Message + USE CRTM_Parameters, ONLY: ZERO, ONE + USE CRTM_SpcCoeff, ONLY: SC, CRTM_SpcCoeff_Load, CRTM_SpcCoeff_Destroy + USE CRTM_AerosolCoeff, ONLY: AeroC, CRTM_AerosolCoeff_Load, CRTM_AerosolCoeff_Destroy + USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type, & + CRTM_Atmosphere_Create, & + CRTM_Atmosphere_Destroy, & + CRTM_Atmosphere_Zero + USE CRTM_AtmOptics_Define, ONLY: CRTM_AtmOptics_type, & + CRTM_AtmOptics_Create, & + CRTM_AtmOptics_Destroy, & + CRTM_AtmOptics_Zero + USE CRTM_AerosolScatter, ONLY: CRTM_Compute_AerosolScatter, & + CRTM_Compute_AerosolScatter_TL, & + CRTM_Compute_AerosolScatter_AD + USE ASvar_Define, ONLY: ASvar_type, & + ASvar_Create, ASvar_Destroy + USE CRTM_ChannelInfo_Define, ONLY: CRTM_ChannelInfo_type + + IMPLICIT NONE + + CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'test_AerosolScatter_AD' + INTEGER :: Error_Status + CHARACTER(256) :: Sensor_Id = 'modis_aqua' + CHARACTER(256) :: File_Path = './testinput/' + + TYPE(CRTM_Atmosphere_type) :: Atm, Atm_TL, Atm_AD + TYPE(CRTM_AtmOptics_type) :: AScat, AScat_TL, AScat_AD + TYPE(ASvar_type) :: ASvar + INTEGER :: SensorIndex = 1 + INTEGER :: ChannelIndex = 1 + INTEGER :: n_Layers = 1 + INTEGER :: n_Aerosols = 1 + INTEGER :: n_Legendre = 4 + INTEGER :: n_Phase = 1 + + REAL(fp) :: sum_tl, sum_ad + INTEGER :: k, n, l, m + + PRINT *, 'Starting test_AerosolScatter_AD...' + + ! Load Coefficients + Error_Status = CRTM_SpcCoeff_Load( (/Sensor_Id/), File_Path=File_Path ) + IF ( Error_Status /= SUCCESS ) STOP 1 + + Error_Status = CRTM_AerosolCoeff_Load('GOCART-GEOS5', 'AerosolCoeff.GOCART-GEOS5.bin', File_Path=File_Path) + IF ( Error_Status /= SUCCESS ) STOP 1 + + ! Create structures + CALL CRTM_Atmosphere_Create( Atm, n_Layers, 1, 0, n_Aerosols ) + CALL CRTM_Atmosphere_Create( Atm_TL, n_Layers, 1, 0, n_Aerosols ) + CALL CRTM_Atmosphere_Create( Atm_AD, n_Layers, 1, 0, n_Aerosols ) + CALL CRTM_AtmOptics_Create( AScat, n_Layers, n_Legendre, n_Phase ) + CALL CRTM_AtmOptics_Create( AScat_TL, n_Layers, n_Legendre, n_Phase ) + CALL CRTM_AtmOptics_Create( AScat_AD, n_Layers, n_Legendre, n_Phase ) + CALL ASvar_Create( ASvar, n_Legendre, n_Phase, n_Layers, n_Aerosols ) + + ! Populate Input + Atm%Absorber_ID(1) = 1 + Atm%Aerosol(1)%Type = 6 ! Sea Salt + Atm%Aerosol(1)%Concentration(1) = 1.0e-4_fp + Atm%Aerosol(1)%Effective_Radius(1) = 1.0_fp + Atm%Relative_Humidity(1) = 0.5_fp + + ! Populate TL Input: Perturb Relative Humidity and Concentration + CALL CRTM_Atmosphere_Zero( Atm_TL ) + Atm_TL%Relative_Humidity(1) = 0.1_fp + Atm_TL%Aerosol(1)%Concentration(1) = 1.0e-5_fp + + AScat%Include_Scattering = .TRUE. + AScat_TL%Include_Scattering = .TRUE. + AScat_AD%Include_Scattering = .TRUE. + + ! Forward Call + Error_Status = CRTM_Compute_AerosolScatter( Atm, SensorIndex, ChannelIndex, AScat, ASvar ) + IF ( Error_Status /= SUCCESS ) STOP 1 + + ! TL Call + CALL CRTM_AtmOptics_Zero( AScat_TL ) + Error_Status = CRTM_Compute_AerosolScatter_TL( Atm, AScat, Atm_TL, SensorIndex, ChannelIndex, AScat_TL, ASvar ) + IF ( Error_Status /= SUCCESS ) STOP 1 + + ! Adjoint Input: Set AScat_AD = AScat_TL + CALL CRTM_AtmOptics_Zero( AScat_AD ) + AScat_AD%Optical_Depth = AScat_TL%Optical_Depth + AScat_AD%Single_Scatter_Albedo = AScat_TL%Single_Scatter_Albedo + AScat_AD%Asymmetry_Factor = AScat_TL%Asymmetry_Factor + AScat_AD%Phase_Coefficient = AScat_TL%Phase_Coefficient + + ! AD Call + CALL CRTM_Atmosphere_Zero( Atm_AD ) + Error_Status = CRTM_Compute_AerosolScatter_AD( Atm, AScat, AScat_AD, SensorIndex, ChannelIndex, Atm_AD, ASvar ) + IF ( Error_Status /= SUCCESS ) STOP 1 + + ! Adjoint Test: sum(AScat_TL * AScat_AD_input) == sum(Atm_TL * Atm_AD_output) + ! sum_tl = sum(AScat_TL**2) + sum_tl = SUM(AScat_TL%Optical_Depth**2) + & + SUM(AScat_TL%Single_Scatter_Albedo**2) + & + SUM(AScat_TL%Asymmetry_Factor**2) + & + SUM(AScat_TL%Phase_Coefficient**2) + + ! sum_ad = sum(Atm_TL * Atm_AD) + sum_ad = SUM(Atm_TL%Relative_Humidity * Atm_AD%Relative_Humidity) + & + SUM(Atm_TL%Aerosol(1)%Concentration * Atm_AD%Aerosol(1)%Concentration) + & + SUM(Atm_TL%Aerosol(1)%Effective_Radius * Atm_AD%Aerosol(1)%Effective_Radius) + & + SUM(Atm_TL%Aerosol(1)%Effective_Variance * Atm_AD%Aerosol(1)%Effective_Variance) + + PRINT *, 'LHS (TL sum): ', sum_tl + PRINT *, 'RHS (AD sum): ', sum_ad + PRINT *, 'Diff: ', ABS(sum_tl - sum_ad) + + IF ( sum_tl > ZERO ) THEN + PRINT *, 'Rel. Diff: ', ABS(sum_tl - sum_ad) / sum_tl + IF ( ABS(sum_tl - sum_ad) / sum_tl < 1.0e-12_fp ) THEN + PRINT *, 'SUCCESS: Adjoint test passed.' + ELSE + PRINT *, 'FAIL: Adjoint test failed.' + STOP 1 + END IF + ELSE + PRINT *, 'FAIL: LHS is ZERO, test is invalid.' + STOP 1 + END IF + + ! Clean up + CALL CRTM_Atmosphere_Destroy( Atm ) + CALL CRTM_Atmosphere_Destroy( Atm_TL ) + CALL CRTM_Atmosphere_Destroy( Atm_AD ) + CALL CRTM_AtmOptics_Destroy( AScat ) + CALL CRTM_AtmOptics_Destroy( AScat_TL ) + CALL CRTM_AtmOptics_Destroy( AScat_AD ) + CALL ASvar_Destroy( ASvar ) + Error_Status = CRTM_SpcCoeff_Destroy() + Error_Status = CRTM_AerosolCoeff_Destroy() + +END PROGRAM test_AerosolScatter_AD diff --git a/test/mains/unit/Unit_Test/test_AerosolScatter_K.f90 b/test/mains/unit/Unit_Test/test_AerosolScatter_K.f90 new file mode 100644 index 00000000..356da66f --- /dev/null +++ b/test/mains/unit/Unit_Test/test_AerosolScatter_K.f90 @@ -0,0 +1,144 @@ + +! +! test_AerosolScatter_K +! +! Unit test for CRTM_K_Matrix focusing on Aerosol sensitivities (Jacobians). +! + +PROGRAM test_AerosolScatter_K + + USE Type_Kinds, ONLY: fp + USE Message_Handler, ONLY: SUCCESS, FAILURE, Display_Message + USE CRTM_Parameters, ONLY: ZERO, ONE + USE CRTM_SpcCoeff, ONLY: SC, CRTM_SpcCoeff_Load, CRTM_SpcCoeff_Destroy + USE CRTM_AerosolCoeff, ONLY: AeroC, CRTM_AerosolCoeff_Load, CRTM_AerosolCoeff_Destroy + USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type, & + CRTM_Atmosphere_Create, & + CRTM_Atmosphere_Destroy, & + CRTM_Atmosphere_Zero, & + H2O_ID, O3_ID, CO2_ID, N2O_ID, CH4_ID, CO_ID, & + MASS_MIXING_RATIO_UNITS + USE CRTM_AtmOptics_Define, ONLY: CRTM_AtmOptics_type, & + CRTM_AtmOptics_Create, & + CRTM_AtmOptics_Destroy, & + CRTM_AtmOptics_Zero + USE CRTM_Surface_Define, ONLY: CRTM_Surface_type, & + CRTM_Surface_Create, & + CRTM_Surface_Destroy, & + CRTM_Surface_Zero + USE CRTM_ChannelInfo_Define, ONLY: CRTM_ChannelInfo_type + USE CRTM_Geometry_Define, ONLY: CRTM_Geometry_type + USE CRTM_RTSolution_Define, ONLY: CRTM_RTSolution_type + USE CRTM_LifeCycle, ONLY: CRTM_Init, CRTM_Destroy + USE CRTM_K_Matrix_Module, ONLY: CRTM_K_Matrix + + IMPLICIT NONE + + CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'test_AerosolScatter_K' + CHARACTER(*), PARAMETER :: COEFFICIENTS_PATH = './testinput/' + + INTEGER :: Error_Status + CHARACTER(256) :: Sensor_Id = 'modis_aqua' + + TYPE(CRTM_ChannelInfo_type) :: ChannelInfo(1) + TYPE(CRTM_Geometry_type) :: Geometry(1) + TYPE(CRTM_Atmosphere_type) :: Atm(1) + TYPE(CRTM_Atmosphere_type), ALLOCATABLE :: Atm_K(:,:) + TYPE(CRTM_Surface_type) :: Sfc(1) + TYPE(CRTM_Surface_type), ALLOCATABLE :: Sfc_K(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution(:,:), RTSolution_K(:,:) + + INTEGER :: n_Channels + INTEGER :: n_Layers = 10 + INTEGER :: n_Absorbers = 6 + INTEGER :: n_Aerosols = 1 + INTEGER :: n_Clouds = 0 + + INTEGER :: l, k + + PRINT *, 'Starting test_AerosolScatter_K...' + + ! Initialize CRTM + Error_Status = CRTM_Init( (/Sensor_Id/), ChannelInfo, File_Path=COEFFICIENTS_PATH ) + IF ( Error_Status /= SUCCESS ) STOP 1 + + ! Load GOCART-GEOS5 + Error_Status = CRTM_AerosolCoeff_Load('GOCART-GEOS5', 'AerosolCoeff.GOCART-GEOS5.bin', File_Path=COEFFICIENTS_PATH) + IF ( Error_Status /= SUCCESS ) STOP 1 + + n_Channels = ChannelInfo(1)%n_Channels + ALLOCATE( RTSolution( n_Channels, 1 ), & + RTSolution_K( n_Channels, 1 ), & + Atm_K( n_Channels, 1 ), & + Sfc_K( n_Channels, 1 ) ) + + ! Create structures + CALL CRTM_Atmosphere_Create( Atm, n_Layers, n_Absorbers, n_Clouds, n_Aerosols ) + CALL CRTM_Atmosphere_Create( Atm_K, n_Layers, n_Absorbers, n_Clouds, n_Aerosols ) + CALL CRTM_Surface_Create( Sfc_K, 0 ) ! Allocate Surface_K too! + + ! Populate Input + Atm(1)%Absorber_ID(1) = H2O_ID + Atm(1)%Absorber_ID(2) = O3_ID + Atm(1)%Absorber_ID(3) = CO2_ID + Atm(1)%Absorber_ID(4) = N2O_ID + Atm(1)%Absorber_ID(5) = CH4_ID + Atm(1)%Absorber_ID(6) = CO_ID + Atm(1)%Aerosol(1)%Type = 6 ! Sea Salt (hygroscopic) + + ! Standard Atmosphere (simplified) + Atm(1)%Pressure(1:n_Layers) = (/ (1000.0_fp - REAL(k-1,fp)*100.0_fp, k=1,n_Layers) /) + Atm(1)%Temperature(1:n_Layers) = 280.0_fp + Atm(1)%Absorber_Units = MASS_MIXING_RATIO_UNITS + Atm(1)%Absorber(1:n_Layers, 1) = 1.0e-3_fp ! H2O + Atm(1)%Absorber(1:n_Layers, 2) = 1.0e-6_fp ! O3 + Atm(1)%Absorber(1:n_Layers, 3:6) = 1.0e-7_fp + Atm(1)%Aerosol(1)%Concentration(1:n_Layers) = 1.0e-5_fp + Atm(1)%Aerosol(1)%Effective_Radius(1:n_Layers) = 1.0_fp + Atm(1)%Relative_Humidity(1:n_Layers) = 0.5_fp ! 50% as fraction + + ! Geometry + Geometry(1)%Sensor_Zenith_Angle = 30.0_fp + + ! Surface + CALL CRTM_Surface_Create( Sfc, 0 ) ! Land + Sfc(1)%Land_Coverage = 1.0_fp + Sfc(1)%Land_Temperature = 290.0_fp + + ! Initialize K-Matrix inputs: Sensitivity to Brightness Temperature + CALL CRTM_Atmosphere_Zero( Atm_K ) + CALL CRTM_Surface_Zero( Sfc_K ) + DO l = 1, n_Channels + RTSolution_K(l,1)%Radiance = ZERO + RTSolution_K(l,1)%Brightness_Temperature = ONE + END DO + + ! Call K-Matrix + Error_Status = CRTM_K_Matrix( Atm, Sfc, RTSolution_K, Geometry, ChannelInfo, Atm_K, Sfc_K, RTSolution ) + IF ( Error_Status /= SUCCESS ) STOP 1 + + ! Inspect Aerosol Jacobians for Channel 1 (VIS) + PRINT *, 'Aerosol Jacobians for Channel 1 (MODIS VIS 0.65um):' + PRINT *, 'Layer, RH_Jacobian, Concentration_Jacobian' + DO k = 1, n_Layers + WRITE(*, '(i3, 2(1x,es12.5))') k, Atm_K(1,1)%Relative_Humidity(k), Atm_K(1,1)%Aerosol(1)%Concentration(k) + END DO + + ! Check if RH Jacobian is non-zero + IF ( ANY(ABS(Atm_K(1,1)%Relative_Humidity) > 1.0e-12_fp) ) THEN + PRINT *, 'SUCCESS: Aerosol RH Jacobians are NON-ZERO.' + ELSE + PRINT *, 'FAIL: Aerosol RH Jacobians are ZERO.' + STOP 1 + END IF + + ! Clean up + CALL CRTM_Atmosphere_Destroy( Atm ) + CALL CRTM_Atmosphere_Destroy( Atm_K ) + CALL CRTM_Surface_Destroy( Sfc ) + CALL CRTM_Surface_Destroy( Sfc_K ) + DEALLOCATE( RTSolution, RTSolution_K ) + Error_Status = CRTM_Destroy( ChannelInfo ) + Error_Status = CRTM_AerosolCoeff_Destroy() + +END PROGRAM test_AerosolScatter_K diff --git a/test/mains/unit/Unit_Test/test_AerosolScatter_TL.f90 b/test/mains/unit/Unit_Test/test_AerosolScatter_TL.f90 new file mode 100644 index 00000000..929457de --- /dev/null +++ b/test/mains/unit/Unit_Test/test_AerosolScatter_TL.f90 @@ -0,0 +1,125 @@ + +! +! test_AerosolScatter_TL +! +! Unit test for CRTM_Compute_AerosolScatter_TL +! Verifies sensitivities to input variables. +! + +PROGRAM test_AerosolScatter_TL + + USE Type_Kinds, ONLY: fp + USE Message_Handler, ONLY: SUCCESS, FAILURE, Display_Message + USE CRTM_Parameters, ONLY: ZERO, ONE + USE CRTM_SpcCoeff, ONLY: SC, CRTM_SpcCoeff_Load, CRTM_SpcCoeff_Destroy + USE CRTM_AerosolCoeff, ONLY: AeroC, CRTM_AerosolCoeff_Load, CRTM_AerosolCoeff_Destroy + USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type, & + CRTM_Atmosphere_Create, & + CRTM_Atmosphere_Destroy, & + CRTM_Atmosphere_Zero + USE CRTM_AtmOptics_Define, ONLY: CRTM_AtmOptics_type, & + CRTM_AtmOptics_Create, & + CRTM_AtmOptics_Destroy, & + CRTM_AtmOptics_Zero + USE CRTM_AerosolScatter, ONLY: CRTM_Compute_AerosolScatter, & + CRTM_Compute_AerosolScatter_TL + USE ASvar_Define, ONLY: ASvar_type, & + ASvar_Create, ASvar_Destroy + USE CRTM_ChannelInfo_Define, ONLY: CRTM_ChannelInfo_type + + IMPLICIT NONE + + CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'test_AerosolScatter_TL' + INTEGER :: Error_Status + CHARACTER(256) :: Sensor_Id = 'modis_aqua' + CHARACTER(256) :: File_Path = './testinput/' + + TYPE(CRTM_Atmosphere_type) :: Atm, Atm_TL + TYPE(CRTM_AtmOptics_type) :: AScat, AScat_TL + TYPE(ASvar_type) :: ASvar + INTEGER :: SensorIndex = 1 + INTEGER :: ChannelIndex = 1 + INTEGER :: n_Layers = 1 + INTEGER :: n_Aerosols = 1 + INTEGER :: n_Legendre = 4 + INTEGER :: n_Phase = 1 + + PRINT *, 'Starting test_AerosolScatter_TL...' + + ! Load Coefficients + Error_Status = CRTM_SpcCoeff_Load( (/Sensor_Id/), File_Path=File_Path ) + IF ( Error_Status /= SUCCESS ) STOP 1 + + ! Using GOCART scheme for RH sensitivity + Error_Status = CRTM_AerosolCoeff_Load('GOCART-GEOS5', 'AerosolCoeff.GOCART-GEOS5.bin', File_Path=File_Path) + IF ( Error_Status /= SUCCESS ) STOP 1 + + ! Create structures + CALL CRTM_Atmosphere_Create( Atm, n_Layers, 1, 0, n_Aerosols ) + CALL CRTM_Atmosphere_Create( Atm_TL, n_Layers, 1, 0, n_Aerosols ) + CALL CRTM_AtmOptics_Create( AScat, n_Layers, n_Legendre, n_Phase ) + CALL CRTM_AtmOptics_Create( AScat_TL, n_Layers, n_Legendre, n_Phase ) + CALL ASvar_Create( ASvar, n_Legendre, n_Phase, n_Layers, n_Aerosols ) + + ! Populate Input + Atm%Absorber_ID(1) = 1 + Atm%Aerosol(1)%Type = 6 + Atm%Aerosol(1)%Concentration(1) = 1.0e-4_fp + Atm%Aerosol(1)%Effective_Radius(1) = 1.0_fp + Atm%Relative_Humidity(1) = 0.5_fp + + ! Populate TL Input: Perturb Relative Humidity + CALL CRTM_Atmosphere_Zero( Atm_TL ) + Atm_TL%Relative_Humidity(1) = 1.0_fp + + AScat%Include_Scattering = .TRUE. + AScat_TL%Include_Scattering = .TRUE. + + ! Forward Call to populate ASvar + Error_Status = CRTM_Compute_AerosolScatter( Atm, SensorIndex, ChannelIndex, AScat, ASvar ) + IF ( Error_Status /= SUCCESS ) STOP 1 + + ! TL Call for RH + CALL CRTM_AtmOptics_Zero( AScat_TL ) + Error_Status = CRTM_Compute_AerosolScatter_TL( Atm, AScat, Atm_TL, SensorIndex, ChannelIndex, AScat_TL, ASvar ) + IF ( Error_Status /= SUCCESS ) STOP 1 + PRINT *, 'AScat_TL%Optical_Depth(1) for RH perturbation: ', AScat_TL%Optical_Depth(1) + IF ( ABS(AScat_TL%Optical_Depth(1)) < 1.0e-12_fp ) THEN + PRINT *, 'FAIL: Sensitivity to RH is ZERO' + ELSE + PRINT *, 'SUCCESS: Sensitivity to RH is NON-ZERO' + END IF + + ! TL Call for Effective Radius (Use a scheme that supports it, e.g., CRTM or CMAQ) + ! We'll stay with GOCART for now but just check if it's zero or non-zero as expected. + ! GOCART-GEOS5 usually has N_RADII=1, so it might be zero. + CALL CRTM_Atmosphere_Zero( Atm_TL ) + Atm_TL%Aerosol(1)%Effective_Radius(1) = 1.0_fp + CALL CRTM_AtmOptics_Zero( AScat_TL ) + Error_Status = CRTM_Compute_AerosolScatter_TL( Atm, AScat, Atm_TL, SensorIndex, ChannelIndex, AScat_TL, ASvar ) + IF ( Error_Status /= SUCCESS ) STOP 1 + PRINT *, 'AScat_TL%Optical_Depth(1) for Reff perturbation: ', AScat_TL%Optical_Depth(1) + + ! TL Call for Concentration + CALL CRTM_Atmosphere_Zero( Atm_TL ) + Atm_TL%Aerosol(1)%Concentration(1) = 1.0_fp + CALL CRTM_AtmOptics_Zero( AScat_TL ) + Error_Status = CRTM_Compute_AerosolScatter_TL( Atm, AScat, Atm_TL, SensorIndex, ChannelIndex, AScat_TL, ASvar ) + IF ( Error_Status /= SUCCESS ) STOP 1 + PRINT *, 'AScat_TL%Optical_Depth(1) for Concentration perturbation: ', AScat_TL%Optical_Depth(1) + IF ( ABS(AScat_TL%Optical_Depth(1)) < 1.0e-12_fp ) THEN + PRINT *, 'FAIL: Sensitivity to Concentration is ZERO' + ELSE + PRINT *, 'SUCCESS: Sensitivity to Concentration is NON-ZERO' + END IF + + ! Clean up + CALL CRTM_Atmosphere_Destroy( Atm ) + CALL CRTM_Atmosphere_Destroy( Atm_TL ) + CALL CRTM_AtmOptics_Destroy( AScat ) + CALL CRTM_AtmOptics_Destroy( AScat_TL ) + CALL ASvar_Destroy( ASvar ) + Error_Status = CRTM_SpcCoeff_Destroy() + Error_Status = CRTM_AerosolCoeff_Destroy() + +END PROGRAM test_AerosolScatter_TL