diff --git a/src/RTSolution/CRTM_Active_Sensor.f90 b/src/RTSolution/CRTM_Active_Sensor.f90 index 5faa24e6..159d26f2 100644 --- a/src/RTSolution/CRTM_Active_Sensor.f90 +++ b/src/RTSolution/CRTM_Active_Sensor.f90 @@ -23,7 +23,7 @@ MODULE CRTM_Active_Sensor SpcCoeff_IsVisibleSensor , & SpcCoeff_IsUltravioletSensor USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type, & - H2O_ID, MASS_MIXING_RATIO_UNITS + H2O_ID USE CRTM_AtmOptics_Define, ONLY: CRTM_AtmOptics_type USE CRTM_RTSolution_Define, ONLY: CRTM_RTSolution_type USE Spectral_Units_Conversion, ONLY: GHz_to_inverse_cm @@ -31,7 +31,6 @@ MODULE CRTM_Active_Sensor USE ODPS_CoordinateMapping, ONLY: Geopotential_Height USE CRTM_GeometryInfo_Define, ONLY: CRTM_GeometryInfo_type USE Message_Handler , ONLY: Display_Message - ! Disable all implicit typing IMPLICIT NONE @@ -54,6 +53,7 @@ MODULE CRTM_Active_Sensor ! 1.0d818 converts from m^3 to mm^6/m^3 (standard radar reflectivity units) REAL(fp), PARAMETER :: M6_MM6 = 1.0d18 REAL(fp), PARAMETER :: REFLECTIVITY_THRESHOLD = TINY(REAL(fp)) + REAL(fp), PARAMETER :: MM6_TO_CM6 = 1.0e-6_fp CONTAINS !-------------------------------------------------------------------------------- @@ -100,6 +100,8 @@ Function Calculate_Height(Atm) RESULT (Height) Atm%Absorber(:, H2O_idx), & ! Input ZERO , & ! Input - surface height Height ) ! Output in km + + Height = Height * ONE_THOUSAND ! convert to meters END FUNCTION Calculate_Height @@ -122,21 +124,22 @@ END FUNCTION Calculate_Height ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN/OUT) ! +! GeometryInfo: Structure containing the view geometry data. +! UNITS: N/A +! TYPE: CRTM_GeometryInfo_type +! DIMENSION: Scalar +! ATTRIBUTES: INTENT(IN) +! !-------------------------------------------------------------------------------- -SUBROUTINE Calculate_Cloud_Water_Density(Atm) +SUBROUTINE Calculate_Cloud_Water_Density(Atm, & + GeometryInfo) TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm + TYPE(CRTM_GeometryInfo_type), OPTIONAL, INTENT(IN) :: GeometryInfo INTEGER :: n_Layers REAL(fp) :: Height(0:Atm%n_Layers), dZ_m(Atm%n_Layers) Integer :: n - ! Function result - INTEGER :: Error_Status - ! Local parameters - CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Calculate_Cloud_Water_Density' - ! Local variables - CHARACTER(256) :: Message - n_Layers = Atm%n_Layers ! Calculate heights if hasn't been set already @@ -144,14 +147,7 @@ SUBROUTINE Calculate_Cloud_Water_Density(Atm) Atm%Height = Calculate_Height(Atm) END IF - dZ_m = (Atm%Height(0:n_Layers-1) - Atm%Height(1:n_Layers)) * ONE_THOUSAND - - IF (ANY(dZ_m <= 0)) THEN - Message = 'Error in calculating cloud water density - layer thickness needs to be greater than zero!' - CALL Display_Message( ROUTINE_NAME, TRIM(Message), Error_Status ) - RETURN - END IF - + dZ_m = (Atm%Height(0:n_Layers-1) - Atm%Height(1:n_Layers)) DO n = 1, Atm%n_Clouds Atm%Cloud(n)%Water_Density = Atm%Cloud(n)%Water_Content / dZ_m @@ -335,18 +331,27 @@ SUBROUTINE CRTM_Compute_Reflectivity(Atm , & ! Input TYPE(CRTM_GeometryInfo_type), INTENT(IN) :: GeometryInfo INTEGER , INTENT(IN) :: SensorIndex INTEGER , INTENT(IN) :: ChannelIndex - TYPE(CRTM_AtmOptics_type) , INTENT(IN) :: AtmOptics + TYPE(CRTM_AtmOptics_type) , INTENT(IN OUT) :: AtmOptics TYPE(CRTM_RTSolution_type), INTENT(IN OUT) :: RTSolution REAL(fp) :: Frequency, Wavenumber, Wavelength_m REAL(fp) :: Reflectivity(AtmOptics%n_Layers) REAL(fp) :: Reflectivity_Attenuated(AtmOptics%n_Layers) REAL(fp) :: Transmittance(AtmOptics%n_Layers) + REAL(fp) :: Optical_Depth(AtmOptics%n_Layers) REAL(fp) :: P1(AtmOptics%n_Layers) - REAL(fp) :: Height(0:AtmOptics%n_Layers), dZ_m(AtmOptics%n_Layers) , Temp_K(AtmOptics%n_Layers) + REAL(fp) :: Height(0:AtmOptics%n_Layers) + REAL(fp) :: dZ_m(AtmOptics%n_Layers) + REAL(fp) :: dS_m(AtmOptics%n_Layers) + REAL(fp) :: Temp_K(AtmOptics%n_Layers) REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Kw_2, perm_re, perm_im COMPLEX :: perm(AtmOptics%n_Layers), kw(AtmOptics%n_Layers) - INTEGER :: k + INTEGER :: k, n_Layers + REAL(fp) :: temp_sum + REAL(fp) :: ext_coef(AtmOptics%n_Layers) ! extinction coefficient + character(len=16) :: ms_env_val + + n_Layers = Atm%n_Layers ! Calculate heights if hasn't been set already IF (ALL(Atm%Height .LT. EPSILON_FP)) THEN @@ -354,8 +359,8 @@ SUBROUTINE CRTM_Compute_Reflectivity(Atm , & ! Input ELSE Height = Atm%Height ENDIF - dZ_m = (Height(0:Atm%n_Layers-1) - Height(1:Atm%n_Layers)) * ONE_THOUSAND - dZ_m = dZ_m / GeometryInfo%Cosine_Sensor_Zenith + dZ_m = (Height(0:n_Layers-1) - Height(1:n_Layers)) + dS_m = dZ_m / GeometryInfo%Cosine_Sensor_Zenith IF ( SpcCoeff_IsMicrowaveSensor(SC(SensorIndex)) ) THEN Frequency = SC(SensorIndex)%Frequency(ChannelIndex) ! GHz @@ -372,201 +377,146 @@ SUBROUTINE CRTM_Compute_Reflectivity(Atm , & ! Input perm_re = REAL(REAL(perm)) ! Double REAL is required to avoud issues with GNU Fortran perm_im = REAL(AIMAG(perm)) ! perm = cmplx(perm_re, perm_im, FP) - kw = (perm - ONE )/(perm + TWO) - Kw_2 = ABS(kw)**TWO - + !kw = (perm**2 - ONE )/(perm**2 + TWO) + kw = (perm - ONE )/(perm + ONE) + Kw_2 = ABS(kw)**TWO P1 = (M6_MM6 * Wavelength_m**4.0_fp) / (PI**5.0_fp * Kw_2) - P1 = P1 / dZ_m ! dZ_m to convert BackScatter coefficient from unitless to 1/cm - ! Calculate transmittance from top to layer k - ! Optical depth is not scaled for zenith angle - DO k = 1, AtmOptics%n_layers - Transmittance(k) = EXP(-TWO * SUM(AtmOptics%optical_depth(1:k) / & - GeometryInfo%Cosine_Sensor_Zenith)) + + ! Calculate transmittance + Optical_Depth = AtmOptics%optical_depth / dS_m + DO k = 1, n_layers + temp_sum = SUM(Optical_Depth(1:k)) + Transmittance(k) = EXP(-TWO * temp_sum) END DO - Reflectivity = P1 * (AtmOptics%Backscat_Coefficient) ! mm^6 m^-3 - Reflectivity_Attenuated = P1 * Transmittance * (AtmOptics%Backscat_Coefficient) ! mm^6 m^-3 - + ! Divide by dZ_m as water conntent calculations are based on + ! vertical thickness of layers + Reflectivity = P1 * (AtmOptics%Backscat_Coefficient / dZ_m) ! mm^6 m^-3 + Reflectivity_Attenuated = Transmittance * Reflectivity ! mm^6 m^-3 ! Convert the unit to dBz WHERE (Reflectivity .GT. REFLECTIVITY_THRESHOLD) - RTSolution%Reflectivity = TEN * LOG10(Reflectivity) ! [dBZ] + RTSolution%Reflectivity = TEN * LOG10(Reflectivity) ! [dBZ] + RTSolution%ReflectivityLinear = MM6_TO_CM6 * Reflectivity ! cm^6m^-3 ELSE WHERE - RTSolution%Reflectivity = MISSING_REFL + RTSolution%Reflectivity = MISSING_REFL + RTSolution%ReflectivityLinear = MISSING_REFL END WHERE ! Convert the unit to dBz ! Note that Reflectivity can be greater than zero but Reflectivity_Attenuated ! can be still zero if Transmittance is zero WHERE (Reflectivity_Attenuated .GT. REFLECTIVITY_THRESHOLD) - RTSolution%Reflectivity_Attenuated = TEN * LOG10(Reflectivity_Attenuated) + RTSolution%Reflectivity_Attenuated = TEN * LOG10(Reflectivity_Attenuated) ! dBZ + RTSolution%Reflectivity_AttenuatedLinear = MM6_TO_CM6 * Reflectivity_attenuated ! cm^6m^-3 ELSE WHERE - RTSolution%Reflectivity_Attenuated = MISSING_REFL + RTSolution%Reflectivity_Attenuated = MISSING_REFL + RTSolution%Reflectivity_AttenuatedLinear = MISSING_REFL END WHERE + END SUBROUTINE CRTM_Compute_Reflectivity +SUBROUTINE CRTM_Compute_Reflectivity_TL(Atm, & + AtmOptics, & + AtmOptics_TL, & + GeometryInfo, & + SensorIndex, & + ChannelIndex, & + RTSolution_TL) + TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm + TYPE(CRTM_AtmOptics_type), INTENT(IN) :: AtmOptics, AtmOptics_TL + TYPE(CRTM_GeometryInfo_type), INTENT(IN) :: GeometryInfo + INTEGER, INTENT(IN) :: SensorIndex, ChannelIndex + TYPE(CRTM_RTSolution_type), INTENT(IN OUT) :: RTSolution_TL + + ! Locals + REAL(fp) :: Frequency, Wavenumber, Wavelength_m + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Reflectivity, Reflectivity_Attenuated + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Reflectivity_TL, Reflectivity_Attenuated_TL + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: ReflectivityLinear, Reflectivity_AttenuatedLinear + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: ReflectivityLinear_TL, Reflectivity_AttenuatedLinear_TL + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Transmittance, Transmittance_TL + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Optical_Depth, Optical_Depth_TL + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: P1, Kw_2, perm_re, perm_im + REAL(fp), DIMENSION(0:AtmOptics%n_Layers) :: Height + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: dZ_m, dS_m + COMPLEX :: perm(AtmOptics%n_Layers) + INTEGER :: k, n_Layers + REAL(fp) :: temp_sum + + n_Layers = Atm%n_Layers + + ! Compute layer interface heights if missing + IF (ALL(Atm%Height .LT. EPSILON_FP)) THEN + Height = Calculate_Height(Atm) + ELSE + Height = Atm%Height + END IF -!-------------------------------------------------------------------------------- -! -! NAME: -! CRTM_Compute_Reflectivity_TL -! -! PURPOSE: -! Subroutine to calculate the tangent-linear of reflectivity for active sensors. -! -! CALLING SEQUENCE: -! CALL CRTM_Compute_Reflectivity_TL(Atm, & ! Input -! AtmOptics , & ! Input -! AtmOptics_TL , & ! Input -! GeometryInfo , & ! Input -! SensorIndex , & ! Input -! ChannelIndex, & ! Input -! RTSolution_TL ) ! Input/Output -! -! INPUT ARGUMENTS: -! INPUT ARGUMENTS: -! Atm: Structure containing the atmospheric state data. -! UNITS: N/A -! TYPE: CRTM_Atmosphere_type -! DIMENSION: Scalar -! ATTRIBUTES: INTENT(IN) -! -! AtmOptics: Structure containing the combined atmospheric -! optical properties for gaseous absorption, clouds, -! and aerosols. -! UNITS: N/A -! TYPE: CRTM_AtmOptics_type -! DIMENSION: Scalar -! ATTRIBUTES: INTENT(IN) -! -! AtmOptics_TL: Structure containing the tangent-linear atmospheric -! optical properties. -! UNITS: N/A -! TYPE: CRTM_AtmOptics_type -! DIMENSION: Scalar -! ATTRIBUTES: INTENT(IN) -! -! GeometryInfo: Structure containing the view geometry data. -! UNITS: N/A -! TYPE: CRTM_GeometryInfo_type -! DIMENSION: Scalar -! ATTRIBUTES: INTENT(IN) -! -! SensorIndex: Sensor index id. This is a unique index associated -! with a (supported) sensor used to access the -! shared coefficient data for a particular sensor. -! See the ChannelIndex argument. -! UNITS: N/A -! TYPE: INTEGER -! DIMENSION: Scalar -! ATTRIBUTES: INTENT(IN) -! -! ChannelIndex: Channel index id. This is a unique index associated -! with a (supported) sensor channel used to access the -! shared coefficient data for a particular sensor's -! channel. -! See the SensorIndex argument. -! UNITS: N/A -! TYPE: INTEGER -! DIMENSION: Scalar -! ATTRIBUTES: INTENT(IN) -! -! OUTPUT ARGUMENTS: -! RTSolution_TL: Structure containing the solution to the tangent-linear -! RT equation for the given inputs. -! UNITS: N/A -! TYPE: CRTM_RTSolution_type -! DIMENSION: Scalar -! ATTRIBUTES: INTENT(IN OUT) -!-------------------------------------------------------------------------------- - - SUBROUTINE CRTM_Compute_Reflectivity_TL(Atm, & ! Input - AtmOptics , & ! Input - AtmOptics_TL , & ! Input - GeometryInfo , & ! Input - SensorIndex , & ! Input - ChannelIndex, & ! Input - RTSolution_TL ) ! Input/Output - ! Arguments - TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm - TYPE(CRTM_GeometryInfo_type), INTENT(IN) :: GeometryInfo - INTEGER , INTENT(IN) :: SensorIndex - INTEGER , INTENT(IN) :: ChannelIndex - TYPE(CRTM_AtmOptics_type) , INTENT(IN) :: AtmOptics, AtmOptics_TL - TYPE(CRTM_RTSolution_type), INTENT(IN OUT) :: RTSolution_TL - - REAL(fp) :: Frequency, Wavenumber, Wavelength_m - REAL(fp) :: Reflectivity(AtmOptics%n_Layers) - REAL(fp) :: Reflectivity_Attenuated(AtmOptics%n_Layers) - REAL(fp) :: Reflectivity_TL(AtmOptics%n_Layers) - REAL(fp) :: Reflectivity_Attenuated_TL(AtmOptics%n_Layers) - REAL(fp) :: Transmittance(AtmOptics%n_Layers) - REAL(fp) :: Transmittance_TL(AtmOptics%n_Layers) - REAL(fp) :: P1(AtmOptics%n_Layers) - REAL(fp) :: Height(0:AtmOptics%n_Layers), dZ_m(AtmOptics%n_Layers) - COMPLEX :: perm(AtmOptics%n_Layers) - REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Kw_2, perm_re, perm_im - INTEGER :: k - - ! Calculate heights if hasn't been set already - IF (ALL(Atm%Height .LT. EPSILON_FP)) THEN - Height = Calculate_Height(Atm) - ELSE - Height = Atm%Height - ENDIF - dZ_m = (Height(0:Atm%n_Layers-1) - Height(1:Atm%n_Layers)) * ONE_THOUSAND - dZ_m = dZ_m / GeometryInfo%Cosine_Sensor_Zenith - - IF ( SpcCoeff_IsMicrowaveSensor(SC(SensorIndex)) ) THEN - Frequency = SC(SensorIndex)%Frequency(ChannelIndex) ! GHz - Wavelength_m = POINT_01 / GHz_to_Inverse_cm( Frequency ) - ELSE IF( SpcCoeff_IsInfraredSensor(SC(SensorIndex)) ) THEN - Wavenumber = SC(SensorIndex)%Wavenumber(ChannelIndex) ! 1/cm - Wavelength_m = POINT_01 / Wavenumber - END IF - - perm = Water_Permittivity_Turner_2016(Frequency * 1.0d9, & ! Input - Atm%Temperature) ! Input - perm_re = REAL(REAL(perm)) ! double REAL is required to avoid problems in GNU Fortran - perm_im = REAL(AIMAG(perm)) - Kw_2 = ((perm_re - ONE )/(perm_re + TWO))**TWO - - ! Calculate transmittance from top to layer k - !Transmittance(1) = ONE - DO k = 1, AtmOptics%n_layers - Transmittance(k) = EXP(-TWO * SUM(AtmOptics%optical_depth(1:k))) - END DO - - P1 = (M6_MM6 * Wavelength_m**4.0_fp) / (PI**5.0_fp * Kw_2) - P1 = P1 / dZ_m ! dZ_m to convert water_content to m/v or cloud water density - Reflectivity = P1 * AtmOptics%Backscat_Coefficient - Reflectivity_Attenuated = Transmittance * Reflectivity - - ! Tanget linear calculations - !Transmittance_TL(1) = ZERO - DO k = 1, AtmOptics%n_layers - Transmittance_TL(k) = - TWO * Transmittance(k) * SUM(AtmOptics_TL%optical_depth(1:k)) - END DO - - Reflectivity_TL = P1 * AtmOptics_TL%Backscat_Coefficient - Reflectivity_Attenuated_TL = Transmittance * Reflectivity_TL + & - Reflectivity * Transmittance_TL - - ! Convert the unit to dBz - WHERE (Reflectivity .GT. REFLECTIVITY_THRESHOLD) - RTSolution_TL%Reflectivity = TEN * Reflectivity_TL / (Reflectivity * LOG(TEN)) - ELSE WHERE - RTSolution_TL%Reflectivity = ZERO - END WHERE + dZ_m = Height(0:n_Layers-1) - Height(1:n_Layers) + dS_m = dZ_m / GeometryInfo%Cosine_Sensor_Zenith - ! Convert the unit to dBz - WHERE (Reflectivity_Attenuated .GT. REFLECTIVITY_THRESHOLD) - RTSolution_TL%Reflectivity_Attenuated = TEN * Reflectivity_Attenuated_TL / & - (Reflectivity_Attenuated * LOG(TEN)) - ELSE WHERE - RTSolution_TL%Reflectivity_Attenuated = ZERO - END WHERE + ! Compute wavelength + IF (SpcCoeff_IsMicrowaveSensor(SC(SensorIndex))) THEN + Frequency = SC(SensorIndex)%Frequency(ChannelIndex) + Wavelength_m = POINT_01 / GHz_to_Inverse_cm(Frequency) + ELSE IF (SpcCoeff_IsInfraredSensor(SC(SensorIndex))) THEN + Wavenumber = SC(SensorIndex)%Wavenumber(ChannelIndex) + Wavelength_m = POINT_01 / Wavenumber + END IF - END SUBROUTINE CRTM_Compute_Reflectivity_TL + ! Compute permittivity and related factor + perm = Water_Permittivity_Turner_2016(Frequency * 1.0d9, Atm%Temperature) + perm_re = REAL(perm) + perm_im = AIMAG(perm) + Kw_2 = ABS((perm_re - ONE) / (perm_re + ONE)) ** TWO + P1 = (M6_MM6 * Wavelength_m**4.0_fp) / (PI**5.0_fp * Kw_2) + + ! Forward (base state) calculations + Optical_Depth = AtmOptics%optical_depth / dS_m + DO k = 1, n_Layers + temp_sum = SUM(Optical_Depth(1:k)) + Transmittance(k) = EXP(-TWO * temp_sum) + END DO + + Reflectivity = P1 * (AtmOptics%Backscat_Coefficient / dZ_m) + Reflectivity_Attenuated = Transmittance * Reflectivity + + ! Tangent-linear calculations + Optical_Depth_TL = AtmOptics_TL%optical_depth / dS_m + DO k = 1, n_Layers + temp_sum = SUM(Optical_Depth_TL(1:k)) + Transmittance_TL(k) = -TWO * Transmittance(k) * temp_sum + END DO + + Reflectivity_TL = P1 * (AtmOptics_TL%Backscat_Coefficient / dZ_m) + Reflectivity_Attenuated_TL = Transmittance * Reflectivity_TL + Reflectivity * Transmittance_TL + + ! Linear reflectivity TL + ReflectivityLinear = MM6_TO_CM6 * Reflectivity + ReflectivityLinear_TL = MM6_TO_CM6 * Reflectivity_TL + Reflectivity_AttenuatedLinear = MM6_TO_CM6 * Reflectivity_Attenuated + Reflectivity_AttenuatedLinear_TL = MM6_TO_CM6 * Reflectivity_Attenuated_TL + + ! Convert to dBZ tangent linear if above threshold + WHERE (Reflectivity > REFLECTIVITY_THRESHOLD) + Reflectivity_TL = TEN * Reflectivity_TL / (Reflectivity * LOG(TEN)) + RTSolution_TL%Reflectivity = Reflectivity_TL + RTSolution_TL%ReflectivityLinear = ReflectivityLinear_TL + ELSEWHERE + RTSolution_TL%Reflectivity = ZERO + RTSolution_TL%ReflectivityLinear = ZERO + END WHERE + + WHERE (Reflectivity_Attenuated > REFLECTIVITY_THRESHOLD) + Reflectivity_Attenuated_TL = TEN * Reflectivity_Attenuated_TL / (Reflectivity_Attenuated * LOG(TEN)) + RTSolution_TL%Reflectivity_Attenuated = Reflectivity_Attenuated_TL + RTSolution_TL%Reflectivity_AttenuatedLinear = Reflectivity_AttenuatedLinear_TL + ELSEWHERE + RTSolution_TL%Reflectivity_Attenuated = ZERO + RTSolution_TL%Reflectivity_AttenuatedLinear = ZERO + END WHERE + +END SUBROUTINE CRTM_Compute_Reflectivity_TL !-------------------------------------------------------------------------------- ! @@ -575,197 +525,222 @@ END SUBROUTINE CRTM_Compute_Reflectivity_TL ! ! PURPOSE: ! Subroutine to calculate the adjoint reflectivity for active instruments. +! Includes both logarithmic (dBZ) and linear reflectivity adjoints: +! ReflectivityLinear_AD and Reflectivity_AttenuatedLinear_AD. ! ! CALLING SEQUENCE: ! CALL CRTM_Compute_Reflectivity_AD(Atm, & ! AtmOptics , & ! Input ! RTSolution , & ! Input -! GeometryInfo , & ! Input +! GeometryInfo , & ! Input ! SensorIndex , & ! Input ! ChannelIndex , & ! Input ! AtmOptics_AD , & ! Input/Output ! RTSolution_AD ) ! Input/Output +! ! INPUT ARGUMENTS: ! Atm: Structure containing the atmospheric state data. -! UNITS: N/A ! TYPE: CRTM_Atmosphere_type -! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! AtmOptics: Structure containing the combined atmospheric ! optical properties for gaseous absorption, clouds, ! and aerosols. -! UNITS: N/A ! TYPE: CRTM_AtmOptics_type -! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! RTSolution: Structure containing the solution to the RT equation ! for the given inputs. -! UNITS: N/A ! TYPE: CRTM_RTSolution_type -! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! GeometryInfo: Structure containing the view geometry data. -! UNITS: N/A ! TYPE: CRTM_GeometryInfo_type -! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! -! SensorIndex: Sensor index id. This is a unique index associated -! with a (supported) sensor used to access the -! shared coefficient data for a particular sensor. -! See the ChannelIndex argument. -! UNITS: N/A +! SensorIndex: Sensor index id used to access shared coefficient data. ! TYPE: INTEGER -! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! -! ChannelIndex: Channel index id. This is a unique index associated -! with a (supported) sensor channel used to access the -! shared coefficient data for a particular sensor's -! channel. -! See the SensorIndex argument. -! UNITS: N/A +! ChannelIndex: Channel index id used to access shared coefficient data. ! TYPE: INTEGER -! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! OUTPUT ARGUMENTS: -! ! RTSolution_AD: Structure containing the RT solution adjoint inputs. -! UNITS: N/A ! TYPE: CRTM_RTSolution_type -! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN OUT) ! ! AtmOptics_AD: Structure containing the adjoint combined atmospheric ! optical properties for gaseous absorption, clouds, -! and aerosols. -! UNITS: N/A +! and aerosols, including: +! - ReflectivityLinear_AD +! - Reflectivity_AttenuatedLinear_AD ! TYPE: CRTM_AtmOptics_type -! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN OUT) ! !-------------------------------------------------------------------------------- +SUBROUTINE CRTM_Compute_Reflectivity_AD(Atm, & + AtmOptics, & + RTSolution, & + GeometryInfo, & + SensorIndex, & + ChannelIndex, & + AtmOptics_AD, & + RTSolution_AD) + !----------------------------------------------------------------------- + ! Compute adjoint sensitivities of reflectivity (linear or dBZ) + ! depending on which adjoint variable is set to one in RTSolution_AD. + ! + ! Supported adjoint targets: + ! - RTSolution_AD%Reflectivity_Attenuated (logarithmic dBZ) + ! - RTSolution_AD%Reflectivity_AttenuatedLinear (linear cm^6/m^3) + ! - RTSolution_AD%Reflectivity (logarithmic dBZ) + ! - RTSolution_AD%ReflectivityLinear (linear cm^6/m^3) + ! + !----------------------------------------------------------------------- + + ! Arguments + TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm + TYPE(CRTM_GeometryInfo_type), INTENT(IN) :: GeometryInfo + INTEGER, INTENT(IN) :: SensorIndex + INTEGER, INTENT(IN) :: ChannelIndex + TYPE(CRTM_AtmOptics_type), INTENT(IN) :: AtmOptics + TYPE(CRTM_RTSolution_type), TARGET, INTENT(IN) :: RTSolution + TYPE(CRTM_AtmOptics_type), INTENT(IN OUT) :: AtmOptics_AD + TYPE(CRTM_RTSolution_type), INTENT(IN OUT) :: RTSolution_AD + + ! Locals + INTEGER :: k, j, n_Layers + REAL(fp) :: Frequency, Wavenumber, Wavelength_m + REAL(fp) :: temp_sum + + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Reflectivity, Reflectivity_Attenuated + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: ReflectivityLinear, Reflectivity_AttenuatedLinear + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Transmittance, Optical_Depth + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: P1, Kw_2, perm_re, perm_im + REAL(fp), DIMENSION(0:AtmOptics%n_Layers) :: Height + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: dZ_m, dS_m + COMPLEX :: perm(AtmOptics%n_Layers) + + ! Adjoint variables + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Transmittance_AD, Reflectivity_AD, Reflectivity_Attenuated_AD, Optical_Depth_AD + REAL(fp), DIMENSION(AtmOptics%n_Layers) :: ReflectivityLinear_AD, Reflectivity_AttenuatedLinear_AD + + !----------------------------------------------------------------------- + ! Initialization + !----------------------------------------------------------------------- + n_Layers = Atm%n_Layers + + Transmittance_AD = ZERO + Reflectivity_AD = ZERO + Reflectivity_Attenuated_AD = ZERO + ReflectivityLinear_AD = ZERO + Reflectivity_AttenuatedLinear_AD = ZERO + Optical_Depth_AD = ZERO + AtmOptics_AD%Backscat_Coefficient = ZERO + AtmOptics_AD%optical_depth = ZERO + + !----------------------------------------------------------------------- + ! Heights and geometry + !----------------------------------------------------------------------- + IF (ALL(Atm%Height .LT. EPSILON_FP)) THEN + Height = Calculate_Height(Atm) + ELSE + Height = Atm%Height + END IF - SUBROUTINE CRTM_Compute_Reflectivity_AD(Atm, & - AtmOptics , & ! Input - RTSolution , & ! Input - GeometryInfo , & ! Input - SensorIndex , & ! Input - ChannelIndex , & ! Input - AtmOptics_AD , & ! Input/Output - RTSolution_AD ) ! Input/Output - ! Arguments - TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm - TYPE(CRTM_GeometryInfo_type), INTENT(IN) :: GeometryInfo - TYPE(CRTM_AtmOptics_type) , INTENT(IN) :: AtmOptics - TYPE(CRTM_RTSolution_type), TARGET, INTENT(IN) :: RTSolution - INTEGER , INTENT(IN) :: SensorIndex - INTEGER , INTENT(IN) :: ChannelIndex - TYPE(CRTM_AtmOptics_type) , INTENT(IN OUT) :: AtmOptics_AD - TYPE(CRTM_RTSolution_type), TARGET, INTENT(IN OUT) :: RTSolution_AD - !TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_AD - - REAL(fp) :: Frequency, Wavenumber, Wavelength_m - REAL(fp) :: Transmittance(AtmOptics%n_Layers) - REAL(fp) :: Transmittance_AD(AtmOptics%n_Layers) - REAL(fp) :: P1(Atm%n_Layers) - REAL(fp) :: Height(0:AtmOptics%n_Layers), dZ_m(AtmOptics%n_Layers) - COMPLEX :: perm(Atm%n_Layers) - REAL(fp), DIMENSION(AtmOptics%n_Layers) :: Kw_2, perm_re, perm_im - INTEGER :: k, j - - ! Note Re and Rea are in dBz and Ra is attenuated reflectivity - ! and R and R_AD needs to be calculated or initilized locally - REAL(fp), DIMENSION(AtmOptics%n_Layers) :: R, Ra, R_AD, Ra_AD - REAL(fp), POINTER, DIMENSION(:) :: Re_AD, Rea_AD - REAL(fp), POINTER, DIMENSION(:) :: Re, Rea - NULLIFY(Re, Rea, Re_AD, Rea_AD) - - ! Calculate heights if hasn't been set already - IF (ALL(Atm%Height .LT. EPSILON_FP)) THEN - Height = Calculate_Height(Atm) - ELSE - Height = Atm%Height - ENDIF - dZ_m = (Height(0:Atm%n_Layers-1) - Height(1:Atm%n_Layers)) * ONE_THOUSAND - dZ_m = dZ_m / GeometryInfo%Cosine_Sensor_Zenith - - IF ( SpcCoeff_IsMicrowaveSensor(SC(SensorIndex)) ) THEN - Frequency = SC(SensorIndex)%Frequency(ChannelIndex) ! GHz - Wavelength_m = POINT_01 / GHz_to_Inverse_cm( Frequency ) - ELSE IF( SpcCoeff_IsInfraredSensor(SC(SensorIndex)) ) THEN - Wavenumber = SC(SensorIndex)%Wavenumber(ChannelIndex) ! 1/cm - Wavelength_m = POINT_01 / Wavenumber - END IF - - perm = Water_Permittivity_Turner_2016(Frequency * 1.0d9, & ! Input - Atm%Temperature) ! Input - - perm_re = REAL(REAL(perm)) - perm_im = REAL(AIMAG(perm)) - Kw_2 = ((perm_re - ONE )/(perm_re + TWO))**TWO + dZ_m = Height(0:n_Layers-1) - Height(1:n_Layers) + dS_m = dZ_m / GeometryInfo%Cosine_Sensor_Zenith + + !----------------------------------------------------------------------- + ! Frequency / wavelength + !----------------------------------------------------------------------- + IF (SpcCoeff_IsMicrowaveSensor(SC(SensorIndex))) THEN + Frequency = SC(SensorIndex)%Frequency(ChannelIndex) + Wavelength_m = POINT_01 / GHz_to_Inverse_cm(Frequency) + ELSE IF (SpcCoeff_IsInfraredSensor(SC(SensorIndex))) THEN + Wavenumber = SC(SensorIndex)%Wavenumber(ChannelIndex) + Wavelength_m = POINT_01 / Wavenumber + END IF - ! Calculate transmittance from top to layer k - DO k = 1, AtmOptics%n_layers - Transmittance(k) = EXP(-TWO * SUM(AtmOptics%optical_depth(1:k))) - END DO + !----------------------------------------------------------------------- + ! Permittivity and forward reflectivity computation + !----------------------------------------------------------------------- + perm = Water_Permittivity_Turner_2016(Frequency * 1.0d9, Atm%Temperature) + perm_re = REAL(perm) + perm_im = AIMAG(perm) + Kw_2 = ABS((perm_re - ONE) / (perm_re + ONE))**TWO + + Optical_Depth = AtmOptics%optical_depth / dS_m + + DO k = 1, n_Layers + temp_sum = SUM(Optical_Depth(1:k)) + Transmittance(k) = EXP(-TWO * temp_sum) + END DO + + P1 = (M6_MM6 * Wavelength_m**4.0_fp) / (PI**5.0_fp * Kw_2) + Reflectivity = P1 * (AtmOptics%Backscat_Coefficient / dZ_m) + Reflectivity_Attenuated = Transmittance * Reflectivity + + ! Linear reflectivities (cm^6/m^3) + ReflectivityLinear = MM6_TO_CM6 * Reflectivity + Reflectivity_AttenuatedLinear = MM6_TO_CM6 * Reflectivity_Attenuated + + !----------------------------------------------------------------------- + ! ADJOINT PROPAGATION: determine which output adjoint is active + !----------------------------------------------------------------------- + + ! === Case 1: Logarithmic attenuated reflectivity (dBZ) + IF (ANY(RTSolution_AD%Reflectivity_Attenuated /= ZERO)) THEN + WHERE (Reflectivity_Attenuated > REFLECTIVITY_THRESHOLD) + Reflectivity_Attenuated_AD = Reflectivity_Attenuated_AD + & + TEN * RTSolution_AD%Reflectivity_Attenuated / (Reflectivity_Attenuated * LOG(TEN)) + END WHERE + + ! === Case 2: Linear attenuated reflectivity (cm^6/m^3) + ELSEIF (ANY(RTSolution_AD%Reflectivity_AttenuatedLinear /= ZERO)) THEN + Reflectivity_AttenuatedLinear_AD = Reflectivity_AttenuatedLinear_AD + RTSolution_AD%Reflectivity_AttenuatedLinear + Reflectivity_Attenuated_AD = Reflectivity_Attenuated_AD + MM6_TO_CM6 * Reflectivity_AttenuatedLinear_AD + END IF - P1 = (M6_MM6 * Wavelength_m**4.0_fp) / (PI**5.0_fp * Kw_2) - P1 = P1 / dZ_m ! dZ_m to convert water_content to m/v or cloud water density - R = P1 * AtmOptics%Backscat_Coefficient - Ra = Transmittance * R - - ! This is just to avoid recalculating the effective reflectivities - Re => RTSolution%Reflectivity ! dBz - Rea => RTSolution%Reflectivity_Attenuated ! dBz - - !============================================================================ - ! Adjoint calcualtions - Re_AD => RTSolution_AD%Reflectivity ! dBz - Rea_AD => RTSolution_AD%Reflectivity_Attenuated ! dBz - R_AD = ZERO ! mm^6 m^-3 - Ra_AD = ZERO ! mm^6 m^-3 - Transmittance_AD = ZERO - - WHERE (R .GT. REFLECTIVITY_THRESHOLD) - R_AD = R_AD + TEN * Re_AD / (R * LOG(TEN)) - ELSE WHERE - R_AD = R_AD + ZERO - END WHERE - ! Note that if transmittance is zero then Ra will be zeroo but not R - WHERE (Ra .GT. REFLECTIVITY_THRESHOLD) - Ra_AD = Ra_AD + TEN * Rea_AD / (Ra * LOG(TEN)) - ELSE WHERE - Ra_AD = Ra_AD + ZERO - END WHERE - Transmittance_AD = Transmittance_AD + R * Ra_AD - ! Note that the follwing two lines are mrged into one line - ! R_AD = R_AD + Transmittance * Ra_AD - ! AtmOptics_AD%Backscat_Coefficient = AtmOptics_AD%Backscat_Coefficient + P1 * R_AD - AtmOptics_AD%Backscat_Coefficient = AtmOptics_AD%Backscat_Coefficient + & - P1 * R_AD + & - P1 * Transmittance * Ra_AD - - ! Calculate transmittance from top (satellite) to layer k - !AtmOptics_AD%optical_depth = ZERO - DO k = 1, AtmOptics%n_layers - Do j = 1, k - AtmOptics_AD%optical_depth(j) = AtmOptics_AD%optical_depth(j) & - - TWO * Transmittance(k) * Transmittance_AD(k) - END DO - END DO + ! === Case 3: Logarithmic reflectivity (dBZ) + IF (ANY(RTSolution_AD%Reflectivity /= ZERO)) THEN + WHERE (Reflectivity > REFLECTIVITY_THRESHOLD) + Reflectivity_AD = Reflectivity_AD + & + TEN * RTSolution_AD%Reflectivity / (Reflectivity * LOG(TEN)) + END WHERE - R_AD = ZERO - Ra_AD = ZERO - Re_AD = ZERO - Rea_AD = ZERO - !============================================================================ + ! === Case 4: Linear reflectivity (cm^6/m^3) + ELSEIF (ANY(RTSolution_AD%ReflectivityLinear /= ZERO)) THEN + ReflectivityLinear_AD = ReflectivityLinear_AD + RTSolution_AD%ReflectivityLinear + Reflectivity_AD = Reflectivity_AD + MM6_TO_CM6 * ReflectivityLinear_AD + END IF - END SUBROUTINE CRTM_Compute_Reflectivity_AD + !----------------------------------------------------------------------- + ! Common adjoint propagation for reflectivity relationships + !----------------------------------------------------------------------- + Transmittance_AD = Transmittance_AD + Reflectivity * Reflectivity_Attenuated_AD + Reflectivity_AD = Reflectivity_AD + Transmittance * Reflectivity_Attenuated_AD + + ! Backscatter coefficient adjoint + AtmOptics_AD%Backscat_Coefficient = AtmOptics_AD%Backscat_Coefficient + Reflectivity_AD * P1 / dZ_m + + !----------------------------------------------------------------------- + ! Adjoint of transmittance → optical depth + !----------------------------------------------------------------------- + DO k = n_Layers, 1, -1 + DO j = 1, k + Optical_Depth_AD(j) = Optical_Depth_AD(j) - TWO * Transmittance(k) * Transmittance_AD(k) + END DO + END DO + + !----------------------------------------------------------------------- + ! Adjoint of optical depth normalization + !----------------------------------------------------------------------- + AtmOptics_AD%optical_depth = AtmOptics_AD%optical_depth + Optical_Depth_AD / dS_m + +END SUBROUTINE CRTM_Compute_Reflectivity_AD END MODULE CRTM_Active_Sensor diff --git a/src/RTSolution/CRTM_RTSolution_Define.f90 b/src/RTSolution/CRTM_RTSolution_Define.f90 index d6c72944..dc4f6b74 100644 --- a/src/RTSolution/CRTM_RTSolution_Define.f90 +++ b/src/RTSolution/CRTM_RTSolution_Define.f90 @@ -155,6 +155,8 @@ MODULE CRTM_RTSolution_Define CHARACTER(*), PARAMETER :: SSA_VARNAME = 'Single_Scatter_Albedo' CHARACTER(*), PARAMETER :: ACREFL_VARNAME = 'Reflectivity' ! Active sensor CHARACTER(*), PARAMETER :: ACRATT_VARNAME = 'Reflectivity_Attenuated' ! Active sensor + CHARACTER(*), PARAMETER :: ACREFLL_VARNAME = 'ReflectivityLinear' ! Active sensor + CHARACTER(*), PARAMETER :: ACRATTL_VARNAME = 'Reflectivity_AttenuatedLinear' ! Active sensor CHARACTER(*), PARAMETER :: BS_VARNAME = 'Backscat_Coefficient' ! Variable units attribute. @@ -166,8 +168,10 @@ MODULE CRTM_RTSolution_Define ! ...Emissivity/Reflectivity CHARACTER(*), PARAMETER :: SEMIS_UNITS = 'fraction (0->1)' CHARACTER(*), PARAMETER :: SREFL_UNITS = 'fraction (0->1)' - CHARACTER(*), PARAMETER :: ACREFL_UNITS = 'Reflectivity' ! Active sensor - CHARACTER(*), PARAMETER :: ACRATT_UNITS = 'Reflectivity_Attenuated' ! Active sensor + CHARACTER(*), PARAMETER :: ACREFL_UNITS = 'dBZ' ! Active sensor + CHARACTER(*), PARAMETER :: ACRATT_UNITS = 'dBZ' ! Active sensor + CHARACTER(*), PARAMETER :: ACREFLL_UNITS = 'cm^6/m^3' ! Active sensor + CHARACTER(*), PARAMETER :: ACRATTL_UNITS = 'cm^6/m^3' ! Active sensor ! ...Cloud CHARACTER(*), PARAMETER :: TCC_UNITS = 'fraction (0->1)' CHARACTER(*), PARAMETER :: RCLEAR_UNITS = 'fraction (0->1)' @@ -246,6 +250,8 @@ MODULE CRTM_RTSolution_Define REAL(fp) :: Reflectance = ZERO REAL(fp), ALLOCATABLE :: Reflectivity(:) ! K REAL(fp), ALLOCATABLE :: Reflectivity_Attenuated(:) ! K + REAL(fp), ALLOCATABLE :: ReflectivityLinear(:) ! K + REAL(fp), ALLOCATABLE :: Reflectivity_AttenuatedLinear(:) ! K END TYPE CRTM_RTSolution_type !:tdoc-: @@ -377,6 +383,8 @@ ELEMENTAL SUBROUTINE CRTM_RTSolution_Create( RTSolution, n_Layers ) RTSolution%Single_Scatter_Albedo(n_Layers), & RTSolution%Reflectivity(n_Layers), & RTSolution%Reflectivity_Attenuated(n_Layers), & + RTSolution%ReflectivityLinear(n_Layers), & + RTSolution%Reflectivity_AttenuatedLinear(n_Layers), & RTSolution%Backscat_Coefficient(n_Layers), & STAT = alloc_stat ) IF ( alloc_stat /= 0 ) RETURN @@ -391,6 +399,8 @@ ELEMENTAL SUBROUTINE CRTM_RTSolution_Create( RTSolution, n_Layers ) RTSolution%Single_Scatter_Albedo = ZERO RTSolution%Reflectivity = ZERO RTSolution%Reflectivity_Attenuated = ZERO + RTSolution%ReflectivityLinear = ZERO + RTSolution%Reflectivity_AttenuatedLinear = ZERO RTSolution%Backscat_Coefficient = ZERO ! Set allocation indicator @@ -457,6 +467,8 @@ ELEMENTAL SUBROUTINE CRTM_RTSolution_Zero( RTSolution ) RTSolution%Single_Scatter_Albedo = ZERO RTSolution%Reflectivity = ZERO RTSolution%Reflectivity_Attenuated = ZERO + RTSolution%ReflectivityLinear = ZERO + RTSolution%Reflectivity_AttenuatedLinear = ZERO RTSolution%Backscat_Coefficient = ZERO END IF @@ -547,6 +559,10 @@ SUBROUTINE Scalar_Inspect( RTSolution, Unit ) WRITE(fid,'(5(1x,es22.15,:))') RTSolution%Reflectivity WRITE(fid,'(3x,"Reflectivity_Attenuated :")') WRITE(fid,'(5(1x,es22.15,:))') RTSolution%Reflectivity_Attenuated + WRITE(fid,'(3x,"ReflectivityLinear :")') + WRITE(fid,'(5(1x,es22.15,:))') RTSolution%ReflectivityLinear + WRITE(fid,'(3x,"Reflectivity_AttenuatedLinear :")') + WRITE(fid,'(5(1x,es22.15,:))') RTSolution%Reflectivity_AttenuatedLinear WRITE(fid,'(3x,"Backscat_Coefficient :")') WRITE(fid,'(5(1x,es22.15,:))') RTSolution%Backscat_Coefficient END IF @@ -700,6 +716,8 @@ ELEMENTAL FUNCTION CRTM_RTSolution_Compare( & (.NOT. ALL(Compares_Within_Tolerance(x%Layer_Optical_Depth , y%Layer_Optical_Depth , n))) .OR. & (.NOT. ALL(Compares_Within_Tolerance(x%Reflectivity , y%Reflectivity , n))) .OR. & (.NOT. ALL(Compares_Within_Tolerance(x%Reflectivity_Attenuated , y%Reflectivity_Attenuated , n))) .OR. & + (.NOT. ALL(Compares_Within_Tolerance(x%ReflectivityLinear , y%ReflectivityLinear , n))) .OR. & + (.NOT. ALL(Compares_Within_Tolerance(x%Reflectivity_AttenuatedLinear , y%Reflectivity_AttenuatedLinear , n))) .OR. & (.NOT. ALL(Compares_Within_Tolerance(x%Backscat_Coefficient , y%Backscat_Coefficient , n)))) RETURN END IF @@ -1634,9 +1652,10 @@ FUNCTION CRTM_RTSolution_ReadFile_NetCDF( & REAL(fp), ALLOCATABLE :: Single_Scatter_Albedo(:,:,:) REAL(fp), ALLOCATABLE :: Reflectivity(:,:,:) REAL(fp), ALLOCATABLE :: Reflectivity_Attenuated(:,:,:) + REAL(fp), ALLOCATABLE :: ReflectivityLinear(:,:,:) + REAL(fp), ALLOCATABLE :: Reflectivity_AttenuatedLinear(:,:,:) REAL(fp), ALLOCATABLE :: Backscat_Coefficient(:,:,:) - ! Set up err_stat = SUCCESS Close_File = .FALSE. @@ -1682,6 +1701,8 @@ FUNCTION CRTM_RTSolution_ReadFile_NetCDF( & Single_Scatter_Albedo( n_Channels, n_Layers, n_Profiles ), & Reflectivity( n_Channels, n_Layers, n_Profiles ), & Reflectivity_Attenuated( n_Channels, n_Layers, n_Profiles ), & + ReflectivityLinear( n_Channels, n_Layers, n_Profiles ), & + Reflectivity_AttenuatedLinear( n_Channels, n_Layers, n_Profiles ), & Backscat_Coefficient( n_Channels, n_Layers, n_Profiles ), & STAT = alloc_stat ) IF ( alloc_stat /= 0 ) THEN @@ -2023,6 +2044,32 @@ FUNCTION CRTM_RTSolution_ReadFile_NetCDF( & ' - '//TRIM(NF90_STRERROR( NF90_Status )) CALL Read_Cleanup(); RETURN END IF + ! ... ReflectivityLinear variable + NF90_Status = NF90_INQ_VARID( FileId,ACREFLL_VARNAME,VarId ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error inquiring '//TRIM(Filename)//' for '//ACREFLL_VARNAME//& + ' variable ID - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Read_Cleanup(); RETURN + END IF + NF90_Status = NF90_GET_VAR( FileId,VarID, ReflectivityLinear) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error writing '//ACREFLL_VARNAME//' to '//TRIM(Filename)//& + ' - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Read_Cleanup(); RETURN + END IF + ! ... Reflectivity_AttenuatedLinear variable + NF90_Status = NF90_INQ_VARID( FileId,ACRATTL_VARNAME,VarId ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error inquiring '//TRIM(Filename)//' for '//ACRATTL_VARNAME//& + ' variable ID - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Read_Cleanup(); RETURN + END IF + NF90_Status = NF90_GET_VAR( FileId,VarID, Reflectivity_AttenuatedLinear) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error writing '//ACRATTL_VARNAME//' to '//TRIM(Filename)//& + ' variable ID - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Read_Cleanup(); RETURN + END IF ! ... Backscat_Coefficient variable NF90_Status = NF90_INQ_VARID( FileId,BS_VARNAME,VarId ) IF ( NF90_Status /= NF90_NOERR ) THEN @@ -2098,6 +2145,8 @@ FUNCTION CRTM_RTSolution_ReadFile_NetCDF( & RTSolution(l,m)%Single_Scatter_Albedo(c) = Single_Scatter_Albedo(l,c,m) RTSolution(l,m)%Reflectivity(c) = Reflectivity(l,c,m) RTSolution(l,m)%Reflectivity_Attenuated(c) = Reflectivity_Attenuated(l,c,m) + RTSolution(l,m)%ReflectivityLinear(c) = ReflectivityLinear(l,c,m) + RTSolution(l,m)%Reflectivity_AttenuatedLinear(c) = Reflectivity_AttenuatedLinear(l,c,m) RTSolution(l,m)%Backscat_Coefficient(c) = Backscat_Coefficient(l,c,m) END DO END DO Channel_Loop @@ -2491,6 +2540,8 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & REAL(fp), ALLOCATABLE :: Single_Scatter_Albedo(:,:,:) REAL(fp), ALLOCATABLE :: Reflectivity(:,:,:) REAL(fp), ALLOCATABLE :: Reflectivity_Attenuated(:,:,:) + REAL(fp), ALLOCATABLE :: ReflectivityLinear(:,:,:) + REAL(fp), ALLOCATABLE :: Reflectivity_AttenuatedLinear(:,:,:) REAL(fp), ALLOCATABLE :: Backscat_Coefficient(:,:,:) ! Set up @@ -2536,6 +2587,8 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & Single_Scatter_Albedo( n_Channels, n_Layers, n_Profiles ), & Reflectivity( n_Channels, n_Layers, n_Profiles ), & Reflectivity_Attenuated( n_Channels, n_Layers, n_Profiles ), & + ReflectivityLinear( n_Channels, n_Layers, n_Profiles ), & + Reflectivity_AttenuatedLinear( n_Channels, n_Layers, n_Profiles ), & Backscat_Coefficient( n_Channels, n_Layers, n_Profiles ), & STAT = alloc_stat ) IF ( alloc_stat /= 0 ) THEN @@ -2576,6 +2629,8 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & Single_Scatter_Albedo(l,c,m) = RTSolution(l,m)%Single_Scatter_Albedo(c) Reflectivity(l,c,m) = RTSolution(l,m)%Reflectivity(c) Reflectivity_Attenuated(l,c,m) = RTSolution(l,m)%Reflectivity_Attenuated(c) + ReflectivityLinear(l,c,m) = RTSolution(l,m)%ReflectivityLinear(c) + Reflectivity_AttenuatedLinear(l,c,m) = RTSolution(l,m)%Reflectivity_AttenuatedLinear(c) Backscat_Coefficient(l,c,m) = RTSolution(l,m)%Backscat_Coefficient(c) END DO END DO Channel_Loop @@ -2930,6 +2985,32 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & ' - '//TRIM(NF90_STRERROR( NF90_Status )) CALL Write_Cleanup(); RETURN END IF + ! ... ReflectivityLinear variable + NF90_Status = NF90_INQ_VARID( FileId,ACREFLL_VARNAME,VarId ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error inquiring '//TRIM(Filename)//' for '//ACREFLL_VARNAME//& + ' variable ID - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Write_Cleanup(); RETURN + END IF + NF90_Status = NF90_PUT_VAR( FileId,VarID, ReflectivityLinear) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error writing '//ACREFLL_VARNAME//' to '//TRIM(Filename)//& + ' - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Write_Cleanup(); RETURN + END IF + ! ... Reflectivity_AttenuatedLinear variable + NF90_Status = NF90_INQ_VARID( FileId,ACRATTL_VARNAME,VarId ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error inquiring '//TRIM(Filename)//' for '//ACRATTL_VARNAME//& + ' variable ID - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Write_Cleanup(); RETURN + END IF + NF90_Status = NF90_PUT_VAR( FileId,VarID, Reflectivity_AttenuatedLinear) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error writing '//ACRATTL_VARNAME//' to '//TRIM(Filename)//& + ' - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Write_Cleanup(); RETURN + END IF ! ... Backscat_Coefficient variable NF90_Status = NF90_INQ_VARID( FileId,BS_VARNAME,VarId ) IF ( NF90_Status /= NF90_NOERR ) THEN @@ -2985,6 +3066,8 @@ FUNCTION CRTM_RTSolution_WriteFile_NetCDF( & Single_Scatter_Albedo, & Reflectivity, & Reflectivity_Attenuated, & + ReflectivityLinear, & + Reflectivity_AttenuatedLinear, & Backscat_Coefficient, & STAT = alloc_stat ) IF ( alloc_stat /= 0 ) THEN @@ -3095,6 +3178,8 @@ ELEMENTAL FUNCTION CRTM_RTSolution_Equal( x, y ) RESULT( is_equal ) ALL(x%Single_Scatter_Albedo .EqualTo. y%Single_Scatter_Albedo ) .AND. & ALL(x%Reflectivity .EqualTo. y%Reflectivity ) .AND. & ALL(x%Reflectivity_Attenuated .EqualTo. y%Reflectivity_Attenuated ) .AND. & + ALL(x%ReflectivityLinear .EqualTo. y%ReflectivityLinear ) .AND. & + ALL(x%Reflectivity_AttenuatedLinear .EqualTo. y%Reflectivity_AttenuatedLinear ) .AND. & ALL(x%Backscat_Coefficient .EqualTo. y%Backscat_Coefficient ) END IF @@ -3188,6 +3273,12 @@ ELEMENTAL FUNCTION CRTM_RTSolution_Add( rts1, rts2 ) RESULT( rtssum ) rtssum%Reflectivity_Attenuated(1:k) = rtssum%Reflectivity_Attenuated(1:k) + & rts2%Reflectivity_Attenuated(1:k) + rtssum%ReflectivityLinear(1:k) = rtssum%ReflectivityLinear(1:k) + & + rts2%ReflectivityLinear(1:k) + + rtssum%Reflectivity_AttenuatedLinear(1:k) = rtssum%Reflectivity_AttenuatedLinear(1:k) + & + rts2%Reflectivity_AttenuatedLinear(1:k) + END IF END FUNCTION CRTM_RTSolution_Add @@ -3279,6 +3370,12 @@ ELEMENTAL FUNCTION CRTM_RTSolution_Subtract( rts1, rts2 ) RESULT( rtsdiff ) rtsdiff%Reflectivity_Attenuated(1:k) = rtsdiff%Reflectivity_Attenuated(1:k) - & rts2%Reflectivity_Attenuated(1:k) + + rtsdiff%ReflectivityLinear(1:k) = rtsdiff%ReflectivityLinear(1:k) - & + rts2%ReflectivityLinear(1:k) + + rtsdiff%Reflectivity_AttenuatedLinear(1:k) = rtsdiff%Reflectivity_AttenuatedLinear(1:k) - & + rts2%Reflectivity_AttenuatedLinear(1:k) END IF END FUNCTION CRTM_RTSolution_Subtract @@ -3358,6 +3455,8 @@ ELEMENTAL FUNCTION CRTM_RTSolution_Exponent( rts, power ) RESULT( rts_power ) rts_power%Layer_Optical_Depth(1:k) = (rts_power%Layer_Optical_Depth(1:k) )**power rts_power%Reflectivity(1:k) = (rts_power%Reflectivity(1:k) )**power rts_power%Reflectivity_Attenuated(1:k) = (rts_power%Reflectivity_Attenuated(1:k) )**power + rts_power%ReflectivityLinear(1:k) = (rts_power%ReflectivityLinear(1:k) )**power + rts_power%Reflectivity_AttenuatedLinear(1:k) = (rts_power%Reflectivity_AttenuatedLinear(1:k) )**power END IF END FUNCTION CRTM_RTSolution_Exponent @@ -3438,6 +3537,8 @@ ELEMENTAL FUNCTION CRTM_RTSolution_Normalise( rts, factor ) RESULT( rts_normal ) rts_normal%Layer_Optical_Depth(1:k) = rts_normal%Layer_Optical_Depth(1:k) /factor rts_normal%Reflectivity(1:k) = rts_normal%Reflectivity(1:k) /factor rts_normal%Reflectivity_Attenuated(1:k) = rts_normal%Reflectivity_Attenuated(1:k) /factor + rts_normal%ReflectivityLinear(1:k) = rts_normal%ReflectivityLinear(1:k) /factor + rts_normal%Reflectivity_AttenuatedLinear(1:k) = rts_normal%Reflectivity_AttenuatedLinear(1:k) /factor END IF END FUNCTION CRTM_RTSolution_Normalise @@ -3509,6 +3610,8 @@ ELEMENTAL FUNCTION CRTM_RTSolution_Sqrt( rts ) RESULT( rts_sqrt ) rts_sqrt%Layer_Optical_Depth(1:k) = SQRT(rts_sqrt%Layer_Optical_Depth(1:k) ) rts_sqrt%Reflectivity(1:k) = SQRT(rts_sqrt%Reflectivity(1:k) ) rts_sqrt%Reflectivity_Attenuated(1:k) = SQRT(rts_sqrt%Reflectivity_Attenuated(1:k) ) + rts_sqrt%ReflectivityLinear(1:k) = SQRT(rts_sqrt%ReflectivityLinear(1:k) ) + rts_sqrt%Reflectivity_AttenuatedLinear(1:k) = SQRT(rts_sqrt%Reflectivity_AttenuatedLinear(1:k) ) END IF END FUNCTION CRTM_RTSolution_Sqrt @@ -3606,6 +3709,8 @@ FUNCTION Read_Record( & rts%Layer_Optical_Depth, & rts%Reflectivity, & rts%Reflectivity_Attenuated, & + rts%ReflectivityLinear, & + rts%Reflectivity_AttenuatedLinear, & rts%Backscat_Coefficient IF ( io_stat /= 0 ) THEN msg = 'Error reading array intermediate results - '//TRIM(io_msg) @@ -3721,6 +3826,8 @@ FUNCTION Write_Record( & rts%Layer_Optical_Depth, & rts%Reflectivity, & rts%Reflectivity_Attenuated, & + rts%ReflectivityLinear, & + rts%Reflectivity_AttenuatedLinear, & rts%Backscat_Coefficient IF ( io_stat /= 0 ) THEN msg = 'Error writing array intermediate results - '//TRIM(io_msg) @@ -4287,7 +4394,7 @@ FUNCTION CreateFile_netCDF( & TRIM(Filename)//' - '//TRIM(NF90_STRERROR( NF90_Status )) CALL Create_Cleanup(); RETURN END IF - Put_Status(1) = NF90_PUT_ATT( FileID,VarID,UNITS_ATTNAME ,ACRATT_UNITS) + Put_Status(1) = NF90_PUT_ATT( FileID,VarID,UNITS_ATTNAME ,ACREFL_UNITS) Put_Status(2) = NF90_PUT_ATT( FileID,VarID,FILLVALUE_ATTNAME ,FILL_FLOAT ) IF ( ANY(Put_Status /= NF90_NOERR) ) THEN msg = 'Error writing '//ACREFL_VARNAME//' variable attributes to '//TRIM(Filename) @@ -4311,6 +4418,41 @@ FUNCTION CreateFile_netCDF( & msg = 'Error writing '//ACRATT_VARNAME//' variable attributes to '//TRIM(Filename) CALL Create_Cleanup(); RETURN END IF + ! ... ReflectivityLinear variable + NF90_Status = NF90_DEF_VAR( FileID, & + ACREFLL_VARNAME, & + FLOAT_TYPE, & + dimIDs=(/n_Channels_DimID, n_Layers_DimID, n_Profiles_DimID/), & + varID=VarID ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error defining '//ACREFLL_VARNAME//' variable in '//& + TRIM(Filename)//' - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Create_Cleanup(); RETURN + END IF + Put_Status(1) = NF90_PUT_ATT( FileID,VarID,UNITS_ATTNAME ,ACREFLL_UNITS) + Put_Status(2) = NF90_PUT_ATT( FileID,VarID,FILLVALUE_ATTNAME ,FILL_FLOAT ) + IF ( ANY(Put_Status /= NF90_NOERR) ) THEN + msg = 'Error writing '//ACREFLL_VARNAME//' variable attributes to '//TRIM(Filename) + CALL Create_Cleanup(); RETURN + END IF + + ! ... Reflectivity_AttenuatedLinear variable + NF90_Status = NF90_DEF_VAR( FileID, & + ACRATTL_VARNAME, & + FLOAT_TYPE, & + dimIDs=(/n_Channels_DimID, n_Layers_DimID, n_Profiles_DimID/), & + varID=VarID ) + IF ( NF90_Status /= NF90_NOERR ) THEN + msg = 'Error defining '//ACRATTL_VARNAME//' variable in '//& + TRIM(Filename)//' - '//TRIM(NF90_STRERROR( NF90_Status )) + CALL Create_Cleanup(); RETURN + END IF + Put_Status(1) = NF90_PUT_ATT( FileID,VarID,UNITS_ATTNAME ,ACRATTL_UNITS) + Put_Status(2) = NF90_PUT_ATT( FileID,VarID,FILLVALUE_ATTNAME ,FILL_FLOAT ) + IF ( ANY(Put_Status /= NF90_NOERR) ) THEN + msg = 'Error writing '//ACRATTL_VARNAME//' variable attributes to '//TRIM(Filename) + CALL Create_Cleanup(); RETURN + END IF ! ... Backscat_Coefficient variable NF90_Status = NF90_DEF_VAR( FileID, & diff --git a/test/mains/unit/Unit_Test/test_active_sensor.f90 b/test/mains/unit/Unit_Test/test_active_sensor.f90 index 7e7f1c7b..3ab4be5c 100644 --- a/test/mains/unit/Unit_Test/test_active_sensor.f90 +++ b/test/mains/unit/Unit_Test/test_active_sensor.f90 @@ -332,11 +332,11 @@ PROGRAM test_active_sensor OPEN(500,FILE='reflectivity.txt',STATUS='UNKNOWN') - write(500,'(4A40)') 'Layer', 'Water Content', 'Reflectivity', 'Attenuated Reflectivity' - write(*,'(4A40)') 'Layer', 'Water Content', 'Reflectivity', 'Attenuated Reflectivity' + write(500,'(4A40)') 'Layer', 'Water Content', 'Reflectivity', 'Attenuated Reflectivity', 'ReflectivityLinear', 'Reflectivity_AttenuatedLinear' + write(*,'(4A40)') 'Layer', 'Water Content', 'Reflectivity', 'Attenuated Reflectivity', 'ReflectivityLinear', 'Reflectivity_AttenuatedLinear' DO ii=1,n_layers - write(500,'(i40,f40.5,f40.5,f40.5)') ii, Atm(1)%Cloud(1)%Water_Content(ii), RTSolution(ichan,iprof)%Reflectivity(ii), RTSolution(ichan,iprof)%Reflectivity_Attenuated(ii) - write(*,'(i40,f40.5,f40.5,f40.5)') ii, Atm(1)%Cloud(1)%Water_Content(ii), RTSolution(ichan,iprof)%Reflectivity(ii), RTSolution(ichan,iprof)%Reflectivity_Attenuated(ii) + write(500,'(i40,f40.5,f40.5,f40.5,f40.5,f40.5)') ii, Atm(1)%Cloud(1)%Water_Content(ii), RTSolution(ichan,iprof)%Reflectivity(ii), RTSolution(ichan,iprof)%Reflectivity_Attenuated(ii), RTSolution(ichan,iprof)%ReflectivityLinear(ii), RTSolution(ichan,iprof)%Reflectivity_AttenuatedLinear(ii) + write(*,'(i40,f40.5,f40.5,f40.5,f40.5,f40.5)') ii, Atm(1)%Cloud(1)%Water_Content(ii), RTSolution(ichan,iprof)%Reflectivity(ii), RTSolution(ichan,iprof)%Reflectivity_Attenuated(ii), RTSolution(ichan,iprof)%ReflectivityLinear(ii), RTSolution(ichan,iprof)%Reflectivity_AttenuatedLinear(ii) ENDDO CLOSE(500)