From 734fc832ee73665f847e7fc9d99d1d02af5b8afb Mon Sep 17 00:00:00 2001 From: chengdang Date: Tue, 25 Mar 2025 11:30:58 -0600 Subject: [PATCH 1/5] Cloudy sky bug fix --- src/RTSolution/Common_RTSolution.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/RTSolution/Common_RTSolution.f90 b/src/RTSolution/Common_RTSolution.f90 index 79486cc7..e224f82f 100644 --- a/src/RTSolution/Common_RTSolution.f90 +++ b/src/RTSolution/Common_RTSolution.f90 @@ -1750,8 +1750,10 @@ FUNCTION Assign_Common_Output_AD( & ELSE IF( RTV%Scattering_RT ) THEN - User_Emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) + !User_Emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) IF( RTV%Diffuse_Surface) THEN + ! bugfix test: CD + User_Emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) IF( RTV%mth_Azi == 0 ) THEN ! Assuming Lambertian surface isn't polarized! DO i = SfcOptics%n_Angles, 1, -1 From 137f7f184dced1e70f33ab21fae95b980cdc4922 Mon Sep 17 00:00:00 2001 From: chengdang Date: Tue, 25 Mar 2025 11:50:39 -0600 Subject: [PATCH 2/5] test for v2.4 implementation --- src/RTSolution/Common_RTSolution.f90 | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/RTSolution/Common_RTSolution.f90 b/src/RTSolution/Common_RTSolution.f90 index e224f82f..d55fcb57 100644 --- a/src/RTSolution/Common_RTSolution.f90 +++ b/src/RTSolution/Common_RTSolution.f90 @@ -1750,16 +1750,20 @@ FUNCTION Assign_Common_Output_AD( & ELSE IF( RTV%Scattering_RT ) THEN + !User_Emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) + !bugfix test: CD option 1, follow CRTM v2.4 + User_Emissivity_AD = ZERO + IF( RTV%Diffuse_Surface) THEN - ! bugfix test: CD - User_Emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) + ! bugfix test: CD, option 2, seems to have some effects on ocean though unclear why + ! User_Emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) IF( RTV%mth_Azi == 0 ) THEN - ! Assuming Lambertian surface isn't polarized! - DO i = SfcOptics%n_Angles, 1, -1 - User_Emissivity_AD = User_Emissivity_AD - & - (SUM(SfcOptics_AD%Reflectivity(1:SfcOptics%n_Angles,1,i,1))*SfcOptics%Weight(i)) - END DO + ! Assuming Lambertian surface isn't polarized! + DO i = SfcOptics%n_Angles, 1, -1 + User_Emissivity_AD = User_Emissivity_AD - & + (SUM(SfcOptics_AD%Reflectivity(1:SfcOptics%n_Angles,1,i,1))*SfcOptics%Weight(i)) + END DO ELSE SfcOptics_AD%Reflectivity(1:SfcOptics%n_Angles, 1, 1:SfcOptics%n_Angles, 1) = ZERO SfcOptics_AD%Direct_Reflectivity(1:SfcOptics%n_Angles,1) = ZERO @@ -1767,9 +1771,9 @@ FUNCTION Assign_Common_Output_AD( & END IF ELSE ! Specular surface - DO i = nZ, 1, -1 - User_Emissivity_AD = User_Emissivity_AD - SfcOptics_AD%Reflectivity(i,1,i,1) - END DO + DO i = nZ, 1, -1 + User_Emissivity_AD = User_Emissivity_AD - SfcOptics_AD%Reflectivity(i,1,i,1) + END DO END IF ! Direct_Reflectivity_AD = SUM(SfcOptics_AD%Direct_Reflectivity(1:nZ,1)) ! SfcOptics_AD%Direct_Reflectivity(1,1) = SfcOptics_AD%Direct_Reflectivity(1,1) + From b6df808a82bd039f7596193f890469ab28d73a3e Mon Sep 17 00:00:00 2001 From: chengdang Date: Tue, 25 Mar 2025 16:05:23 -0600 Subject: [PATCH 3/5] Bug fix for surface emissivity jacobian Minor modification to CRTM test: including amusa test over water covered surface --- src/RTSolution/Common_RTSolution.f90 | 19 ++++---------- test/CMakeLists.txt | 1 + .../k_matrix/test_Simple/Load_Sfc_Data.inc | 26 +++++++++---------- .../k_matrix/test_Simple/test_Simple.f90 | 3 ++- 4 files changed, 21 insertions(+), 28 deletions(-) diff --git a/src/RTSolution/Common_RTSolution.f90 b/src/RTSolution/Common_RTSolution.f90 index d55fcb57..a6fdcf45 100644 --- a/src/RTSolution/Common_RTSolution.f90 +++ b/src/RTSolution/Common_RTSolution.f90 @@ -59,8 +59,8 @@ MODULE Common_RTSolution ! Module parameters ! ----------------- ! Version Id for the module - CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = & - '$Id: $' + ! CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = & + ! '$Id: $' CONTAINS @@ -1723,12 +1723,9 @@ FUNCTION Assign_Common_Output_AD( & ! ----------------------------------------- ! Populate the SfcOptics_AD structure based ! on FORWARD model SfcOptics Compute_Switch - ! ----------------------------------------- - + ! ----------------------------------------- IF ( SfcOptics%Compute ) THEN - RTSolution_AD%Surface_Emissivity = SfcOptics_AD%Emissivity(SfcOptics_AD%Index_Sat_Ang,1) - Error_Status = CRTM_Compute_SfcOptics_AD( & Surface , & ! Input SfcOptics , & ! Input @@ -1750,14 +1747,9 @@ FUNCTION Assign_Common_Output_AD( & ELSE IF( RTV%Scattering_RT ) THEN - - !User_Emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) - !bugfix test: CD option 1, follow CRTM v2.4 - User_Emissivity_AD = ZERO - + ! User_Emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) ! V3.0 implementation + User_Emissivity_AD = ZERO ! V2.4 implementation IF( RTV%Diffuse_Surface) THEN - ! bugfix test: CD, option 2, seems to have some effects on ocean though unclear why - ! User_Emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) IF( RTV%mth_Azi == 0 ) THEN ! Assuming Lambertian surface isn't polarized! DO i = SfcOptics%n_Angles, 1, -1 @@ -1769,7 +1761,6 @@ FUNCTION Assign_Common_Output_AD( & SfcOptics_AD%Direct_Reflectivity(1:SfcOptics%n_Angles,1) = ZERO SfcOptics_AD%Emissivity(1:SfcOptics%n_Angles,1) = ZERO END IF - ELSE ! Specular surface DO i = nZ, 1, -1 User_Emissivity_AD = User_Emissivity_AD - SfcOptics_AD%Reflectivity(i,1,i,1) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c1a190a0..a7e58005 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -125,6 +125,7 @@ list( APPEND Simple_Sensor_Ids v.abi_gr abi_g18 modis_aqua + amsua_n19 ) list( APPEND ScatteringSwitch_Sensor_Ids diff --git a/test/mains/regression/k_matrix/test_Simple/Load_Sfc_Data.inc b/test/mains/regression/k_matrix/test_Simple/Load_Sfc_Data.inc index 3b2aec4a..f969f884 100644 --- a/test/mains/regression/k_matrix/test_Simple/Load_Sfc_Data.inc +++ b/test/mains/regression/k_matrix/test_Simple/Load_Sfc_Data.inc @@ -21,24 +21,24 @@ ! 4a.1 Profile #1 ! --------------- ! ...Land surface characteristics - sfc(1)%Land_Coverage = 0.1_fp - sfc(1)%Land_Type = TUNDRA_SURFACE_TYPE - sfc(1)%Land_Temperature = 272.0_fp - sfc(1)%Lai = 0.17_fp - sfc(1)%Soil_Type = COARSE_SOIL_TYPE - sfc(1)%Vegetation_Type = GROUNDCOVER_VEGETATION_TYPE + !sfc(1)%Land_Coverage = 0.1_fp + !sfc(1)%Land_Type = TUNDRA_SURFACE_TYPE + !sfc(1)%Land_Temperature = 272.0_fp + !sfc(1)%Lai = 0.17_fp + !sfc(1)%Soil_Type = COARSE_SOIL_TYPE + !sfc(1)%Vegetation_Type = GROUNDCOVER_VEGETATION_TYPE ! ...Water surface characteristics - sfc(1)%Water_Coverage = 0.5_fp + sfc(1)%Water_Coverage = 1.0_fp sfc(1)%Water_Type = SEA_WATER_TYPE sfc(1)%Water_Temperature = 275.0_fp ! ...Snow coverage characteristics - sfc(1)%Snow_Coverage = 0.25_fp - sfc(1)%Snow_Type = FRESH_SNOW_TYPE - sfc(1)%Snow_Temperature = 265.0_fp + !sfc(1)%Snow_Coverage = 0.25_fp + !sfc(1)%Snow_Type = FRESH_SNOW_TYPE + !sfc(1)%Snow_Temperature = 265.0_fp ! ...Ice surface characteristics - sfc(1)%Ice_Coverage = 0.15_fp - sfc(1)%Ice_Type = FRESH_ICE_TYPE - sfc(1)%Ice_Temperature = 269.0_fp + !sfc(1)%Ice_Coverage = 0.15_fp + !sfc(1)%Ice_Type = FRESH_ICE_TYPE + !sfc(1)%Ice_Temperature = 269.0_fp diff --git a/test/mains/regression/k_matrix/test_Simple/test_Simple.f90 b/test/mains/regression/k_matrix/test_Simple/test_Simple.f90 index e79beb13..1ee2942a 100644 --- a/test/mains/regression/k_matrix/test_Simple/test_Simple.f90 +++ b/test/mains/regression/k_matrix/test_Simple/test_Simple.f90 @@ -160,7 +160,7 @@ PROGRAM test_Simple Message = 'Error allocating CRTM Atmosphere_K structure' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) STOP 1 - END IF + END IF ! Deleted in V2.4.1 ! The comparative K-MATRIX structure inside the results file !CALL CRTM_Atmosphere_Create( atm_K, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) @@ -257,6 +257,7 @@ PROGRAM test_Simple CALL CRTM_RTSolution_Inspect(RTSolution(l,m)) ! K-MATRIX output WRITE( *, '(/3x,"K-MATRIX OUTPUT")') + CALL CRTM_RTSolution_Inspect(RTSolution_K(l,m)) CALL CRTM_Surface_Inspect(Surface_K(l,m)) CALL CRTM_Atmosphere_Inspect(Atmosphere_K(l,m)) END DO From 964e34ef6f69466621c8e6fdae7fe8bd74cdd185 Mon Sep 17 00:00:00 2001 From: chengdang Date: Wed, 2 Apr 2025 15:37:20 -0600 Subject: [PATCH 4/5] Revert changes to Load_Sfc_Data.inc --- .../k_matrix/test_Simple/Load_Sfc_Data.inc | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/test/mains/regression/k_matrix/test_Simple/Load_Sfc_Data.inc b/test/mains/regression/k_matrix/test_Simple/Load_Sfc_Data.inc index f969f884..3b2aec4a 100644 --- a/test/mains/regression/k_matrix/test_Simple/Load_Sfc_Data.inc +++ b/test/mains/regression/k_matrix/test_Simple/Load_Sfc_Data.inc @@ -21,24 +21,24 @@ ! 4a.1 Profile #1 ! --------------- ! ...Land surface characteristics - !sfc(1)%Land_Coverage = 0.1_fp - !sfc(1)%Land_Type = TUNDRA_SURFACE_TYPE - !sfc(1)%Land_Temperature = 272.0_fp - !sfc(1)%Lai = 0.17_fp - !sfc(1)%Soil_Type = COARSE_SOIL_TYPE - !sfc(1)%Vegetation_Type = GROUNDCOVER_VEGETATION_TYPE + sfc(1)%Land_Coverage = 0.1_fp + sfc(1)%Land_Type = TUNDRA_SURFACE_TYPE + sfc(1)%Land_Temperature = 272.0_fp + sfc(1)%Lai = 0.17_fp + sfc(1)%Soil_Type = COARSE_SOIL_TYPE + sfc(1)%Vegetation_Type = GROUNDCOVER_VEGETATION_TYPE ! ...Water surface characteristics - sfc(1)%Water_Coverage = 1.0_fp + sfc(1)%Water_Coverage = 0.5_fp sfc(1)%Water_Type = SEA_WATER_TYPE sfc(1)%Water_Temperature = 275.0_fp ! ...Snow coverage characteristics - !sfc(1)%Snow_Coverage = 0.25_fp - !sfc(1)%Snow_Type = FRESH_SNOW_TYPE - !sfc(1)%Snow_Temperature = 265.0_fp + sfc(1)%Snow_Coverage = 0.25_fp + sfc(1)%Snow_Type = FRESH_SNOW_TYPE + sfc(1)%Snow_Temperature = 265.0_fp ! ...Ice surface characteristics - !sfc(1)%Ice_Coverage = 0.15_fp - !sfc(1)%Ice_Type = FRESH_ICE_TYPE - !sfc(1)%Ice_Temperature = 269.0_fp + sfc(1)%Ice_Coverage = 0.15_fp + sfc(1)%Ice_Type = FRESH_ICE_TYPE + sfc(1)%Ice_Temperature = 269.0_fp From 677afc8457e8bd4c502a61371db99d6a88c23f77 Mon Sep 17 00:00:00 2001 From: chengdang Date: Mon, 7 Apr 2025 15:14:55 -0600 Subject: [PATCH 5/5] Fix code format for CRTM_K_Matrix_Module.f90 --- src/CRTM_K_Matrix_Module.f90 | 588 +++++++++++++++++------------------ 1 file changed, 293 insertions(+), 295 deletions(-) diff --git a/src/CRTM_K_Matrix_Module.f90 b/src/CRTM_K_Matrix_Module.f90 index 86be339d..a71b45bd 100644 --- a/src/CRTM_K_Matrix_Module.f90 +++ b/src/CRTM_K_Matrix_Module.f90 @@ -130,7 +130,7 @@ MODULE CRTM_K_Matrix_Module USE CRTM_Planck_Functions, ONLY: CRTM_Planck_Temperature , & CRTM_Planck_Temperature_AD USE CRTM_CloudCover_Define, ONLY: CRTM_CloudCover_type - USE CRTM_Active_Sensor, ONLY: CRTM_Compute_Reflectivity, & + USE CRTM_Active_Sensor, ONLY: CRTM_Compute_Reflectivity , & CRTM_Compute_Reflectivity_AD, & Calculate_Cloud_Water_Density @@ -449,26 +449,25 @@ FUNCTION CRTM_K_Matrix( & !$OMP PARALLEL DO PRIVATE (nc,Message) NUM_THREADS(n_profile_threads) Profile_Loop1: DO m = 1, n_Profiles ! Check the cloud and aerosol coeff. data for cases with clouds and aerosol - IF ( Atmosphere(m)%n_Clouds > 0) THEN - !** clear clouds where cloud_fraction < threshold - DO nc = 1, Atmosphere(m)%n_clouds - WHERE (Atmosphere(m)%Cloud_Fraction(:) < MIN_COVERAGE_THRESHOLD) - Atmosphere(m)%Cloud_Fraction(:) = ZERO - Atmosphere(m)%Cloud(nc)%Water_Content(:) = ZERO - Atmosphere(m)%Cloud(nc)%Effective_Radius(:) = ZERO - END WHERE - END DO - - - IF(.NOT. CRTM_CloudCoeff_IsLoaded() )THEN - - Error_Status = FAILURE - WRITE( Message,'("The CloudCoeff data must be loaded (with CRTM_Init routine) ", & - &"for the cloudy case profile #",i0)' ) m - CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) - CYCLE Profile_Loop1 - END IF - END IF + ! ...cloud + IF ( Atmosphere(m)%n_Clouds > 0 ) THEN + !** clear clouds where cloud_fraction < threshold + DO nc = 1, Atmosphere(m)%n_clouds + WHERE (Atmosphere(m)%Cloud_Fraction(:) < MIN_COVERAGE_THRESHOLD) + Atmosphere(m)%Cloud_Fraction(:) = ZERO + Atmosphere(m)%Cloud(nc)%Water_Content(:) = ZERO + Atmosphere(m)%Cloud(nc)%Effective_Radius(:) = ZERO + END WHERE + END DO + IF( .NOT. CRTM_CloudCoeff_IsLoaded() )THEN + Error_Status = FAILURE + WRITE( Message,'("The CloudCoeff data must be loaded (with CRTM_Init routine) ", & + &"for the cloudy case profile #",i0)' ) m + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + CYCLE Profile_Loop1 + END IF + END IF + ! ...aerosol IF( Atmosphere(m)%n_Aerosols > 0 .AND. .NOT. CRTM_AerosolCoeff_IsLoaded() )THEN Error_Status = FAILURE WRITE( Message,'("The AerosolCoeff data must be loaded (with CRTM_Init routine) ", & @@ -599,7 +598,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) DO nt = 1, n_channel_threads CALL CRTM_SfcOptics_Create( SfcOptics(nt) , MAX_N_ANGLES, MAX_N_STOKES ) CALL CRTM_SfcOptics_Create( SfcOptics_K(nt), MAX_N_ANGLES, MAX_N_STOKES ) - IF ( (.NOT. CRTM_SfcOptics_Associated(SfcOptics(nt) )) .OR. & + IF ( (.NOT. CRTM_SfcOptics_Associated(SfcOptics(nt)) ) .OR. & (.NOT. CRTM_SfcOptics_Associated(SfcOptics_K(nt))) ) THEN Error_Status = FAILURE Message = 'Error allocating SfcOptics data structures' @@ -610,7 +609,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) SfcOptics_K(nt)%n_Stokes = RTV(nt)%n_Stokes END DO !$OMP END PARALLEL DO - IF ( Error_Status == FAILURE) RETURN + IF ( Error_Status == FAILURE ) RETURN ! Copy over forward "non-variable" inputs to K-matrix outputs !$OMP PARALLEL DO NUM_THREADS(n_channel_threads) @@ -710,28 +709,28 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Check n_Stokes and number of phase elements IF ( CRTM_CloudCoeff_IsLoaded() .AND. & (RTV(1)%n_Stokes > 1 .AND. CloudC%N_PHASE_ELEMENTS < 6 )) THEN - Error_Status = FAILURE - WRITE( Message,'("N_PHASE_ELEMENTS OF CLOUD LUT NOT RIGHT ",i0)' ) CloudC%N_PHASE_ELEMENTS - CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) - RETURN + Error_Status = FAILURE + WRITE( Message,'("N_PHASE_ELEMENTS OF CLOUD LUT NOT RIGHT ",i0)' ) CloudC%N_PHASE_ELEMENTS + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN END IF IF ( CRTM_AerosolCoeff_IsLoaded() .AND. & (RTV(1)%n_Stokes > 1 .AND. AeroC%N_PHASE_ELEMENTS < 6 )) THEN - Error_Status = FAILURE - WRITE( Message,'("N_PHASE_ELEMENTS OF AEROSOL LUT NOT RIGHT ",i0)' ) AeroC%N_PHASE_ELEMENTS - CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) - RETURN + Error_Status = FAILURE + WRITE( Message,'("N_PHASE_ELEMENTS OF AEROSOL LUT NOT RIGHT ",i0)' ) AeroC%N_PHASE_ELEMENTS + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN END IF IF ( CRTM_CloudCoeff_IsLoaded() .AND. CRTM_AerosolCoeff_IsLoaded() .AND. & (CloudC%N_PHASE_ELEMENTS /= AeroC%N_PHASE_ELEMENTS) ) THEN - Error_Status = FAILURE - WRITE( Message,'("N_PHASE_ELEMENTS OF CLOUD AND AEROSOL LUTS DO NOT MATCH")' ) - CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) - RETURN + Error_Status = FAILURE + WRITE( Message,'("N_PHASE_ELEMENTS OF CLOUD AND AEROSOL LUTS DO NOT MATCH")' ) + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN END IF - + ! Calculate cloud water density CALL Calculate_Cloud_Water_Density(Atm) @@ -739,24 +738,24 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! ...Allocate the atmospheric optics structures based on Atm extension !$OMP PARALLEL DO NUM_THREADS(n_channel_threads) PRIVATE(Message) DO nt = 1, n_channel_threads - CALL CRTM_AtmOptics_Create( AtmOptics(nt), & + CALL CRTM_AtmOptics_Create( AtmOptics(nt) , & Atm%n_Layers , & MAX_N_LEGENDRE_TERMS, & CloudC%N_PHASE_ELEMENTS ) - CALL CRTM_AtmOptics_Create( AtmOptics_K(nt), & + CALL CRTM_AtmOptics_Create( AtmOptics_K(nt) , & Atm%n_Layers , & MAX_N_LEGENDRE_TERMS, & CloudC%N_PHASE_ELEMENTS ) - IF ( Options_Present ) THEN - AtmOptics(nt)%depolarization = Opt%depolarization - AtmOptics_K(nt)%depolarization = Opt%depolarization - IF( Opt%n_Stokes > 0 ) RTV(nt)%n_Stokes = Opt%n_Stokes - AtmOptics(nt)%n_Stokes = RTV(nt)%n_Stokes - AtmOptics_K(nt)%n_Stokes = RTV(nt)%n_Stokes - END IF + IF ( Options_Present ) THEN + AtmOptics(nt)%depolarization = Opt%depolarization + AtmOptics_K(nt)%depolarization = Opt%depolarization + IF( Opt%n_Stokes > 0 ) RTV(nt)%n_Stokes = Opt%n_Stokes + AtmOptics(nt)%n_Stokes = RTV(nt)%n_Stokes + AtmOptics_K(nt)%n_Stokes = RTV(nt)%n_Stokes + END IF - IF ( .NOT. CRTM_AtmOptics_Associated( Atmoptics(nt) ) .OR. & + IF ( .NOT. CRTM_AtmOptics_Associated( Atmoptics(nt) ) .OR. & .NOT. CRTM_AtmOptics_Associated( Atmoptics_K(nt) ) ) THEN Error_Status = FAILURE WRITE( Message,'("Error allocating AtmOptics data structures for profile #",i0)' ) m @@ -773,23 +772,23 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Allocate the scattering internal variables if necessary ! ...Cloud IF ( Atm%n_Clouds > 0 ) THEN - CALL CSvar_Create( CSvar(nt), & - MAX_N_LEGENDRE_TERMS, & + CALL CSvar_Create( CSvar(nt) , & + MAX_N_LEGENDRE_TERMS , & CloudC%N_PHASE_ELEMENTS, & - Atm%n_Layers , & - Atm%n_Clouds ) + Atm%n_Layers , & + Atm%n_Clouds ) END IF ! ...Aerosol IF ( Atm%n_Aerosols > 0 ) THEN - CALL ASvar_Create( ASvar(nt), & - MAX_N_LEGENDRE_TERMS, & + CALL ASvar_Create( ASvar(nt) , & + MAX_N_LEGENDRE_TERMS , & AeroC%N_PHASE_ELEMENTS, & - Atm%n_Layers , & - Atm%n_Aerosols ) + Atm%n_Layers , & + Atm%n_Aerosols ) END IF - END DO + END DO ! Channel threads !$OMP END PARALLEL DO - IF ( Error_Status == FAILURE) RETURN + IF ( Error_Status == FAILURE ) RETURN ! Determine the type of cloud coverage cloud_coverage_flag = CRTM_Atmosphere_Coverage( atm ) @@ -845,9 +844,9 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) CALL CRTM_Compute_SurfaceT( Surface(m), SfcOptics_Clear(nt) ) END DO !$OMP END PARALLEL DO - IF ( Error_Status == FAILURE) RETURN + IF ( Error_Status == FAILURE ) RETURN - END IF + END IF ! If ractional cloud coverage ! Average surface skin temperature for multi-surface types @@ -904,14 +903,14 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) Atm , & ! Input GeometryInfo , & ! Input AncillaryInput, & ! Input - Predictor(nt) , & ! Output - PVar(nt) ) ! Internal variable output + Predictor(nt) , & ! Output + PVar(nt) ) ! Internal variable output ! Allocate the AtmAbsorption predictor structures CALL CRTM_Predictor_Create( & - Predictor_K(nt) , & + Predictor_K(nt), & atm%n_Layers, & SensorIndex ) IF ( .NOT. CRTM_Predictor_Associated(Predictor_K(nt)) ) THEN @@ -925,7 +924,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) IF( ( Atm%n_Clouds > 0 .OR. & Atm%n_Aerosols > 0 .OR. & SpcCoeff_IsVisibleSensor(SC(SensorIndex)).OR.SpcCoeff_IsUltravioletSensor(SC(SensorIndex)) ) .AND. & - AtmOptics(nt)%Include_Scattering ) THEN + AtmOptics(nt)%Include_Scattering ) THEN RTV(nt)%RT_Algorithm_Id = Opt%RT_Algorithm_Id CALL RTV_Create( RTV(nt), MAX_N_ANGLES, MAX_N_LEGENDRE_TERMS, Atm%n_Layers ) IF ( .NOT. RTV_Associated(RTV(nt)) ) THEN @@ -938,9 +937,9 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Assign algorithm selector END IF - END DO + END DO ! n_channel_threads !$OMP END PARALLEL DO - IF ( Error_Status == FAILURE) RETURN + IF ( Error_Status == FAILURE ) RETURN ! Compute NLTE correction predictors IF ( Opt%Apply_NLTE_Correction ) THEN @@ -988,7 +987,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) Thread_Loop: DO nt = 1, n_channel_threads start_ch = (nt - 1) * chunk_ch + 1 - IF ( nt == n_channel_threads) THEN + IF ( nt == n_channel_threads ) THEN end_ch = n_sensor_channels ELSE end_ch = start_ch + chunk_ch - 1 @@ -1080,20 +1079,22 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) CALL CRTM_Compute_AtmAbsorption( SensorIndex , & ! Input ChannelIndex , & ! Input AncillaryInput, & ! Input - Predictor(nt) , & ! Input - AtmOptics(nt) , & ! Output - AAvar(nt) ) ! Internal variable output + Predictor(nt) , & ! Input + AtmOptics(nt) , & ! Output + AAvar(nt) ) ! Internal variable output ! Compute the molecular scattering properties ! ...Solar radiation IF( SC(SensorIndex)%Solar_Irradiance(ChannelIndex) > ZERO .AND. & - Source_ZA < MAX_SOURCE_ZENITH_ANGLE) THEN - RTV(nt)%Solar_Flag_true = .TRUE. + Source_ZA < MAX_SOURCE_ZENITH_ANGLE ) THEN + RTV(nt)%Solar_Flag_true = .TRUE. END IF + ! ...Visible channel with solar radiation IF( (SpcCoeff_IsVisibleSensor(SC(SensorIndex)).OR. & - SpcCoeff_IsUltravioletSensor(SC(SensorIndex))) .AND. RTV(nt)%Solar_Flag_true ) THEN + SpcCoeff_IsUltravioletSensor(SC(SensorIndex))) .AND. & + RTV(nt)%Solar_Flag_true ) THEN RTV(nt)%Visible_Flag_true = .TRUE. ! Two cases ! (1) If clear sky, AtmOptics(nt)%n_Legendre_Terms == 0, compute Rayleigh scattering @@ -1132,6 +1133,8 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) END IF + + ! Copy the clear-sky AtmOptics IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN Status_FWD = CRTM_AtmOptics_NoScatterCopy( AtmOptics(nt), AtmOptics_Clear(nt) ) @@ -1151,12 +1154,12 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Compute the cloud particle absorption/scattering properties IF( Atm%n_Clouds > 0 ) THEN - Error_Status = CRTM_Compute_CloudScatter( Atm , & ! Input - GeometryInfo, & ! Input - SensorIndex , & ! Input - ChannelIndex, & ! Input - AtmOptics(nt) , & ! Output - CSvar(nt) ) ! Internal variable output + Error_Status = CRTM_Compute_CloudScatter( Atm , & ! Input + GeometryInfo , & ! Input + SensorIndex , & ! Input + ChannelIndex , & ! Input + AtmOptics(nt), & ! Output + CSvar(nt) ) ! Internal variable output IF (Error_Status /= SUCCESS) THEN WRITE( Message,'("Error computing CloudScatter for ",a,& &", channel ",i0,", profile #",i0)' ) & @@ -1169,11 +1172,11 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Compute the aerosol absorption/scattering properties IF ( Atm%n_Aerosols > 0 ) THEN - Error_Status = CRTM_Compute_AerosolScatter( Atm , & ! Input - SensorIndex , & ! Input - ChannelIndex, & ! Input - AtmOptics(nt) , & ! In/Output - ASvar(nt) ) ! Internal variable output + Error_Status = CRTM_Compute_AerosolScatter( Atm , & ! Input + SensorIndex , & ! Input + ChannelIndex , & ! Input + AtmOptics(nt), & ! In/Output + ASvar(nt) ) ! Internal variable output IF ( Error_Status /= SUCCESS ) THEN WRITE( Message,'("Error computing AerosolScatter for ",a,& &", channel ",i0,", profile #",i0)' ) & @@ -1232,249 +1235,245 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! mth_Azi = 0 is for an azimuth-averaged value (IR, MW) ! ...Initialise radiance RTSolution(ln,m)%Radiance = ZERO + IF ( SpcCoeff_IsVisibleSensor( SC(SensorIndex) ) .OR. & + SpcCoeff_IsUltravioletSensor( SC(SensorIndex) ) ) THEN + IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN + RTSolution(ln,m)%Total_Cloud_Cover = CloudCover%Total_Cloud_Cover + IF(RTV(nt)%n_Stokes == 1) THEN + RTSolution_Clear_K(nt)%Radiance = (ONE - CloudCover%Total_Cloud_Cover) * RTSolution_K(ln,m)%Radiance + RTSolution_K(ln,m)%Radiance = CloudCover%Total_Cloud_Cover * RTSolution_K(ln,m)%Radiance + ELSE + DO ks = 1, RTV(nt)%n_Stokes + ! RTSolution_Clear_K(nt)%Stokes(ks)=(ONE-CloudCover%Total_Cloud_Cover)*RTSolution_K(ln,m)%Stokes(ks) + ! Clear_sky case uses Emission_Module.f90 without atmospheric scatterings right now. + IF( ks == 1 ) & + RTSolution_Clear_K(nt)%Radiance=(ONE-CloudCover%Total_Cloud_Cover)*RTSolution_K(ln,m)%Stokes(ks) + RTSolution_K(ln,m)%Stokes(ks) = CloudCover%Total_Cloud_Cover*RTSolution_K(ln,m)%Stokes(ks) + END DO + END IF + END IF + END IF - IF ( SpcCoeff_IsVisibleSensor( SC(SensorIndex) ) .OR. & - SpcCoeff_IsUltravioletSensor( SC(SensorIndex) ) ) THEN - IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN - RTSolution(ln,m)%Total_Cloud_Cover = CloudCover%Total_Cloud_Cover - IF(RTV(nt)%n_Stokes == 1) THEN - RTSolution_Clear_K(nt)%Radiance = (ONE - CloudCover%Total_Cloud_Cover) * RTSolution_K(ln,m)%Radiance - RTSolution_K(ln,m)%Radiance = CloudCover%Total_Cloud_Cover * RTSolution_K(ln,m)%Radiance - ELSE - DO ks = 1, RTV(nt)%n_Stokes -! RTSolution_Clear_K(nt)%Stokes(ks)=(ONE-CloudCover%Total_Cloud_Cover)*RTSolution_K(ln,m)%Stokes(ks) -! Clear_sky case uses Emission_Module.f90 without atmospheric scatterings right now. - IF( ks == 1 ) & - RTSolution_Clear_K(nt)%Radiance=(ONE-CloudCover%Total_Cloud_Cover)*RTSolution_K(ln,m)%Stokes(ks) - RTSolution_K(ln,m)%Stokes(ks) = CloudCover%Total_Cloud_Cover*RTSolution_K(ln,m)%Stokes(ks) - END DO - END IF - END IF - END IF - - ! -------------- - ! ...Fourier expansion over azimuth angle - - Azimuth_Fourier_Loop: DO mth_Azi = 0, RTV(nt)%n_Azi + ! -------------- + ! ...Fourier expansion over azimuth angle + + Azimuth_Fourier_Loop: DO mth_Azi = 0, RTV(nt)%n_Azi + + ! Set dependent component counters + RTV(nt)%mth_Azi = mth_Azi + SfcOptics(nt)%mth_Azi = mth_Azi + + ! Solve the forward radiative transfer problem + Error_Status = CRTM_Compute_RTSolution( & + Atm , & ! Input + Surface(m) , & ! Input + AtmOptics(nt) , & ! Input + SfcOptics(nt) , & ! Input + GeometryInfo , & ! Input + SensorIndex , & ! Input + ChannelIndex , & ! Input + RTSolution(ln,m), & ! Output + RTV(nt) ) ! Internal variable output + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'( "Error computing RTSolution for ", a, & + &", channel ", i0,", profile #",i0)' ) & + TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + CYCLE Thread_Loop + END IF - ! Set dependent component counters - RTV(nt)%mth_Azi = mth_Azi - SfcOptics(nt)%mth_Azi = mth_Azi - ! Solve the forward radiative transfer problem + ! Repeat clear sky for fractionally cloudy atmospheres + IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag).and.RTV(nt)%mth_Azi==0 ) THEN + RTV_Clear(nt)%mth_Azi = RTV(nt)%mth_Azi + SfcOptics_Clear(nt)%mth_Azi = SfcOptics(nt)%mth_Azi Error_Status = CRTM_Compute_RTSolution( & - Atm , & ! Input - Surface(m) , & ! Input - AtmOptics(nt) , & ! Input - SfcOptics(nt) , & ! Input - GeometryInfo , & ! Input - SensorIndex , & ! Input - ChannelIndex , & ! Input - RTSolution(ln,m), & ! Output - RTV(nt) ) ! Internal variable output + Atm_Clear , & ! Input + Surface(m) , & ! Input + AtmOptics_Clear(nt) , & ! Input + SfcOptics_Clear(nt) , & ! Input + GeometryInfo , & ! Input + SensorIndex , & ! Input + ChannelIndex , & ! Input + RTSolution_Clear(nt), & ! Output + RTV_Clear(nt) ) ! Internal variable output IF ( Error_Status /= SUCCESS ) THEN - WRITE( Message,'( "Error computing RTSolution for ", a, & + WRITE( Message,'( "Error computing CLEAR SKY RTSolution for ", a, & &", channel ", i0,", profile #",i0)' ) & TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) CYCLE Thread_Loop END IF + END IF - - ! Repeat clear sky for fractionally cloudy atmospheres - IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag).and.RTV(nt)%mth_Azi==0 ) THEN - RTV_Clear(nt)%mth_Azi = RTV(nt)%mth_Azi - SfcOptics_Clear(nt)%mth_Azi = SfcOptics(nt)%mth_Azi - Error_Status = CRTM_Compute_RTSolution( & - Atm_Clear , & ! Input - Surface(m) , & ! Input - AtmOptics_Clear(nt) , & ! Input - SfcOptics_Clear(nt) , & ! Input - GeometryInfo , & ! Input - SensorIndex , & ! Input - ChannelIndex , & ! Input - RTSolution_Clear(nt), & ! Output - RTV_Clear(nt) ) ! Internal variable output - IF ( Error_Status /= SUCCESS ) THEN - WRITE( Message,'( "Error computing CLEAR SKY RTSolution for ", a, & - &", channel ", i0,", profile #",i0)' ) & - TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m - CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) - CYCLE Thread_Loop + IF( mth_Azi == RTV(nt)%n_Azi ) THEN + + ! Combine cloudy and clear radiances for fractional cloud coverage + IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN + IF( RTV(nt)%n_Stokes == 1 ) THEN + r_cloudy(1) = RTSolution(ln,m)%Radiance ! Save the 100% cloudy radiance + RTSolution(ln,m)%Radiance = & + ((ONE - CloudCover%Total_Cloud_Cover) * RTSolution_Clear(nt)%Radiance) + & + (CloudCover%Total_Cloud_Cover * RTSolution(ln,m)%Radiance) + RTSolution(ln,m)%Stokes(1) = RTSolution(ln,m)%Radiance + ELSE + DO ks = 1, RTV(nt)%n_Stokes + r_cloudy(ks) = RTSolution(ln,m)%Stokes(ks) ! Save the 100% cloudy radiance + RTSolution(ln,m)%Stokes(ks) = & + ((ONE - CloudCover%Total_Cloud_Cover) * RTSolution_Clear(nt)%Stokes(ks)) + & + (CloudCover%Total_Cloud_Cover * RTSolution(ln,m)%Stokes(ks)) + END DO + RTSolution(ln,m)%Radiance = RTSolution(ln,m)%Stokes(1) END IF - END IF - - IF(mth_Azi == RTV(nt)%n_Azi) THEN - ! Combine cloudy and clear radiances for fractional cloud coverage - IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN - IF( RTV(nt)%n_Stokes == 1 ) THEN - r_cloudy(1) = RTSolution(ln,m)%Radiance ! Save the 100% cloudy radiance - RTSolution(ln,m)%Radiance = & - ((ONE - CloudCover%Total_Cloud_Cover) * RTSolution_Clear(nt)%Radiance) + & - (CloudCover%Total_Cloud_Cover * RTSolution(ln,m)%Radiance) - RTSolution(ln,m)%Stokes(1) = RTSolution(ln,m)%Radiance - ELSE - DO ks = 1, RTV(nt)%n_Stokes - r_cloudy(ks) = RTSolution(ln,m)%Stokes(ks) ! Save the 100% cloudy radiance - RTSolution(ln,m)%Stokes(ks) = & - ((ONE - CloudCover%Total_Cloud_Cover) * RTSolution_Clear(nt)%Stokes(ks)) + & - (CloudCover%Total_Cloud_Cover * RTSolution(ln,m)%Stokes(ks)) - END DO - RTSolution(ln,m)%Radiance = RTSolution(ln,m)%Stokes(1) - END IF - ! ...Save the cloud cover in the output structure - RTSolution(ln,m)%Total_Cloud_Cover = CloudCover%Total_Cloud_Cover - END IF + ! ...Save the cloud cover in the output structure + RTSolution(ln,m)%Total_Cloud_Cover = CloudCover%Total_Cloud_Cover + END IF - ! The radiance post-processing - CALL Post_Process_RTSolution(Opt, RTSolution(ln,m), & - NLTE_Predictor, & - ChannelIndex, SensorIndex, & - compute_antenna_correction, GeometryInfo) - - !======= Active sensor ======= - ! Calculate reflectivity for active instruments - IF ( SC(SensorIndex)%Is_Active_Sensor .AND. AtmOptics(nt)%Include_Scattering) THEN - CALL CRTM_Compute_Reflectivity(Atm , & ! Input - AtmOptics(nt) , & ! Input - GeometryInfo , & ! Input - SensorIndex , & ! Input - ChannelIndex , & ! Input - RTSolution(ln,m)) ! Input/Output - ENDIF - !============================= - - - IF ( SpcCoeff_IsInfraredSensor( SC(SensorIndex) ) .OR. & - SpcCoeff_IsMicrowaveSensor( SC(SensorIndex) ) ) THEN - ! Perform clear-sky post and pre-processing - IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN - ! Radiance post-processing - CALL Post_Process_RTSolution(Opt, RTSolution_Clear(nt), & + ! The radiance post-processing + CALL Post_Process_RTSolution(Opt, RTSolution(ln,m), & NLTE_Predictor, & ChannelIndex, SensorIndex, & compute_antenna_correction, GeometryInfo) - RTSolution(ln,m)%R_Clear = RTSolution_Clear(nt)%Radiance - RTSolution(ln,m)%Tb_Clear = RTSolution_Clear(nt)%Brightness_Temperature - END IF - - - ! The adjoint radiance pre-processing - CALL Pre_Process_RTSolution_K(Opt, RTSolution(ln,m), RTSolution_K(ln,m), & - NLTE_Predictor, NLTE_Predictor_K(nt), & - ChannelIndex, SensorIndex, & - compute_antenna_correction, GeometryInfo) - - ! More fractionally cloudy atmospheres processing - IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN - ! The adjoint of the clear and cloudy radiance combination - CloudCover_K(nt)%Total_Cloud_Cover = RTSolution_K(ln,m)%Total_Cloud_Cover - RTSolution_K(ln,m)%Total_Cloud_Cover = ZERO - RTSolution_Clear_K(nt)%Radiance = (ONE - CloudCover%Total_Cloud_Cover) * RTSolution_K(ln,m)%Radiance - CloudCover_K(nt)%Total_Cloud_Cover = CloudCover_K(nt)%Total_Cloud_Cover + & - ((r_cloudy(1) - RTSolution_Clear(nt)%Radiance) * RTSolution_K(ln,m)%Radiance) - RTSolution_K(ln,m)%Radiance = CloudCover%Total_Cloud_Cover * RTSolution_K(ln,m)%Radiance - END IF - - END IF + !======= Active sensor ======= + ! Calculate reflectivity for active instruments + IF ( SC(SensorIndex)%Is_Active_Sensor .AND. AtmOptics(nt)%Include_Scattering ) THEN + CALL CRTM_Compute_Reflectivity(Atm , & ! Input + AtmOptics(nt) , & ! Input + GeometryInfo , & ! Input + SensorIndex , & ! Input + ChannelIndex , & ! Input + RTSolution(ln,m)) ! Input/Output + ENDIF + !============================= + + + IF ( SpcCoeff_IsInfraredSensor( SC(SensorIndex) ) .OR. & + SpcCoeff_IsMicrowaveSensor( SC(SensorIndex) ) ) THEN + ! Perform clear-sky post and pre-processing + IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN + ! Radiance post-processing + CALL Post_Process_RTSolution(Opt, RTSolution_Clear(nt), & + NLTE_Predictor, & + ChannelIndex, SensorIndex, & + compute_antenna_correction, GeometryInfo) + RTSolution(ln,m)%R_Clear = RTSolution_Clear(nt)%Radiance + RTSolution(ln,m)%Tb_Clear = RTSolution_Clear(nt)%Brightness_Temperature + END IF + ! The adjoint radiance pre-processing + CALL Pre_Process_RTSolution_K(Opt, RTSolution(ln,m), RTSolution_K(ln,m), & + NLTE_Predictor, NLTE_Predictor_K(nt), & + ChannelIndex, SensorIndex, & + compute_antenna_correction, GeometryInfo) + ! More fractionally cloudy atmospheres processing + IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN + ! The adjoint of the clear and cloudy radiance combination + CloudCover_K(nt)%Total_Cloud_Cover = RTSolution_K(ln,m)%Total_Cloud_Cover + RTSolution_K(ln,m)%Total_Cloud_Cover = ZERO + RTSolution_Clear_K(nt)%Radiance = (ONE - CloudCover%Total_Cloud_Cover) * RTSolution_K(ln,m)%Radiance + CloudCover_K(nt)%Total_Cloud_Cover = CloudCover_K(nt)%Total_Cloud_Cover + & + ((r_cloudy(1) - RTSolution_Clear(nt)%Radiance) * RTSolution_K(ln,m)%Radiance) + RTSolution_K(ln,m)%Radiance = CloudCover%Total_Cloud_Cover * RTSolution_K(ln,m)%Radiance + END IF + END IF - END IF + END IF ! mth_Azi == RTV(nt)%n_Azi - IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag).and.RTV(nt)%mth_Azi==0 ) THEN + IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag).and.RTV(nt)%mth_Azi==0 ) THEN ! The adjoint of the clear sky radiative transfer for fractionally cloudy atmospheres RTV_Clear(nt)%mth_Azi = RTV(nt)%mth_Azi SfcOptics_Clear(nt)%mth_Azi = SfcOptics(nt)%mth_Azi + Error_Status = CRTM_Compute_RTSolution_AD( & + Atm_Clear , & ! FWD Input + Surface(m) , & ! FWD Input + AtmOptics_Clear(nt) , & ! FWD Input + SfcOptics_Clear(nt) , & ! FWD Input + RTSolution_Clear(nt) , & ! FWD Input + RTSolution_Clear_K(nt), & ! K Input + GeometryInfo , & ! Input + SensorIndex , & ! Input + ChannelIndex , & ! Input + Atm_Clear_K(nt) , & ! K Output + Surface_K(ln,m) , & ! K Output + AtmOptics_Clear_K(nt) , & ! K Output + SfcOptics_Clear_K(nt) , & ! K Output + RTV_Clear(nt) ) ! Internal variable input + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'( "Error computing CLEAR SKY RTSolution_K for ", a, & + &", channel ", i0,", profile #",i0)' ) & + TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + CYCLE Thread_Loop + END IF + END IF + + + ! The adjoint of the radiative transfer Error_Status = CRTM_Compute_RTSolution_AD( & - Atm_Clear , & ! FWD Input + Atm , & ! FWD Input Surface(m) , & ! FWD Input - AtmOptics_Clear(nt) , & ! FWD Input - SfcOptics_Clear(nt) , & ! FWD Input - RTSolution_Clear(nt) , & ! FWD Input - RTSolution_Clear_K(nt), & ! K Input + AtmOptics(nt) , & ! FWD Input + SfcOptics(nt) , & ! FWD Input + RTSolution(ln,m) , & ! FWD Input + RTSolution_K(ln,m), & ! K Input GeometryInfo , & ! Input SensorIndex , & ! Input ChannelIndex , & ! Input - Atm_Clear_K(nt) , & ! K Output + Atm_K(nt) , & ! K Output Surface_K(ln,m) , & ! K Output - AtmOptics_Clear_K(nt) , & ! K Output - SfcOptics_Clear_K(nt) , & ! K Output - RTV_Clear(nt) ) ! Internal variable input + AtmOptics_K(nt) , & ! K Output + SfcOptics_K(nt) , & ! K Output + RTV(nt) ) ! Internal variable input IF ( Error_Status /= SUCCESS ) THEN - WRITE( Message,'( "Error computing CLEAR SKY RTSolution_K for ", a, & + WRITE( Message,'( "Error computing RTSolution_K for ", a, & &", channel ", i0,", profile #",i0)' ) & TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) CYCLE Thread_Loop END IF - END IF - - - ! The adjoint of the radiative transfer - Error_Status = CRTM_Compute_RTSolution_AD( & - Atm , & ! FWD Input - Surface(m) , & ! FWD Input - AtmOptics(nt) , & ! FWD Input - SfcOptics(nt) , & ! FWD Input - RTSolution(ln,m) , & ! FWD Input - RTSolution_K(ln,m), & ! K Input - GeometryInfo , & ! Input - SensorIndex , & ! Input - ChannelIndex , & ! Input - Atm_K(nt) , & ! K Output - Surface_K(ln,m) , & ! K Output - AtmOptics_K(nt) , & ! K Output - SfcOptics_K(nt) , & ! K Output - RTV(nt) ) ! Internal variable input - IF ( Error_Status /= SUCCESS ) THEN - WRITE( Message,'( "Error computing RTSolution_K for ", a, & - &", channel ", i0,", profile #",i0)' ) & - TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m - CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) - CYCLE Thread_Loop - END IF - ! Calculate the adjoint for the active sensor reflectivity - IF ( SC(SensorIndex)%Is_Active_Sensor .AND. AtmOptics(nt)%Include_Scattering) THEN - CALL CRTM_Compute_Reflectivity_AD(Atm , & ! Input - AtmOptics(nt) , & ! Input - RTSolution(ln,m), & ! Input - GeometryInfo , & ! Input - SensorIndex , & ! Input - ChannelIndex , & ! Input - AtmOptics_K(nt) , & ! Input/Output - RTSolution_K(ln,m)) ! Input/Output - ENDIF + ! Calculate the adjoint for the active sensor reflectivity + IF ( SC(SensorIndex)%Is_Active_Sensor .AND. AtmOptics(nt)%Include_Scattering ) THEN + CALL CRTM_Compute_Reflectivity_AD(Atm , & ! Input + AtmOptics(nt) , & ! Input + RTSolution(ln,m), & ! Input + GeometryInfo , & ! Input + SensorIndex , & ! Input + ChannelIndex , & ! Input + AtmOptics_K(nt) , & ! Input/Output + RTSolution_K(ln,m)) ! Input/Output + END IF - END DO Azimuth_Fourier_Loop + END DO Azimuth_Fourier_Loop -! The following part corresponding TL part is moved from above CRTM_Compute_RTSolution -! to here because forward radiance is not available above calling CRTM_Compute_RTSolution. - IF ( SpcCoeff_IsVisibleSensor( SC(SensorIndex) ) .OR. & - SpcCoeff_IsUltravioletSensor( SC(SensorIndex) ) ) THEN - IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN - RTSolution(ln,m)%R_Clear = RTSolution_Clear(nt)%Radiance - RTSolution(ln,m)%Tb_Clear = RTSolution_Clear(nt)%Brightness_Temperature + ! The following part corresponding TL part is moved from above CRTM_Compute_RTSolution + ! to here because forward radiance is not available above calling CRTM_Compute_RTSolution. + IF ( SpcCoeff_IsVisibleSensor( SC(SensorIndex) ) .OR. & + SpcCoeff_IsUltravioletSensor( SC(SensorIndex) ) ) THEN - IF( CloudCover%Total_Cloud_Cover > ZERO) THEN - IF( RTV(nt)%n_Stokes == 1 ) THEN - CloudCover_K(nt)%Total_Cloud_Cover = (r_cloudy(1) - RTSolution_Clear(nt)%Radiance) & - * RTSolution_K(ln,m)%Radiance/CloudCover%Total_Cloud_Cover - ELSE - DO ks = 1, RTV(nt)%n_Stokes - CloudCover_K(nt)%Total_Cloud_Cover = CloudCover_K(nt)%Total_Cloud_Cover & - + (r_cloudy(ks) - RTSolution_Clear(nt)%Stokes(ks)) & - * RTSolution_K(ln,m)%Stokes(ks)/CloudCover%Total_Cloud_Cover - END DO + IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) THEN + RTSolution(ln,m)%R_Clear = RTSolution_Clear(nt)%Radiance + RTSolution(ln,m)%Tb_Clear = RTSolution_Clear(nt)%Brightness_Temperature + + IF( CloudCover%Total_Cloud_Cover > ZERO ) THEN + IF( RTV(nt)%n_Stokes == 1 ) THEN + CloudCover_K(nt)%Total_Cloud_Cover = (r_cloudy(1) - RTSolution_Clear(nt)%Radiance) & + * RTSolution_K(ln,m)%Radiance/CloudCover%Total_Cloud_Cover + ELSE + DO ks = 1, RTV(nt)%n_Stokes + CloudCover_K(nt)%Total_Cloud_Cover = CloudCover_K(nt)%Total_Cloud_Cover & + + (r_cloudy(ks) - RTSolution_Clear(nt)%Stokes(ks)) & + * RTSolution_K(ln,m)%Stokes(ks)/CloudCover%Total_Cloud_Cover + END DO + END IF END IF END IF END IF - END IF + ! ################################################### ! TEMPORARY FIX : SENSOR-DEPENDENT AZIMUTH ANGLE LOOP ! ################################################### - ! Compute the adjoint of the all-sky atmospheric transmittance ! for use in FASTEM-X reflection correction transmittance_K = SfcOptics_K(nt)%transmittance @@ -1491,7 +1490,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Compute the adjoint of the combined atmospheric optical properties AtmOptics_K(nt)%Scattering_Optical_Depth = AtmOptics_K(nt)%Scattering_Optical_Depth + & RTSolution_K(ln,m)%SOD - RTSolution_K(ln,m)%SOD = ZERO + RTSolution_K(ln,m)%SOD = ZERO IF( AtmOptics(nt)%Include_Scattering ) THEN CALL CRTM_AtmOptics_Combine_AD( AtmOptics(nt), AtmOptics_K(nt), AOvar(nt) ) END IF @@ -1499,11 +1498,11 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Compute the adjoint aerosol absorption/scattering properties IF ( Atm%n_Aerosols > 0 ) THEN - Error_Status = CRTM_Compute_AerosolScatter_AD( Atm , & ! FWD Input + Error_Status = CRTM_Compute_AerosolScatter_AD( Atm , & ! FWD Input AtmOptics(nt) , & ! FWD Input AtmOptics_K(nt) , & ! K Input - SensorIndex , & ! Input - ChannelIndex, & ! Input + SensorIndex , & ! Input + ChannelIndex , & ! Input Atm_K(nt) , & ! K Output ASvar(nt) ) ! Internal variable input IF ( Error_Status /= SUCCESS ) THEN @@ -1518,12 +1517,12 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Compute the adjoint cloud absorption/scattering properties IF ( Atm%n_Clouds > 0 ) THEN - Error_Status = CRTM_Compute_CloudScatter_AD( Atm , & ! FWD Input + Error_Status = CRTM_Compute_CloudScatter_AD( Atm , & ! FWD Input AtmOptics(nt) , & ! FWD Input AtmOptics_K(nt) , & ! K Input GeometryInfo , & ! Input - SensorIndex , & ! Input - ChannelIndex, & ! Input + SensorIndex , & ! Input + ChannelIndex , & ! Input Atm_K(nt) , & ! K Output CSvar(nt) ) ! Internal variable input IF ( Error_Status /= SUCCESS ) THEN @@ -1569,8 +1568,8 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Compute the adjoint gaseous absorption - CALL CRTM_Compute_AtmAbsorption_AD( SensorIndex , & ! Input - ChannelIndex, & ! Input + CALL CRTM_Compute_AtmAbsorption_AD( SensorIndex , & ! Input + ChannelIndex , & ! Input Predictor(nt) , & ! FWD Input AtmOptics_K(nt) , & ! K Input Predictor_K(nt) , & ! K Output @@ -1580,18 +1579,18 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! K-matrix of the NLTE correction predictor calculations IF ( Opt%Apply_NLTE_Correction ) THEN CALL Compute_NLTE_Predictor_AD( & - NLTE_Predictor , & ! Input + NLTE_Predictor , & ! Input NLTE_Predictor_K(nt) , & ! Input Atm_K(nt) ) ! Output END IF ! K-matrix of the predictor calculations - CALL CRTM_Compute_Predictors_AD( SensorIndex , & ! Input - Atm , & ! FWD Input + CALL CRTM_Compute_Predictors_AD( SensorIndex , & ! Input + Atm , & ! FWD Input Predictor(nt) , & ! FWD Input Predictor_K(nt) , & ! K Input - AncillaryInput, & ! Input + AncillaryInput , & ! Input Atm_K(nt) , & ! K Output PVar(nt) ) ! Internal variable input @@ -1667,7 +1666,6 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) !! CALL CRTM_Predictor_Destroy( Predictor ) DO nt = 1, n_channel_threads - CALL CRTM_Predictor_Destroy( Predictor_K(nt) ) CALL CRTM_AtmOptics_Destroy( AtmOptics(nt) ) CALL CRTM_AtmOptics_Destroy( AtmOptics_K(nt) ) @@ -1747,7 +1745,7 @@ SUBROUTINE Pre_Process_RTSolution_K(Opt, rts, rts_K, & NLTE_Predictor, NLTE_Predictor_K, & ChannelIndex, SensorIndex, & compute_antenna_correction, GeometryInfo) - TYPE(CRTM_Options_type), INTENT(IN) :: Opt + TYPE(CRTM_Options_type), INTENT(IN) :: Opt TYPE(CRTM_RTSolution_type), INTENT(IN OUT) :: rts, rts_K TYPE(NLTE_Predictor_type), INTENT(IN) :: NLTE_Predictor TYPE(NLTE_Predictor_type), INTENT(IN OUT) :: NLTE_Predictor_K