Feature/cd ssp input

.gitignore

@@ -14,3 +14,4 @@ build.bash
 rebuild.bash
 *~
 conductor/
+.vscode/

src/CMakeLists.txt

@@ -143,6 +143,7 @@ list( APPEND crtm_src_files
       Options/CRTM_Options_Define.f90
       Options/SSU_Input/SSU_Input_Define.f90
       Options/Zeeman_Input/Zeeman_Input_Define.f90
+      Options/OP_Input/OP_Input_Define.f90
       RTSolution/ADA/ADA_Module.f90
       RTSolution/Common_RTSolution.f90
       RTSolution/CRTM_RTSolution_Define.f90

src/CRTM_Forward_Module.f90

@@ -100,6 +100,7 @@ MODULE CRTM_Forward_Module
   USE CRTM_CloudCover_Define,     ONLY: CRTM_CloudCover_type
   USE CRTM_Active_Sensor,         ONLY: CRTM_Compute_Reflectivity, &
                                         Calculate_Cloud_Water_Density
+  USE OP_Input_Define,            ONLY: OP_Input_Channel_Position
 
   ! Internal variable definition modules
   ! ...AtmOptics
@@ -455,6 +456,8 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status )
       LOGICAL :: Atmosphere_Invalid, Surface_Invalid, Geometry_Invalid, Options_Invalid
       INTEGER :: iFOV
       INTEGER :: n, l    ! sensor index, channel index
+      INTEGER :: ilay, iphas, ileg, ileg1
+      INTEGER :: pos  ! column position of ChannelIndex within opt%TOP/COP/AOP
       INTEGER :: SensorIndex
       INTEGER :: ChannelIndex
       INTEGER :: ln, nc, ks
@@ -1022,39 +1025,101 @@ FUNCTION profile_solution (m, Opt, AncillaryInput) RESULT( Error_Status )
                   END IF
                END IF
 
-               ! Compute the cloud particle absorption/scattering properties
-               IF( Atm%n_Clouds > 0 ) THEN
-                  Err_Thread = CRTM_Compute_CloudScatter( Atm          , &  ! Input
-                                                            GeometryInfo , &  ! Input
-                                                            SensorIndex  , &  ! Input
-                                                            ChannelIndex , &  ! Input
-                                                            AtmOptics(nt), &  ! Output
-                                                            CSvar(nt)       ) ! Internal variable output
-                  IF ( Err_Thread /= SUCCESS ) THEN
-                     WRITE( Message,'("Error computing CloudScatter for ",a,&
-                          &", channel ",i0,", profile #",i0)' ) &
-                          TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m
-                     CALL Display_Message( ROUTINE_NAME, Message, Err_Thread )
-                     thread_error = MAX(thread_error, Err_Thread)
-                  END IF
-               END IF
+               ! Get or compute cloud particle absorption/scattering properties
+               !
+               ! pos maps ChannelIndex to the column of opt%TOP/COP/AOP holding this
+               ! channel's data. It is 0 (not found) whenever the corresponding
+               ! Use_*_OP flag is off, or when it's on but this particular channel
+               ! isn't covered by the supplied OP_Input (e.g. a single-channel
+               ! OP_Input in a multi-channel run) - either way, fall through to the
+               ! usual internal computation for this channel.
+               pos = 0
+               IF ( Options_Present .AND. opt%Use_Total_OP ) pos = OP_Input_Channel_Position( opt%TOP, ChannelIndex )
+               IF ( pos >= 1 ) THEN ! use total optical properties from user-defined OP_Input
+                  AtmOptics(nt)%Include_Scattering = .TRUE.
+                  DO ilay = 1, Atm%n_Layers
+                     AtmOptics(nt)%Optical_Depth(ilay)         = AtmOptics(nt)%Optical_Depth(ilay)         + opt%TOP%tau(pos, ilay)
+                     AtmOptics(nt)%Single_Scatter_Albedo(ilay) = AtmOptics(nt)%Single_Scatter_Albedo(ilay) + opt%TOP%bs(pos, ilay)
+                     AtmOptics(nt)%Backscat_Coefficient(ilay)  = AtmOptics(nt)%Backscat_Coefficient(ilay)  + opt%TOP%kb(pos, ilay)
+                     DO iphas = 1, 1
+                        DO ileg = 0, AtmOptics(nt)%n_Legendre_Terms
+                        ileg1 = ileg + 1
+                        AtmOptics(nt)%Phase_Coefficient(ileg,iphas,ilay) = AtmOptics(nt)%Phase_Coefficient(ileg,iphas,ilay) + &
+                                                                           opt%TOP%pcoeff(pos,ilay,iphas,ileg1)
+                        END DO
+                     END DO
+                  END DO
 
-               ! Compute the aerosol absorption/scattering properties
-               IF ( Atm%n_Aerosols > 0 ) THEN
-                  Err_Thread = CRTM_Compute_AerosolScatter( Atm          , &  ! Input
-                                                              SensorIndex  , &  ! Input
-                                                              ChannelIndex , &  ! Input
-                                                              AtmOptics(nt), &  ! In/Output
-                                                              ASvar(nt)       ) ! Internal variable output
+               ELSE
 
-                  IF ( Err_Thread /= SUCCESS ) THEN
-                     WRITE( Message,'("Error computing AerosolScatter for ",a,&
-                          &", channel ",i0,", profile #",i0)' ) &
-                          TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m
-                     CALL Display_Message( ROUTINE_NAME, Message, Err_Thread )
-                     thread_error = MAX(thread_error, Err_Thread)
+                  ! ...Clouds
+                  pos = 0
+                  IF ( Options_Present .AND. opt%Use_Cloud_OP ) pos = OP_Input_Channel_Position( opt%COP, ChannelIndex )
+                  IF ( pos >= 1 ) THEN ! use cloud optical properties from user-defined OP_Input
+                     ! Use user-defined cloud optical profiles
+                     AtmOptics(nt)%Include_Scattering = .TRUE.
+                     DO ilay = 1, Atm%n_Layers
+                        AtmOptics(nt)%Optical_Depth(ilay)         = AtmOptics(nt)%Optical_Depth(ilay)         + opt%COP%tau(pos, ilay)
+                        AtmOptics(nt)%Single_Scatter_Albedo(ilay) = AtmOptics(nt)%Single_Scatter_Albedo(ilay) + opt%COP%bs(pos, ilay)
+                        AtmOptics(nt)%Backscat_Coefficient(ilay)  = AtmOptics(nt)%Backscat_Coefficient(ilay)  + opt%COP%kb(pos, ilay)
+                        DO iphas = 1, 1
+                        DO ileg = 0, AtmOptics(nt)%n_Legendre_Terms
+                           ileg1 = ileg + 1
+                           AtmOptics(nt)%Phase_Coefficient(ileg,iphas,ilay) = AtmOptics(nt)%Phase_Coefficient(ileg,iphas,ilay) + &
+                                                                              opt%COP%pcoeff(pos,ilay,iphas,ileg1)
+                        END DO
+                        END DO
+                     END DO
+                  ELSEIF( Atm%n_Clouds > 0 ) THEN ! else compute, CRTM default
+                     ! Compute the cloud particle absorption/scattering properties
+                     Err_Thread = CRTM_Compute_CloudScatter( Atm          , &  ! Input
+                                                             GeometryInfo , &  ! Input
+                                                             SensorIndex  , &  ! Input
+                                                             ChannelIndex , &  ! Input
+                                                             AtmOptics(nt), &  ! Output
+                                                             CSvar(nt)       ) ! Internal variable output
+                     IF ( Err_Thread /= SUCCESS ) THEN
+                        WRITE( Message,'("Error computing CloudScatter for ",a,&
+                           &", channel ",i0,", profile #",i0)' ) &
+                           TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m
+                           CALL Display_Message( ROUTINE_NAME, Message, Err_Thread )
+                        thread_error = MAX(thread_error, Err_Thread)
+                     END IF
                   END IF
-               END IF
+
+                  ! ...Aerosols
+                  pos = 0
+                  IF ( Options_Present .AND. opt%Use_Aerosol_OP ) pos = OP_Input_Channel_Position( opt%AOP, ChannelIndex )
+                  IF ( pos >= 1 ) THEN  ! use aerosol optical properties from user-defined OP_Input
+                     AtmOptics(nt)%Include_Scattering = .TRUE.
+                     DO ilay = 1, Atm%n_Layers
+                        AtmOptics(nt)%Optical_Depth(ilay)         = AtmOptics(nt)%Optical_Depth(ilay)         + opt%AOP%tau(pos, ilay)
+                        AtmOptics(nt)%Single_Scatter_Albedo(ilay) = AtmOptics(nt)%Single_Scatter_Albedo(ilay) + opt%AOP%bs(pos, ilay)
+                        AtmOptics(nt)%Backscat_Coefficient(ilay)  = AtmOptics(nt)%Backscat_Coefficient(ilay)  + opt%AOP%kb(pos, ilay)
+                        DO iphas = 1, 1
+                        DO ileg = 0, AtmOptics(nt)%n_Legendre_Terms
+                           ileg1 = ileg + 1
+                           AtmOptics(nt)%Phase_Coefficient(ileg,iphas,ilay) = AtmOptics(nt)%Phase_Coefficient(ileg,iphas,ilay) + &
+                                                                              opt%AOP%pcoeff(pos,ilay,iphas,ileg1)
+                        END DO
+                        END DO
+                     END DO
+                  ELSEIF ( Atm%n_Aerosols > 0 ) THEN ! else, CRTM default
+                     ! Compute the aerosol absorption/scattering properties
+                     Err_Thread = CRTM_Compute_AerosolScatter( Atm          , &  ! Input
+                                                               SensorIndex  , &  ! Input
+                                                               ChannelIndex , &  ! Input
+                                                               AtmOptics(nt), &  ! In/Output
+                                                               ASvar(nt)      )  ! Internal variable output
+                     IF ( Err_Thread /= SUCCESS ) THEN
+                        WRITE( Message,'("Error computing AerosolScatter for ",a,&
+                              &", channel ",i0,", profile #",i0)' ) &
+                              TRIM(ChannelInfo(n)%Sensor_ID), ChannelInfo(n)%Sensor_Channel(l), m
+                        CALL Display_Message( ROUTINE_NAME, Message, Err_Thread )
+                        thread_error = MAX(thread_error, Err_Thread)
+                     END IF
+                  END IF
+               END IF  !IF ( Options_Present .AND. opt%Use_Total_OP )
 
                ! Compute the combined atmospheric optical properties
                IF( AtmOptics(nt)%Include_Scattering ) THEN

src/CRTM_Module.F90

@@ -21,6 +21,7 @@ MODULE CRTM_Module
   USE CRTM_RTSolution_Define
   USE CRTM_Options_Define
   USE CRTM_AncillaryInput_Define
+  USE OP_Input_Define
   USE CRTM_IRlandCoeff          , ONLY: CRTM_IRlandCoeff_Classification
 
   ! Parameter definition module

src/Options/CRTM_Options_Define.f90

@@ -51,6 +51,10 @@ MODULE CRTM_Options_Define
                                     CloudCover_Overcast_Overlap, &
                                     CloudCover_Overlap_IsValid, &
                                     CloudCover_Overlap_Name
+  USE OP_Input_Define      , ONLY: OP_Input_type, &
+                                   OPERATOR(==), &
+                                   OP_Input_IsValid, &
+                                   OP_Input_Inspect
   ! Disable implicit typing
   IMPLICIT NONE
 
@@ -196,6 +200,15 @@ MODULE CRTM_Options_Define
 
     ! Whether to skip this profile
     LOGICAL :: Skip_Profile = .FALSE.
+
+    ! User defined optical profiles
+    LOGICAL :: Use_Aerosol_OP = .FALSE.  ! Aerosol
+    LOGICAL :: Use_Cloud_OP   = .FALSE.  ! Cloud
+    LOGICAL :: Use_Total_OP   = .FALSE.  ! Aerosol + Cloud
+    TYPE(OP_Input_Type) :: AOP           ! Aerosol optical profiles
+    TYPE(OP_Input_Type) :: COP           ! Cloud optical profiles
+    TYPE(OP_Input_Type) :: TOP           ! Total optical profiles
+
   END TYPE CRTM_Options_type
   !:tdoc-:
 
@@ -420,8 +433,11 @@ ELEMENTAL SUBROUTINE CRTM_Options_SetValue( &
     Set_Overcast_Overlap         , &
     Use_Emissivity               , &
     Use_Direct_Reflectivity      , &
+    Use_Aerosol_OP               , &
+    Use_Cloud_OP                 , &
+    Use_Total_OP                 , &
     n_Streams                    , &
-    Aircraft_Pressure            , &
+    Aircraft_Pressure            , &  
     Compute_Down_Radiance        , &
     Compute_Down_Radiance_Profile, &
     Compute_Up_Radiance_Profile    )
@@ -441,6 +457,9 @@ ELEMENTAL SUBROUTINE CRTM_Options_SetValue( &
     LOGICAL ,      OPTIONAL, INTENT(IN)     :: Set_Overcast_Overlap
     LOGICAL ,      OPTIONAL, INTENT(IN)     :: Use_Emissivity
     LOGICAL ,      OPTIONAL, INTENT(IN)     :: Use_Direct_Reflectivity
+    LOGICAL ,      OPTIONAL, INTENT(IN)     :: Use_Aerosol_OP
+    LOGICAL ,      OPTIONAL, INTENT(IN)     :: Use_Cloud_OP
+    LOGICAL ,      OPTIONAL, INTENT(IN)     :: Use_Total_OP
     INTEGER ,      OPTIONAL, INTENT(IN)     :: n_Streams
     REAL(fp),      OPTIONAL, INTENT(IN)     :: Aircraft_Pressure
     LOGICAL ,      OPTIONAL, INTENT(IN)     :: Compute_Down_Radiance
@@ -494,6 +513,18 @@ ELEMENTAL SUBROUTINE CRTM_Options_SetValue( &
     IF ( PRESENT(Use_Direct_Reflectivity) ) &
       self%Use_Direct_Reflectivity = Use_Direct_Reflectivity .AND. self%Is_Allocated
 
+    ! Aerosol optical profiles
+    IF ( PRESENT(Use_Aerosol_OP) ) &
+      self%Use_Aerosol_OP = Use_Aerosol_OP .AND. self%Is_Allocated
+
+    ! Cloud optical profiles
+    IF ( PRESENT(Use_Cloud_OP) ) &
+      self%Use_Cloud_OP = Use_Cloud_OP .AND. self%Is_Allocated
+
+    ! Total optical profiles
+    IF ( PRESENT(Use_Total_OP) ) &
+      self%Use_Total_OP = Use_Total_OP .AND. self%Is_Allocated
+
   END SUBROUTINE CRTM_Options_SetValue
 
 
@@ -849,6 +880,19 @@ FUNCTION CRTM_Options_IsValid( self ) RESULT( IsValid )
     ! Check cloud overlap option validity
     IsValid = CloudCover_Overlap_IsValid( self%Overlap_Id ) .AND. IsValid
 
+    ! Check OP input Options
+    IF ( self%Use_Aerosol_OP ) THEN
+      IsValid = OP_Input_IsValid( self%AOP ) .AND. IsValid
+    END IF
+
+    IF ( self%Use_Cloud_OP ) THEN
+      IsValid = OP_Input_IsValid( self%COP ) .AND. IsValid
+    END IF
+
+    IF ( self%Use_Total_OP ) THEN
+      IsValid = OP_Input_IsValid( self%TOP ) .AND. IsValid
+    END IF
+
   END FUNCTION CRTM_Options_IsValid
 
 
@@ -909,6 +953,10 @@ SUBROUTINE CRTM_Options_Inspect( self )
     CALL SSU_Input_Inspect( self%SSU )
     ! ...Zeeman input
     CALL Zeeman_Input_Inspect( self%Zeeman )
+    ! ...OP input
+    CALL OP_Input_Inspect( self%AOP )
+    CALL OP_Input_Inspect( self%COP )
+    CALL OP_Input_Inspect( self%TOP )
 
   END SUBROUTINE CRTM_Options_Inspect