diff --git a/CMakeLists.txt b/CMakeLists.txt index 7d845b22..1bbd59ae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ cmake_minimum_required(VERSION 3.20) -project(crtm VERSION 3.1.2 LANGUAGES Fortran) +project(crtm VERSION 3.2.0 LANGUAGES Fortran) option(OPENMP "Build crtm with OpenMP support" ON) option(FIX_FILE_PATH "Path to fix files (default: fix/)" OFF) @@ -150,4 +150,3 @@ else() endif() #export the CRTM_SOURCE_DIR for use in the test framework (unchanged from original) - diff --git a/README.md b/README.md index 4ea14182..590106d0 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -CRTM REL-3.1.2 +CRTM REL-3.2.0 ==================== [![Build Status](https://app.travis-ci.com/JCSDA/CRTMv3.svg?branch=develop)](https://app.travis-ci.com/JCSDA/CRTMv3) @@ -6,8 +6,9 @@ CRTM REL-3.1.2 Preamble -------- -CRTM v3.1.2 release (`REL-3.1.2`) +CRTM v3.2.0 release (`REL-3.2.0`) +v3.2.0 released TBD v3.1.2 released July 11, 2025 v3.1.1 released August 12, 2024 v3.1.0 (alpha) Released October 31, 2023 @@ -57,7 +58,7 @@ Contents Configuration, building, and testing the library ================================================ -JCSDA CRTM v3.1.2 Build Instructions +JCSDA CRTM v3.2.0 Build Instructions The CRTM repository directory structure looks (something) like: @@ -196,4 +197,3 @@ Known Issues - diff --git a/VERSION.cmake b/VERSION.cmake index 0d30ae20..c7674450 100644 --- a/VERSION.cmake +++ b/VERSION.cmake @@ -3,5 +3,4 @@ # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. -set( ${PROJECT_NAME}_VERSION_STR "3.1.2" ) - +set( ${PROJECT_NAME}_VERSION_STR "3.2.0" ) diff --git a/src/AtmScatter/CRTM_CloudScatter.f90 b/src/AtmScatter/CRTM_CloudScatter.f90 index 36a19eee..9f0f3251 100644 --- a/src/AtmScatter/CRTM_CloudScatter.f90 +++ b/src/AtmScatter/CRTM_CloudScatter.f90 @@ -609,9 +609,11 @@ FUNCTION CRTM_Compute_CloudScatter_TL( & (CSV%ke(kc,n) * Atm_TL%Cloud(n)%Water_Content(kc)) ! Compute the Backscatering coefficient - CScat_TL%Backscat_Coefficient(kc) = CScat_TL%Backscat_Coefficient(kc) + & - (kb_TL * Atm%Cloud(n)%Water_Content(kc)) + & - (CSV%kb(kc,n) * Atm_TL%Cloud(n)%Water_Content(kc)) + IF ( SC(SensorIndex)%Is_Active_Sensor .and. CScat%Include_Scattering ) THEN + CScat_TL%Backscat_Coefficient(kc) = CScat_TL%Backscat_Coefficient(kc) + & + (kb_TL * Atm%Cloud(n)%Water_Content(kc)) + & + (CSV%kb(kc,n) * Atm_TL%Cloud(n)%Water_Content(kc)) + END IF ! Compute the phase matrix coefficients IF( n_Phase_Elements > 0 .and. CScat%Include_Scattering ) THEN diff --git a/src/CRTM_Adjoint_Module.f90 b/src/CRTM_Adjoint_Module.f90 index f62bdb31..0d802470 100644 --- a/src/CRTM_Adjoint_Module.f90 +++ b/src/CRTM_Adjoint_Module.f90 @@ -33,6 +33,7 @@ MODULE CRTM_Adjoint_Module MAX_N_AZIMUTH_FOURIER , & MAX_SOURCE_ZENITH_ANGLE, & MAX_N_STREAMS , & + AIRCRAFT_PRESSURE_THRESHOLD, & MIN_COVERAGE_THRESHOLD , & SCATTERING_ALBEDO_THRESHOLD USE CRTM_SpcCoeff, ONLY: SC, & @@ -327,7 +328,7 @@ FUNCTION CRTM_Adjoint( & ! Local variables CHARACTER(256) :: Message LOGICAL :: Options_Present - INTEGER :: n_Sensors, nc + INTEGER :: n_Sensors INTEGER :: n_Channels INTEGER :: m, n_Profiles ! Local ancillary input structure @@ -400,8 +401,7 @@ FUNCTION CRTM_Adjoint( & ! ------------ ! PROFILE LOOPS ! ------------ - -!JR First loop just checks validity of Atmosphere(m) contents +!JR Loop handles per-profile solution !$OMP PARALLEL DO PRIVATE (m, Opt, AncillaryInput) SCHEDULE (runtime) Profile_Loop2: DO m = 1, n_Profiles ! Check the optional Options structure argument @@ -448,6 +448,7 @@ FUNCTION CRTM_Adjoint( & ! and Post_Process_RTSolution_K also access CRTM_K_Matrix data, but multi-level function ! "contain" clauses cause compiler errors so arguments to these functions were needed. FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) + USE ADA_Module, ONLY: CRTM_SurfRef ! INTEGER, INTENT(in) :: m ! profile index TYPE(CRTM_Options_type), INTENT(IN) :: Opt @@ -463,7 +464,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) INTEGER :: n, l ! sensor index, channel index INTEGER :: SensorIndex INTEGER :: ChannelIndex - INTEGER :: ln, ks + INTEGER :: ln, ks, nc INTEGER :: n_Full_Streams, mth_Azi INTEGER :: cloud_coverage_flag REAL(fp) :: Source_ZA @@ -474,6 +475,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Local atmosphere structure for extra layering TYPE(CRTM_Atmosphere_type) :: Atm, Atm_AD + TYPE(CRTM_Atmosphere_type) :: Atm_Input ! Clear sky structures TYPE(CRTM_Atmosphere_type) :: Atm_Clear , Atm_Clear_AD @@ -500,7 +502,23 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) Error_Status = SUCCESS - ! Allocate the profile independent surface optics local structure + ! Reinitialise the output RTSolution + CALL CRTM_RTSolution_Zero(RTSolution(:,m)) + + Atm_Input = Atmosphere(m) + ! Fix for cloud_Fraction < MIN_COVERAGE_THRESHOLD + IF ( Atm_Input%n_Clouds > 0 ) THEN + !** clear clouds where cloud_fraction < threshold + DO nc = 1, Atm_Input%n_Clouds + WHERE (Atm_Input%Cloud_Fraction(:) < MIN_COVERAGE_THRESHOLD) + Atm_Input%Cloud_Fraction(:) = ZERO + Atm_Input%Cloud(nc)%Water_Content(:) = ZERO + Atm_Input%Cloud(nc)%Effective_Radius(:) = ZERO + END WHERE + END DO + END IF + + ! Allocate the profile independent surface optics local structure CALL CRTM_SfcOptics_Create( SfcOptics , MAX_N_ANGLES, MAX_N_STOKES ) CALL CRTM_SfcOptics_Create( SfcOptics_AD, MAX_N_ANGLES, MAX_N_STOKES ) IF ( (.NOT. CRTM_SfcOptics_Associated(SfcOptics )) .OR. & @@ -512,14 +530,14 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) END IF ! Check the cloud and aerosol coeff. data for cases with clouds and aerosol - IF( Atmosphere(m)%n_Clouds > 0 .AND. .NOT. CRTM_CloudCoeff_IsLoaded() )THEN + IF( Atm_Input%n_Clouds > 0 .AND. .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 ) RETURN END IF - IF( Atmosphere(m)%n_Aerosols > 0 .AND. .NOT. CRTM_AerosolCoeff_IsLoaded() )THEN + IF( Atm_Input%n_Aerosols > 0 .AND. .NOT. CRTM_AerosolCoeff_IsLoaded() )THEN Error_Status = FAILURE WRITE( Message,'("The AerosolCoeff data must be loaded (with CRTM_Init routine) ", & &"for the aerosol case profile #",i0)' ) m @@ -529,10 +547,15 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Copy over forward "non-variable" inputs to adjoint outputs - CALL CRTM_Atmosphere_NonVariableCopy( Atmosphere(m), Atmosphere_AD(m) ) + CALL CRTM_Atmosphere_NonVariableCopy( Atm_Input, Atmosphere_AD(m) ) CALL CRTM_Surface_NonVariableCopy( Surface(m), Surface_AD(m) ) - IF( Opt%n_Stokes > 0 ) RTV%n_Stokes = Opt%n_Stokes + IF ( Options_Present ) THEN + IF( Opt%n_Stokes > 0 ) THEN + RTV%n_Stokes = Opt%n_Stokes + RTV_Clear%n_Stokes = Opt%n_Stokes + END IF + END IF AtmOptics%n_Stokes = RTV%n_Stokes ! ...Assign the option specific SfcOptics input SfcOptics%n_Stokes = RTV%n_Stokes @@ -546,7 +569,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Check the input data if required IF ( Opt%Check_Input ) THEN ! ...Mandatory inputs - Atmosphere_Invalid = .NOT. CRTM_Atmosphere_IsValid( Atmosphere(m) ) + Atmosphere_Invalid = .NOT. CRTM_Atmosphere_IsValid( Atm_Input ) Surface_Invalid = .NOT. CRTM_Surface_IsValid( Surface(m) ) Geometry_Invalid = .NOT. CRTM_Geometry_IsValid( Geometry(m) ) IF ( Atmosphere_Invalid .OR. Surface_Invalid .OR. Geometry_Invalid ) THEN @@ -602,7 +625,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Add extra layers to current atmosphere profile ! if necessary to handle upper atmosphere - Error_Status = CRTM_Atmosphere_AddLayers( Atmosphere(m), Atm ) + Error_Status = CRTM_Atmosphere_AddLayers( Atm_Input, Atm ) IF ( Error_Status /= SUCCESS ) THEN Error_Status = FAILURE WRITE( Message,'("Error adding FWD extra layers to profile #",i0)' ) m @@ -618,9 +641,6 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) RETURN END IF - ! Calculate cloud water density - CALL Calculate_Cloud_Water_Density(Atm) - ! ...Similarly extend a copy of the input adjoint atmosphere Atm_AD = CRTM_Atmosphere_AddLayerCopy( Atmosphere_AD(m), Atm%n_Added_Layers ) @@ -652,6 +672,9 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) RETURN END IF + + ! Calculate cloud water density + CALL Calculate_Cloud_Water_Density(Atm) ! Prepare the atmospheric optics structures ! ...Allocate the atmospheric optics structures based on Atm extension CALL CRTM_AtmOptics_Create( AtmOptics, & @@ -663,13 +686,10 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) MAX_N_LEGENDRE_TERMS, & CloudC%N_PHASE_ELEMENTS ) - IF ( Options_Present ) THEN - AtmOptics%depolarization = Opt%depolarization - AtmOptics_AD%depolarization = Opt%depolarization - IF( Opt%n_Stokes > 0 ) RTV%n_Stokes = Opt%n_Stokes - AtmOptics%n_Stokes = RTV%n_Stokes - AtmOptics_AD%n_Stokes = RTV%n_Stokes - END IF + AtmOptics%depolarization = Opt%depolarization + AtmOptics_AD%depolarization = Opt%depolarization + AtmOptics%n_Stokes = RTV%n_Stokes + AtmOptics_AD%n_Stokes = RTV%n_Stokes IF ( .NOT. CRTM_AtmOptics_Associated( Atmoptics ) .OR. & .NOT. CRTM_AtmOptics_Associated( Atmoptics_AD ) ) THEN @@ -748,8 +768,8 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! ...Copy over surface optics input SfcOptics_Clear%Use_New_MWSSEM = .NOT. Opt%Use_Old_MWSSEM SfcOptics_Clear_AD%Use_New_MWSSEM = .NOT. Opt%Use_Old_MWSSEM - SfcOptics_Clear%n_Stokes = RTV%n_Stokes - SfcOptics_Clear_AD%n_Stokes = RTV%n_Stokes + SfcOptics_Clear%n_Stokes = RTV_Clear%n_Stokes + SfcOptics_Clear_AD%n_Stokes = RTV_Clear%n_Stokes ! ...CLEAR SKY average surface skin temperature for multi-surface types CALL CRTM_Compute_SurfaceT( Surface(m), SfcOptics_Clear ) END IF @@ -776,6 +796,41 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ACCoeff_Associated( SC(SensorIndex)%AC ) .AND. & iFOV /= 0 ) + ! Process aircraft pressure altitude + IF ( Opt%Aircraft_Pressure > ZERO ) THEN + RTV%aircraft%rt = .TRUE. + RTV%aircraft%idx = CRTM_Get_PressureLevelIdx(Atm, Opt%Aircraft_Pressure) + ! ...Issue warning if profile level is TOO different from flight level + IF ( ABS(Atm%Level_Pressure(RTV%aircraft%idx)-Opt%Aircraft_Pressure) > AIRCRAFT_PRESSURE_THRESHOLD ) THEN + WRITE( Message,'("Difference between aircraft pressure level (",es22.15,& + &"hPa) and closest input profile level (",es22.15,& + &"hPa) is larger than recommended (",f4.1,"hPa) for profile #",i0)') & + Opt%Aircraft_Pressure, Atm%Level_Pressure(RTV%aircraft%idx), & + AIRCRAFT_PRESSURE_THRESHOLD, m + CALL Display_Message( ROUTINE_NAME, Message, WARNING ) + END IF + ELSE + RTV%aircraft%rt = .FALSE. + END IF + + ! Process observing downward radiance, Obs_4_downward_P = ZERO means at surface + ! Obs_4_downward_P > ZERO, sensor at the pressure + IF ( Opt%Obs_4_downward_P > ZERO ) THEN + RTV%Obs_4_downward%rt = .TRUE. + RTV%Obs_4_downward%idx = CRTM_Get_PressureLevelIdx(Atm, Opt%Obs_4_downward_P) + ! ...Issue warning if profile level is TOO different from flight level + IF ( ABS(Atm%Level_Pressure(RTV%Obs_4_downward%idx)-Opt%Obs_4_downward_P) > AIRCRAFT_PRESSURE_THRESHOLD ) THEN + WRITE( Message,'("Difference between Obs pressure level (",es22.15,& + &"hPa) and closest input profile level (",es22.15,& + &"hPa) is larger than recommended (",f4.1,"hPa) for profile #",i0)') & + Opt%Obs_4_downward_P, Atm%Level_Pressure(RTV%Obs_4_downward%idx), & + AIRCRAFT_PRESSURE_THRESHOLD, m + CALL Display_Message( ROUTINE_NAME, Message, WARNING ) + END IF + ELSE + RTV%Obs_4_downward%rt = .FALSE. + END IF + ! Allocate the AtmAbsorption predictor structures CALL CRTM_Predictor_Create( & @@ -811,7 +866,9 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) SpcCoeff_IsVisibleSensor(SC(SensorIndex)) ) .AND. & AtmOptics%Include_Scattering ) THEN ! Assign algorithm selector - RTV%RT_Algorithm_Id = Opt%RT_Algorithm_Id + IF ( Options_Present ) THEN + RTV%RT_Algorithm_Id = Opt%RT_Algorithm_Id + END IF CALL RTV_Create( RTV, MAX_N_ANGLES, MAX_N_LEGENDRE_TERMS, Atm%n_Layers ) IF ( .NOT. RTV_Associated(RTV) ) THEN Error_Status=FAILURE @@ -878,7 +935,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Determine the number of streams (n_Full_Streams) in up+downward directions - IF ( Opt%Use_N_Streams ) THEN + IF ( Options_Present .AND. Opt%Use_N_Streams ) THEN n_Full_Streams = Options(m)%n_Streams RTSolution(ln,m)%n_Full_Streams = n_Full_Streams + 2 RTSolution(ln,m)%Scattering_Flag = .TRUE. @@ -909,6 +966,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) IF( SC(SensorIndex)%Solar_Irradiance(ChannelIndex) > ZERO .AND. & Source_ZA < MAX_SOURCE_ZENITH_ANGLE) THEN RTV%Solar_Flag_true = .TRUE. + IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) RTV_Clear%Solar_Flag_true = .TRUE. END IF ! ...Visible channel with solar radiation IF( (SpcCoeff_IsVisibleSensor(SC(SensorIndex)).or.SpcCoeff_IsUltravioletSensor(SC(SensorIndex))) & @@ -1025,7 +1083,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Fill the SfcOptics structure for the optional emissivity input case. SfcOptics%Compute = .TRUE. SfcOptics_Clear%Compute = .TRUE. - IF ( Opt%Use_Emissivity ) THEN + IF ( Options_Present .AND. Opt%Use_Emissivity ) THEN SfcOptics%Compute = .FALSE. SfcOptics%Emissivity(1,1) = Opt%Emissivity(ln) SfcOptics%Reflectivity(1,1,1,1) = ONE - Opt%Emissivity(ln) @@ -1109,6 +1167,15 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) RETURN END IF + IF( Options_Present ) THEN + IF(opt%Derive_Surface_Refl.AND.RTV%mth_Azi==0.AND.RTV%COS_SUN>ZERO) THEN + CALL CRTM_SurfRef(Atm%n_Layers,SUM( AtmOptics%Optical_Depth(:)), & ! Input layer optical depth + SfcOptics%Direct_Reflectivity(SfcOptics%Index_Sat_Ang,1), & + SfcOptics%Index_Sat_Ang, RTSolution(ln,m)%Surface_Planck_Radiance, & + RTSolution(ln,m)%Up_Radiance, RTSolution(ln,m)%Down_Radiance,RTV, Error_Status) + END IF + END IF + ! Repeat clear sky for fractionally cloudy atmospheres IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag).and.RTV%mth_Azi==0 ) THEN RTV_Clear%mth_Azi = RTV%mth_Azi @@ -1300,6 +1367,12 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) END IF END IF + !** output Tb_clear in the case of n_clouds = 0 (note this is NOT aerosol cleared) + IF (Atm%n_Clouds == 0 .OR. CloudCover%Total_Cloud_Cover < MIN_COVERAGE_THRESHOLD) THEN + RTSolution(ln,m)%Tb_clear = RTSolution(ln,m)%Brightness_Temperature + RTSolution(ln,m)%R_clear = RTSolution(ln,m)%Radiance + END IF + ! ################################################### ! TEMPORARY FIX : SENSOR-DEPENDENT AZIMUTH ANGLE LOOP ! ################################################### @@ -1457,7 +1530,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) END IF ! Adjoint of the atmosphere layer addition - Error_Status = CRTM_Atmosphere_AddLayers_AD( Atmosphere(m), Atm_AD, Atmosphere_AD(m) ) + Error_Status = CRTM_Atmosphere_AddLayers_AD( Atm_Input, Atm_AD, Atmosphere_AD(m) ) IF ( Error_Status /= SUCCESS ) THEN Error_Status = FAILURE diff --git a/src/CRTM_Forward_Module.f90 b/src/CRTM_Forward_Module.f90 index f0141d01..41c8bba1 100644 --- a/src/CRTM_Forward_Module.f90 +++ b/src/CRTM_Forward_Module.f90 @@ -802,7 +802,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) WRITE( Message,'("Difference between Obs pressure level (",es22.15,& &"hPa) and closest input profile level (",es22.15,& &"hPa) is larger than recommended (",f4.1,"hPa) for profile #",i0)') & - Opt%Obs_4_downward_P, Atm%Level_Pressure(RTV%Obs_4_downward%idx), & + Opt%Obs_4_downward_P, Atm%Level_Pressure(RTV(nt)%Obs_4_downward%idx), & AIRCRAFT_PRESSURE_THRESHOLD, m CALL Display_Message( ROUTINE_NAME, Message, WARNING ) END IF @@ -912,7 +912,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) END IF ! Determine the number of streams (n_Full_Streams) in up+downward directions - IF ( Opt%Use_N_Streams ) THEN + IF ( Options_Present .AND. Opt%Use_N_Streams ) THEN n_Full_Streams = Opt%n_Streams RTSolution(ln,m)%n_Full_Streams = n_Full_Streams + 2 RTSolution(ln,m)%Scattering_Flag = .TRUE. @@ -1048,7 +1048,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Fill the SfcOptics structures for the optional emissivity input case. SfcOptics(nt)%Compute = .TRUE. - IF ( Opt%Use_Emissivity ) THEN + IF ( Options_Present .AND. Opt%Use_Emissivity ) THEN ! ...Cloudy/all-sky case SfcOptics(nt)%Compute = .FALSE. SfcOptics(nt)%Emissivity(1,1) = Opt%Emissivity(ln) diff --git a/src/CRTM_K_Matrix_Module.f90 b/src/CRTM_K_Matrix_Module.f90 index 66f6eebf..04e00fe2 100644 --- a/src/CRTM_K_Matrix_Module.f90 +++ b/src/CRTM_K_Matrix_Module.f90 @@ -31,6 +31,7 @@ MODULE CRTM_K_Matrix_Module MAX_N_AZIMUTH_FOURIER , & MAX_SOURCE_ZENITH_ANGLE, & MAX_N_STREAMS , & + AIRCRAFT_PRESSURE_THRESHOLD, & MIN_COVERAGE_THRESHOLD , & SCATTERING_ALBEDO_THRESHOLD USE CRTM_SpcCoeff, ONLY: SC, & @@ -429,14 +430,7 @@ FUNCTION CRTM_K_Matrix( & ELSE n_profile_threads = n_Profiles -!** BTJ: temporary preprocessor directive for openMP over channels bypass, permitting modern ifort / ifx versions to run properly -!** https://github.com/JCSDA/CRTMv3/issues/231 - -#if 1 - n_channel_threads = 1 -#else n_channel_threads = MIN(n_Channels, n_omp_threads / n_Profiles) -#endif IF(n_channel_threads > 1) THEN CALL OMP_SET_MAX_ACTIVE_LEVELS(2) ELSE @@ -541,6 +535,7 @@ FUNCTION CRTM_K_Matrix( & ! and Post_Process_RTSolution_K also access CRTM_K_Matrix data, but multi-level function ! "contain" clauses cause compiler errors so arguments to these functions were needed. FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) + USE ADA_Module, ONLY: CRTM_SurfRef INTEGER, INTENT(in) :: m ! profile index TYPE(CRTM_Options_type), INTENT(IN) :: Opt TYPE(CRTM_AncillaryInput_type), INTENT(IN) :: AncillaryInput @@ -599,6 +594,13 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Reinitialise the output RTSolution CALL CRTM_RTSolution_Zero(RTSolution(:,m)) + IF ( Options_Present ) THEN + IF( Opt%n_Stokes > 0 ) THEN + RTV(:)%n_Stokes = Opt%n_Stokes + RTV_Clear(:)%n_Stokes = Opt%n_Stokes + END IF + RTV(:)%RT_Algorithm_Id = Opt%RT_Algorithm_Id + END IF ! Allocate the profile independent surface optics local structure @@ -619,18 +621,11 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) !$OMP END PARALLEL DO IF ( Error_Status == FAILURE ) RETURN - ! Copy over forward "non-variable" inputs to K-matrix outputs -!$OMP PARALLEL DO NUM_THREADS(n_channel_threads) - DO l = 1, n_Channels - CALL CRTM_Atmosphere_NonVariableCopy( Atmosphere(m), Atmosphere_K(l,m) ) - CALL CRTM_Surface_NonVariableCopy( Surface(m), Surface_K(l,m) ) - END DO -!$OMP END PARALLEL DO - ! ...Assign the option specific SfcOptics input !$OMP PARALLEL DO NUM_THREADS(n_channel_threads) DO nt = 1, n_channel_threads SfcOptics(nt)%Use_New_MWSSEM = .NOT. Opt%Use_Old_MWSSEM + SfcOptics_K(nt)%Use_New_MWSSEM = .NOT. Opt%Use_Old_MWSSEM END DO !$OMP END PARALLEL DO @@ -683,6 +678,14 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) END IF END IF + ! Copy over forward "non-variable" inputs to K-matrix outputs +!$OMP PARALLEL DO NUM_THREADS(n_channel_threads) + DO l = 1, n_Channels + CALL CRTM_Atmosphere_NonVariableCopy( Atmosphere(m), Atmosphere_K(l,m) ) + CALL CRTM_Surface_NonVariableCopy( Surface(m), Surface_K(l,m) ) + END DO +!$OMP END PARALLEL DO + ! Process geometry ! ...Compute derived geometry @@ -755,13 +758,10 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) 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 + AtmOptics(nt)%depolarization = Opt%depolarization + AtmOptics_K(nt)%depolarization = Opt%depolarization + AtmOptics(nt)%n_Stokes = RTV(nt)%n_Stokes + AtmOptics_K(nt)%n_Stokes = RTV(nt)%n_Stokes IF ( .NOT. CRTM_AtmOptics_Associated( Atmoptics(nt) ) .OR. & .NOT. CRTM_AtmOptics_Associated( Atmoptics_K(nt) ) ) THEN @@ -798,6 +798,14 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) !$OMP END PARALLEL DO IF ( Error_Status == FAILURE ) RETURN + ! Ensure SfcOptics n_Stokes matches any option overrides +!$OMP PARALLEL DO NUM_THREADS(n_channel_threads) + DO nt = 1, n_channel_threads + SfcOptics(nt)%n_Stokes = RTV(nt)%n_Stokes + SfcOptics_K(nt)%n_Stokes = RTV(nt)%n_Stokes + END DO +!$OMP END PARALLEL DO + ! Determine the type of cloud coverage cloud_coverage_flag = CRTM_Atmosphere_Coverage( atm ) @@ -845,8 +853,9 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) END IF ! ...Copy over surface optics input SfcOptics_Clear(nt)%Use_New_MWSSEM = .NOT. Opt%Use_Old_MWSSEM - SfcOptics_Clear(nt)%n_Stokes = RTV(nt)%n_Stokes - SfcOptics_Clear_K(nt)%n_Stokes = RTV(nt)%n_Stokes + SfcOptics_Clear_K(nt)%Use_New_MWSSEM = .NOT. Opt%Use_Old_MWSSEM + SfcOptics_Clear(nt)%n_Stokes = RTV_Clear(nt)%n_Stokes + SfcOptics_Clear_K(nt)%n_Stokes = RTV_Clear(nt)%n_Stokes ! ...CLEAR SKY average surface skin temperature for multi-surface types CALL CRTM_Compute_SurfaceT( Surface(m), SfcOptics_Clear(nt) ) @@ -889,6 +898,41 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) !$OMP PARALLEL DO NUM_THREADS(n_channel_threads) PRIVATE(Message) DO nt = 1, n_channel_threads + ! Process aircraft pressure altitude + IF ( Opt%Aircraft_Pressure > ZERO ) THEN + RTV(nt)%aircraft%rt = .TRUE. + RTV(nt)%aircraft%idx = CRTM_Get_PressureLevelIdx(Atm, Opt%Aircraft_Pressure) + ! ...Issue warning if profile level is TOO different from flight level + IF ( ABS(Atm%Level_Pressure(RTV(nt)%aircraft%idx)-Opt%Aircraft_Pressure) > AIRCRAFT_PRESSURE_THRESHOLD ) THEN + WRITE( Message,'("Difference between aircraft pressure level (",es22.15,& + &"hPa) and closest input profile level (",es22.15,& + &"hPa) is larger than recommended (",f4.1,"hPa) for profile #",i0)') & + Opt%Aircraft_Pressure, Atm%Level_Pressure(RTV(nt)%aircraft%idx), & + AIRCRAFT_PRESSURE_THRESHOLD, m + CALL Display_Message( ROUTINE_NAME, Message, WARNING ) + END IF + ELSE + RTV(nt)%aircraft%rt = .FALSE. + END IF + + ! Process observing downward radiance, Obs_4_downward_P = ZERO means at surface + ! Obs_4_downward_P > ZERO, sensor at the pressure + IF ( Opt%Obs_4_downward_P > ZERO ) THEN + RTV(nt)%Obs_4_downward%rt = .TRUE. + RTV(nt)%Obs_4_downward%idx = CRTM_Get_PressureLevelIdx(Atm, Opt%Obs_4_downward_P) + ! ...Issue warning if profile level is TOO different from flight level + IF ( ABS(Atm%Level_Pressure(RTV(nt)%Obs_4_downward%idx)-Opt%Obs_4_downward_P) > AIRCRAFT_PRESSURE_THRESHOLD ) THEN + WRITE( Message,'("Difference between Obs pressure level (",es22.15,& + &"hPa) and closest input profile level (",es22.15,& + &"hPa) is larger than recommended (",f4.1,"hPa) for profile #",i0)') & + Opt%Obs_4_downward_P, Atm%Level_Pressure(RTV(nt)%Obs_4_downward%idx), & + AIRCRAFT_PRESSURE_THRESHOLD, m + CALL Display_Message( ROUTINE_NAME, Message, WARNING ) + END IF + ELSE + RTV(nt)%Obs_4_downward%rt = .FALSE. + END IF + ! Allocate the AtmAbsorption predictor structures CALL CRTM_Predictor_Create( & @@ -933,7 +977,9 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) Atm%n_Aerosols > 0 .OR. & SpcCoeff_IsVisibleSensor(SC(SensorIndex)).OR.SpcCoeff_IsUltravioletSensor(SC(SensorIndex)) ) .AND. & AtmOptics(nt)%Include_Scattering ) THEN - RTV(nt)%RT_Algorithm_Id = Opt%RT_Algorithm_Id + IF ( Options_Present ) THEN + RTV(nt)%RT_Algorithm_Id = Opt%RT_Algorithm_Id + END IF CALL RTV_Create( RTV(nt), MAX_N_ANGLES, MAX_N_LEGENDRE_TERMS, Atm%n_Layers ) IF ( .NOT. RTV_Associated(RTV(nt)) ) THEN Error_Status=FAILURE @@ -1077,7 +1123,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) Atm_K(nt)%Height = Atm%Height ! Determine the number of streams (n_Full_Streams) in up+downward directions - IF ( Opt%Use_N_Streams ) THEN + IF ( Options_Present .AND. Opt%Use_N_Streams ) THEN n_Full_Streams = Options(m)%n_Streams RTSolution(ln,m)%n_Full_Streams = n_Full_Streams + 2 RTSolution(ln,m)%Scattering_Flag = .TRUE. @@ -1108,6 +1154,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) IF( SC(SensorIndex)%Solar_Irradiance(ChannelIndex) > ZERO .AND. & Source_ZA < MAX_SOURCE_ZENITH_ANGLE ) THEN RTV(nt)%Solar_Flag_true = .TRUE. + IF ( CRTM_Atmosphere_IsFractional(cloud_coverage_flag) ) RTV_Clear(nt)%Solar_Flag_true = .TRUE. END IF ! ...Visible channel with solar radiation @@ -1228,7 +1275,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Fill the SfcOptics structure for the optional emissivity input case. SfcOptics(nt)%Compute = .TRUE. SfcOptics_Clear(nt)%Compute = .TRUE. - IF ( Opt%Use_Emissivity ) THEN + IF ( Options_Present .AND. Opt%Use_Emissivity ) THEN SfcOptics(nt)%Compute = .FALSE. SfcOptics(nt)%Emissivity(1,1) = Opt%Emissivity(ln) SfcOptics(nt)%Reflectivity(1,1,1,1) = ONE - Opt%Emissivity(ln) @@ -1301,6 +1348,15 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) CYCLE Thread_Loop END IF + IF( Options_Present ) THEN + IF(opt%Derive_Surface_Refl.AND.RTV(nt)%mth_Azi==0.AND.RTV(nt)%COS_SUN>ZERO) THEN + CALL CRTM_SurfRef(Atm%n_Layers,SUM( AtmOptics(nt)%Optical_Depth(:)), & ! Input layer optical depth + SfcOptics(nt)%Direct_Reflectivity(SfcOptics(nt)%Index_Sat_Ang,1), & + SfcOptics(nt)%Index_Sat_Ang, RTSolution(ln,m)%Surface_Planck_Radiance, & + RTSolution(ln,m)%Up_Radiance, RTSolution(ln,m)%Down_Radiance,RTV(nt), Error_Status) + 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 @@ -1380,6 +1436,11 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) RTSolution(ln,m)%Tb_Clear = RTSolution_Clear(nt)%Brightness_Temperature END IF ! The adjoint radiance pre-processing + IF ( RTV(nt)%obs_4_downward%rt ) THEN + IF ( RTSolution_K(ln,m)%Radiance == ZERO ) THEN + RTSolution_K(ln,m)%Radiance = RTSolution_K(ln,m)%Down_Radiance + END IF + END IF CALL Pre_Process_RTSolution_K(Opt, RTSolution(ln,m), RTSolution_K(ln,m), & NLTE_Predictor, NLTE_Predictor_K(nt), & ChannelIndex, SensorIndex, & @@ -1489,6 +1550,12 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) END IF END IF + !** output Tb_clear in the case of n_clouds = 0 (note this is NOT aerosol cleared) + IF (Atm%n_Clouds == 0 .OR. CloudCover%Total_Cloud_Cover < MIN_COVERAGE_THRESHOLD) THEN + RTSolution(ln,m)%Tb_clear = RTSolution(ln,m)%Brightness_Temperature + RTSolution(ln,m)%R_clear = RTSolution(ln,m)%Radiance + END IF + ! ################################################### ! TEMPORARY FIX : SENSOR-DEPENDENT AZIMUTH ANGLE LOOP ! ################################################### diff --git a/src/CRTM_Tangent_Linear_Module.f90 b/src/CRTM_Tangent_Linear_Module.f90 index b16cc514..233039d1 100644 --- a/src/CRTM_Tangent_Linear_Module.f90 +++ b/src/CRTM_Tangent_Linear_Module.f90 @@ -295,7 +295,7 @@ FUNCTION CRTM_Tangent_Linear( & ! Function result INTEGER :: Error_Status ! Local parameters - CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Forward' + CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Tangent_Linear' ! Local variables CHARACTER(256) :: Message LOGICAL :: Options_Present @@ -471,12 +471,12 @@ FUNCTION CRTM_Tangent_Linear( & CALL SYSTEM_CLOCK (count=count_end) elapsed = REAL (count_end - count_start) / REAL (count_rate) elapsed_running = elapsed_running + elapsed - WRITE(6,*) 'CRTM_Forward elapsed =',elapsed - WRITE(6,*) 'CRTM_Forward elapsed running total=',elapsed_running + WRITE(6,*) 'CRTM_Tangent_Linear elapsed =',elapsed + WRITE(6,*) 'CRTM_Tangent_Linear elapsed running total=',elapsed_running END IF IF (output_verification) THEN - WRITE(6,*)'CRTM_Forward inspecting RTSolution...' + WRITE(6,*)'CRTM_Tangent_Linear inspecting RTSolution...' CALL CRTM_RTSolution_Inspect (RTSolution(:,:)) END IF RETURN @@ -544,16 +544,18 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Reinitialise the output RTSolution CALL CRTM_RTSolution_Zero(RTSolution(:,m)) + CALL CRTM_RTSolution_Zero(RTSolution_TL(:,m)) + IF ( Options_Present ) THEN + IF( Opt%n_Stokes > 0 ) THEN + RTV(:)%n_Stokes = Opt%n_Stokes + RTV_Clear(:)%n_Stokes = Opt%n_Stokes + END IF + RTV(:)%RT_Algorithm_Id = Opt%RT_Algorithm_Id + END IF ! Allocate the profile independent surface opticss local structure !$OMP PARALLEL DO NUM_THREADS(n_channel_threads) PRIVATE(Message) DO nt = 1, n_channel_threads - IF ( Options_Present ) THEN - ! ...Assign the option specific SfcOptics input - IF( Opt%n_Stokes > 0 ) RTV(nt)%n_Stokes = Opt%n_Stokes - RTV(nt)%RT_Algorithm_Id = Opt%RT_Algorithm_Id - END IF - CALL CRTM_SfcOptics_Create( SfcOptics(nt) , MAX_N_ANGLES, MAX_N_STOKES ) CALL CRTM_SfcOptics_Create( SfcOptics_TL(nt), MAX_N_ANGLES, MAX_N_STOKES ) IF ( (.NOT. CRTM_SfcOptics_Associated(SfcOptics(nt) )) .OR. & @@ -647,9 +649,6 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) RETURN END IF - ! Calculate cloud water density - CALL Calculate_Cloud_Water_Density(Atm) - Error_Status = CRTM_Atmosphere_AddLayers_TL( Atmosphere(m), Atmosphere_TL(m), Atm_TL ) IF ( Error_Status /= SUCCESS ) THEN Error_Status = FAILURE @@ -695,6 +694,9 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) RETURN END IF + ! Calculate cloud water density + CALL Calculate_Cloud_Water_Density(Atm) + !$OMP PARALLEL DO NUM_THREADS(n_channel_threads) PRIVATE(Message) DO nt = 1, n_channel_threads ! Prepare the atmospheric optics structures @@ -708,15 +710,12 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) MAX_N_LEGENDRE_TERMS, & CloudC%N_PHASE_ELEMENTS ) - IF ( Options_Present ) THEN AtmOptics(nt)%depolarization = Opt%depolarization AtmOptics_TL(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_TL(nt)%n_Stokes = RTV(nt)%n_Stokes AtmOptics(nt)%Include_Scattering = Opt%Include_Scattering AtmOptics_TL(nt)%Include_Scattering = Opt%Include_Scattering - END IF IF ( (.NOT. CRTM_AtmOptics_Associated( Atmoptics(nt) )) .OR. & (.NOT. CRTM_AtmOptics_Associated( Atmoptics_TL(nt) )) ) THEN @@ -798,8 +797,8 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! ...Copy over surface optics input SfcOptics_Clear(nt)%Use_New_MWSSEM = .NOT. Opt%Use_Old_MWSSEM SfcOptics_Clear_TL(nt)%Use_New_MWSSEM = .NOT. Opt%Use_Old_MWSSEM - SfcOptics_Clear(nt)%n_Stokes = RTV(nt)%n_Stokes ! It may be changed for CSEM. - SfcOptics_Clear_TL(nt)%n_Stokes = RTV(nt)%n_Stokes ! It may be changed for CSEM. + SfcOptics_Clear(nt)%n_Stokes = RTV_Clear(nt)%n_Stokes ! It may be changed for CSEM. + SfcOptics_Clear_TL(nt)%n_Stokes = RTV_Clear(nt)%n_Stokes ! It may be changed for CSEM. ! ...CLEAR SKY average surface skin temperature for multi-surface types CALL CRTM_Compute_SurfaceT( Surface(m), SfcOptics_Clear(nt) ) @@ -875,6 +874,41 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) Predictor_TL(nt) , & ! Output PVar(nt) ) ! Internal variable input + ! Process aircraft pressure altitude + IF ( Opt%Aircraft_Pressure > ZERO ) THEN + RTV(nt)%aircraft%rt = .TRUE. + RTV(nt)%aircraft%idx = CRTM_Get_PressureLevelIdx(Atm, Opt%Aircraft_Pressure) + ! ...Issue warning if profile level is TOO different from flight level + IF ( ABS(Atm%Level_Pressure(RTV(nt)%aircraft%idx)-Opt%Aircraft_Pressure) > AIRCRAFT_PRESSURE_THRESHOLD ) THEN + WRITE( Message,'("Difference between aircraft pressure level (",es22.15,& + &"hPa) and closest input profile level (",es22.15,& + &"hPa) is larger than recommended (",f4.1,"hPa) for profile #",i0)') & + Opt%Aircraft_Pressure, Atm%Level_Pressure(RTV(nt)%aircraft%idx), & + AIRCRAFT_PRESSURE_THRESHOLD, m + CALL Display_Message( ROUTINE_NAME, Message, WARNING ) + END IF + ELSE + RTV(nt)%aircraft%rt = .FALSE. + END IF + + ! Process observing downward radiance, Obs_4_downward_P = ZERO means at surface + ! Obs_4_downward_P > ZERO, sensor at the pressure + IF ( Opt%Obs_4_downward_P > ZERO ) THEN + RTV(nt)%Obs_4_downward%rt = .TRUE. + RTV(nt)%Obs_4_downward%idx = CRTM_Get_PressureLevelIdx(Atm, Opt%Obs_4_downward_P) + ! ...Issue warning if profile level is TOO different from flight level + IF ( ABS(Atm%Level_Pressure(RTV(nt)%Obs_4_downward%idx)-Opt%Obs_4_downward_P) > AIRCRAFT_PRESSURE_THRESHOLD ) THEN + WRITE( Message,'("Difference between Obs pressure level (",es22.15,& + &"hPa) and closest input profile level (",es22.15,& + &"hPa) is larger than recommended (",f4.1,"hPa) for profile #",i0)') & + Opt%Obs_4_downward_P, Atm%Level_Pressure(RTV(nt)%Obs_4_downward%idx), & + AIRCRAFT_PRESSURE_THRESHOLD, m + CALL Display_Message( ROUTINE_NAME, Message, WARNING ) + END IF + ELSE + RTV(nt)%Obs_4_downward%rt = .FALSE. + END IF + ! Compute predictors for AtmAbsorption calcs ! ...Allocate the predictor structure @@ -886,7 +920,9 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) SpcCoeff_IsVisibleSensor(SC(SensorIndex)) ) .AND. & AtmOptics(nt)%Include_Scattering ) THEN ! Assign algorithm selector - RTV(nt)%RT_Algorithm_Id = Opt%RT_Algorithm_Id + IF ( Options_Present ) THEN + RTV(nt)%RT_Algorithm_Id = Opt%RT_Algorithm_Id + END IF CALL RTV_Create( RTV(nt), MAX_N_ANGLES, MAX_N_LEGENDRE_TERMS, Atm%n_Layers ) IF ( .NOT. RTV_Associated(RTV(nt)) ) THEN @@ -992,7 +1028,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! CALL CRTM_SfcOptics_Zero( SfcOptics(nt) ) ! CALL CRTM_SfcOptics_Zero( SfcOptics_TL(nt) ) ! Determine the number of streams (n_Full_Streams) in up+downward directions - IF ( Opt%Use_N_Streams ) THEN + IF ( Options_Present .AND. Opt%Use_N_Streams ) THEN n_Full_Streams = Opt%n_Streams RTSolution(ln,m)%n_Full_Streams = n_Full_Streams + 2 RTSolution(ln,m)%Scattering_Flag = .TRUE. @@ -1168,7 +1204,7 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) ! Fill the SfcOptics structures for the optional emissivity input case. SfcOptics(nt)%Compute = .TRUE. SfcOptics_Clear(nt)%Compute = .TRUE. - IF ( Opt%Use_Emissivity ) THEN + IF ( Options_Present .AND. Opt%Use_Emissivity ) THEN ! ...Cloudy/all-sky case SfcOptics(nt)%Compute = .FALSE. SfcOptics(nt)%Emissivity(1,1) = Opt%Emissivity(ln) @@ -1226,6 +1262,15 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) END IF + IF( Options_Present ) THEN + IF(opt%Derive_Surface_Refl.AND.RTV(nt)%mth_Azi==0.AND.RTV(nt)%COS_SUN>ZERO) THEN + CALL CRTM_SurfRef(Atm%n_Layers,SUM( AtmOptics(nt)%Optical_Depth(:)), & ! Input layer optical depth + SfcOptics(nt)%Direct_Reflectivity(SfcOptics(nt)%Index_Sat_Ang,1), & + SfcOptics(nt)%Index_Sat_Ang, RTSolution(ln,m)%Surface_Planck_Radiance, & + RTSolution(ln,m)%Up_Radiance, RTSolution(ln,m)%Down_Radiance,RTV(nt), Error_Status) + END IF + END IF + ! ...Tangent-linear model Error_Status = CRTM_Compute_RTSolution_TL( & @@ -1347,6 +1392,14 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status ) RTSolution_TL(ln,m)%Tb_Clear = RTSolution_Clear_TL(nt)%Brightness_Temperature END IF + !** output Tb_clear in the case of n_clouds = 0 (note this is NOT aerosol cleared) + IF (Atm%n_Clouds == 0 .OR. CloudCover%Total_Cloud_Cover < MIN_COVERAGE_THRESHOLD) THEN + RTSolution(ln,m)%Tb_clear = RTSolution(ln,m)%Brightness_Temperature + RTSolution(ln,m)%R_clear = RTSolution(ln,m)%Radiance + RTSolution_TL(ln,m)%Tb_clear = RTSolution_TL(ln,m)%Brightness_Temperature + RTSolution_TL(ln,m)%R_clear = RTSolution_TL(ln,m)%Radiance + END IF + ! Calculate reflectivity for active instruments IF ( SC(SensorIndex)%Is_Active_Sensor .AND. AtmOptics(nt)%Include_Scattering) THEN CALL CRTM_Compute_Reflectivity(Atm , & ! Input diff --git a/src/CRTM_Version.inc b/src/CRTM_Version.inc index a0fd4b03..87b761c1 100644 --- a/src/CRTM_Version.inc +++ b/src/CRTM_Version.inc @@ -1 +1 @@ -#define CRTM_VERSION 'v3.1.2' +#define CRTM_VERSION 'v3.2.0' diff --git a/src/RTSolution/ADA/ADA_Module.f90 b/src/RTSolution/ADA/ADA_Module.f90 index 42553096..da9ffbca 100644 --- a/src/RTSolution/ADA/ADA_Module.f90 +++ b/src/RTSolution/ADA/ADA_Module.f90 @@ -135,9 +135,9 @@ SUBROUTINE CRTM_ADA(n_Layers, & ! Input number of atmospheric layers CHARACTER(256) :: Message total_opt(0) = ZERO - DO k = 1, n_Layers - total_opt(k) = total_opt(k-1) + T_OD(k) - END DO + DO k = 1, n_Layers + total_opt(k) = total_opt(k-1) + T_OD(k) + END DO nZ = RTV%n_Angles * RTV%n_Stokes RTV%s_Layer_Trans = ZERO @@ -1248,7 +1248,8 @@ SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers direct_reflectivity_TL, & ! Input TL direct reflectivity Pff_TL, & ! Input TL forward phase matrix Pbb_TL, & ! Input TL backward phase matrix - s_rad_up_TL) ! Output TL upward radiance + s_rad_up_TL, & ! Output TL upward radiance + down_rad_TL_out) ! Optional output TL downwelling radiance ! ------------------------------------------------------------------------- ! ! FUNCTION: ! ! This subroutine calculates IR/MW tangent-linear radiance at the top of ! @@ -1277,6 +1278,7 @@ SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers REAL (fp),INTENT(IN),DIMENSION( :,: ) :: reflectivity_TL REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_TL REAL (fp),INTENT(INOUT),DIMENSION( : ) :: direct_reflectivity_TL + REAL (fp),INTENT(OUT), OPTIONAL, DIMENSION( : ) :: down_rad_TL_out ! -------------- internal variables --------------------------------- ! ! Abbreviations: ! ! s: scattering, rad: radiance, trans: transmission, ! @@ -1290,9 +1292,29 @@ SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers REAL (fp), DIMENSION( RTV%n_Angles*RTV%n_Stokes, RTV%n_Angles*RTV%n_Stokes ) :: s_trans_TL,s_refl_TL,Refl_Trans_TL REAL (fp), DIMENSION( RTV%n_Angles*RTV%n_Stokes, RTV%n_Angles*RTV%n_Stokes ) :: s_refl_up_TL,Inv_Gamma_TL,Inv_GammaT_TL REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_TL + LOGICAL :: compute_down + INTEGER :: obs_idx + REAL(fp), ALLOCATABLE :: rad_up_TL(:,:) + REAL(fp), ALLOCATABLE :: refl_up_TL(:,:,:) + REAL(fp), ALLOCATABLE :: layer_trans_TL(:,:,:) + REAL(fp), ALLOCATABLE :: layer_refl_TL(:,:,:) + REAL(fp), ALLOCATABLE :: layer_source_up_TL(:,:) + REAL(fp), ALLOCATABLE :: layer_source_down_TL(:,:) + REAL(fp), ALLOCATABLE :: rad_down_fwd(:,:) + REAL(fp), ALLOCATABLE :: rad_down_TL(:,:) + REAL(fp), ALLOCATABLE :: refl_down_TL_mat(:,:,:) + REAL(fp), ALLOCATABLE :: inv_gamma2_TL(:,:,:) + REAL(fp), ALLOCATABLE :: inv_gamma2T_TL(:,:,:) + REAL(fp), ALLOCATABLE :: inv_gamma3_TL(:,:,:) + REAL(fp), ALLOCATABLE :: temp_vec(:) + REAL(fp), ALLOCATABLE :: temp_vec_TL(:) + REAL(fp), ALLOCATABLE :: rad_downt_TL(:) + REAL(fp), ALLOCATABLE :: temporal_matrix_TL2(:,:) INTEGER :: i, j, k, nZ ! nZ = RTV%n_Angles*RTV%n_Stokes + compute_down = PRESENT(down_rad_TL_out) .AND. RTV%obs_4_downward%rt + obs_idx = RTV%obs_4_downward%idx total_opt(0) = ZERO total_opt_TL(0) = ZERO DO k = 1, n_Layers @@ -1304,16 +1326,49 @@ SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers s_rad_up_TL = ZERO s_refl_up_TL = reflectivity_TL(1:nZ,1:nZ) + IF ( compute_down ) THEN + ALLOCATE(rad_up_TL(nZ,0:n_Layers)) + ALLOCATE(refl_up_TL(nZ,nZ,0:n_Layers)) + ALLOCATE(layer_trans_TL(nZ,nZ,n_Layers)) + ALLOCATE(layer_refl_TL(nZ,nZ,n_Layers)) + ALLOCATE(layer_source_up_TL(nZ,n_Layers)) + ALLOCATE(layer_source_down_TL(nZ,n_Layers)) + ALLOCATE(rad_down_fwd(nZ,0:obs_idx)) + ALLOCATE(rad_down_TL(nZ,0:obs_idx)) + ALLOCATE(refl_down_TL_mat(nZ,nZ,0:obs_idx)) + ALLOCATE(inv_gamma2_TL(nZ,nZ,obs_idx)) + ALLOCATE(inv_gamma2T_TL(nZ,nZ,obs_idx)) + ALLOCATE(inv_gamma3_TL(nZ,nZ,obs_idx)) + ALLOCATE(temp_vec(nZ)) + ALLOCATE(temp_vec_TL(nZ)) + ALLOCATE(rad_downt_TL(nZ)) + ALLOCATE(temporal_matrix_TL2(nZ,nZ)) + rad_up_TL = ZERO + refl_up_TL = ZERO + layer_trans_TL = ZERO + layer_refl_TL = ZERO + layer_source_up_TL = ZERO + layer_source_down_TL = ZERO + rad_down_fwd = ZERO + rad_down_TL = ZERO + refl_down_TL_mat = ZERO + inv_gamma2_TL = ZERO + inv_gamma2T_TL = ZERO + inv_gamma3_TL = ZERO + refl_up_TL(:,:,n_Layers) = s_refl_up_TL + END IF IF( RTV%mth_Azi == 0 ) THEN s_rad_up_TL = emissivity_TL(1:nZ)* RTV%Planck_Surface + emissivity(1:nZ) * Planck_Surface_TL END IF - IF( RTV%Solar_Flag_true ) THEN s_rad_up_TL = s_rad_up_TL+direct_reflectivity_TL(1:nZ)*RTV%COS_SUN*RTV%Solar_irradiance/PI & * exp(-total_opt(n_Layers)/RTV%COS_SUN) & - direct_reflectivity(1:nZ) * RTV%Solar_irradiance/PI & * total_opt_TL(n_Layers) * exp(-total_opt(n_Layers)/RTV%COS_SUN) END IF + IF ( compute_down ) THEN + rad_up_TL(:,n_Layers) = s_rad_up_TL + END IF DO 10 k = n_Layers, 1, -1 s_source_up_TL = ZERO @@ -1395,6 +1450,14 @@ SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers ENDDO ENDIF + IF ( compute_down ) THEN + layer_trans_TL(:,:,k) = s_trans_TL + layer_refl_TL(:,:,k) = s_refl_TL + layer_source_up_TL(:,k) = s_source_up_TL + layer_source_down_TL(:,k) = s_source_down_TL + rad_up_TL(:,k-1) = s_rad_up_TL + refl_up_TL(:,:,k-1) = s_refl_up_TL + END IF 10 CONTINUE ! ! Adding reflected cosmic background radiation @@ -1404,6 +1467,83 @@ SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers ENDDO END IF + IF ( compute_down ) THEN + IF ( RTV%mth_Azi == 0 ) THEN + DO i = 1, nZ, RTV%n_Stokes + rad_down_fwd(i,0) = cosmic_background + END DO + END IF + + DO k = 1, obs_idx + IF(w(k) > SCATTERING_ALBEDO_THRESHOLD .and. maxval(abs(RTV%Pff(1:nZ,1:nZ,k))) > ZERO) THEN + temporal_matrix_TL2 = -matmul(refl_down_TL_mat(:,:,k-1), RTV%s_Layer_Refl(1:nZ,1:nZ,k)) & + - matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), layer_refl_TL(:,:,k)) + inv_gamma2_TL(:,:,k) = -matmul(RTV%Inv_Gamma2(1:nZ,1:nZ,k), & + matmul(temporal_matrix_TL2, RTV%Inv_Gamma2(1:nZ,1:nZ,k))) + inv_gamma2T_TL(:,:,k) = matmul(layer_trans_TL(:,:,k), RTV%Inv_Gamma2(1:nZ,1:nZ,k)) + & + matmul(RTV%s_Layer_Trans(1:nZ,1:nZ,k), inv_gamma2_TL(:,:,k)) + + temp_vec = matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), RTV%s_Layer_Source_UP(1:nZ,k)) + temp_vec_TL = matmul(refl_down_TL_mat(:,:,k-1), RTV%s_Layer_Source_UP(1:nZ,k)) + & + matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), layer_source_up_TL(:,k)) + + rad_down_fwd(:,k) = RTV%s_Layer_Source_DOWN(1:nZ,k) + & + matmul(RTV%Inv_Gamma2T(1:nZ,1:nZ,k), temp_vec + rad_down_fwd(:,k-1)) + rad_down_TL(:,k) = layer_source_down_TL(:,k) + & + matmul(inv_gamma2T_TL(:,:,k), temp_vec + rad_down_fwd(:,k-1)) + & + matmul(RTV%Inv_Gamma2T(1:nZ,1:nZ,k), temp_vec_TL + rad_down_TL(:,k-1)) + + temporal_matrix_TL2 = matmul(refl_down_TL_mat(:,:,k-1), RTV%s_Layer_Trans(1:nZ,1:nZ,k)) + & + matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), layer_trans_TL(:,:,k)) + refl_down_TL_mat(:,:,k) = layer_refl_TL(:,:,k) + & + matmul(inv_gamma2T_TL(:,:,k), RTV%Refl_Trans_DOWN(1:nZ,1:nZ,k)) + & + matmul(RTV%Inv_Gamma2T(1:nZ,1:nZ,k), temporal_matrix_TL2) + ELSE + temp_vec = matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), RTV%s_Layer_Source_UP(1:nZ,k)) + temp_vec_TL = matmul(refl_down_TL_mat(:,:,k-1), RTV%s_Layer_Source_UP(1:nZ,k)) + & + matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), layer_source_up_TL(:,k)) + DO i = 1, nZ + rad_down_fwd(i,k) = RTV%s_Layer_Source_DOWN(i,k) + & + RTV%s_Layer_Trans(i,i,k) * (temp_vec(i) + rad_down_fwd(i,k-1)) + rad_down_TL(i,k) = layer_source_down_TL(i,k) + & + layer_trans_TL(i,i,k) * (temp_vec(i) + rad_down_fwd(i,k-1)) + & + RTV%s_Layer_Trans(i,i,k) * (temp_vec_TL(i) + rad_down_TL(i,k-1)) + END DO + DO i = 1, nZ + DO j = 1, nZ + refl_down_TL_mat(i,j,k) = layer_trans_TL(i,i,k) * RTV%s_Level_Refl_DOWN(i,j,k-1) * & + RTV%s_Layer_Trans(j,j,k) + & + RTV%s_Layer_Trans(i,i,k) * refl_down_TL_mat(i,j,k-1) * & + RTV%s_Layer_Trans(j,j,k) + & + RTV%s_Layer_Trans(i,i,k) * RTV%s_Level_Refl_DOWN(i,j,k-1) * & + layer_trans_TL(j,j,k) + END DO + END DO + END IF + + IF (maxval(abs(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k))) > ZERO) THEN + temporal_matrix_TL2 = -matmul(refl_down_TL_mat(:,:,k), RTV%s_Level_Refl_UP(1:nZ,1:nZ,k)) - & + matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k), refl_up_TL(:,:,k)) + inv_gamma3_TL(:,:,k) = -matmul(RTV%Inv_Gamma3(1:nZ,1:nZ,k), & + matmul(temporal_matrix_TL2, RTV%Inv_Gamma3(1:nZ,1:nZ,k))) + temp_vec = matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k), RTV%s_Level_Rad_UP(1:nZ,k)) + & + rad_down_fwd(:,k) + temp_vec_TL = matmul(refl_down_TL_mat(:,:,k), RTV%s_Level_Rad_UP(1:nZ,k)) + & + matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k), rad_up_TL(:,k)) + & + rad_down_TL(:,k) + rad_downt_TL(:) = matmul(inv_gamma3_TL(:,:,k), temp_vec) + & + matmul(RTV%Inv_Gamma3(1:nZ,1:nZ,k), temp_vec_TL) + ELSE + rad_downt_TL(:) = rad_down_TL(:,k) + END IF + END DO + + down_rad_TL_out(:) = ZERO + IF ( obs_idx > 0 ) THEN + down_rad_TL_out(1:nZ) = rad_downt_TL(1:nZ) + END IF + END IF + RETURN END SUBROUTINE CRTM_ADA_TL ! @@ -1683,7 +1823,8 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers reflectivity_AD, & ! Output AD surface reflectivity direct_reflectivity_AD, & ! Output AD surface direct reflectivity Pff_AD, & ! Output AD forward phase matrix - Pbb_AD) ! Output AD backward phase matrix + Pbb_AD, & ! Output AD backward phase matrix + down_rad_AD_in) ! Optional input AD downwelling radiance ! ------------------------------------------------------------------------- ! ! FUNCTION: ! ! This subroutine calculates IR/MW adjoint radiance at the top of ! @@ -1711,6 +1852,7 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers REAL (fp),INTENT(INOUT),DIMENSION( : ) :: emissivity_AD,direct_reflectivity_AD REAL (fp),INTENT(INOUT),DIMENSION( :,: ) :: reflectivity_AD REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_AD + REAL (fp),INTENT(IN), OPTIONAL, DIMENSION( : ) :: down_rad_AD_in ! -------------- internal variables --------------------------------- ! ! Abbreviations: ! ! s: scattering, rad: radiance, trans: transmission, ! @@ -1727,7 +1869,27 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers Inv_Gamma_AD,Inv_GammaT_AD REAL (fp) :: sum_s_AD, sums_AD, xx REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_AD - INTEGER :: i, j, k,nZ + INTEGER :: i, j, k, nZ + LOGICAL :: compute_down + INTEGER :: obs_idx + REAL(fp), ALLOCATABLE :: rad_down_fwd(:,:) + REAL(fp), ALLOCATABLE :: rad_down_AD(:,:) + REAL(fp), ALLOCATABLE :: refl_down_AD_mat(:,:,:) + REAL(fp), ALLOCATABLE :: layer_trans_AD_dn(:,:,:) + REAL(fp), ALLOCATABLE :: layer_refl_AD_dn(:,:,:) + REAL(fp), ALLOCATABLE :: layer_source_up_AD_dn(:,:) + REAL(fp), ALLOCATABLE :: layer_source_down_AD_dn(:,:) + REAL(fp), ALLOCATABLE :: rad_up_AD_src(:,:) + REAL(fp), ALLOCATABLE :: refl_up_AD_src(:,:,:) + REAL(fp), ALLOCATABLE :: rprev_ad(:,:) + REAL(fp), ALLOCATABLE :: inv_gamma2T_AD_dn(:,:) + REAL(fp), ALLOCATABLE :: inv_gamma2_AD_dn(:,:) + REAL(fp), ALLOCATABLE :: refl_trans_AD_dn(:,:) + REAL(fp), ALLOCATABLE :: ltrans_ad_diag(:) + REAL(fp), ALLOCATABLE :: temp_vec(:) + REAL(fp), ALLOCATABLE :: temp_vec_AD(:) + REAL(fp), ALLOCATABLE :: temporal_matrix(:,:) + REAL(fp), ALLOCATABLE :: temporal_matrix_AD_dn(:,:) ! s_trans_AD = ZERO @@ -1736,12 +1898,179 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers s_refl_up_AD = ZERO nZ = RTV%n_Angles*RTV%n_Stokes + compute_down = PRESENT(down_rad_AD_in) .AND. RTV%obs_4_downward%rt + obs_idx = RTV%obs_4_downward%idx + IF ( compute_down ) THEN + ALLOCATE(rad_down_fwd(nZ,0:obs_idx)) + ALLOCATE(rad_down_AD(nZ,0:obs_idx)) + ALLOCATE(refl_down_AD_mat(nZ,nZ,0:obs_idx)) + ALLOCATE(layer_trans_AD_dn(nZ,nZ,n_Layers)) + ALLOCATE(layer_refl_AD_dn(nZ,nZ,n_Layers)) + ALLOCATE(layer_source_up_AD_dn(nZ,n_Layers)) + ALLOCATE(layer_source_down_AD_dn(nZ,n_Layers)) + ALLOCATE(rad_up_AD_src(nZ,0:n_Layers)) + ALLOCATE(refl_up_AD_src(nZ,nZ,0:n_Layers)) + ALLOCATE(rprev_ad(nZ,nZ)) + ALLOCATE(inv_gamma2T_AD_dn(nZ,nZ)) + ALLOCATE(inv_gamma2_AD_dn(nZ,nZ)) + ALLOCATE(refl_trans_AD_dn(nZ,nZ)) + ALLOCATE(ltrans_ad_diag(nZ)) + ALLOCATE(temp_vec(nZ)) + ALLOCATE(temp_vec_AD(nZ)) + ALLOCATE(temporal_matrix(nZ,nZ)) + ALLOCATE(temporal_matrix_AD_dn(nZ,nZ)) + rad_down_fwd = ZERO + rad_down_AD = ZERO + refl_down_AD_mat = ZERO + layer_trans_AD_dn = ZERO + layer_refl_AD_dn = ZERO + layer_source_up_AD_dn = ZERO + layer_source_down_AD_dn = ZERO + rad_up_AD_src = ZERO + refl_up_AD_src = ZERO + rprev_ad = ZERO + inv_gamma2T_AD_dn = ZERO + inv_gamma2_AD_dn = ZERO + refl_trans_AD_dn = ZERO + ltrans_ad_diag = ZERO + temp_vec = ZERO + temp_vec_AD = ZERO + temporal_matrix = ZERO + temporal_matrix_AD_dn = ZERO + END IF total_opt_AD = ZERO total_opt(0) = ZERO DO k = 1, n_Layers total_opt(k) = total_opt(k-1) + T_OD(k) END DO + IF ( compute_down ) THEN + IF ( RTV%mth_Azi == 0 ) THEN + DO i = 1, nZ, RTV%n_Stokes + rad_down_fwd(i,0) = cosmic_background + END DO + END IF + DO k = 1, obs_idx + IF(w(k) > SCATTERING_ALBEDO_THRESHOLD .and. maxval(abs(RTV%Pff(1:nZ,1:nZ,k))) > ZERO) THEN + temp_vec = matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), RTV%s_Layer_Source_UP(1:nZ,k)) + rad_down_fwd(:,k) = RTV%s_Layer_Source_DOWN(1:nZ,k) + & + matmul(RTV%Inv_Gamma2T(1:nZ,1:nZ,k), temp_vec + rad_down_fwd(:,k-1)) + ELSE + temp_vec = matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), RTV%s_Layer_Source_UP(1:nZ,k)) + DO i = 1, nZ + rad_down_fwd(i,k) = RTV%s_Layer_Source_DOWN(i,k) + & + RTV%s_Layer_Trans(i,i,k) * (temp_vec(i) + rad_down_fwd(i,k-1)) + END DO + END IF + END DO + + rad_down_AD(:,obs_idx) = down_rad_AD_in(1:nZ) + IF (maxval(abs(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,obs_idx))) > ZERO) THEN + temp_vec = matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,obs_idx), RTV%s_Level_Rad_UP(1:nZ,obs_idx)) + & + rad_down_fwd(:,obs_idx) + temp_vec_AD = matmul(transpose(RTV%Inv_Gamma3(1:nZ,1:nZ,obs_idx)), rad_down_AD(:,obs_idx)) + temporal_matrix = matmul(reshape(rad_down_AD(:,obs_idx), (/nZ,1/)), & + reshape(temp_vec, (/1,nZ/))) + temporal_matrix_AD_dn = -matmul(transpose(RTV%Inv_Gamma3(1:nZ,1:nZ,obs_idx)), & + matmul(temporal_matrix, transpose(RTV%Inv_Gamma3(1:nZ,1:nZ,obs_idx)))) + + refl_down_AD_mat(:,:,obs_idx) = refl_down_AD_mat(:,:,obs_idx) + & + matmul(reshape(temp_vec_AD, (/nZ,1/)), & + reshape(RTV%s_Level_Rad_UP(1:nZ,obs_idx), (/1,nZ/))) + rad_up_AD_src(:,obs_idx) = rad_up_AD_src(:,obs_idx) + & + matmul(transpose(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,obs_idx)), temp_vec_AD) + rad_down_AD(:,obs_idx) = rad_down_AD(:,obs_idx) + temp_vec_AD + + refl_down_AD_mat(:,:,obs_idx) = refl_down_AD_mat(:,:,obs_idx) - & + matmul(temporal_matrix_AD_dn, transpose(RTV%s_Level_Refl_UP(1:nZ,1:nZ,obs_idx))) + refl_up_AD_src(:,:,obs_idx) = refl_up_AD_src(:,:,obs_idx) - & + matmul(transpose(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,obs_idx)), temporal_matrix_AD_dn) + END IF + + DO k = obs_idx, 1, -1 + rprev_ad = ZERO + inv_gamma2T_AD_dn = ZERO + inv_gamma2_AD_dn = ZERO + refl_trans_AD_dn = ZERO + ltrans_ad_diag = ZERO + temp_vec_AD = ZERO + + IF(w(k) > SCATTERING_ALBEDO_THRESHOLD .and. maxval(abs(RTV%Pff(1:nZ,1:nZ,k))) > ZERO) THEN + temp_vec = matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), RTV%s_Layer_Source_UP(1:nZ,k)) + + layer_source_down_AD_dn(:,k) = layer_source_down_AD_dn(:,k) + rad_down_AD(:,k) + inv_gamma2T_AD_dn = inv_gamma2T_AD_dn + & + matmul(reshape(rad_down_AD(:,k), (/nZ,1/)), & + reshape(temp_vec + rad_down_fwd(:,k-1), (/1,nZ/))) + temp_vec_AD = temp_vec_AD + matmul(transpose(RTV%Inv_Gamma2T(1:nZ,1:nZ,k)), rad_down_AD(:,k)) + rad_down_AD(:,k-1) = rad_down_AD(:,k-1) + temp_vec_AD + + layer_refl_AD_dn(:,:,k) = layer_refl_AD_dn(:,:,k) + refl_down_AD_mat(:,:,k) + inv_gamma2T_AD_dn = inv_gamma2T_AD_dn + & + matmul(refl_down_AD_mat(:,:,k), transpose(RTV%Refl_Trans_DOWN(1:nZ,1:nZ,k))) + refl_trans_AD_dn = refl_trans_AD_dn + & + matmul(transpose(RTV%Inv_Gamma2T(1:nZ,1:nZ,k)), refl_down_AD_mat(:,:,k)) + + rprev_ad = rprev_ad + matmul(reshape(temp_vec_AD, (/nZ,1/)), & + reshape(RTV%s_Layer_Source_UP(1:nZ,k), (/1,nZ/))) + layer_source_up_AD_dn(:,k) = layer_source_up_AD_dn(:,k) + & + matmul(transpose(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1)), temp_vec_AD) + + rprev_ad = rprev_ad + matmul(refl_trans_AD_dn, transpose(RTV%s_Layer_Trans(1:nZ,1:nZ,k))) + layer_trans_AD_dn(:,:,k) = layer_trans_AD_dn(:,:,k) + & + matmul(transpose(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1)), refl_trans_AD_dn) + + layer_trans_AD_dn(:,:,k) = layer_trans_AD_dn(:,:,k) + & + matmul(inv_gamma2T_AD_dn, transpose(RTV%Inv_Gamma2(1:nZ,1:nZ,k))) + inv_gamma2_AD_dn = matmul(transpose(RTV%s_Layer_Trans(1:nZ,1:nZ,k)), inv_gamma2T_AD_dn) + + temporal_matrix_AD_dn = -matmul(transpose(RTV%Inv_Gamma2(1:nZ,1:nZ,k)), & + matmul(inv_gamma2_AD_dn, transpose(RTV%Inv_Gamma2(1:nZ,1:nZ,k)))) + rprev_ad = rprev_ad - matmul(temporal_matrix_AD_dn, transpose(RTV%s_Layer_Refl(1:nZ,1:nZ,k))) + layer_refl_AD_dn(:,:,k) = layer_refl_AD_dn(:,:,k) - & + matmul(transpose(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1)), temporal_matrix_AD_dn) + + refl_down_AD_mat(:,:,k-1) = refl_down_AD_mat(:,:,k-1) + rprev_ad + ELSE + temp_vec = matmul(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1), RTV%s_Layer_Source_UP(1:nZ,k)) + + DO i = 1, nZ + layer_source_down_AD_dn(i,k) = layer_source_down_AD_dn(i,k) + rad_down_AD(i,k) + temp_vec_AD(i) = temp_vec_AD(i) + RTV%s_Layer_Trans(i,i,k) * rad_down_AD(i,k) + ltrans_ad_diag(i) = ltrans_ad_diag(i) + rad_down_AD(i,k) * & + (temp_vec(i) + rad_down_fwd(i,k-1)) + rad_down_AD(i,k-1) = rad_down_AD(i,k-1) + RTV%s_Layer_Trans(i,i,k) * rad_down_AD(i,k) + END DO + + DO i = 1, nZ + DO j = 1, nZ + rprev_ad(i,j) = rprev_ad(i,j) + RTV%s_Layer_Trans(i,i,k) * refl_down_AD_mat(i,j,k) * & + RTV%s_Layer_Trans(j,j,k) + ltrans_ad_diag(i) = ltrans_ad_diag(i) + refl_down_AD_mat(i,j,k) * & + RTV%s_Level_Refl_DOWN(i,j,k-1) * RTV%s_Layer_Trans(j,j,k) + ltrans_ad_diag(j) = ltrans_ad_diag(j) + refl_down_AD_mat(i,j,k) * & + RTV%s_Layer_Trans(i,i,k) * RTV%s_Level_Refl_DOWN(i,j,k-1) + END DO + END DO + + DO i = 1, nZ + layer_trans_AD_dn(i,i,k) = layer_trans_AD_dn(i,i,k) + ltrans_ad_diag(i) + END DO + + rprev_ad = rprev_ad + matmul(reshape(temp_vec_AD, (/nZ,1/)), & + reshape(RTV%s_Layer_Source_UP(1:nZ,k), (/1,nZ/))) + layer_source_up_AD_dn(:,k) = layer_source_up_AD_dn(:,k) + & + matmul(transpose(RTV%s_Level_Refl_DOWN(1:nZ,1:nZ,k-1)), temp_vec_AD) + refl_down_AD_mat(:,:,k-1) = refl_down_AD_mat(:,:,k-1) + rprev_ad + END IF + END DO + END IF + + IF ( compute_down .AND. obs_idx == 0 ) THEN + s_rad_up_AD = s_rad_up_AD + rad_up_AD_src(:,0) + s_refl_up_AD = s_refl_up_AD + refl_up_AD_src(:,:,0) + END IF + ! Adding reflected cosmic background radiation DO i = 1, nZ, RTV%n_Stokes @@ -1753,9 +2082,20 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers ! DO 10 k = 1, n_Layers + IF ( compute_down .AND. obs_idx == (k-1) ) THEN + s_rad_up_AD = s_rad_up_AD + rad_up_AD_src(:,k-1) + s_refl_up_AD = s_refl_up_AD + refl_up_AD_src(:,:,k-1) + END IF s_source_up_AD = ZERO s_source_down_AD = ZERO s_trans_AD = ZERO + s_refl_AD = ZERO + IF ( compute_down ) THEN + s_source_up_AD = s_source_up_AD + layer_source_up_AD_dn(:,k) + s_source_down_AD = s_source_down_AD + layer_source_down_AD_dn(:,k) + s_trans_AD = s_trans_AD + layer_trans_AD_dn(:,:,k) + s_refl_AD = s_refl_AD + layer_refl_AD_dn(:,:,k) + END IF ! ! Compute tranmission and reflection matrices for a layer IF(w(k) > SCATTERING_ALBEDO_THRESHOLD .and. maxval(abs(RTV%Pff(1:nZ,1:nZ,k))) > ZERO) THEN @@ -1763,7 +2103,7 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers refl_down(1:nZ,k) = matmul(RTV%s_Level_Refl_UP(1:nZ,1:nZ,k), & RTV%s_Layer_Source_DOWN(1:nZ,k)) - s_refl_AD = s_refl_up_AD + s_refl_AD = s_refl_AD + s_refl_up_AD Inv_GammaT_AD = matmul(s_refl_up_AD,transpose(RTV%Refl_Trans(1:nZ,1:nZ,k))) Refl_Trans_AD = matmul(transpose(RTV%Inv_GammaT(1:nZ,1:nZ,k)),s_refl_up_AD) @@ -1851,6 +2191,11 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers ENDIF 10 CONTINUE + IF ( compute_down .AND. obs_idx == n_Layers ) THEN + s_rad_up_AD = s_rad_up_AD + rad_up_AD_src(:,n_Layers) + s_refl_up_AD = s_refl_up_AD + refl_up_AD_src(:,:,n_Layers) + END IF + ! IF( RTV%Solar_Flag_true ) THEN xx = RTV%Solar_irradiance/PI * exp(-total_opt(n_Layers)/RTV%COS_SUN) diff --git a/src/RTSolution/CRTM_RTSolution.f90 b/src/RTSolution/CRTM_RTSolution.f90 index b6654eca..697bfcfe 100644 --- a/src/RTSolution/CRTM_RTSolution.f90 +++ b/src/RTSolution/CRTM_RTSolution.f90 @@ -558,6 +558,7 @@ FUNCTION CRTM_Compute_RTSolution_TL( & ! Local variables CHARACTER(256) :: Message INTEGER :: nZ + INTEGER :: i REAL(fp) :: User_Emissivity_TL, Direct_Reflectivity_TL REAL(fp) :: Planck_Surface_TL ! Surface TL radiance REAL(fp), DIMENSION( 0:Atmosphere%n_Layers ) :: Planck_Atmosphere_TL ! *LAYER* TL radiances @@ -570,12 +571,16 @@ FUNCTION CRTM_Compute_RTSolution_TL( & (RTV%n_Angles+1) * RTV%n_Stokes, & Atmosphere%n_Layers ) :: Pbb_TL ! Backward scattering TL phase matrix REAL(fp), DIMENSION( RTV%n_Angles * RTV%n_Stokes ) :: Scattering_Radiance_TL + REAL(fp), DIMENSION( RTV%n_Angles * RTV%n_Stokes ) :: Down_Scatter_TL REAL(fp) :: Radiance_TL + REAL(fp) :: Down_Radiance_TL ! ------ ! Set up ! ------ Error_Status = SUCCESS + Down_Radiance_TL = ZERO + Down_Scatter_TL = ZERO RTSolution_TL%RT_Algorithm_Name = RTSolution%RT_Algorithm_Name @@ -607,6 +612,13 @@ FUNCTION CRTM_Compute_RTSolution_TL( & END IF nZ = RTV%n_Angles * RTV%n_Stokes + IF ( RTV%Visible_Flag_true ) THEN + DO i = 1, nZ + IF ( SfcOptics%Direct_Reflectivity(i,1) >= ONE ) THEN + SfcOptics_TL%Direct_Reflectivity(i,1) = ZERO + END IF + END DO + END IF IF( RTV%n_Stokes > 1 ) THEN CALL Reshape_Surf_Opt(RTV%n_Angles, RTV%n_Stokes, SfcOptics_TL%Emissivity, SfcOptics_TL%Direct_Reflectivity, & SfcOptics_TL%Reflectivity, SfcOptics_TL%S_Emissivity, SfcOptics_TL%S_Direct_Ref, SfcOptics_TL%S_Reflectivity) @@ -615,6 +627,28 @@ FUNCTION CRTM_Compute_RTSolution_TL( & ! --------------------------------------------- ! Select the RT model IF( RTV%Scattering_RT ) THEN + IF ( RTV%obs_4_downward%rt ) THEN + CALL CRTM_ADA_TL( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo + AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + RTV%Cosmic_Background_Radiance, & ! cosmic background radiation + SfcOptics%S_Emissivity(1:nZ), & ! Input, FWD surface emissivity + SfcOptics%S_Direct_Ref(1:nZ), & ! Input, surface direct reflectivity + RTV, & ! Input, structure containing forward results + Planck_Atmosphere_TL, & ! Input, TL layer radiances + Planck_Surface_TL, & ! Input, TL surface radiance + AtmOptics_TL%Single_Scatter_Albedo, & ! Input, TL layer single scattering albedo + AtmOptics_TL%Optical_Depth, & ! Input, TL layer optical depth + SfcOptics_TL%S_Emissivity(1:nZ), & ! Input, TL surface emissivity + SfcOptics_TL%S_Reflectivity(1:nZ,1:nZ), & ! Input, TL surface reflectivity + SfcOptics_TL%S_Direct_Ref(1:nZ), & ! Input, TL surface direct reflectivity + Pff_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer forward phase matrix + Pbb_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer backward phase matrix + Scattering_Radiance_TL(1:nZ), & ! Output, TL radiances + down_rad_TL_out=Down_Scatter_TL(1:nZ) ) ! Optional output, TL downwelling radiance + Scattering_Radiance_TL(1:nZ) = Down_Scatter_TL(1:nZ) + ELSE CALL CRTM_ADA_TL( & Atmosphere%n_Layers, & ! Input, number of atmospheric layers AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo @@ -633,7 +667,9 @@ FUNCTION CRTM_Compute_RTSolution_TL( & Pff_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer forward phase matrix Pbb_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer backward phase matrix Scattering_Radiance_TL(1:nZ) ) ! Output, TL radiances + END IF ELSE + Down_Radiance_TL = ZERO CALL CRTM_Emission_TL( & Atmosphere%n_Layers, & ! Input, number of atmospheric layers RTV%n_Angles, & ! Input, number of discrete zenith angles @@ -653,7 +689,8 @@ FUNCTION CRTM_Compute_RTSolution_TL( & SfcOptics_TL%S_Emissivity(1:nZ), & ! Input, TL surface emissivity SfcOptics_TL%S_Reflectivity(1:nZ,1:nZ), & ! Input, TL surface reflectivity SfcOptics_TL%S_Direct_Ref(1:nZ), & ! Input, TL surface reflectivity for a point source - Radiance_TL ) ! Output, TL radiances + Radiance_TL, & ! Output, TL radiances + Down_Radiance_TL ) ! Optional output, TL downwelling radiance END IF ELSE IF( RTV%Scattering_RT ) THEN @@ -661,44 +698,89 @@ FUNCTION CRTM_Compute_RTSolution_TL( & SELECT CASE(RTV%RT_Algorithm_Id) CASE (RT_ADA, RT_VMOM) ! NESDIS advanced adding-doubling method - CALL CRTM_ADA_TL( & - Atmosphere%n_Layers, & ! Input, number of atmospheric layers - AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo - AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth - RTV%Cosmic_Background_Radiance, & ! cosmic background radiation - SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity - SfcOptics%Direct_Reflectivity(1:nZ,1), & ! Input, surface direct reflectivity - RTV, & ! Input, structure containing forward results - Planck_Atmosphere_TL, & ! Input, TL layer radiances - Planck_Surface_TL, & ! Input, TL surface radiance - AtmOptics_TL%Single_Scatter_Albedo, & ! Input, TL layer single scattering albedo - AtmOptics_TL%Optical_Depth, & ! Input, TL layer optical depth - SfcOptics_TL%Emissivity(1:nZ,1), & ! Input, TL surface emissivity - SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, TL surface reflectivity - SfcOptics_TL%Direct_Reflectivity(1:nZ,1), & ! Input, TL surface direct reflectivity - Pff_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer forward phase matrix - Pbb_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer backward phase matrix - Scattering_Radiance_TL(1:nZ) ) ! Output, TL radiances + IF ( RTV%obs_4_downward%rt ) THEN + CALL CRTM_ADA_TL( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo + AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + RTV%Cosmic_Background_Radiance, & ! cosmic background radiation + SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity + SfcOptics%Direct_Reflectivity(1:nZ,1), & ! Input, surface direct reflectivity + RTV, & ! Input, structure containing forward results + Planck_Atmosphere_TL, & ! Input, TL layer radiances + Planck_Surface_TL, & ! Input, TL surface radiance + AtmOptics_TL%Single_Scatter_Albedo, & ! Input, TL layer single scattering albedo + AtmOptics_TL%Optical_Depth, & ! Input, TL layer optical depth + SfcOptics_TL%Emissivity(1:nZ,1), & ! Input, TL surface emissivity + SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, TL surface reflectivity + SfcOptics_TL%Direct_Reflectivity(1:nZ,1), & ! Input, TL surface direct reflectivity + Pff_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer forward phase matrix + Pbb_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer backward phase matrix + Scattering_Radiance_TL(1:nZ), & ! Output, TL radiances + down_rad_TL_out=Down_Scatter_TL(1:nZ) ) ! Optional output, TL downwelling radiance + Scattering_Radiance_TL(1:nZ) = Down_Scatter_TL(1:nZ) + ELSE + CALL CRTM_ADA_TL( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo + AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + RTV%Cosmic_Background_Radiance, & ! cosmic background radiation + SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity + SfcOptics%Direct_Reflectivity(1:nZ,1), & ! Input, surface direct reflectivity + RTV, & ! Input, structure containing forward results + Planck_Atmosphere_TL, & ! Input, TL layer radiances + Planck_Surface_TL, & ! Input, TL surface radiance + AtmOptics_TL%Single_Scatter_Albedo, & ! Input, TL layer single scattering albedo + AtmOptics_TL%Optical_Depth, & ! Input, TL layer optical depth + SfcOptics_TL%Emissivity(1:nZ,1), & ! Input, TL surface emissivity + SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, TL surface reflectivity + SfcOptics_TL%Direct_Reflectivity(1:nZ,1), & ! Input, TL surface direct reflectivity + Pff_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer forward phase matrix + Pbb_TL(1:nZ,1:(nZ+1),:), & ! Input, TL layer backward phase matrix + Scattering_Radiance_TL(1:nZ) ) ! Output, TL radiances + END IF CASE (RT_SOI) ! UW SOI RT solver - CALL CRTM_SOI_TL( & - Atmosphere%n_Layers, & ! Input, number of atmospheric layers - AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo - AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth - SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity - SfcOptics%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, surface reflectivity - SfcOptics%Index_Sat_Ang, & ! Input, Satellite angle index - RTV, & ! Input, structure containing forward results - Planck_Atmosphere_TL, & ! Input, TL layer radiances - Planck_Surface_TL, & ! Input, TL surface radiance - AtmOptics_TL%Single_Scatter_Albedo, & ! Input, TL layer single scattering albedo - AtmOptics_TL%Optical_Depth, & ! Input, TL layer optical depth - SfcOptics_TL%Emissivity(1:nZ,1), & ! Input, TL surface emissivity - SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, TL surface reflectivity - Pff_TL(1:nZ,1:nZ,:), & ! Input, TL layer forward phase matrix - Pbb_TL(1:nZ,1:nZ,:), & ! Input, TL layer backward phase matrix - Scattering_Radiance_TL(1:nZ) ) ! Output, TL radiances + IF ( RTV%obs_4_downward%rt ) THEN + CALL CRTM_SOI_TL( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo + AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity + SfcOptics%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, surface reflectivity + SfcOptics%Index_Sat_Ang, & ! Input, Satellite angle index + RTV, & ! Input, structure containing forward results + Planck_Atmosphere_TL, & ! Input, TL layer radiances + Planck_Surface_TL, & ! Input, TL surface radiance + AtmOptics_TL%Single_Scatter_Albedo, & ! Input, TL layer single scattering albedo + AtmOptics_TL%Optical_Depth, & ! Input, TL layer optical depth + SfcOptics_TL%Emissivity(1:nZ,1), & ! Input, TL surface emissivity + SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, TL surface reflectivity + Pff_TL(1:nZ,1:nZ,:), & ! Input, TL layer forward phase matrix + Pbb_TL(1:nZ,1:nZ,:), & ! Input, TL layer backward phase matrix + Scattering_Radiance_TL(1:nZ), & ! Output, TL radiances + down_rad_TL_out=Down_Scatter_TL(1:nZ) ) ! Optional output, TL downwelling radiance + Scattering_Radiance_TL(1:nZ) = Down_Scatter_TL(1:nZ) + ELSE + CALL CRTM_SOI_TL( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo + AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity + SfcOptics%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, surface reflectivity + SfcOptics%Index_Sat_Ang, & ! Input, Satellite angle index + RTV, & ! Input, structure containing forward results + Planck_Atmosphere_TL, & ! Input, TL layer radiances + Planck_Surface_TL, & ! Input, TL surface radiance + AtmOptics_TL%Single_Scatter_Albedo, & ! Input, TL layer single scattering albedo + AtmOptics_TL%Optical_Depth, & ! Input, TL layer optical depth + SfcOptics_TL%Emissivity(1:nZ,1), & ! Input, TL surface emissivity + SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, TL surface reflectivity + Pff_TL(1:nZ,1:nZ,:), & ! Input, TL layer forward phase matrix + Pbb_TL(1:nZ,1:nZ,:), & ! Input, TL layer backward phase matrix + Scattering_Radiance_TL(1:nZ) ) ! Output, TL radiances + END IF CASE DEFAULT Error_Status = FAILURE WRITE(Message,'("Incorrect TL RT_Algorithm_ID, ",i0,", do not fit model")') & @@ -710,6 +792,7 @@ FUNCTION CRTM_Compute_RTSolution_TL( & ! ----------------- ! Emission model RT ! ----------------- + Down_Radiance_TL = ZERO CALL CRTM_Emission_TL( & Atmosphere%n_Layers, & ! Input, number of atmospheric layers RTV%n_Angles, & ! Input, number of discrete zenith angles @@ -729,7 +812,8 @@ FUNCTION CRTM_Compute_RTSolution_TL( & SfcOptics_TL%Emissivity(1:nZ,1), & ! Input, TL surface emissivity SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, TL surface reflectivity SfcOptics_TL%Direct_Reflectivity(1:nZ,1), & ! Input, TL surface reflectivity for a point source - Radiance_TL ) ! Output, TL radiances + Radiance_TL, & ! Output, TL radiances + Down_Radiance_TL ) ! Optional output, TL downwelling radiance END IF Error_Status = Assign_Common_Output_TL( SfcOptics , & @@ -740,7 +824,8 @@ FUNCTION CRTM_Compute_RTSolution_TL( & SensorIndex , & ChannelIndex , & RTV , & - RTSolution_TL ) + RTSolution_TL , & + Down_Radiance_TL ) IF ( Error_Status /= SUCCESS ) THEN Message = 'Error assigning output for TL RTSolution algorithms' CALL Display_Message( ROUTINE_NAME, TRIM(Message), Error_Status ) @@ -943,6 +1028,7 @@ FUNCTION CRTM_Compute_RTSolution_AD( & (RTV%n_Angles+1) * RTV%n_Stokes, & Atmosphere%n_Layers ) :: Pbb_AD ! Backward scattering AD phase matrix REAL (fp),DIMENSION( RTV%n_Angles * RTV%n_Stokes ) :: Scattering_Radiance_AD + REAL (fp),DIMENSION( RTV%n_Angles * RTV%n_Stokes ) :: Down_Scatter_AD REAL (fp) :: Radiance_AD(MAX_N_STOKES) @@ -950,6 +1036,7 @@ FUNCTION CRTM_Compute_RTSolution_AD( & ! Setup ! ----- Error_Status = SUCCESS + Down_Scatter_AD = ZERO Pff_AD = ZERO Pbb_AD = ZERO @@ -986,32 +1073,81 @@ FUNCTION CRTM_Compute_RTSolution_AD( & ! Initialise the input adjoint radiance Scattering_Radiance_AD = ZERO - Scattering_Radiance_AD(n1:n1-1+RTV%n_Stokes) = Radiance_AD(1:RTV%n_Stokes) ! qliu 1/11 + IF ( RTV%obs_4_downward%rt ) THEN + Down_Scatter_AD(n1:n1-1+RTV%n_Stokes) = Radiance_AD(1:RTV%n_Stokes) + ELSE + Scattering_Radiance_AD(n1:n1-1+RTV%n_Stokes) = Radiance_AD(1:RTV%n_Stokes) ! qliu 1/11 + END IF END IF IF( RTV%n_Stokes > 1 ) THEN ! Select the RT model IF( RTV%Scattering_RT ) THEN - CALL CRTM_ADA_AD( & + IF ( RTV%obs_4_downward%rt ) THEN + CALL CRTM_ADA_AD( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo + AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + RTV%Cosmic_Background_Radiance, & ! Input, cosmic background radiation + SfcOptics%S_Emissivity(1:nZ), & ! Input, FWD surface emissivity + SfcOptics%S_Direct_Ref(1:nZ), & ! Input, FWD surface reflectivity for a point source + RTV, & ! In/Output, internal variables + Scattering_Radiance_AD(1:nZ), & ! Input, AD radiances + Planck_Atmosphere_AD, & ! Output, AD layer radiances + Planck_Surface_AD, & ! Output, AD surface radiance + AtmOptics_AD%Single_Scatter_Albedo, & ! Output, AD layer single scattering albedo + AtmOptics_AD%Optical_Depth, & ! Output, AD layer optical depth + SfcOptics_AD%S_Emissivity(1:nZ), & ! Output, AD surface emissivity + SfcOptics_AD%S_Reflectivity(1:nZ,1:nZ), & ! Output, AD surface reflectivity + SfcOptics_AD%S_Direct_Ref(1:nZ), & ! Output, AD surface reflectivity for a point source + Pff_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer forward phase matrix + Pbb_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer backward phase matrix + down_rad_AD_in=Down_Scatter_AD(1:nZ) ) ! Optional input, AD downwelling radiance + ELSE + CALL CRTM_ADA_AD( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo + AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + RTV%Cosmic_Background_Radiance, & ! Input, cosmic background radiation + SfcOptics%S_Emissivity(1:nZ), & ! Input, FWD surface emissivity + SfcOptics%S_Direct_Ref(1:nZ), & ! Input, FWD surface reflectivity for a point source + RTV, & ! In/Output, internal variables + Scattering_Radiance_AD(1:nZ), & ! Input, AD radiances + Planck_Atmosphere_AD, & ! Output, AD layer radiances + Planck_Surface_AD, & ! Output, AD surface radiance + AtmOptics_AD%Single_Scatter_Albedo, & ! Output, AD layer single scattering albedo + AtmOptics_AD%Optical_Depth, & ! Output, AD layer optical depth + SfcOptics_AD%S_Emissivity(1:nZ), & ! Output, AD surface emissivity + SfcOptics_AD%S_Reflectivity(1:nZ,1:nZ), & ! Output, AD surface reflectivity + SfcOptics_AD%S_Direct_Ref(1:nZ), & ! Output, AD surface reflectivity for a point source + Pff_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer forward phase matrix + Pbb_AD(1:nZ,1:(nZ+1),:) ) ! Output, AD layer backward phase matrix + END IF + ELSE + IF ( RTV%obs_4_downward%rt ) THEN + CALL CRTM_Emission_AD( & Atmosphere%n_Layers, & ! Input, number of atmospheric layers - AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo - AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth - RTV%Cosmic_Background_Radiance, & ! Input, cosmic background radiation - SfcOptics%S_Emissivity(1:nZ), & ! Input, FWD surface emissivity + RTV%n_Angles, & ! Input, number of discrete zenith angles + GeometryInfo%Cosine_Sensor_Zenith, & ! Input, cosine of sensor zenith angle + RTV%Planck_Atmosphere, & ! Input, FWD layer radiances + RTV%Planck_Surface, & ! Input, FWD surface radiance + SfcOptics%S_Emissivity(1:nZ), & ! Input, FWD surface emissivity + SfcOptics%S_Reflectivity(1:nZ,1:nZ), & ! Input, FWD surface reflectivity SfcOptics%S_Direct_Ref(1:nZ), & ! Input, FWD surface reflectivity for a point source - RTV, & ! In/Output, internal variables - Scattering_Radiance_AD(1:nZ), & ! Input, AD radiances + RTV%Solar_Irradiance, & ! Input, Source irradiance at TOA + RTV%Is_Solar_Channel, & ! Input, Source sensitive channel info. + GeometryInfo%Source_Zenith_Radian, & ! Input, Source zenith angle + RTV, & ! Input, internal variables + Radiance_AD(1), & ! Input, AD radiance + AtmOptics_AD%Optical_Depth, & ! Output, AD layer optical depth Planck_Atmosphere_AD, & ! Output, AD layer radiances Planck_Surface_AD, & ! Output, AD surface radiance - AtmOptics_AD%Single_Scatter_Albedo, & ! Output, AD layer single scattering albedo - AtmOptics_AD%Optical_Depth, & ! Output, AD layer optical depth SfcOptics_AD%S_Emissivity(1:nZ), & ! Output, AD surface emissivity SfcOptics_AD%S_Reflectivity(1:nZ,1:nZ), & ! Output, AD surface reflectivity SfcOptics_AD%S_Direct_Ref(1:nZ), & ! Output, AD surface reflectivity for a point source - Pff_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer forward phase matrix - Pbb_AD(1:nZ,1:(nZ+1),:) ) ! Output, AD layer backward phase matrix - ELSE - CALL CRTM_Emission_AD( & + down_rad_AD_in=Radiance_AD(1) ) ! Optional input, AD downwelling radiance + ELSE + CALL CRTM_Emission_AD( & Atmosphere%n_Layers, & ! Input, number of atmospheric layers RTV%n_Angles, & ! Input, number of discrete zenith angles GeometryInfo%Cosine_Sensor_Zenith, & ! Input, cosine of sensor zenith angle @@ -1031,6 +1167,7 @@ FUNCTION CRTM_Compute_RTSolution_AD( & SfcOptics_AD%S_Emissivity(1:nZ), & ! Output, AD surface emissivity SfcOptics_AD%S_Reflectivity(1:nZ,1:nZ), & ! Output, AD surface reflectivity SfcOptics_AD%S_Direct_Ref(1:nZ) ) ! Output, AD surface reflectivity for a point source + END IF END IF CALL Reshape_Surf_Opt_AD(RTV%n_Angles, RTV%n_Stokes, SfcOptics_AD%Emissivity, SfcOptics_AD%Direct_Reflectivity, & SfcOptics_AD%Reflectivity, SfcOptics_AD%S_Emissivity, SfcOptics_AD%S_Direct_Ref, SfcOptics_AD%S_Reflectivity) @@ -1042,7 +1179,8 @@ FUNCTION CRTM_Compute_RTSolution_AD( & CASE (RT_ADA, RT_VMOM) ! NESDIS advanced adding-doubling method - CALL CRTM_ADA_AD( & + IF ( RTV%obs_4_downward%rt ) THEN + CALL CRTM_ADA_AD( & Atmosphere%n_Layers, & ! Input, number of atmospheric layers AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth @@ -1059,27 +1197,69 @@ FUNCTION CRTM_Compute_RTSolution_AD( & SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1), & ! Output, AD surface reflectivity SfcOptics_AD%Direct_Reflectivity(1:nZ,1), & ! Output, AD surface reflectivity for a point source Pff_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer forward phase matrix - Pbb_AD(1:nZ,1:(nZ+1),:) ) ! Output, AD layer backward phase matrix - - CASE (RT_SOI) - ! UW SOI RT solver - CALL CRTM_SOI_AD( & + Pbb_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer backward phase matrix + down_rad_AD_in=Down_Scatter_AD(1:nZ) ) ! Optional input, AD downwelling radiance + ELSE + CALL CRTM_ADA_AD( & Atmosphere%n_Layers, & ! Input, number of atmospheric layers AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + RTV%Cosmic_Background_Radiance, & ! Input, cosmic background radiation SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity - SfcOptics%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, FWD surface reflectivity - SfcOptics%Index_Sat_Ang, & ! Input, Satellite angle index + SfcOptics%Direct_Reflectivity(1:nZ,1), & ! Input, FWD surface reflectivity for a point source RTV, & ! In/Output, internal variables Scattering_Radiance_AD(1:nZ), & ! Input, AD radiances - Planck_Atmosphere_AD, & ! Output AD atmospheric layer Planck radiance - Planck_Surface_AD, & ! Output AD surface Planck radiance + Planck_Atmosphere_AD, & ! Output, AD layer radiances + Planck_Surface_AD, & ! Output, AD surface radiance AtmOptics_AD%Single_Scatter_Albedo, & ! Output, AD layer single scattering albedo AtmOptics_AD%Optical_Depth, & ! Output, AD layer optical depth SfcOptics_AD%Emissivity(1:nZ,1), & ! Output, AD surface emissivity SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1), & ! Output, AD surface reflectivity + SfcOptics_AD%Direct_Reflectivity(1:nZ,1), & ! Output, AD surface reflectivity for a point source Pff_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer forward phase matrix Pbb_AD(1:nZ,1:(nZ+1),:) ) ! Output, AD layer backward phase matrix + END IF + + CASE (RT_SOI) + ! UW SOI RT solver + IF ( RTV%obs_4_downward%rt ) THEN + CALL CRTM_SOI_AD( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo + AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity + SfcOptics%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, FWD surface reflectivity + SfcOptics%Index_Sat_Ang, & ! Input, Satellite angle index + RTV, & ! In/Output, internal variables + Scattering_Radiance_AD(1:nZ), & ! Input, AD radiances + Planck_Atmosphere_AD, & ! Output AD atmospheric layer Planck radiance + Planck_Surface_AD, & ! Output AD surface Planck radiance + AtmOptics_AD%Single_Scatter_Albedo, & ! Output, AD layer single scattering albedo + AtmOptics_AD%Optical_Depth, & ! Output, AD layer optical depth + SfcOptics_AD%Emissivity(1:nZ,1), & ! Output, AD surface emissivity + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1), & ! Output, AD surface reflectivity + Pff_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer forward phase matrix + Pbb_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer backward phase matrix + down_rad_AD_in=Down_Scatter_AD(1:nZ) ) ! Optional input, AD downwelling radiance + ELSE + CALL CRTM_SOI_AD( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + AtmOptics%Single_Scatter_Albedo, & ! Input, FWD layer single scattering albedo + AtmOptics%Optical_Depth, & ! Input, FWD layer optical depth + SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity + SfcOptics%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, FWD surface reflectivity + SfcOptics%Index_Sat_Ang, & ! Input, Satellite angle index + RTV, & ! In/Output, internal variables + Scattering_Radiance_AD(1:nZ), & ! Input, AD radiances + Planck_Atmosphere_AD, & ! Output AD atmospheric layer Planck radiance + Planck_Surface_AD, & ! Output AD surface Planck radiance + AtmOptics_AD%Single_Scatter_Albedo, & ! Output, AD layer single scattering albedo + AtmOptics_AD%Optical_Depth, & ! Output, AD layer optical depth + SfcOptics_AD%Emissivity(1:nZ,1), & ! Output, AD surface emissivity + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1), & ! Output, AD surface reflectivity + Pff_AD(1:nZ,1:(nZ+1),:), & ! Output, AD layer forward phase matrix + Pbb_AD(1:nZ,1:(nZ+1),:) ) ! Output, AD layer backward phase matrix + END IF CASE DEFAULT Error_Status = FAILURE WRITE(Message,'("Incorrect AD RT_Algorithm_ID, ",i0,", do not fit model")') & @@ -1093,7 +1273,30 @@ FUNCTION CRTM_Compute_RTSolution_AD( & ! ----------------- ! Emission model RT ! ----------------- - CALL CRTM_Emission_AD( & + IF ( RTV%obs_4_downward%rt ) THEN + CALL CRTM_Emission_AD( & + Atmosphere%n_Layers, & ! Input, number of atmospheric layers + RTV%n_Angles, & ! Input, number of discrete zenith angles + GeometryInfo%Cosine_Sensor_Zenith, & ! Input, cosine of sensor zenith angle + RTV%Planck_Atmosphere, & ! Input, FWD layer radiances + RTV%Planck_Surface, & ! Input, FWD surface radiance + SfcOptics%Emissivity(1:nZ,1), & ! Input, FWD surface emissivity + SfcOptics%Reflectivity(1:nZ,1,1:nZ,1), & ! Input, FWD surface reflectivity + SfcOptics%Direct_Reflectivity(1:nZ,1), & ! Input, FWD surface reflectivity for a point source + RTV%Solar_Irradiance, & ! Input, Source irradiance at TOA + RTV%Is_Solar_Channel, & ! Input, Source sensitive channel info. + GeometryInfo%Source_Zenith_Radian, & ! Input, Source zenith angle + RTV, & ! Input, internal variables + Radiance_AD(1), & ! Input, AD radiance + AtmOptics_AD%Optical_Depth, & ! Output, AD layer optical depth + Planck_Atmosphere_AD, & ! Output, AD layer radiances + Planck_Surface_AD, & ! Output, AD surface radiance + SfcOptics_AD%Emissivity(1:nZ,1), & ! Output, AD surface emissivity + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1), & ! Output, AD surface reflectivity + SfcOptics_AD%Direct_Reflectivity(1:nZ,1), & ! Output, AD surface reflectivity for a point source + down_rad_AD_in=Radiance_AD(1) ) ! Optional input, AD downwelling radiance + ELSE + CALL CRTM_Emission_AD( & Atmosphere%n_Layers, & ! Input, number of atmospheric layers RTV%n_Angles, & ! Input, number of discrete zenith angles GeometryInfo%Cosine_Sensor_Zenith, & ! Input, cosine of sensor zenith angle @@ -1113,6 +1316,7 @@ FUNCTION CRTM_Compute_RTSolution_AD( & SfcOptics_AD%Emissivity(1:nZ,1), & ! Output, AD surface emissivity SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1), & ! Output, AD surface reflectivity SfcOptics_AD%Direct_Reflectivity(1:nZ,1) ) ! Output, AD surface reflectivity for a point source + END IF END IF Error_Status = Assign_Common_Output_AD( Atmosphere , & ! Input diff --git a/src/RTSolution/Common_RTSolution.f90 b/src/RTSolution/Common_RTSolution.f90 index a6fdcf45..673fc65c 100644 --- a/src/RTSolution/Common_RTSolution.f90 +++ b/src/RTSolution/Common_RTSolution.f90 @@ -769,7 +769,7 @@ FUNCTION Assign_Common_Input_TL( & SfcOptics_TL%Reflectivity = ZERO User_Emissivity_TL = SfcOptics_TL%Emissivity(1,1) SfcOptics_TL%Emissivity(1,1) = ZERO - Direct_Reflectivity_TL = SfcOptics_TL%Direct_Reflectivity(1,1)/PI + Direct_Reflectivity_TL = SfcOptics_TL%Direct_Reflectivity(1,1) SfcOptics_TL%Emissivity(1:nZ,1) = User_Emissivity_TL ! Replicate the user reflectivities for all angles SfcOptics_TL%Direct_Reflectivity(1:nZ,1) = Direct_Reflectivity_TL @@ -1372,7 +1372,8 @@ FUNCTION Assign_Common_Output_TL( & SensorIndex , & ! Input ChannelIndex , & ! Input RTV , & ! Input - RTSolution_TL ) & ! Output + RTSolution_TL , & ! Output + Down_Radiance_TL ) & ! Optional input RESULT( Error_Status ) ! Arguments TYPE(CRTM_SfcOptics_type) , INTENT(IN) :: SfcOptics @@ -1384,6 +1385,7 @@ FUNCTION Assign_Common_Output_TL( & INTEGER , INTENT(IN) :: ChannelIndex TYPE(RTV_type) , INTENT(IN) :: RTV TYPE(CRTM_RTSolution_type) , INTENT(IN OUT) :: RTSolution_TL + REAL(fp) , INTENT(IN), OPTIONAL :: Down_Radiance_TL ! Function Result INTEGER :: Error_Status,n1 @@ -1399,7 +1401,11 @@ FUNCTION Assign_Common_Output_TL( & SRadiance_TL(:) = Scattering_Radiance_TL(n1:n1-1+RTV%n_Stokes) ! Emission specific assignments ELSE - SRadiance_TL(1) = Radiance_TL + IF ( RTV%obs_4_downward%rt .AND. PRESENT(Down_Radiance_TL) ) THEN + SRadiance_TL(1) = Down_Radiance_TL + ELSE + SRadiance_TL(1) = Radiance_TL + END IF END IF ! accumulate Fourier component @@ -1662,6 +1668,7 @@ FUNCTION Assign_Common_Output_AD( & CHARACTER(256) :: Message INTEGER :: no, na, nt INTEGER :: k, i + REAL(fp) :: Direct_Reflectivity_AD Error_Status = SUCCESS @@ -1709,6 +1716,14 @@ FUNCTION Assign_Common_Output_AD( & SfcOptics_AD%Angle = SfcOptics%Angle SfcOptics_AD%Weight = SfcOptics%Weight + IF ( RTV%Visible_Flag_true ) THEN + DO i = 1, nZ + IF ( SfcOptics%Direct_Reflectivity(i,1) >= ONE ) THEN + SfcOptics_AD%Direct_Reflectivity(i,1) = ZERO + END IF + END DO + END IF + ! -------------------------------------------------------------------------------- ! This part is designed for specific requirement for the total sensitivity @@ -1766,9 +1781,9 @@ FUNCTION Assign_Common_Output_AD( & 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) + -! (Direct_Reflectivity_AD/PI) + Direct_Reflectivity_AD = SUM(SfcOptics_AD%Direct_Reflectivity(1:nZ,1)) + SfcOptics_AD%Direct_Reflectivity(1:nZ,1) = ZERO + SfcOptics_AD%Direct_Reflectivity(1,1) = Direct_Reflectivity_AD RTSolution_AD%Surface_Emissivity = User_Emissivity_AD ELSE !! need to check !! diff --git a/src/RTSolution/Emission/Emission_Module.f90 b/src/RTSolution/Emission/Emission_Module.f90 index 8e7589de..3c501578 100644 --- a/src/RTSolution/Emission/Emission_Module.f90 +++ b/src/RTSolution/Emission/Emission_Module.f90 @@ -173,7 +173,7 @@ SUBROUTINE CRTM_Emission(n_Layers, & ! Input number of atmospheric layers END SUBROUTINE CRTM_Emission - SUBROUTINE CRTM_Emission_TL(n_Layers, & ! Input number of atmospheric layers + SUBROUTINE CRTM_Emission_TL(n_Layers, & ! Input number of atmospheric layers n_Angles, & ! number angles used in SfcOptics u, & ! Input cosine of local viewing angle Planck_Atmosphere, & ! Input atmospheric layer Planck radiance @@ -188,10 +188,11 @@ SUBROUTINE CRTM_Emission_TL(n_Layers, & ! Input number of atmospheric layers T_OD_TL, & ! Input tangent-linear of layer optical depth Planck_Atmosphere_TL, & ! Input TL atmospheric layer Planck radiance Planck_Surface_TL, & ! Input TL surface Planck radiance - emissivity_TL, & ! Input TL surface emissivity - reflectivity_TL, & ! Input TL surface reflectivity matrix + emissivity_TL, & ! Input TL surface emissivity + reflectivity_TL, & ! Input TL surface reflectivity matrix direct_reflectivity_TL, & ! Input TL surface ditrct reflectivity - up_rad_TL) ! Output TL TOA radiance + up_rad_TL, & ! Output TL TOA radiance + down_rad_TL_out) ! Optional output TL downwelling radiance ! --------------------------------------------------------------------------- ! ! FUNCTION: Compute tangent-linear upward radiance at the top of the ! ! atmosphere using carried results in RTV structure from forward ! @@ -208,13 +209,16 @@ SUBROUTINE CRTM_Emission_TL(n_Layers, & ! Input number of atmospheric layers REAL (fp), INTENT(IN), DIMENSION( 0: ) :: Planck_Atmosphere,Planck_Atmosphere_TL REAL (fp), INTENT(IN) :: Planck_Surface,u,Planck_Surface_TL REAL (fp), INTENT(INOUT) :: up_rad_TL + REAL (fp), INTENT(OUT), OPTIONAL :: down_rad_TL_out ! Structure RTV carried in variables from forward calculation. TYPE(RTV_type), INTENT( IN) :: RTV ! internal variables REAL (fp) :: layer_source_up_TL, layer_source_down_TL,a_TL,down_rad_TL REAL (fp) :: Total_OD, Total_OD_TL + REAL (fp) :: down_rad_TL_level INTEGER :: k + INTEGER :: obs_idx REAL( fp) :: cosine_u0 !#--------------------------------------------------------------------------# @@ -225,6 +229,8 @@ SUBROUTINE CRTM_Emission_TL(n_Layers, & ! Input number of atmospheric layers Total_OD_TL = ZERO Total_OD = RTV%Total_OD + down_rad_TL_level = ZERO + obs_idx = RTV%obs_4_downward%idx DO k = 1, n_Layers ! accumulate tangent-linear optical depth @@ -238,6 +244,7 @@ SUBROUTINE CRTM_Emission_TL(n_Layers, & ! Input number of atmospheric layers ! down_rad(k) = down_rad(k-1) * layer_trans(k) + layer_source_down down_rad_TL = down_rad_TL*RTV%e_Layer_Trans_DOWN(k) & +RTV%e_Level_Rad_DOWN(k-1)*RTV%e_Layer_Trans_DOWN(k)*a_TL+layer_source_down_TL + IF ( RTV%obs_4_downward%rt .AND. k == obs_idx ) down_rad_TL_level = down_rad_TL ENDDO !#--------------------------------------------------------------------------# @@ -273,6 +280,14 @@ SUBROUTINE CRTM_Emission_TL(n_Layers, & ! Input number of atmospheric layers +RTV%e_Level_Rad_UP(k)*RTV%e_Layer_Trans_UP(k)*a_TL+layer_source_up_TL ENDDO ! + IF ( PRESENT(down_rad_TL_out) ) THEN + IF ( RTV%obs_4_downward%rt ) THEN + down_rad_TL_out = down_rad_TL_level + ELSE + down_rad_TL_out = ZERO + END IF + END IF + RETURN END SUBROUTINE CRTM_Emission_TL ! @@ -295,7 +310,8 @@ SUBROUTINE CRTM_Emission_AD(n_Layers, & ! Input number of atmospheric layers Planck_Surface_AD, & ! Output AD surface Planck radiance emissivity_AD, & ! Output AD surface emissivity reflectivity_AD, & ! Output AD surface reflectivity matrix - direct_reflectivity_AD) ! Output AD surface direct reflectivity + direct_reflectivity_AD, & ! Output AD surface direct reflectivity + down_rad_AD_in) ! Optional input AD downwelling radiance ! --------------------------------------------------------------------------- ! ! FUNCTION: Compute adjoint upward radiance at the top of the ! ! atmosphere using carried results in RTV structure from forward ! @@ -312,6 +328,7 @@ SUBROUTINE CRTM_Emission_AD(n_Layers, & ! Input number of atmospheric layers REAL (fp), INTENT(IN), DIMENSION( 0: ) :: Planck_Atmosphere REAL (fp), INTENT(IN) :: Planck_Surface,u REAL (fp), INTENT(IN) :: up_rad_AD_in + REAL (fp), INTENT(IN), OPTIONAL :: down_rad_AD_in REAL (fp), INTENT(IN OUT), DIMENSION( : ) :: T_OD_AD,emissivity_AD REAL (fp), INTENT(IN OUT), DIMENSION( :,: ) :: reflectivity_AD REAL (fp), INTENT(IN OUT), DIMENSION( : ) :: direct_reflectivity_AD @@ -322,6 +339,7 @@ SUBROUTINE CRTM_Emission_AD(n_Layers, & ! Input number of atmospheric layers REAL (fp) :: layer_source_up_AD, layer_source_down_AD,a_AD,down_rad_AD REAL (fp) :: cosine_u0, up_rad_AD, Total_OD, Total_OD_AD INTEGER :: k + INTEGER :: obs_idx ! ! Initialize variables Total_OD_AD = ZERO @@ -335,6 +353,24 @@ SUBROUTINE CRTM_Emission_AD(n_Layers, & ! Input number of atmospheric layers ! Total column optical depth carried from forward part Total_OD = RTV%Total_OD + obs_idx = RTV%obs_4_downward%idx + + IF ( PRESENT(down_rad_AD_in) ) THEN + down_rad_AD = down_rad_AD_in + DO k = obs_idx, 1, -1 + a_AD = RTV%e_Level_Rad_DOWN(k-1)*RTV%e_Layer_Trans_DOWN(k)*down_rad_AD + layer_source_down_AD = down_rad_AD + down_rad_AD = down_rad_AD*RTV%e_Layer_Trans_DOWN(k) + + Planck_Atmosphere_AD(k) = Planck_Atmosphere_AD(k) + layer_source_down_AD * & + (ONE - RTV%e_Layer_Trans_DOWN(k)) + a_AD = a_AD - Planck_Atmosphere(k) * RTV%e_Layer_Trans_DOWN(k)* layer_source_down_AD + + T_OD_AD(k) = T_OD_AD(k) - a_AD * RTV%Secant_Down_Angle + ENDDO + down_rad_AD = ZERO + RETURN + END IF !#--------------------------------------------------------------------------# !# -- Upwelling adjoint radiance -- # diff --git a/src/RTSolution/SOI/SOI_Module.f90 b/src/RTSolution/SOI/SOI_Module.f90 index 8cca3043..d59dff79 100644 --- a/src/RTSolution/SOI/SOI_Module.f90 +++ b/src/RTSolution/SOI/SOI_Module.f90 @@ -90,7 +90,7 @@ SUBROUTINE CRTM_SOI(n_Layers, & ! Input number of atmospheric layers ! -------------- internal variables --------------------------------- ! - INTEGER :: i, k, niter + INTEGER :: i, k, niter, iter REAL(fp) :: rad, rad_change REAL(fp), PARAMETER :: initial_error = 1.E10 REAL(fp), PARAMETER :: SMALL = 1.E-15 @@ -277,6 +277,17 @@ SUBROUTINE CRTM_SOI(n_Layers, & ! Input number of atmospheric layers END DO soi_loop + IF ( RTV%obs_4_downward%rt ) THEN + RTV%s_Level_Rad_DOWN( 1 : RTV%n_Angles, 0 : n_Layers ) = ZERO + DO iter = 1, niter + RTV%s_Level_Rad_DOWN( 1 : RTV%n_Angles, 0 : n_Layers ) = & + RTV%s_Level_Rad_DOWN( 1 : RTV%n_Angles, 0 : n_Layers ) + & + RTV%s_Level_IterRad_DOWN( 1 : RTV%n_Angles, 0 : n_Layers, iter ) + END DO + RTV%s_Level_Rad_DOWNT( 1 : RTV%n_Angles, 0 : n_Layers ) = & + RTV%s_Level_Rad_DOWN( 1 : RTV%n_Angles, 0 : n_Layers ) + END IF + RTV%Number_SOI_Iter = niter RETURN @@ -298,7 +309,8 @@ SUBROUTINE CRTM_SOI_TL(n_Layers, & ! Input number of atmospheric layers reflectivity_TL, & ! Input TL reflectivity Pff_TL, & ! Input TL forward phase matrix Pbb_TL, & ! Input TL backward phase matrix - s_rad_up_TL) ! Output TL upward radiance + s_rad_up_TL, & ! Output TL upward radiance + down_rad_TL_out) ! Optional output, TL downwelling radiance ! ------------------------------------------------------------------------- ! ! ! ! FUNCTION: ! @@ -323,11 +335,14 @@ SUBROUTINE CRTM_SOI_TL(n_Layers, & ! Input number of atmospheric layers REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity_TL REAL (fp), INTENT(IN), DIMENSION( :, : ) :: reflectivity_TL REAL (fp), INTENT(IN), DIMENSION( :, :, : ) :: Pff_TL, Pbb_TL - REAL (fp), INTENT(INOUT), DIMENSION( : ) :: s_rad_up_TL + REAL (fp), INTENT(INOUT), DIMENSION( : ) :: s_rad_up_TL + REAL (fp), INTENT(OUT), OPTIONAL, DIMENSION( : ) :: down_rad_TL_out ! -------------- internal variables --------------------------------- ! INTEGER :: i, k, iter + INTEGER :: obs_idx + LOGICAL :: compute_down REAL(fp), PARAMETER :: SNGL_SCAT_ALB_THRESH = 0.8 REAL(fp), PARAMETER :: OPT_DEPTH_THRESH = 4.0 REAL(fp), DIMENSION( RTV%n_Angles ) :: source_TL @@ -335,6 +350,8 @@ SUBROUTINE CRTM_SOI_TL(n_Layers, & ! Input number of atmospheric layers REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: s_source_up_TL, s_source_down_TL REAL(fp), DIMENSION( RTV%n_Angles, RTV%n_Angles, n_Layers ) :: s_trans_TL, s_refl_TL REAL(fp), DIMENSION( RTV%n_Angles, 0:n_Layers, RTV%Number_SOI_Iter ) :: s_IterRad_UP_TL, s_IterRad_DOWN_TL + compute_down = PRESENT(down_rad_TL_out) .AND. RTV%obs_4_downward%rt + obs_idx = RTV%obs_4_downward%idx DO k = 1, n_Layers @@ -467,10 +484,20 @@ SUBROUTINE CRTM_SOI_TL(n_Layers, & ! Input number of atmospheric layers !---------------- ! Add up terms !---------------- - s_Rad_UP_TL( Index_Sat_Angle ) = s_Rad_UP_TL( Index_Sat_Angle ) + s_IterRad_UP_TL( Index_Sat_Angle, 0, iter ) + s_Rad_UP_TL( Index_Sat_Angle ) = s_Rad_UP_TL( Index_Sat_Angle ) + s_IterRad_UP_TL( Index_Sat_Angle, 0, iter ) END DO + IF ( compute_down ) THEN + down_rad_TL_out = ZERO + IF ( obs_idx > 0 ) THEN + DO iter = 1, RTV%Number_SOI_Iter + down_rad_TL_out( 1 : RTV%n_Angles ) = down_rad_TL_out( 1 : RTV%n_Angles ) + & + s_IterRad_DOWN_TL( 1 : RTV%n_Angles, obs_idx, iter ) + END DO + END IF + END IF + RETURN END SUBROUTINE CRTM_SOI_TL @@ -489,7 +516,8 @@ SUBROUTINE CRTM_SOI_AD(n_Layers, & ! Input number of atmospheric layers emissivity_AD, & ! Output AD surface emissivity reflectivity_AD, & ! Output AD surface reflectivity Pff_AD, & ! Output AD forward phase matrix - Pbb_AD) ! Output AD backward phase matrix + Pbb_AD, & ! Output AD backward phase matrix + down_rad_AD_in) ! Optional input, AD downwelling radiance ! ------------------------------------------------------------------------- ! ! FUNCTION: ! ! This subroutine calculates IR/MW adjoint radiance at the top of ! @@ -511,12 +539,15 @@ SUBROUTINE CRTM_SOI_AD(n_Layers, & ! Input number of atmospheric layers REAL (fp),INTENT(INOUT) :: Planck_Surface_AD REAL (fp),INTENT(INOUT),DIMENSION( : ) :: emissivity_AD REAL (fp),INTENT(INOUT),DIMENSION( :, : ) :: reflectivity_AD - REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_AD + REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_AD + REAL (fp), INTENT(IN), OPTIONAL, DIMENSION( : ) :: down_rad_AD_in ! Local variables REAL(fp), PARAMETER :: SNGL_SCAT_ALB_THRESH = 0.8 REAL(fp), PARAMETER :: OPT_DEPTH_THRESH = 4.0 INTEGER :: iter, k, i, j + INTEGER :: obs_idx + LOGICAL :: compute_down REAL(fp), DIMENSION( RTV%n_Angles, 0:n_Layers, RTV%Number_SOI_Iter ) :: s_IterRad_UP_AD REAL(fp), DIMENSION( RTV%n_Angles, 0:n_Layers, RTV%Number_SOI_Iter ) :: s_IterRad_DOWN_AD REAL(fp), DIMENSION( RTV%n_Angles ) :: source_AD @@ -538,6 +569,18 @@ SUBROUTINE CRTM_SOI_AD(n_Layers, & ! Input number of atmospheric layers Pff_AD = ZERO Pbb_AD = ZERO + compute_down = PRESENT(down_rad_AD_in) .AND. RTV%obs_4_downward%rt + obs_idx = RTV%obs_4_downward%idx + + IF ( compute_down ) THEN + IF ( obs_idx > 0 ) THEN + DO iter = 1, RTV%Number_SOI_Iter + s_IterRad_DOWN_AD( 1 : RTV%n_Angles, obs_idx, iter ) = & + s_IterRad_DOWN_AD( 1 : RTV%n_Angles, obs_idx, iter ) + down_rad_AD_in( 1 : RTV%n_Angles ) + END DO + END IF + END IF + !-------------------------------------------------------------------- ! Loop through each successive order of interaction in reverse order @@ -1482,4 +1525,3 @@ SUBROUTINE CRTM_Doubling_layer_AD(n_streams, & ! Input, number of streams END SUBROUTINE CRTM_Doubling_layer_AD END MODULE SOI_Module - diff --git a/src/SfcOptics/CRTM_SfcOptics.f90 b/src/SfcOptics/CRTM_SfcOptics.f90 index ea3565c2..d72f4310 100644 --- a/src/SfcOptics/CRTM_SfcOptics.f90 +++ b/src/SfcOptics/CRTM_SfcOptics.f90 @@ -1712,13 +1712,126 @@ FUNCTION CRTM_Compute_SfcOptics_TL( & ELSE IF ( SpcCoeff_IsVisibleSensor( SC(SensorIndex) ) ) THEN + mth_Azi_Test: IF( SfcOptics%mth_Azi == 0 ) THEN + + ! ================== + ! Lambertian surface + ! ================== + + ! ------------------------------------- + ! Visible LAND emissivity/reflectivity + ! ------------------------------------- + Visible_Land: IF( Surface%Land_Coverage > ZERO ) THEN + + ! Compute the surface optics + Error_Status = Compute_VIS_Land_SfcOptics_TL( SfcOptics_TL ) + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'("Error computing VIS land SfcOptics_TL at ",& + &"channel index ",i0)' ) ChannelIndex + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN + END IF + + ! Accumulate the surface optics properties + ! based on land coverage fraction + Emissivity_TL(1:nZ,1) = & + SfcOptics_TL%Emissivity(1:nZ,1) * Surface%Land_Coverage + Reflectivity_TL(1:nZ,1,1:nZ,1) = & + SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1) * Surface%Land_Coverage + Direct_Reflectivity_TL(1:nZ,1) = & + SfcOptics_TL%Direct_Reflectivity(1:nZ,1) * Surface%Land_Coverage + END IF Visible_Land + + + ! ------------------------------------- + ! Visible WATER emissivity/reflectivity + ! ------------------------------------- + Visible_Water: IF( Surface%Water_Coverage > ZERO ) THEN + + ! Compute the surface optics + Error_Status = Compute_VIS_Water_SfcOptics_TL( SfcOptics_TL ) + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'("Error computing VIS water SfcOptics_TL at ",& + &"channel index ",i0)' ) ChannelIndex + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN + END IF + + ! Accumulate the surface optics properties + ! based on water coverage fraction + Emissivity_TL(1:nZ,1) = Emissivity_TL(1:nZ,1) + & + ( SfcOptics_TL%Emissivity(1:nZ,1) * Surface%Water_Coverage ) + Reflectivity_TL(1:nZ,1,1:nZ,1) = Reflectivity_TL(1:nZ,1,1:nZ,1) + & + ( SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1) * Surface%Water_Coverage ) + Direct_Reflectivity_TL(1:nZ,1) = Direct_Reflectivity_TL(1:nZ,1) + & + ( SfcOptics_TL%Direct_Reflectivity(1:nZ,1) * Surface%Water_Coverage ) + END IF Visible_Water + + + ! ------------------------------------ + ! Visible SNOW emissivity/reflectivity + ! ------------------------------------ + Visible_Snow: IF( Surface%Snow_Coverage > ZERO ) THEN + + ! Compute the surface optics + Error_Status = Compute_VIS_Snow_SfcOptics_TL( SfcOptics_TL ) + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'("Error computing VIS snow SfcOptics_TL at ",& + &"channel index ",i0)' ) ChannelIndex + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN + END IF + + ! Accumulate the surface optics properties + ! based on snow coverage fraction + Emissivity_TL(1:nZ,1) = Emissivity_TL(1:nZ,1) + & + ( SfcOptics_TL%Emissivity(1:nZ,1) * Surface%Snow_Coverage ) + Reflectivity_TL(1:nZ,1,1:nZ,1) = Reflectivity_TL(1:nZ,1,1:nZ,1) + & + ( SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1) * Surface%Snow_Coverage ) + Direct_Reflectivity_TL(1:nZ,1) = Direct_Reflectivity_TL(1:nZ,1) + & + ( SfcOptics_TL%Direct_Reflectivity(1:nZ,1) * Surface%Snow_Coverage ) + END IF Visible_Snow + + + ! ----------------------------------- + ! Visible ICE emissivity/reflectivity + ! ----------------------------------- + Visible_Ice: IF( Surface%Ice_Coverage > ZERO ) THEN + + ! Compute the surface optics + Error_Status = Compute_VIS_Ice_SfcOptics_TL( SfcOptics_TL ) + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'("Error computing VIS ice SfcOptics_TL at ",& + &"channel index ",i0)' ) ChannelIndex + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN + END IF - ! ------------------- - ! Default values only - ! ------------------- - SfcOptics_TL%Emissivity(1:nZ,1) = ZERO - SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1) = ZERO - SfcOptics_TL%Direct_Reflectivity = ZERO + ! Accumulate the surface optics properties + ! based on ice coverage fraction + Emissivity_TL(1:nZ,1) = Emissivity_TL(1:nZ,1) + & + ( SfcOptics_TL%Emissivity(1:nZ,1) * Surface%Ice_Coverage ) + Reflectivity_TL(1:nZ,1,1:nZ,1) = Reflectivity_TL(1:nZ,1,1:nZ,1) + & + ( SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1) * Surface%Ice_Coverage ) + Direct_Reflectivity_TL(1:nZ,1) = Direct_Reflectivity_TL(1:nZ,1) + & + ( SfcOptics_TL%Direct_Reflectivity(1:nZ,1) * Surface%Ice_Coverage ) + END IF Visible_Ice + + + ! ----------------------- + ! Assign the final result + ! ----------------------- + SfcOptics_TL%Emissivity(1:nZ,1) = Emissivity_TL(1:nZ,1) + SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1) = Reflectivity_TL(1:nZ,1,1:nZ,1) + SfcOptics_TL%Direct_Reflectivity(1:nZ,1) = Direct_Reflectivity_TL(1:nZ,1) + + ELSE + + SfcOptics_TL%Emissivity(1:nZ,1) = ZERO + SfcOptics_TL%Reflectivity(1:nZ,1,1:nZ,1) = ZERO + SfcOptics_TL%Direct_Reflectivity = ZERO + + END IF mth_Azi_Test !########################################################################## @@ -1939,22 +2052,22 @@ FUNCTION CRTM_Compute_SfcOptics_AD( & ! r = (rV + rH)/2 ! Note: INTENSITY == UNPOLARIZED == FIRST_STOKES_COMPONENT CASE( INTENSITY ) - Emissivity_AD(1:nZ,1) = SfcOptics_AD%Emissivity(1:nZ,1) - Emissivity_AD(1:nZ,2) = SfcOptics_AD%Emissivity(1:nZ,1) + Emissivity_AD(1:nZ,1) = POINT_5 * SfcOptics_AD%Emissivity(1:nZ,1) + Emissivity_AD(1:nZ,2) = POINT_5 * SfcOptics_AD%Emissivity(1:nZ,1) SfcOptics_AD%Emissivity = ZERO - Reflectivity_AD(1:nZ,1,1:nZ,1) = SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) - Reflectivity_AD(1:nZ,2,1:nZ,2) = SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) + Reflectivity_AD(1:nZ,1,1:nZ,1) = POINT_5 * SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) + Reflectivity_AD(1:nZ,2,1:nZ,2) = POINT_5 * SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) SfcOptics_AD%Reflectivity = ZERO ! The second Stokes component, Q, the polarisation difference. ! e = (eV - eH)/2 ! r = (rV - rH)/2 CASE( SECOND_STOKES_COMPONENT ) - Emissivity_AD(1:nZ,1) = SfcOptics_AD%Emissivity(1:nZ,1) - Emissivity_AD(1:nZ,2) = -SfcOptics_AD%Emissivity(1:nZ,1) + Emissivity_AD(1:nZ,1) = POINT_5 * SfcOptics_AD%Emissivity(1:nZ,1) + Emissivity_AD(1:nZ,2) = -POINT_5 * SfcOptics_AD%Emissivity(1:nZ,1) SfcOptics_AD%Emissivity = ZERO - Reflectivity_AD(1:nZ,1,1:nZ,1) = SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) - Reflectivity_AD(1:nZ,2,1:nZ,2) = -SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) + Reflectivity_AD(1:nZ,1,1:nZ,1) = POINT_5 * SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) + Reflectivity_AD(1:nZ,2,1:nZ,2) = -POINT_5 * SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) SfcOptics_AD%Reflectivity = ZERO ! The third Stokes component, U. @@ -2395,13 +2508,98 @@ FUNCTION CRTM_Compute_SfcOptics_AD( & ELSE IF ( SpcCoeff_IsVisibleSensor( SC(SensorIndex) ) ) THEN + mth_Azi_Test: IF( SfcOptics%mth_Azi == 0 ) THEN + + Emissivity_AD(1:nZ,1) = SfcOptics_AD%Emissivity(1:nZ,1) + Reflectivity_AD(1:nZ,1,1:nZ,1) = SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) + Direct_Reflectivity_AD(1:nZ,1) = SfcOptics_AD%Direct_Reflectivity(1:nZ,1) + SfcOptics_AD%Emissivity = ZERO + SfcOptics_AD%Reflectivity = ZERO + SfcOptics_AD%Direct_Reflectivity = ZERO + + ! ----------------------------------- + ! Visible ICE emissivity/reflectivity + ! ----------------------------------- + Visible_Ice: IF( Surface%Ice_Coverage > ZERO ) THEN + SfcOptics_AD%Emissivity(1:nZ,1) = & + SfcOptics_AD%Emissivity(1:nZ,1) + (Emissivity_AD(1:nZ,1) * Surface%Ice_Coverage) + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) = & + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) + (Reflectivity_AD(1:nZ,1,1:nZ,1) * Surface%Ice_Coverage) + SfcOptics_AD%Direct_Reflectivity(1:nZ,1) = & + SfcOptics_AD%Direct_Reflectivity(1:nZ,1) + (Direct_Reflectivity_AD(1:nZ,1) * Surface%Ice_Coverage) + Error_Status = Compute_VIS_Ice_SfcOptics_AD( SfcOptics_AD ) + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'("Error computing VIS ice SfcOptics_AD at ",& + &"channel index ",i0)' ) ChannelIndex + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN + END IF + END IF Visible_Ice + + ! ------------------------------------ + ! Visible SNOW emissivity/reflectivity + ! ------------------------------------ + Visible_Snow: IF( Surface%Snow_Coverage > ZERO ) THEN + SfcOptics_AD%Emissivity(1:nZ,1) = & + SfcOptics_AD%Emissivity(1:nZ,1) + (Emissivity_AD(1:nZ,1) * Surface%Snow_Coverage) + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) = & + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) + (Reflectivity_AD(1:nZ,1,1:nZ,1) * Surface%Snow_Coverage) + SfcOptics_AD%Direct_Reflectivity(1:nZ,1) = & + SfcOptics_AD%Direct_Reflectivity(1:nZ,1) + (Direct_Reflectivity_AD(1:nZ,1) * Surface%Snow_Coverage) + Error_Status = Compute_VIS_Snow_SfcOptics_AD( SfcOptics_AD ) + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'("Error computing VIS snow SfcOptics_AD at ",& + &"channel index ",i0)' ) ChannelIndex + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN + END IF + END IF Visible_Snow + + ! ------------------------------------- + ! Visible WATER emissivity/reflectivity + ! ------------------------------------- + Visible_Water: IF( Surface%Water_Coverage > ZERO ) THEN + SfcOptics_AD%Emissivity(1:nZ,1) = & + SfcOptics_AD%Emissivity(1:nZ,1) + (Emissivity_AD(1:nZ,1) * Surface%Water_Coverage) + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) = & + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) + (Reflectivity_AD(1:nZ,1,1:nZ,1) * Surface%Water_Coverage) + SfcOptics_AD%Direct_Reflectivity(1:nZ,1) = & + SfcOptics_AD%Direct_Reflectivity(1:nZ,1) + (Direct_Reflectivity_AD(1:nZ,1) * Surface%Water_Coverage) + Error_Status = Compute_VIS_Water_SfcOptics_AD( SfcOptics_AD ) + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'("Error computing VIS water SfcOptics_AD at ",& + &"channel index ",i0)' ) ChannelIndex + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN + END IF + END IF Visible_Water + + ! ------------------------------------- + ! Visible LAND emissivity/reflectivity + ! ------------------------------------- + Visible_Land: IF( Surface%Land_Coverage > ZERO ) THEN + SfcOptics_AD%Emissivity(1:nZ,1) = & + SfcOptics_AD%Emissivity(1:nZ,1) + (Emissivity_AD(1:nZ,1) * Surface%Land_Coverage) + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) = & + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) + (Reflectivity_AD(1:nZ,1,1:nZ,1) * Surface%Land_Coverage) + SfcOptics_AD%Direct_Reflectivity(1:nZ,1) = & + SfcOptics_AD%Direct_Reflectivity(1:nZ,1) + (Direct_Reflectivity_AD(1:nZ,1) * Surface%Land_Coverage) + Error_Status = Compute_VIS_Land_SfcOptics_AD( SfcOptics_AD ) + IF ( Error_Status /= SUCCESS ) THEN + WRITE( Message,'("Error computing VIS land SfcOptics_AD at ",& + &"channel index ",i0)' ) ChannelIndex + CALL Display_Message( ROUTINE_NAME, Message, Error_Status ) + RETURN + END IF + END IF Visible_Land - ! ------------------- - ! Default values only - ! ------------------- - SfcOptics_AD%Emissivity(1:nZ,1) = ZERO - SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) = ZERO - SfcOptics_AD%Direct_Reflectivity = ZERO + ELSE + + SfcOptics_AD%Emissivity(1:nZ,1) = ZERO + SfcOptics_AD%Reflectivity(1:nZ,1,1:nZ,1) = ZERO + SfcOptics_AD%Direct_Reflectivity = ZERO + + END IF mth_Azi_Test !########################################################################## diff --git a/src/SfcOptics/CRTM_VIS_Ice_SfcOptics.f90 b/src/SfcOptics/CRTM_VIS_Ice_SfcOptics.f90 index a1eaaf4d..393f80bf 100644 --- a/src/SfcOptics/CRTM_VIS_Ice_SfcOptics.f90 +++ b/src/SfcOptics/CRTM_VIS_Ice_SfcOptics.f90 @@ -247,17 +247,22 @@ FUNCTION Compute_VIS_Ice_SfcOptics_TL( & ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_VIS_Ice_SfcOptics_TL' ! Local variables + INTEGER :: j + INTEGER :: nZ + REAL(fp) :: emissivity_TL ! Set up err_stat = SUCCESS - ! Compute the tangent-linear surface optical parameters - ! ***No TL models yet, so default TL output is zero*** - SfcOptics_TL%Reflectivity = ZERO - SfcOptics_TL%Direct_Reflectivity = ZERO - SfcOptics_TL%Emissivity = ZERO + nZ = SfcOptics_TL%n_Angles + emissivity_TL = SfcOptics_TL%Emissivity(1,1) + SfcOptics_TL%Direct_Reflectivity(1:nZ,1) = -emissivity_TL + SfcOptics_TL%Emissivity(1:nZ,1) = emissivity_TL + DO j = 1, nZ + SfcOptics_TL%Reflectivity(1:nZ,1,j,1) = -emissivity_TL * SfcOptics_TL%Weight(j) + END DO END FUNCTION Compute_VIS_Ice_SfcOptics_TL @@ -318,17 +323,27 @@ FUNCTION Compute_VIS_Ice_SfcOptics_AD( & ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_VIS_Ice_SfcOptics_AD' ! Local variables + INTEGER :: j + INTEGER :: nZ + REAL(fp) :: emissivity_AD ! Set up err_stat = SUCCESS - ! Compute the adjoint surface optical parameters - ! ***No AD models yet, so there is no impact on AD result*** + nZ = SfcOptics_AD%n_Angles + emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) - & + SUM(SfcOptics_AD%Direct_Reflectivity(1:nZ,1)) + DO j = 1, nZ + emissivity_AD = emissivity_AD - & + ( SfcOptics_AD%Weight(j) * SUM(SfcOptics_AD%Reflectivity(1:nZ,1,j,1)) ) + END DO + SfcOptics_AD%Reflectivity = ZERO SfcOptics_AD%Direct_Reflectivity = ZERO SfcOptics_AD%Emissivity = ZERO + SfcOptics_AD%Emissivity(1,1) = emissivity_AD END FUNCTION Compute_VIS_Ice_SfcOptics_AD diff --git a/src/SfcOptics/CRTM_VIS_Land_SfcOptics.f90 b/src/SfcOptics/CRTM_VIS_Land_SfcOptics.f90 index fa180453..e0a20759 100644 --- a/src/SfcOptics/CRTM_VIS_Land_SfcOptics.f90 +++ b/src/SfcOptics/CRTM_VIS_Land_SfcOptics.f90 @@ -248,17 +248,22 @@ FUNCTION Compute_VIS_Land_SfcOptics_TL( & ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_VIS_Land_SfcOptics_TL' ! Local variables + INTEGER :: j + INTEGER :: nZ + REAL(fp) :: emissivity_TL ! Set up err_stat = SUCCESS - ! Compute the tangent-linear surface optical parameters - ! ***No TL models yet, so default TL output is zero*** - SfcOptics_TL%Reflectivity = ZERO - SfcOptics_TL%Direct_Reflectivity = ZERO - SfcOptics_TL%Emissivity = ZERO + nZ = SfcOptics_TL%n_Angles + emissivity_TL = SfcOptics_TL%Emissivity(1,1) + SfcOptics_TL%Direct_Reflectivity(1:nZ,1) = -emissivity_TL + SfcOptics_TL%Emissivity(1:nZ,1) = emissivity_TL + DO j = 1, nZ + SfcOptics_TL%Reflectivity(1:nZ,1,j,1) = -emissivity_TL * SfcOptics_TL%Weight(j) + END DO END FUNCTION Compute_VIS_Land_SfcOptics_TL @@ -319,17 +324,27 @@ FUNCTION Compute_VIS_Land_SfcOptics_AD( & ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_VIS_Land_SfcOptics_AD' ! Local variables + INTEGER :: j + INTEGER :: nZ + REAL(fp) :: emissivity_AD ! Set up err_stat = SUCCESS - ! Compute the adjoint surface optical parameters - ! ***No AD models yet, so there is no impact on AD result*** + nZ = SfcOptics_AD%n_Angles + emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) - & + SUM(SfcOptics_AD%Direct_Reflectivity(1:nZ,1)) + DO j = 1, nZ + emissivity_AD = emissivity_AD - & + ( SfcOptics_AD%Weight(j) * SUM(SfcOptics_AD%Reflectivity(1:nZ,1,j,1)) ) + END DO + SfcOptics_AD%Reflectivity = ZERO SfcOptics_AD%Direct_Reflectivity = ZERO SfcOptics_AD%Emissivity = ZERO + SfcOptics_AD%Emissivity(1,1) = emissivity_AD END FUNCTION Compute_VIS_Land_SfcOptics_AD diff --git a/src/SfcOptics/CRTM_VIS_Snow_SfcOptics.f90 b/src/SfcOptics/CRTM_VIS_Snow_SfcOptics.f90 index 0496672c..fa62a5b8 100644 --- a/src/SfcOptics/CRTM_VIS_Snow_SfcOptics.f90 +++ b/src/SfcOptics/CRTM_VIS_Snow_SfcOptics.f90 @@ -248,17 +248,22 @@ FUNCTION Compute_VIS_Snow_SfcOptics_TL( & ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_VIS_Snow_SfcOptics_TL' ! Local variables + INTEGER :: j + INTEGER :: nZ + REAL(fp) :: emissivity_TL ! Set up err_stat = SUCCESS - ! Compute the tangent-linear surface optical parameters - ! ***No TL models yet, so default TL output is zero*** - SfcOptics_TL%Reflectivity = ZERO - SfcOptics_TL%Direct_Reflectivity = ZERO - SfcOptics_TL%Emissivity = ZERO + nZ = SfcOptics_TL%n_Angles + emissivity_TL = SfcOptics_TL%Emissivity(1,1) + SfcOptics_TL%Direct_Reflectivity(1:nZ,1) = -emissivity_TL + SfcOptics_TL%Emissivity(1:nZ,1) = emissivity_TL + DO j = 1, nZ + SfcOptics_TL%Reflectivity(1:nZ,1,j,1) = -emissivity_TL * SfcOptics_TL%Weight(j) + END DO END FUNCTION Compute_VIS_Snow_SfcOptics_TL @@ -319,17 +324,27 @@ FUNCTION Compute_VIS_Snow_SfcOptics_AD( & ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_VIS_Snow_SfcOptics_AD' ! Local variables + INTEGER :: j + INTEGER :: nZ + REAL(fp) :: emissivity_AD ! Set up err_stat = SUCCESS - ! Compute the adjoint surface optical parameters - ! ***No AD models yet, so there is no impact on AD result*** + nZ = SfcOptics_AD%n_Angles + emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) - & + SUM(SfcOptics_AD%Direct_Reflectivity(1:nZ,1)) + DO j = 1, nZ + emissivity_AD = emissivity_AD - & + ( SfcOptics_AD%Weight(j) * SUM(SfcOptics_AD%Reflectivity(1:nZ,1,j,1)) ) + END DO + SfcOptics_AD%Reflectivity = ZERO SfcOptics_AD%Direct_Reflectivity = ZERO SfcOptics_AD%Emissivity = ZERO + SfcOptics_AD%Emissivity(1,1) = emissivity_AD END FUNCTION Compute_VIS_Snow_SfcOptics_AD diff --git a/src/SfcOptics/CRTM_VIS_Water_SfcOptics.f90 b/src/SfcOptics/CRTM_VIS_Water_SfcOptics.f90 index 3b676f1d..969b742a 100644 --- a/src/SfcOptics/CRTM_VIS_Water_SfcOptics.f90 +++ b/src/SfcOptics/CRTM_VIS_Water_SfcOptics.f90 @@ -248,17 +248,22 @@ FUNCTION Compute_VIS_Water_SfcOptics_TL( & ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_VIS_Water_SfcOptics_TL' ! Local variables + INTEGER :: j + INTEGER :: nZ + REAL(fp) :: emissivity_TL ! Set up err_stat = SUCCESS - ! Compute the tangent-linear surface optical parameters - ! ***No TL models yet, so default TL output is zero*** - SfcOptics_TL%Reflectivity = ZERO - SfcOptics_TL%Direct_Reflectivity = ZERO - SfcOptics_TL%Emissivity = ZERO + nZ = SfcOptics_TL%n_Angles + emissivity_TL = SfcOptics_TL%Emissivity(1,1) + SfcOptics_TL%Direct_Reflectivity(1:nZ,1) = -emissivity_TL + SfcOptics_TL%Emissivity(1:nZ,1) = emissivity_TL + DO j = 1, nZ + SfcOptics_TL%Reflectivity(1:nZ,1,j,1) = -emissivity_TL * SfcOptics_TL%Weight(j) + END DO END FUNCTION Compute_VIS_Water_SfcOptics_TL @@ -319,17 +324,27 @@ FUNCTION Compute_VIS_Water_SfcOptics_AD( & ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Compute_VIS_Water_SfcOptics_AD' ! Local variables + INTEGER :: j + INTEGER :: nZ + REAL(fp) :: emissivity_AD ! Set up err_stat = SUCCESS - ! Compute the adjoint surface optical parameters - ! ***No AD models yet, so there is no impact on AD result*** + nZ = SfcOptics_AD%n_Angles + emissivity_AD = SUM(SfcOptics_AD%Emissivity(1:nZ,1)) - & + SUM(SfcOptics_AD%Direct_Reflectivity(1:nZ,1)) + DO j = 1, nZ + emissivity_AD = emissivity_AD - & + ( SfcOptics_AD%Weight(j) * SUM(SfcOptics_AD%Reflectivity(1:nZ,1,j,1)) ) + END DO + SfcOptics_AD%Reflectivity = ZERO SfcOptics_AD%Direct_Reflectivity = ZERO SfcOptics_AD%Emissivity = ZERO + SfcOptics_AD%Emissivity(1,1) = emissivity_AD END FUNCTION Compute_VIS_Water_SfcOptics_AD diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 692dad4c..a297fbbc 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -4,7 +4,7 @@ # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. cmake_minimum_required (VERSION 3.12) -project("CRTM_Tests" VERSION 3.1.2 LANGUAGES Fortran C) +project("CRTM_Tests" VERSION 3.2.0 LANGUAGES Fortran C) enable_testing () @@ -220,13 +220,13 @@ list( APPEND Aircraft_Sensor_Ids list( APPEND Downwelling_Radiance_Sensor_Ids atms_n21 - cris-fsr_n21 - v.abi_g18 - atms_npp - cris399_npp - v.abi_gr - abi_g18 - modis_aqua +# cris-fsr_n21 +# v.abi_g18 +# atms_npp +# cris399_npp +# v.abi_gr +# abi_g18 +# modis_aqua ) @@ -281,13 +281,6 @@ add_test(NAME test_check_crtm COMMAND test_check_crtm) set_tests_properties(test_check_crtm PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") -add_executable(test_check_crtm_random mains/application/check_crtm_random_profiles.F90) -target_link_libraries(test_check_crtm_random PRIVATE crtm) -add_test(NAME test_check_crtm_random - COMMAND test_check_crtm_random) -set_tests_properties(test_check_crtm_random PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") - - add_executable(Unit_TL_TEST mains/unit/Unit_Test/test_TL.f90) target_link_libraries(Unit_TL_TEST PRIVATE crtm) add_test(NAME test_Unit_TL_TEST @@ -318,6 +311,12 @@ add_test(NAME test_Unit_Aerosol_Bypass_k_matrix COMMAND $) set_tests_properties(test_Unit_Aerosol_Bypass_k_matrix PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") +add_executable(test_Options_Overrides mains/unit/Unit_Test/test_Options_Overrides.f90) +target_link_libraries(test_Options_Overrides PRIVATE crtm) +add_test(NAME test_Options_Overrides + COMMAND $) +set_tests_properties(test_Options_Overrides PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=$ENV{OMP_NUM_THREADS}") + #SpcCoeff utilities list (APPEND SCoeff_Utils SpcCoeff_Edit diff --git a/test/mains/regression/adjoint/test_Simple/test_Simple.f90 b/test/mains/regression/adjoint/test_Simple/test_Simple.f90 index 4cd13f8e..18731c0f 100644 --- a/test/mains/regression/adjoint/test_Simple/test_Simple.f90 +++ b/test/mains/regression/adjoint/test_Simple/test_Simple.f90 @@ -78,6 +78,7 @@ PROGRAM test_Simple INTEGER :: Allocate_Status INTEGER :: n_Channels INTEGER :: l, m + INTEGER :: obs_level ! Declarations for Adjoint comparisons INTEGER :: n_la, n_ma INTEGER :: n_ls, n_ms @@ -104,6 +105,7 @@ PROGRAM test_Simple TYPE(CRTM_Atmosphere_type) :: Atmosphere_AD(N_PROFILES) TYPE(CRTM_Surface_type) :: Surface_AD(N_PROFILES) TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_AD(:,:) + TYPE(CRTM_Options_type) :: Options(N_PROFILES) ! ============================================================================ @@ -244,6 +246,10 @@ PROGRAM test_Simple ! -------------------------------- CALL Load_Atm_Data() CALL Load_Sfc_Data() + obs_level = N_LAYERS / 2 + DO m = 1, N_PROFILES + Options(m)%Obs_4_downward_P = Atm(m)%Level_Pressure(obs_level) + END DO ! 4b. GeometryInfo input @@ -271,8 +277,11 @@ PROGRAM test_Simple ! and the visible results are dR/dx ! ------------------------------------- DO l = 1, n_Channels - IF ( ChannelInfo(1)%Sensor_Type == INFRARED_SENSOR .OR. & - ChannelInfo(1)%Sensor_Type == MICROWAVE_SENSOR ) THEN + IF ( Options(1)%Obs_4_downward_P > ZERO ) THEN + RTSolution_AD(l,:)%Radiance = ONE + RTSolution_AD(l,:)%Brightness_Temperature = ZERO + ELSE IF ( ChannelInfo(1)%Sensor_Type == INFRARED_SENSOR .OR. & + ChannelInfo(1)%Sensor_Type == MICROWAVE_SENSOR ) THEN RTSolution_AD(l,:)%Radiance = ZERO RTSolution_AD(l,:)%Brightness_Temperature = ONE ELSE @@ -296,7 +305,8 @@ PROGRAM test_Simple ChannelInfo , & Atmosphere_AD, & Surface_AD , & - RTSolution ) + RTSolution , & + Options=Options ) IF ( Error_Status /= SUCCESS ) THEN Message = 'Error in CRTM Adjoint Model' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) 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 1615f2cb..72ddd13b 100644 --- a/test/mains/regression/k_matrix/test_Simple/test_Simple.f90 +++ b/test/mains/regression/k_matrix/test_Simple/test_Simple.f90 @@ -54,6 +54,7 @@ PROGRAM test_Simple INTEGER :: Allocate_Status INTEGER :: n_Channels INTEGER :: l, m + INTEGER :: obs_level ! Declarations for Jacobian comparisons INTEGER :: n_la, n_ma INTEGER :: n_ls, n_ms @@ -80,6 +81,7 @@ PROGRAM test_Simple TYPE(CRTM_Atmosphere_type), ALLOCATABLE :: Atmosphere_K(:,:) TYPE(CRTM_Surface_type) , ALLOCATABLE :: Surface_K(:,:) TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_K(:,:) + TYPE(CRTM_Options_type) :: Options(N_PROFILES) ! ============================================================================ !First, make sure the right number of inputs have been provided @@ -184,6 +186,10 @@ PROGRAM test_Simple ! -------------------------------- CALL Load_Atm_Data() CALL Load_Sfc_Data() + obs_level = N_LAYERS / 2 + DO m = 1, N_PROFILES + Options(m)%Obs_4_downward_P = Atm(m)%Level_Pressure(obs_level) + END DO ! 4b. GeometryInfo input @@ -212,8 +218,12 @@ PROGRAM test_Simple ! and the visible results are dR/dx ! ------------------------------------- DO l = 1, n_Channels - IF ( ChannelInfo(1)%Sensor_Type == INFRARED_SENSOR .OR. & - ChannelInfo(1)%Sensor_Type == MICROWAVE_SENSOR ) THEN + IF ( Options(1)%Obs_4_downward_P > ZERO ) THEN + RTSolution_K(l,:)%Radiance = ZERO + RTSolution_K(l,:)%Down_Radiance = ONE + RTSolution_K(l,:)%Brightness_Temperature = ZERO + ELSE IF ( ChannelInfo(1)%Sensor_Type == INFRARED_SENSOR .OR. & + ChannelInfo(1)%Sensor_Type == MICROWAVE_SENSOR ) THEN RTSolution_K(l,:)%Radiance = ZERO RTSolution_K(l,:)%Brightness_Temperature = ONE ELSE @@ -236,7 +246,8 @@ PROGRAM test_Simple ChannelInfo , & Atmosphere_K, & Surface_K , & - RTSolution ) + RTSolution , & + Options=Options ) IF ( Error_Status /= SUCCESS ) THEN Message = 'Error in CRTM K_Matrix Model' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) diff --git a/test/mains/regression/tangent_linear/test_Simple/test_Simple.f90 b/test/mains/regression/tangent_linear/test_Simple/test_Simple.f90 index b329e346..f0a259fb 100644 --- a/test/mains/regression/tangent_linear/test_Simple/test_Simple.f90 +++ b/test/mains/regression/tangent_linear/test_Simple/test_Simple.f90 @@ -53,6 +53,7 @@ PROGRAM test_Simple INTEGER :: Allocate_Status INTEGER :: n_Channels INTEGER :: l, m + INTEGER :: obs_level ! Declarations for RTSolution comparison INTEGER :: n_l, n_m CHARACTER(256) :: rts_File @@ -67,6 +68,7 @@ PROGRAM test_Simple TYPE(CRTM_Atmosphere_type) :: Atm(N_PROFILES), Atm_TL(N_PROFILES) TYPE(CRTM_Surface_type) :: Sfc(N_PROFILES), Sfc_TL(N_PROFILES) TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution(:,:), RTSolution_TL(:,:) + TYPE(CRTM_Options_type) :: Options(N_PROFILES) ! ============================================================================ !First, make sure the right number of inputs have been provided @@ -159,6 +161,10 @@ PROGRAM test_Simple ! ------------------------------------ CALL Load_Atm_Data() CALL Load_Sfc_Data() + obs_level = N_LAYERS / 2 + DO m = 1, N_PROFILES + Options(m)%Obs_4_downward_P = Atm(m)%Level_Pressure(obs_level) + END DO ! 4b. GeometryInfo input @@ -206,7 +212,8 @@ PROGRAM test_Simple Geometry , & ChannelInfo , & RTSolution , & - RTSolution_TL ) + RTSolution_TL , & + Options=Options ) IF ( Error_Status /= SUCCESS ) THEN Message = 'Error in CRTM Tangent-Linear Model' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) diff --git a/test/mains/unit/Unit_Test/test_Options_Overrides.f90 b/test/mains/unit/Unit_Test/test_Options_Overrides.f90 new file mode 100644 index 00000000..a027c17a --- /dev/null +++ b/test/mains/unit/Unit_Test/test_Options_Overrides.f90 @@ -0,0 +1,344 @@ +! +! test_Options_Overrides.f90 +! +! Unit test to verify option overrides affect Forward/TL/AD/K outputs +! only when the Options argument is present. +! +PROGRAM test_Options_Overrides + + ! ============================================================================ + ! **** ENVIRONMENT SETUP FOR RTM USAGE **** + ! + USE CRTM_Module + IMPLICIT NONE + ! ============================================================================ + + ! ---------- + ! Parameters + ! ---------- + CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'test_Options_Overrides' + CHARACTER(*), PARAMETER :: COEFFICIENTS_PATH = './testinput/' + INTEGER, PARAMETER :: N_PROFILES = 2 + INTEGER, PARAMETER :: N_LAYERS = 92 + INTEGER, PARAMETER :: N_ABSORBERS = 2 + INTEGER, PARAMETER :: N_CLOUDS = 0 + INTEGER, PARAMETER :: N_AEROSOLS = 0 + INTEGER, PARAMETER :: N_SENSORS = 1 + INTEGER, PARAMETER :: N_STREAMS = 8 + REAL(fp), PARAMETER :: ZENITH_ANGLE = 30.0_fp + REAL(fp), PARAMETER :: SCAN_ANGLE = 26.37293341421_fp + REAL(fp), PARAMETER :: DIFF_TOL = 1.0e-6_fp + + ! --------- + ! Variables + ! --------- + CHARACTER(256) :: Message + CHARACTER(256) :: Version + CHARACTER(256) :: Sensor_Id + INTEGER :: Error_Status + INTEGER :: Allocate_Status + INTEGER :: n_Channels + INTEGER :: m + REAL(fp) :: r_noopt, r_opt + INTEGER :: n_full_opt + + REAL(fp), ALLOCATABLE :: Emissivity(:) + + TYPE(CRTM_ChannelInfo_type) :: ChannelInfo(N_SENSORS) + TYPE(CRTM_Geometry_type) :: Geometry(N_PROFILES) + TYPE(CRTM_Atmosphere_type) :: Atm(N_PROFILES) + TYPE(CRTM_Surface_type) :: Sfc(N_PROFILES) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution(:,:) + TYPE(CRTM_Options_type) :: Opt(N_PROFILES) + + TYPE(CRTM_Atmosphere_type) :: Atmosphere_TL(N_PROFILES) + TYPE(CRTM_Surface_type) :: Surface_TL(N_PROFILES) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_TL(:,:) + + TYPE(CRTM_Atmosphere_type) :: Atmosphere_AD(N_PROFILES) + TYPE(CRTM_Surface_type) :: Surface_AD(N_PROFILES) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_AD(:,:) + + TYPE(CRTM_Atmosphere_type), ALLOCATABLE :: Atmosphere_K(:,:) + TYPE(CRTM_Surface_type) , ALLOCATABLE :: Surface_K(:,:) + TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_K(:,:) + + ! Program header + ! -------------- + CALL CRTM_Version( Version ) + CALL Program_Message( PROGRAM_NAME, '', & + 'CRTM Version: '//TRIM(Version) ) + + Sensor_Id = 'atms_n21' + Sensor_Id = ADJUSTL(Sensor_Id) + WRITE( *,'(//5x,"Running CRTM for ",a," sensor...")' ) TRIM(Sensor_Id) + + ! ============================================================================ + ! 1. **** INITIALIZE THE CRTM **** + ! + WRITE( *,'(/5x,"Initializing the CRTM...")' ) + Error_Status = CRTM_Init( (/Sensor_Id/), & + ChannelInfo, & + File_Path=COEFFICIENTS_PATH) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error initializing CRTM' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + n_Channels = SUM(CRTM_ChannelInfo_n_Channels(ChannelInfo)) + + ! ============================================================================ + ! 2. **** ALLOCATE STRUCTURE ARRAYS **** + ! + ALLOCATE( RTSolution( n_Channels, N_PROFILES ), & + RTSolution_TL( n_Channels, N_PROFILES ), & + RTSolution_AD( n_Channels, N_PROFILES ), & + RTSolution_K( n_Channels, N_PROFILES ), & + Atmosphere_K( n_Channels, N_PROFILES ), & + Surface_K( n_Channels, N_PROFILES ), & + STAT = Allocate_Status ) + IF ( Allocate_Status /= 0 ) THEN + Message = 'Error allocating structure arrays' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + + CALL CRTM_Atmosphere_Create( Atm, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atm)) ) THEN + Message = 'Error allocating CRTM Atmosphere structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + CALL CRTM_Atmosphere_Create( Atmosphere_TL, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + CALL CRTM_Atmosphere_Create( Atmosphere_AD, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + CALL CRTM_Atmosphere_Create( Atmosphere_K, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS ) + + ! ============================================================================ + ! 3. **** ASSIGN INPUT DATA **** + ! + CALL Load_Atm_Data() + CALL Load_Sfc_Data() + + CALL CRTM_Geometry_SetValue( Geometry, & + Sensor_Zenith_Angle = ZENITH_ANGLE, & + Sensor_Scan_Angle = SCAN_ANGLE ) + + ! ============================================================================ + ! 4. **** SET UP OPTIONS **** + ! + ALLOCATE( Emissivity( n_Channels ), STAT=Allocate_Status ) + IF ( Allocate_Status /= 0 ) THEN + Message = 'Error allocating Emissivity' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + Emissivity(:) = 0.5_fp + + CALL CRTM_Options_Create( Opt, n_Channels ) + IF ( .NOT. ALL(CRTM_Options_Associated( Opt )) ) THEN + Message = 'Error allocating Options structure' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + DO m = 1, N_PROFILES + CALL CRTM_Options_SetValue( Opt(m), n_Streams = N_STREAMS ) + CALL CRTM_Options_SetEmissivity( Opt(m), Emissivity ) + END DO + + ! ============================================================================ + ! 5. **** FORWARD WITHOUT OPTIONS **** + ! + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo, & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Forward Model (no options)' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + r_noopt = RTSolution(1,1)%Radiance + + ! ============================================================================ + ! 6. **** FORWARD WITH OPTIONS **** + ! + Error_Status = CRTM_Forward( Atm , & + Sfc , & + Geometry , & + ChannelInfo, & + RTSolution , & + Options = Opt ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Forward Model (options)' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + r_opt = RTSolution(1,1)%Radiance + n_full_opt = RTSolution(1,1)%n_Full_Streams + IF ( n_full_opt /= (Opt(1)%n_Streams + 2) ) THEN + WRITE(*,'(a,2i6)') 'n_Full_Streams mismatch: ', n_full_opt, Opt(1)%n_Streams + 2 + STOP 1 + END IF + IF ( ABS(r_opt - r_noopt) <= DIFF_TOL ) THEN + WRITE(*,'(a,2es14.6)') 'Radiance not affected by options: ', r_opt, r_noopt + STOP 1 + END IF + + ! ============================================================================ + ! 7. **** TL WITH/WITHOUT OPTIONS **** + ! + CALL CRTM_Atmosphere_Zero( Atmosphere_TL ) + CALL CRTM_Surface_Zero( Surface_TL ) + CALL CRTM_RTSolution_Zero( RTSolution_TL ) + Error_Status = CRTM_Tangent_Linear( Atm , & + Sfc , & + Atmosphere_TL , & + Surface_TL , & + Geometry , & + ChannelInfo , & + RTSolution , & + RTSolution_TL ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Tangent-Linear Model (no options)' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + IF ( ABS(RTSolution(1,1)%Radiance - r_noopt) > DIFF_TOL ) THEN + WRITE(*,'(a,2es14.6)') 'TL forward mismatch (no options): ', RTSolution(1,1)%Radiance, r_noopt + STOP 1 + END IF + + CALL CRTM_RTSolution_Zero( RTSolution_TL ) + Error_Status = CRTM_Tangent_Linear( Atm , & + Sfc , & + Atmosphere_TL , & + Surface_TL , & + Geometry , & + ChannelInfo , & + RTSolution , & + RTSolution_TL , & + Options = Opt ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Tangent-Linear Model (options)' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + IF ( RTSolution(1,1)%n_Full_Streams /= (Opt(1)%n_Streams + 2) ) THEN + WRITE(*,'(a,2i6)') 'TL n_Full_Streams mismatch: ', RTSolution(1,1)%n_Full_Streams, Opt(1)%n_Streams + 2 + STOP 1 + END IF + IF ( ABS(RTSolution(1,1)%Radiance - r_opt) > DIFF_TOL ) THEN + WRITE(*,'(a,2es14.6)') 'TL forward mismatch (options): ', RTSolution(1,1)%Radiance, r_opt + STOP 1 + END IF + + ! ============================================================================ + ! 8. **** AD WITH/WITHOUT OPTIONS **** + ! + CALL CRTM_Atmosphere_Zero( Atmosphere_AD ) + CALL CRTM_Surface_Zero( Surface_AD ) + CALL CRTM_RTSolution_Zero( RTSolution_AD ) + Error_Status = CRTM_Adjoint( Atm , & + Sfc , & + RTSolution_AD, & + Geometry , & + ChannelInfo , & + Atmosphere_AD, & + Surface_AD , & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Adjoint Model (no options)' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + IF ( ABS(RTSolution(1,1)%Radiance - r_noopt) > DIFF_TOL ) THEN + WRITE(*,'(a,2es14.6)') 'AD forward mismatch (no options): ', RTSolution(1,1)%Radiance, r_noopt + STOP 1 + END IF + + CALL CRTM_Atmosphere_Zero( Atmosphere_AD ) + CALL CRTM_Surface_Zero( Surface_AD ) + CALL CRTM_RTSolution_Zero( RTSolution_AD ) + Error_Status = CRTM_Adjoint( Atm , & + Sfc , & + RTSolution_AD, & + Geometry , & + ChannelInfo , & + Atmosphere_AD, & + Surface_AD , & + RTSolution , & + Options = Opt ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM Adjoint Model (options)' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + IF ( RTSolution(1,1)%n_Full_Streams /= (Opt(1)%n_Streams + 2) ) THEN + WRITE(*,'(a,2i6)') 'AD n_Full_Streams mismatch: ', RTSolution(1,1)%n_Full_Streams, Opt(1)%n_Streams + 2 + STOP 1 + END IF + IF ( ABS(RTSolution(1,1)%Radiance - r_opt) > DIFF_TOL ) THEN + WRITE(*,'(a,2es14.6)') 'AD forward mismatch (options): ', RTSolution(1,1)%Radiance, r_opt + STOP 1 + END IF + + ! ============================================================================ + ! 9. **** K-MATRIX WITH/WITHOUT OPTIONS **** + ! + CALL CRTM_Atmosphere_Zero( Atmosphere_K ) + CALL CRTM_Surface_Zero( Surface_K ) + CALL CRTM_RTSolution_Zero( RTSolution_K ) + Error_Status = CRTM_K_Matrix( Atm , & + Sfc , & + RTSolution_K, & + Geometry , & + ChannelInfo , & + Atmosphere_K, & + Surface_K , & + RTSolution ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM K-Matrix Model (no options)' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + IF ( ABS(RTSolution(1,1)%Radiance - r_noopt) > DIFF_TOL ) THEN + WRITE(*,'(a,2es14.6)') 'K forward mismatch (no options): ', RTSolution(1,1)%Radiance, r_noopt + STOP 1 + END IF + + CALL CRTM_Atmosphere_Zero( Atmosphere_K ) + CALL CRTM_Surface_Zero( Surface_K ) + CALL CRTM_RTSolution_Zero( RTSolution_K ) + Error_Status = CRTM_K_Matrix( Atm , & + Sfc , & + RTSolution_K, & + Geometry , & + ChannelInfo , & + Atmosphere_K, & + Surface_K , & + RTSolution , & + Options = Opt ) + IF ( Error_Status /= SUCCESS ) THEN + Message = 'Error in CRTM K-Matrix Model (options)' + CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + STOP 1 + END IF + IF ( RTSolution(1,1)%n_Full_Streams /= (Opt(1)%n_Streams + 2) ) THEN + WRITE(*,'(a,2i6)') 'K n_Full_Streams mismatch: ', RTSolution(1,1)%n_Full_Streams, Opt(1)%n_Streams + 2 + STOP 1 + END IF + IF ( ABS(RTSolution(1,1)%Radiance - r_opt) > DIFF_TOL ) THEN + WRITE(*,'(a,2es14.6)') 'K forward mismatch (options): ', RTSolution(1,1)%Radiance, r_opt + STOP 1 + END IF + + STOP 0 + +CONTAINS + + INCLUDE 'Load_Atm_Data.inc' + INCLUDE 'Load_Sfc_Data.inc' + +END PROGRAM test_Options_Overrides