aircraft & downwelling RT support, surface-reflectance derivation, and clear-sky bookkeeping fixes

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)
-

README.md

@@ -1,13 +1,14 @@
-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)
 
 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
 
 
 
-

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" )

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

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

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)