diff --git a/src/AtmOptics/CRTM_AtmOptics.f90 b/src/AtmOptics/CRTM_AtmOptics.f90 index 83ba94fe..b1ded782 100644 --- a/src/AtmOptics/CRTM_AtmOptics.f90 +++ b/src/AtmOptics/CRTM_AtmOptics.f90 @@ -72,6 +72,8 @@ MODULE CRTM_AtmOptics ! --------------- ! Message string length INTEGER, PARAMETER :: ML = 256 + REAL(fp), PARAMETER :: MIN_TRANSMITTANCE = TINY(ONE) + REAL(fp), PARAMETER :: MAX_OPTICAL_DEPTH = -LOG(MIN_TRANSMITTANCE) CONTAINS @@ -211,8 +213,14 @@ SUBROUTINE CRTM_Compute_Transmittance( atmoptics, transmittance ) TYPE(CRTM_AtmOptics_type), INTENT(IN) :: atmoptics REAL(fp) , INTENT(OUT) :: transmittance INTEGER :: k + REAL(fp) :: optical_depth k = atmoptics%n_layers - transmittance = EXP(-ONE*SUM(atmoptics%optical_depth(1:k))) + optical_depth = SUM(atmoptics%optical_depth(1:k)) + IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN + transmittance = MIN_TRANSMITTANCE + ELSE + transmittance = EXP(-optical_depth) + END IF END SUBROUTINE CRTM_Compute_Transmittance !-------------------------------------------------------------------------------- @@ -265,10 +273,16 @@ SUBROUTINE CRTM_Compute_Transmittance_TL( & ! Local variables INTEGER :: k REAL(fp) :: transmittance + REAL(fp) :: optical_depth k = atmoptics%n_layers - transmittance = EXP(-ONE*SUM(atmoptics%optical_depth(1:k))) - transmittance_TL = -transmittance * SUM(atmoptics_TL%optical_depth(1:k)) + optical_depth = SUM(atmoptics%optical_depth(1:k)) + IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN + transmittance_TL = ZERO + ELSE + transmittance = EXP(-optical_depth) + transmittance_TL = -transmittance * SUM(atmoptics_TL%optical_depth(1:k)) + END IF END SUBROUTINE CRTM_Compute_Transmittance_TL @@ -322,9 +336,15 @@ SUBROUTINE CRTM_Compute_Transmittance_AD( & ! Local variables INTEGER :: k REAL(fp) :: transmittance, t_delstar_t + REAL(fp) :: optical_depth k = atmoptics%n_layers - transmittance = EXP(-ONE*SUM(atmoptics%optical_depth(1:k))) + optical_depth = SUM(atmoptics%optical_depth(1:k)) + IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN + transmittance_AD = ZERO + RETURN + END IF + transmittance = EXP(-optical_depth) t_delstar_t = transmittance*transmittance_AD DO k = 1, atmoptics%n_layers atmoptics_AD%optical_depth(k) = atmoptics_AD%optical_depth(k) - t_delstar_t diff --git a/src/CRTM_Forward_Module.f90 b/src/CRTM_Forward_Module.f90 index f0141d01..307439de 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 diff --git a/src/CRTM_K_Matrix_Module.f90 b/src/CRTM_K_Matrix_Module.f90 index 66f6eebf..d7c4424e 100644 --- a/src/CRTM_K_Matrix_Module.f90 +++ b/src/CRTM_K_Matrix_Module.f90 @@ -530,11 +530,293 @@ FUNCTION CRTM_K_Matrix( & WRITE(6,*)'CRTM_K_Matrix inspecting Surface_K...' CALL CRTM_Surface_Inspect (Surface_K(:,:)) END IF + + ! Zero subnormal K-matrix fields introduced by safe exp floors. + CALL Zero_Subnormal_Atmosphere_K(Atmosphere_K) + CALL Zero_Subnormal_Surface_K(Surface_K) + CALL Zero_Subnormal_RTSolution_K(RTSolution_K) RETURN CONTAINS + SUBROUTINE Zero_Subnormal_Atmosphere_K(Atm_K) + TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_K(:,:) + INTEGER :: i_channel, i_profile, i_layer, i_cloud, i_aero, i_abs + INTEGER :: i0, i1, j0, j1 + + DO i_profile = 1, SIZE(Atm_K, DIM=2) + DO i_channel = 1, SIZE(Atm_K, DIM=1) + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Level_Pressure) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Level_Pressure, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Level_Pressure, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Level_Pressure(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Level_Pressure(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Height) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Height, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Height, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Height(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Height(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Pressure) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Pressure, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Pressure, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Pressure(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Pressure(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Temperature) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Temperature, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Temperature, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Temperature(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Temperature(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Absorber) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Absorber, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Absorber, 1) + j0 = LBOUND(Atm_K(i_channel,i_profile)%Absorber, 2) + j1 = UBOUND(Atm_K(i_channel,i_profile)%Absorber, 2) + DO i_abs = j0, j1 + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Absorber(i_layer,i_abs)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Absorber(i_layer,i_abs) = ZERO + END IF + END DO + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Cloud_Fraction) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Cloud_Fraction, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Cloud_Fraction, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Cloud_Fraction(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Cloud_Fraction(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Relative_Humidity) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Relative_Humidity, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Relative_Humidity, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Relative_Humidity(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Relative_Humidity(i_layer) = ZERO + END IF + END DO + END IF + IF ( Atm_K(i_channel,i_profile)%n_Clouds > 0 .AND. ALLOCATED(Atm_K(i_channel,i_profile)%Cloud) ) THEN + DO i_cloud = 1, Atm_K(i_channel,i_profile)%n_Clouds + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Radius) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Radius, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Radius, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Radius(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Radius(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Variance) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Variance, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Variance, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Variance(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Effective_Variance(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Content) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Content, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Content, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Content(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Content(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Density) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Density, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Density, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Density(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Cloud(i_cloud)%Water_Density(i_layer) = ZERO + END IF + END DO + END IF + END DO + END IF + IF ( Atm_K(i_channel,i_profile)%n_Aerosols > 0 .AND. ALLOCATED(Atm_K(i_channel,i_profile)%Aerosol) ) THEN + DO i_aero = 1, Atm_K(i_channel,i_profile)%n_Aerosols + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Radius) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Radius, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Radius, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Radius(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Radius(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Variance) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Variance, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Variance, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Variance(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Effective_Variance(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Concentration) ) THEN + i0 = LBOUND(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Concentration, 1) + i1 = UBOUND(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Concentration, 1) + DO i_layer = i0, i1 + IF ( ABS(Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Concentration(i_layer)) < TINY(ONE) ) THEN + Atm_K(i_channel,i_profile)%Aerosol(i_aero)%Concentration(i_layer) = ZERO + END IF + END DO + END IF + END DO + END IF + END DO + END DO + END SUBROUTINE Zero_Subnormal_Atmosphere_K + + + SUBROUTINE Zero_Subnormal_Surface_K(Sfc_K) + TYPE(CRTM_Surface_type), INTENT(IN OUT) :: Sfc_K(:,:) + INTEGER :: i_channel, i_profile, i_tb + + DO i_profile = 1, SIZE(Sfc_K, DIM=2) + DO i_channel = 1, SIZE(Sfc_K, DIM=1) + IF ( ABS(Sfc_K(i_channel,i_profile)%Land_Coverage) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Land_Coverage = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Water_Coverage) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Water_Coverage = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Snow_Coverage) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Snow_Coverage = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Ice_Coverage) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Ice_Coverage = ZERO + + IF ( ABS(Sfc_K(i_channel,i_profile)%Land_Temperature) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Land_Temperature = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Soil_Moisture_Content) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Soil_Moisture_Content = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Canopy_Water_Content) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Canopy_Water_Content = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Vegetation_Fraction) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Vegetation_Fraction = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Soil_Temperature) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Soil_Temperature = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%LAI) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%LAI = ZERO + + IF ( ABS(Sfc_K(i_channel,i_profile)%Water_Temperature) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Water_Temperature = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Wind_Speed) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Wind_Speed = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Wind_Direction) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Wind_Direction = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Salinity) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Salinity = ZERO + + IF ( ABS(Sfc_K(i_channel,i_profile)%Snow_Temperature) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Snow_Temperature = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Snow_Depth) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Snow_Depth = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Snow_Density) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Snow_Density = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Snow_Grain_Size) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Snow_Grain_Size = ZERO + + IF ( ABS(Sfc_K(i_channel,i_profile)%Ice_Temperature) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Ice_Temperature = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Ice_Thickness) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Ice_Thickness = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Ice_Density) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Ice_Density = ZERO + IF ( ABS(Sfc_K(i_channel,i_profile)%Ice_Roughness) < TINY(ONE) ) Sfc_K(i_channel,i_profile)%Ice_Roughness = ZERO + + IF ( ALLOCATED(Sfc_K(i_channel,i_profile)%SensorData%Tb) ) THEN + DO i_tb = LBOUND(Sfc_K(i_channel,i_profile)%SensorData%Tb, 1), & + UBOUND(Sfc_K(i_channel,i_profile)%SensorData%Tb, 1) + IF ( ABS(Sfc_K(i_channel,i_profile)%SensorData%Tb(i_tb)) < TINY(ONE) ) THEN + Sfc_K(i_channel,i_profile)%SensorData%Tb(i_tb) = ZERO + END IF + END DO + END IF + END DO + END DO + END SUBROUTINE Zero_Subnormal_Surface_K + + + SUBROUTINE Zero_Subnormal_RTSolution_K(RTS_K) + TYPE(CRTM_RTSolution_type), INTENT(IN OUT) :: RTS_K(:,:) + INTEGER :: i_channel, i_profile, i_layer, i_stokes + + DO i_profile = 1, SIZE(RTS_K, DIM=2) + DO i_channel = 1, SIZE(RTS_K, DIM=1) + IF ( ABS(RTS_K(i_channel,i_profile)%SSA_Max) < TINY(ONE) ) RTS_K(i_channel,i_profile)%SSA_Max = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%SOD) < TINY(ONE) ) RTS_K(i_channel,i_profile)%SOD = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Surface_Emissivity) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Surface_Emissivity = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Surface_Reflectivity) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Surface_Reflectivity = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Up_Radiance) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Up_Radiance = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Down_Radiance) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Down_Radiance = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Down_Solar_Radiance) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Down_Solar_Radiance = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Surface_Planck_Radiance) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Surface_Planck_Radiance = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Total_Cloud_Cover) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Total_Cloud_Cover = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%R_clear) < TINY(ONE) ) RTS_K(i_channel,i_profile)%R_clear = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Tb_clear) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Tb_clear = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Reflectance_clear) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Reflectance_clear = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Radiance) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Radiance = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Brightness_Temperature) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Brightness_Temperature = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Solar_Irradiance) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Solar_Irradiance = ZERO + IF ( ABS(RTS_K(i_channel,i_profile)%Reflectance) < TINY(ONE) ) RTS_K(i_channel,i_profile)%Reflectance = ZERO + + DO i_stokes = 1, SIZE(RTS_K(i_channel,i_profile)%Stokes) + IF ( ABS(RTS_K(i_channel,i_profile)%Stokes(i_stokes)) < TINY(ONE) ) THEN + RTS_K(i_channel,i_profile)%Stokes(i_stokes) = ZERO + END IF + END DO + + IF ( ALLOCATED(RTS_K(i_channel,i_profile)%Upwelling_Overcast_Radiance) ) THEN + DO i_layer = LBOUND(RTS_K(i_channel,i_profile)%Upwelling_Overcast_Radiance, 1), & + UBOUND(RTS_K(i_channel,i_profile)%Upwelling_Overcast_Radiance, 1) + IF ( ABS(RTS_K(i_channel,i_profile)%Upwelling_Overcast_Radiance(i_layer)) < TINY(ONE) ) THEN + RTS_K(i_channel,i_profile)%Upwelling_Overcast_Radiance(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(RTS_K(i_channel,i_profile)%Upwelling_Radiance) ) THEN + DO i_layer = LBOUND(RTS_K(i_channel,i_profile)%Upwelling_Radiance, 1), & + UBOUND(RTS_K(i_channel,i_profile)%Upwelling_Radiance, 1) + IF ( ABS(RTS_K(i_channel,i_profile)%Upwelling_Radiance(i_layer)) < TINY(ONE) ) THEN + RTS_K(i_channel,i_profile)%Upwelling_Radiance(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(RTS_K(i_channel,i_profile)%Layer_Optical_Depth) ) THEN + DO i_layer = LBOUND(RTS_K(i_channel,i_profile)%Layer_Optical_Depth, 1), & + UBOUND(RTS_K(i_channel,i_profile)%Layer_Optical_Depth, 1) + IF ( ABS(RTS_K(i_channel,i_profile)%Layer_Optical_Depth(i_layer)) < TINY(ONE) ) THEN + RTS_K(i_channel,i_profile)%Layer_Optical_Depth(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(RTS_K(i_channel,i_profile)%Single_Scatter_Albedo) ) THEN + DO i_layer = LBOUND(RTS_K(i_channel,i_profile)%Single_Scatter_Albedo, 1), & + UBOUND(RTS_K(i_channel,i_profile)%Single_Scatter_Albedo, 1) + IF ( ABS(RTS_K(i_channel,i_profile)%Single_Scatter_Albedo(i_layer)) < TINY(ONE) ) THEN + RTS_K(i_channel,i_profile)%Single_Scatter_Albedo(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(RTS_K(i_channel,i_profile)%Reflectivity) ) THEN + DO i_layer = LBOUND(RTS_K(i_channel,i_profile)%Reflectivity, 1), & + UBOUND(RTS_K(i_channel,i_profile)%Reflectivity, 1) + IF ( ABS(RTS_K(i_channel,i_profile)%Reflectivity(i_layer)) < TINY(ONE) ) THEN + RTS_K(i_channel,i_profile)%Reflectivity(i_layer) = ZERO + END IF + END DO + END IF + IF ( ALLOCATED(RTS_K(i_channel,i_profile)%Reflectivity_Attenuated) ) THEN + DO i_layer = LBOUND(RTS_K(i_channel,i_profile)%Reflectivity_Attenuated, 1), & + UBOUND(RTS_K(i_channel,i_profile)%Reflectivity_Attenuated, 1) + IF ( ABS(RTS_K(i_channel,i_profile)%Reflectivity_Attenuated(i_layer)) < TINY(ONE) ) THEN + RTS_K(i_channel,i_profile)%Reflectivity_Attenuated(i_layer) = ZERO + END IF + END DO + END IF + END DO + END DO + END SUBROUTINE Zero_Subnormal_RTSolution_K + ! Function profile_solution contains all the computational code inside of CRTM_K_Matrix that ! is contained inside an OMP loop. It is "contain"ed inside the function mainly so it can ! access arrays which are arguments to CRTM_K_Matrix. Subroutines Pre_Process_RTSolution_K diff --git a/src/Coefficients/FitCoeff/FitCoeff_SetValue.inc b/src/Coefficients/FitCoeff/FitCoeff_SetValue.inc index be6d1d55..fa2bb394 100644 --- a/src/Coefficients/FitCoeff/FitCoeff_SetValue.inc +++ b/src/Coefficients/FitCoeff/FitCoeff_SetValue.inc @@ -2,11 +2,13 @@ CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'FitCoeff_SetValue' ! Local variables CHARACTER(ML) :: msg + INTEGER :: dimensions(SIZE(self%Dimensions)) INTEGER :: i ! Allocate structure if necessary IF ( .NOT. FitCoeff_Associated(self) ) THEN - CALL FitCoeff_Create( self, SHAPE(C) ) + dimensions = SHAPE(C) + CALL FitCoeff_Create( self, dimensions ) IF ( .NOT. FitCoeff_Associated(self) ) THEN msg = 'Allocation of FitCoeff structure failed' CALL Display_Message( ROUTINE_NAME, msg, FAILURE ); RETURN diff --git a/src/RTSolution/ADA/ADA_Module.f90 b/src/RTSolution/ADA/ADA_Module.f90 index 42553096..fa79560d 100644 --- a/src/RTSolution/ADA/ADA_Module.f90 +++ b/src/RTSolution/ADA/ADA_Module.f90 @@ -39,9 +39,41 @@ MODULE ADA_Module ! ----------------- ! Module parameters ! ----------------- + REAL(fp), PARAMETER :: MIN_TRANSMITTANCE = TINY(ONE) + REAL(fp), PARAMETER :: MAX_OPTICAL_DEPTH = -LOG(MIN_TRANSMITTANCE) CONTAINS + SUBROUTINE Safe_Exp_Neg(x, exp_x, is_clamped) + REAL(fp), INTENT(IN) :: x + REAL(fp), INTENT(OUT) :: exp_x + LOGICAL , INTENT(OUT) :: is_clamped + + IF ( x >= MAX_OPTICAL_DEPTH ) THEN + exp_x = MIN_TRANSMITTANCE + is_clamped = .TRUE. + ELSE + exp_x = EXP(-x) + is_clamped = .FALSE. + END IF + END SUBROUTINE Safe_Exp_Neg + + + SUBROUTINE Safe_Exp_Neg_TL(x, x_TL, exp_x, exp_x_TL) + REAL(fp), INTENT(IN) :: x + REAL(fp), INTENT(IN) :: x_TL + REAL(fp), INTENT(OUT) :: exp_x + REAL(fp), INTENT(OUT) :: exp_x_TL + LOGICAL :: is_clamped + + CALL Safe_Exp_Neg(x, exp_x, is_clamped) + IF ( is_clamped ) THEN + exp_x_TL = ZERO + ELSE + exp_x_TL = -exp_x * x_TL + END IF + END SUBROUTINE Safe_Exp_Neg_TL + !################################################################################ !################################################################################ !## ## @@ -56,9 +88,11 @@ SUBROUTINE CRTM_SurfRef(n_Layers,total_od,dire,Index_Sat_Angle,Rsphere,Ref_0,Ref TYPE(RTV_type), INTENT( INOUT ) :: RTV REAL (fp), INTENT(IN) :: total_od REAL (fp) :: Rsphere, Ref_0, Ref_1, Rfac, dire, SRef(3), S0a(3) + REAL (fp) :: optical_depth, exp_term REAL (fp), DIMENSION(RTV%n_Angles*RTV%n_Stokes, RTV%n_Angles*RTV%n_Stokes) :: temporal_matrix,R0, T0, R2, R1 REAL (fp), DIMENSION( RTV%n_Angles*RTV%n_Stokes) :: refl_down, S0, Sun0 INTEGER :: i, j, k, L, Error_Status, nZ,Index_Sat_Angle + LOGICAL :: is_clamped REAL(fp), PARAMETER :: rFactor(3) = (/ 0.0_fp, 0.5_fp, 1.0_fp/) Sun0(:) = ZERO DO i = 1, RTV%n_Angles @@ -70,7 +104,9 @@ SUBROUTINE CRTM_SurfRef(n_Layers,total_od,dire,Index_Sat_Angle,Rsphere,Ref_0,Ref DO 101 L = 1, 2 R0(:,:) = RTV%s_Level_Refl_UP(1:nZ,1:nZ,n_Layers)*rFactor(L) SRef(L) = dire*rFactor(L) - S0(:)=(ONE-dire)*RTV%Planck_Surface+SRef(L)*RTV%COS_SUN*RTV%Solar_irradiance/PI*exp(-total_od/RTV%COS_SUN)*Sun0(:) + optical_depth = total_od / RTV%COS_SUN + CALL Safe_Exp_Neg(optical_depth, exp_term, is_clamped) + S0(:)=(ONE-dire)*RTV%Planck_Surface+SRef(L)*RTV%COS_SUN*RTV%Solar_irradiance/PI*exp_term*Sun0(:) T0(:,:) = RTV%s_Layer_Trans(1:nZ,1:nZ,n_Layers) DO k = n_Layers, 1, -1 temporal_matrix = -matmul(R0,RTV%s_Layer_Refl(1:nZ,1:nZ,k)) @@ -131,6 +167,8 @@ SUBROUTINE CRTM_ADA(n_Layers, & ! Input number of atmospheric layers REAL (fp), DIMENSION(RTV%n_Angles*RTV%n_Stokes) :: temporal_vector REAL (fp), DIMENSION(0:n_Layers) :: total_opt INTEGER :: i, j, k, Error_Status + REAL(fp) :: optical_depth, exp_term + LOGICAL :: is_clamped CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_ADA' CHARACTER(256) :: Message @@ -156,8 +194,10 @@ SUBROUTINE CRTM_ADA(n_Layers, & ! Input number of atmospheric layers END IF IF( RTV%Solar_Flag_true ) THEN + optical_depth = total_opt(n_Layers) / RTV%COS_SUN + CALL Safe_Exp_Neg(optical_depth, exp_term, is_clamped) RTV%s_Level_Rad_UP(1:nZ,n_Layers ) = RTV%s_Level_Rad_UP(1:nZ,n_Layers )+direct_reflectivity(1:nZ)* & - RTV%COS_SUN*RTV%Solar_irradiance/PI*exp(-total_opt(n_Layers)/RTV%COS_SUN) + RTV%COS_SUN*RTV%Solar_irradiance/PI*exp_term END IF ! UPWARD ADDING LOOP STARTS FROM BOTTOM LAYER TO ATMOSPHERIC TOP LAYER. @@ -226,7 +266,8 @@ SUBROUTINE CRTM_ADA(n_Layers, & ! Input number of atmospheric layers matmul(RTV%Inv_GammaT(1:nZ,1:nZ,k),RTV%Refl_Trans(1:nZ,1:nZ,k)) ELSE DO i = 1, nZ - RTV%s_Layer_Trans(i,i,k) = exp(-T_OD(k)/RTV%COS_AngleS(i)) + optical_depth = T_OD(k)/RTV%COS_AngleS(i) + CALL Safe_Exp_Neg(optical_depth, RTV%s_Layer_Trans(i,i,k), is_clamped) END DO DO i = 1, nZ, RTV%n_Stokes RTV%s_Layer_Source_UP(i,k) = RTV%Planck_Atmosphere(k) * (ONE - RTV%s_Layer_Trans(i,i,k) ) @@ -393,6 +434,7 @@ SUBROUTINE MOM(KL,nZ,optical_depth,trans,refl,RTV, Error_Status) CHARACTER(256) :: Message CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MOM' REAL(fp) :: xx + LOGICAL :: is_clamped REAL(fp), DIMENSION(nZ,nZ) :: tempo RTV%PPM(1:nZ,1:nZ,KL) = RTV%PP(1:nZ,1:nZ,KL) - RTV%PM(1:nZ,1:nZ,KL) @@ -444,7 +486,7 @@ SUBROUTINE MOM(KL,nZ,optical_depth,trans,refl,RTV, Error_Status) DO i = 1, nZ xx = RTV%EigValue(i,KL)*optical_depth - RTV%Exp_x(i,KL) = exp(-xx) + CALL Safe_Exp_Neg(xx, RTV%Exp_x(i,KL), is_clamped) END DO DO i = 1, nZ @@ -489,7 +531,7 @@ SUBROUTINE MOM_TL(KL,nZ,optical_depth,RTV,optical_depth_TL,PP_TL,PM_TL,trans_TL, INTEGER :: i, j CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MOM_TL' - REAL(fp) :: xx_TL + REAL(fp) :: xx_TL, xx REAL(fp), DIMENSION(nZ,nZ) :: PPM_TL, i_PPM_TL, PPP_TL, HH_TL REAL(fp), DIMENSION(nZ,nZ) :: EigVeF_TL,Gp_TL,Gm_TL,A1_TL,A2_TL,A3_TL,A4_TL,A5_TL,A6_TL REAL(fp), DIMENSION(nZ,nZ) :: Gm_A5_TL,i_Gm_A5_TL,EigVe_TL,EigVeVa_TL,i_Gm_TL @@ -529,7 +571,12 @@ SUBROUTINE MOM_TL(KL,nZ,optical_depth,RTV,optical_depth_TL,PP_TL,PM_TL,trans_TL, i_Gm_TL = -matmul( RTV%i_Gm(1:nZ,1:nZ,KL), matmul(Gm_TL,RTV%i_Gm(1:nZ,1:nZ,KL)) ) DO i = 1, nZ xx_TL = EigValue_TL(i)*optical_depth+RTV%EigValue(i,KL)*optical_depth_TL - Exp_x_TL(i) = -xx_TL*RTV%Exp_x(i,KL) + xx = RTV%EigValue(i,KL)*optical_depth + IF ( xx >= MAX_OPTICAL_DEPTH ) THEN + Exp_x_TL(i) = ZERO + ELSE + Exp_x_TL(i) = -xx_TL*RTV%Exp_x(i,KL) + END IF END DO DO i = 1, nZ @@ -573,7 +620,7 @@ SUBROUTINE MOM_AD(KL,nZ,optical_depth,RTV, optical_depth_AD,PP_AD,PM_AD,trans_AD INTEGER :: i, j CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MOM_TL' - REAL(fp) :: xx_AD + REAL(fp) :: xx_AD, xx REAL(fp), DIMENSION(nZ) :: Exp_x_AD,EigValue_AD,EigVa_AD @@ -617,8 +664,14 @@ SUBROUTINE MOM_AD(KL,nZ,optical_depth,RTV, optical_depth_AD,PP_AD,PM_AD,trans_AD END DO DO i = nZ, 1, -1 - xx_AD = -Exp_x_AD(i)*RTV%Exp_x(i,KL) - Exp_x_AD(i) = ZERO + xx = RTV%EigValue(i,KL)*optical_depth + IF ( xx >= MAX_OPTICAL_DEPTH ) THEN + Exp_x_AD(i) = ZERO + xx_AD = ZERO + ELSE + xx_AD = -Exp_x_AD(i)*RTV%Exp_x(i,KL) + Exp_x_AD(i) = ZERO + END IF EigValue_AD(i) = xx_AD*optical_depth optical_depth_AD = optical_depth_AD + RTV%EigValue(i,KL)*xx_AD END DO @@ -714,6 +767,8 @@ SUBROUTINE CRTM_AMOM_layer( n_streams, & ! Input, number of streams INTEGER :: Error_Status REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ) REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ) + REAL(fp) :: optical_depth_sun, total_opt_sun + LOGICAL :: is_clamped CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer' CHARACTER(256) :: Message @@ -807,8 +862,10 @@ SUBROUTINE CRTM_AMOM_layer( n_streams, & ! Input, number of streams ! Solar source Sfactor = single_albedo*RTV%Solar_irradiance/PI IF( RTV%mth_Azi == 0 ) Sfactor = Sfactor/TWO - EXPfactor = exp(-optical_depth/RTV%COS_SUN) - s_transmittance = exp(-total_opt/RTV%COS_SUN) + optical_depth_sun = optical_depth / RTV%COS_SUN + CALL Safe_Exp_Neg(optical_depth_sun, EXPfactor, is_clamped) + total_opt_sun = total_opt / RTV%COS_SUN + CALL Safe_Exp_Neg(total_opt_sun, s_transmittance, is_clamped) DO i = 1, nZ Solar(i) = -bb(i,nZ+1)*Sfactor @@ -1290,6 +1347,7 @@ 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 + REAL (fp) :: optical_depth, optical_depth_TL, exp_term, exp_term_TL INTEGER :: i, j, k, nZ ! nZ = RTV%n_Angles*RTV%n_Stokes @@ -1309,10 +1367,12 @@ SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers 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) + optical_depth = total_opt(n_Layers) / RTV%COS_SUN + optical_depth_TL = total_opt_TL(n_Layers) / RTV%COS_SUN + CALL Safe_Exp_Neg_TL(optical_depth, optical_depth_TL, exp_term, exp_term_TL) + s_rad_up_TL = s_rad_up_TL + direct_reflectivity_TL(1:nZ)*RTV%COS_SUN*RTV%Solar_irradiance/PI & + * exp_term + direct_reflectivity(1:nZ) * RTV%Solar_irradiance/PI & + * RTV%COS_SUN * exp_term_TL END IF DO 10 k = n_Layers, 1, -1 @@ -1465,6 +1525,8 @@ SUBROUTINE CRTM_AMOM_layer_TL( n_streams, & ! Input, number of streams REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ) REAL(fp) :: EXPfactor_TL,Sfactor_TL,s_transmittance_TL,Solar_TL(2*nZ),V0_TL(2*nZ,2*nZ),Solar1_TL(2*nZ) REAL(fp) :: Sfac2_TL + REAL(fp) :: optical_depth_sun, total_opt_sun + REAL(fp) :: optical_depth_sun_TL, total_opt_sun_TL REAL(fp), DIMENSION( nZ ) :: thermal_up_TL,thermal_down_TL REAL(fp), DIMENSION(nZ,nZ) :: trans, refl INTEGER :: N2, N2_1 @@ -1581,11 +1643,13 @@ SUBROUTINE CRTM_AMOM_layer_TL( n_streams, & ! Input, number of streams Sfactor = Sfactor/TWO Sfactor_TL = Sfactor_TL/TWO END IF - EXPfactor = exp(-optical_depth/RTV%COS_SUN) - EXPfactor_TL = -optical_depth_TL/RTV%COS_SUN*EXPfactor + optical_depth_sun = optical_depth / RTV%COS_SUN + optical_depth_sun_TL = optical_depth_TL / RTV%COS_SUN + CALL Safe_Exp_Neg_TL(optical_depth_sun, optical_depth_sun_TL, EXPfactor, EXPfactor_TL) - s_transmittance = exp(-total_opt/RTV%COS_SUN) - s_transmittance_TL = -total_opt_TL/RTV%COS_SUN*s_transmittance + total_opt_sun = total_opt / RTV%COS_SUN + total_opt_sun_TL = total_opt_TL / RTV%COS_SUN + CALL Safe_Exp_Neg_TL(total_opt_sun, total_opt_sun_TL, s_transmittance, s_transmittance_TL) DO i = 1, nZ Solar(i) = -bb(i,nZ+1)*Sfactor @@ -1727,6 +1791,8 @@ 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 + REAL (fp) :: optical_depth, exp_term + LOGICAL :: is_clamped INTEGER :: i, j, k,nZ ! @@ -1853,9 +1919,13 @@ SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers ! IF( RTV%Solar_Flag_true ) THEN - xx = RTV%Solar_irradiance/PI * exp(-total_opt(n_Layers)/RTV%COS_SUN) - total_opt_AD(n_Layers) = total_opt_AD(n_Layers) & - - xx*sum(direct_reflectivity(1:nZ)*s_rad_up_AD(1:nZ)) + optical_depth = total_opt(n_Layers) / RTV%COS_SUN + CALL Safe_Exp_Neg(optical_depth, exp_term, is_clamped) + xx = RTV%Solar_irradiance/PI * exp_term + IF ( .NOT. is_clamped ) THEN + total_opt_AD(n_Layers) = total_opt_AD(n_Layers) & + - xx*sum(direct_reflectivity(1:nZ)*s_rad_up_AD(1:nZ)) + END IF direct_reflectivity_AD(1:nZ) = direct_reflectivity_AD(1:nZ) & + s_rad_up_AD(1:nZ) * RTV%COS_SUN * xx END IF @@ -1940,6 +2010,8 @@ SUBROUTINE CRTM_AMOM_layer_AD( n_streams, & ! Input, number of streams REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ) REAL(fp) :: EXPfactor_AD,Sfactor_AD,s_transmittance_AD,Solar_AD(2*nZ),V0_AD(2*nZ,2*nZ),Solar1_AD(2*nZ) REAL(fp) :: V1_AD(2*nZ,2*nZ),Sfac2_AD + REAL(fp) :: optical_depth_sun, total_opt_sun + LOGICAL :: is_clamped_exp, is_clamped_total REAL(fp), DIMENSION(nZ,nZ) :: trans, refl INTEGER :: N2, N2_1 REAL(fp) :: Thermal_C_AD @@ -2011,8 +2083,10 @@ SUBROUTINE CRTM_AMOM_layer_AD( n_streams, & ! Input, number of streams ! Solar source Sfactor = single_albedo*RTV%Solar_irradiance/PI IF( RTV%mth_Azi == 0 ) Sfactor = Sfactor/TWO - EXPfactor = exp(-optical_depth/RTV%COS_SUN) - s_transmittance = exp(-total_opt/RTV%COS_SUN) + optical_depth_sun = optical_depth / RTV%COS_SUN + CALL Safe_Exp_Neg(optical_depth_sun, EXPfactor, is_clamped_exp) + total_opt_sun = total_opt / RTV%COS_SUN + CALL Safe_Exp_Neg(total_opt_sun, s_transmittance, is_clamped_total) DO i = 1, nZ Solar(i) = -bb(i,nZ+1)*Sfactor @@ -2126,8 +2200,12 @@ SUBROUTINE CRTM_AMOM_layer_AD( n_streams, & ! Input, number of streams bb_AD(i,nZ+1)=bb_AD(i,nZ+1) - Solar_AD(i)*Sfactor ENDDO - total_opt_AD = total_opt_AD -s_transmittance_AD/RTV%COS_SUN*s_transmittance - optical_depth_AD = optical_depth_AD -EXPfactor_AD/RTV%COS_SUN*EXPfactor + IF ( .NOT. is_clamped_total ) THEN + total_opt_AD = total_opt_AD - s_transmittance_AD/RTV%COS_SUN*s_transmittance + END IF + IF ( .NOT. is_clamped_exp ) THEN + optical_depth_AD = optical_depth_AD -EXPfactor_AD/RTV%COS_SUN*EXPfactor + END IF IF( RTV%mth_Azi == 0 ) THEN Sfactor_AD = Sfactor_AD/TWO diff --git a/src/RTSolution/Emission/Emission_Module.f90 b/src/RTSolution/Emission/Emission_Module.f90 index 8e7589de..8a0aa59a 100644 --- a/src/RTSolution/Emission/Emission_Module.f90 +++ b/src/RTSolution/Emission/Emission_Module.f90 @@ -39,6 +39,8 @@ MODULE Emission_Module ! Version Id for the module CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = & '$Id: $' + REAL(fp), PARAMETER :: MIN_TRANSMITTANCE = TINY(ONE) + REAL(fp), PARAMETER :: MAX_OPTICAL_DEPTH = -LOG(MIN_TRANSMITTANCE) CONTAINS @@ -99,7 +101,8 @@ SUBROUTINE CRTM_Emission(n_Layers, & ! Input number of atmospheric layers REAL(fp), INTENT(IN) :: Source_Zenith_Radian TYPE(RTV_type), INTENT(IN OUT) :: RTV ! Local variables - REAL(fp) :: layer_source_up, cosine_u0 + REAL(fp) :: layer_source_up, cosine_u0 + REAL(fp) :: layer_optical_depth, optical_depth, exp_term INTEGER :: k ! -------------------- @@ -121,12 +124,22 @@ SUBROUTINE CRTM_Emission(n_Layers, & ! Input number of atmospheric layers ! Accumulate optical depth RTV%Total_OD = RTV%Total_OD + T_OD(k) ! Layer downward transmittance - RTV%e_Layer_Trans_DOWN(k) = EXP(-T_OD(k)*RTV%Secant_Down_Angle) + layer_optical_depth = T_OD(k) * RTV%Secant_Down_Angle + IF ( layer_optical_depth >= MAX_OPTICAL_DEPTH ) THEN + RTV%e_Layer_Trans_DOWN(k) = MIN_TRANSMITTANCE + ELSE + RTV%e_Layer_Trans_DOWN(k) = EXP(-layer_optical_depth) + END IF ! Downward radiance RTV%e_Level_Rad_DOWN(k) = (RTV%e_Level_Rad_DOWN(k-1)*RTV%e_Layer_Trans_DOWN(k)) + & (Planck_Atmosphere(k)*(ONE-RTV%e_Layer_Trans_DOWN(k))) - RTV%e_Layer_Trans_UP(k) = EXP(-T_OD(k)/u) + layer_optical_depth = T_OD(k) / u + IF ( layer_optical_depth >= MAX_OPTICAL_DEPTH ) THEN + RTV%e_Layer_Trans_UP(k) = MIN_TRANSMITTANCE + ELSE + RTV%e_Layer_Trans_UP(k) = EXP(-layer_optical_depth) + END IF ! GSI cloud detection RTV%e_Cloud_Radiance_UP(k) = RTV%e_Source_UP(k-1) + Planck_Atmosphere(k)*RTV%e_Level_Trans_UP(k-1) @@ -146,7 +159,13 @@ SUBROUTINE CRTM_Emission(n_Layers, & ! Input number of atmospheric layers IF( Is_Solar_Channel ) THEN cosine_u0 = COS(Source_Zenith_Radian) IF( cosine_u0 > ZERO) THEN - RTV%Down_Solar_Radiance = cosine_u0*EXP(-RTV%Total_OD/cosine_u0)*Solar_Irradiance/PI + optical_depth = RTV%Total_OD / cosine_u0 + IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN + exp_term = MIN_TRANSMITTANCE + ELSE + exp_term = EXP(-optical_depth) + END IF + RTV%Down_Solar_Radiance = cosine_u0 * exp_term * Solar_Irradiance / PI RTV%e_Level_Rad_UP(n_Layers) = RTV%e_Level_Rad_UP(n_Layers) + & (RTV%Down_Solar_Radiance*direct_reflectivity(1)) END IF @@ -215,7 +234,7 @@ SUBROUTINE CRTM_Emission_TL(n_Layers, & ! Input number of atmospheric layers REAL (fp) :: layer_source_up_TL, layer_source_down_TL,a_TL,down_rad_TL REAL (fp) :: Total_OD, Total_OD_TL INTEGER :: k - REAL( fp) :: cosine_u0 + REAL( fp) :: cosine_u0, optical_depth, exp_term !#--------------------------------------------------------------------------# !# -- Downwelling TL radiance -- # @@ -252,10 +271,18 @@ SUBROUTINE CRTM_Emission_TL(n_Layers, & ! Input number of atmospheric layers IF( Is_Solar_Channel ) THEN cosine_u0 = cos(Source_Zenith_Radian) IF( cosine_u0 > ZERO) THEN - up_rad_TL = up_rad_TL + cosine_u0*Solar_Irradiance/PI & - * direct_reflectivity_TL(1) * exp(-Total_OD/cosine_u0) & - - Solar_Irradiance/PI * direct_reflectivity(1) & - * Total_OD_TL * exp(-Total_OD/cosine_u0) + optical_depth = Total_OD / cosine_u0 + IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN + exp_term = MIN_TRANSMITTANCE + up_rad_TL = up_rad_TL + cosine_u0*Solar_Irradiance/PI & + * direct_reflectivity_TL(1) * exp_term + ELSE + exp_term = EXP(-optical_depth) + up_rad_TL = up_rad_TL + cosine_u0*Solar_Irradiance/PI & + * direct_reflectivity_TL(1) * exp_term & + - Solar_Irradiance/PI * direct_reflectivity(1) & + * Total_OD_TL * exp_term + ENDIF ENDIF ENDIF @@ -321,6 +348,7 @@ SUBROUTINE CRTM_Emission_AD(n_Layers, & ! Input number of atmospheric layers ! internal variables 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 + REAL (fp) :: optical_depth, exp_term INTEGER :: k ! ! Initialize variables @@ -358,10 +386,17 @@ SUBROUTINE CRTM_Emission_AD(n_Layers, & ! Input number of atmospheric layers IF( Is_Solar_Channel ) THEN cosine_u0 = cos(Source_Zenith_Radian) IF( cosine_u0 > ZERO) THEN - Total_OD_AD = -Solar_Irradiance/PI * direct_reflectivity(1) & - * up_rad_AD * exp(-Total_OD/cosine_u0) + optical_depth = Total_OD / cosine_u0 + IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN + exp_term = MIN_TRANSMITTANCE + Total_OD_AD = ZERO + ELSE + exp_term = EXP(-optical_depth) + Total_OD_AD = -Solar_Irradiance/PI * direct_reflectivity(1) & + * up_rad_AD * exp_term + END IF direct_reflectivity_AD(1) = cosine_u0 * Solar_Irradiance/PI & - * up_rad_AD* exp(-Total_OD/cosine_u0) + * up_rad_AD * exp_term ENDIF ENDIF diff --git a/src/RTSolution/SOI/SOI_Module.f90 b/src/RTSolution/SOI/SOI_Module.f90 index 8cca3043..4d41fbf2 100644 --- a/src/RTSolution/SOI/SOI_Module.f90 +++ b/src/RTSolution/SOI/SOI_Module.f90 @@ -47,6 +47,8 @@ MODULE SOI_Module ! ----------------- ! Module parameters ! ----------------- + REAL(fp), PARAMETER :: MIN_TRANSMITTANCE = TINY(ONE) + REAL(fp), PARAMETER :: MAX_OPTICAL_DEPTH = -LOG(MIN_TRANSMITTANCE) CONTAINS @@ -98,11 +100,19 @@ SUBROUTINE CRTM_SOI(n_Layers, & ! Input number of atmospheric layers REAL(fp), PARAMETER :: OPT_DEPTH_THRESH = 4.0 REAL(fp) :: radiance_thresh REAL(fp), DIMENSION( MAX_N_ANGLES ) :: source + REAL(fp) :: optical_depth ! Precompute layer R/T matrices and thermal sources DO k = 1, n_Layers ! Precompute simple layer properties - RTV%e_Layer_Trans( 1 : RTV%n_Angles, k ) = EXP( -T_OD( k ) / RTV%COS_Angle( 1 : RTV%n_Angles ) ) + DO i = 1, RTV%n_Angles + optical_depth = T_OD(k) / RTV%COS_Angle(i) + IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN + RTV%e_Layer_Trans(i, k) = MIN_TRANSMITTANCE + ELSE + RTV%e_Layer_Trans(i, k) = EXP(-optical_depth) + END IF + END DO IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD ) THEN IF ( w( k ) < SNGL_SCAT_ALB_THRESH .AND. T_OD( k ) < OPT_DEPTH_THRESH ) THEN @@ -335,12 +345,19 @@ 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 + REAL(fp) :: optical_depth DO k = 1, n_Layers - e_Trans_TL( 1 : RTV%n_Angles, k ) = -T_OD_TL(k) * (EXP( -T_OD( k ) / RTV%COS_Angle( 1 : RTV%n_Angles ) ) ) / & - RTV%COS_Angle( 1 : RTV%n_Angles ) + DO i = 1, RTV%n_Angles + optical_depth = T_OD(k) / RTV%COS_Angle(i) + IF ( optical_depth >= MAX_OPTICAL_DEPTH ) THEN + e_Trans_TL(i, k) = ZERO + ELSE + e_Trans_TL(i, k) = -T_OD_TL(k) * RTV%e_Layer_Trans(i, k) / RTV%COS_Angle(i) + END IF + END DO IF ( w( k ) > SCATTERING_ALBEDO_THRESHOLD ) THEN IF ( w( k ) < SNGL_SCAT_ALB_THRESH .AND. T_OD( k ) < OPT_DEPTH_THRESH ) THEN @@ -525,6 +542,7 @@ SUBROUTINE CRTM_SOI_AD(n_Layers, & ! Input number of atmospheric layers REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: e_Trans_AD REAL(fp), DIMENSION( RTV%n_Angles, RTV%n_Angles, n_Layers ) :: s_Refl_AD REAL(fp), DIMENSION( RTV%n_Angles, RTV%n_Angles, n_Layers ) :: s_Trans_AD + REAL(fp) :: optical_depth ! Zero out all local variables s_IterRad_UP_AD = ZERO @@ -678,8 +696,12 @@ SUBROUTINE CRTM_SOI_AD(n_Layers, & ! Input number of atmospheric layers END DO END IF - T_OD_AD( k ) = T_OD_AD( k ) - SUM( e_Trans_AD( :, k ) * RTV%e_Layer_Trans( 1:RTV%n_Angles, k ) / & - RTV%COS_Angle(1:RTV%n_Angles) ) + DO i = 1, RTV%n_Angles + optical_depth = T_OD(k) / RTV%COS_Angle(i) + IF ( optical_depth < MAX_OPTICAL_DEPTH ) THEN + T_OD_AD( k ) = T_OD_AD( k ) - e_Trans_AD(i, k) * RTV%e_Layer_Trans(i, k) / RTV%COS_Angle(i) + END IF + END DO e_Trans_AD( 1 : RTV%n_Angles, k ) = ZERO END DO @@ -1182,7 +1204,11 @@ SUBROUTINE CRTM_Doubling_layer(n_streams, & ! Input, number of streams print *,' error at matinv in CRTM_Doubling_layer ' RTV%s_Layer_Trans(1:NANG,1:NANG,KL) = ZERO DO i = 1, NANG - RTV%s_Layer_Trans(i,i,KL) = exp(-optical_depth/COS_Angle(i)) + IF ( optical_depth/COS_Angle(i) >= MAX_OPTICAL_DEPTH ) THEN + RTV%s_Layer_Trans(i,i,KL) = MIN_TRANSMITTANCE + ELSE + RTV%s_Layer_Trans(i,i,KL) = exp(-optical_depth/COS_Angle(i)) + END IF ENDDO RTV%s_Layer_Refl(1:NANG,1:NANG,KL) = ZERO RTV%s_Layer_Source_DOWN(1:NANG,KL) = ZERO @@ -1482,4 +1508,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/MW_Water/FASTEM_MWSSEM/Azimuth_Emissivity_F6_Module.f90 b/src/SfcOptics/MW_Water/FASTEM_MWSSEM/Azimuth_Emissivity_F6_Module.f90 index 2101ff47..cc9679f0 100644 --- a/src/SfcOptics/MW_Water/FASTEM_MWSSEM/Azimuth_Emissivity_F6_Module.f90 +++ b/src/SfcOptics/MW_Water/FASTEM_MWSSEM/Azimuth_Emissivity_F6_Module.f90 @@ -47,6 +47,8 @@ MODULE Azimuth_Emissivity_F6_Module REAL(fp), PARAMETER :: FOUR = 4.0_fp REAL(fp), PARAMETER :: PI = 3.141592653589793238462643383279_fp REAL(fp), PARAMETER :: DEGREES_TO_RADIANS = PI / 180.0_fp + REAL(fp), PARAMETER :: MIN_EXP = TINY(ONE) + REAL(fp), PARAMETER :: MAX_EXP_ARG = -LOG(MIN_EXP) ! Dimensions ! ...Number of Stokes parameters handled @@ -131,6 +133,8 @@ SUBROUTINE Azimuth_Emissivity_F6( & TYPE(iVar_type) , INTENT(IN OUT) :: iVar ! Local variables INTEGER :: j + REAL(fp) :: exp_arg_v, exp_arg_h + REAL(fp) :: exp_term_v, exp_term_h ! Initialise output e_Azimuth = ZERO @@ -152,14 +156,27 @@ SUBROUTINE Azimuth_Emissivity_F6( & ! Loop over frequencies to compute the intermediate terms Frequency_Loop: DO j = 1, N_FREQUENCIES - iVar%A1v(j) = AZCoeff%C(1,j,IVPOL) * ( EXP(-AZCoeff%C(5,j,IVPOL) * iVar%w18**2 ) - ONE ) * & + exp_arg_v = AZCoeff%C(5,j,IVPOL) * iVar%w18**2 + IF ( exp_arg_v >= MAX_EXP_ARG ) THEN + exp_term_v = MIN_EXP + ELSE + exp_term_v = EXP(-exp_arg_v) + END IF + exp_arg_h = AZCoeff%C(6,j,IHPOL) * iVar%w18**2 + IF ( exp_arg_h >= MAX_EXP_ARG ) THEN + exp_term_h = MIN_EXP + ELSE + exp_term_h = EXP(-exp_arg_h) + END IF + + iVar%A1v(j) = AZCoeff%C(1,j,IVPOL) * ( exp_term_v - ONE ) * & ( AZCoeff%C(2,j,IVPOL) * iVar%w18 + & AZCoeff%C(3,j,IVPOL) * iVar%w18**2 + & AZCoeff%C(4,j,IVPOL) * iVar%w18**3 ) iVar%A2v(j) = AZCoeff%C(6,j,IVPOL) * iVar%w18 iVar%A1h(j) = AZCoeff%C(1,j,IHPOL) * iVar%w18 - iVar%A2h(j) = AZCoeff%C(2,j,IHPOL) * ( EXP(-AZCoeff%C(6,j,IHPOL) * iVar%w18**2 ) - ONE ) * & + iVar%A2h(j) = AZCoeff%C(2,j,IHPOL) * ( exp_term_h - ONE ) * & ( AZCoeff%C(3,j,IHPOL) * iVar%w18 + & AZCoeff%C(4,j,IHPOL) * iVar%w18**2 + & AZCoeff%C(5,j,IHPOL) * iVar%w18**3 ) @@ -239,6 +256,9 @@ SUBROUTINE Azimuth_Emissivity_F6_TL( & REAL(fp), DIMENSION(N_FREQUENCIES) :: A1s1_theta_TL, A1s2_theta_TL, A2s1_theta_TL, A2s2_theta_TL REAL(fp), DIMENSION(N_FREQUENCIES) :: A1v_theta_TL , A1h_theta_TL , A2v_theta_TL , A2h_theta_TL REAL(fp) :: azimuth_component_TL(N_FREQUENCIES, N_STOKES) + REAL(fp) :: exp_arg_v, exp_arg_h + REAL(fp) :: exp_term_v, exp_term_h + LOGICAL :: exp_clamped_v, exp_clamped_h ! Initialise output e_Azimuth_TL = ZERO @@ -256,28 +276,59 @@ SUBROUTINE Azimuth_Emissivity_F6_TL( & A1h_TL(j) = ZERO A2h_TL(j) = ZERO ELSE - A1v_TL(j) = ( AZCoeff%C(1,j,IVPOL) * ( EXP(-AZCoeff%C(5,j,IVPOL) * iVar%w18**2 ) - ONE ) * & - ( AZCoeff%C(2,j,IVPOL) + & - TWO * AZCoeff%C(3,j,IVPOL) * iVar%w18 + & - THREE * AZCoeff%C(4,j,IVPOL) * iVar%w18**2 ) - & - TWO * AZCoeff%C(1,j,IVPOL) * AZCoeff%C(5,j,IVPOL) * iVar%w18 * & - EXP(-AZCoeff%C(5,j,IVPOL) * iVar%w18**2 ) * & - ( AZCoeff%C(2,j,IVPOL) * iVar%w18 + & - AZCoeff%C(3,j,IVPOL) * iVar%w18**2 + & - AZCoeff%C(4,j,IVPOL) * iVar%w18**3 ) ) * Wind_Speed_TL + exp_arg_v = AZCoeff%C(5,j,IVPOL) * iVar%w18**2 + IF ( exp_arg_v >= MAX_EXP_ARG ) THEN + exp_term_v = MIN_EXP + exp_clamped_v = .TRUE. + ELSE + exp_term_v = EXP(-exp_arg_v) + exp_clamped_v = .FALSE. + END IF + exp_arg_h = AZCoeff%C(6,j,IHPOL) * iVar%w18**2 + IF ( exp_arg_h >= MAX_EXP_ARG ) THEN + exp_term_h = MIN_EXP + exp_clamped_h = .TRUE. + ELSE + exp_term_h = EXP(-exp_arg_h) + exp_clamped_h = .FALSE. + END IF + + IF ( exp_clamped_v ) THEN + A1v_TL(j) = ( AZCoeff%C(1,j,IVPOL) * ( exp_term_v - ONE ) * & + ( AZCoeff%C(2,j,IVPOL) + & + TWO * AZCoeff%C(3,j,IVPOL) * iVar%w18 + & + THREE * AZCoeff%C(4,j,IVPOL) * iVar%w18**2 ) ) * Wind_Speed_TL + ELSE + A1v_TL(j) = ( AZCoeff%C(1,j,IVPOL) * ( exp_term_v - ONE ) * & + ( AZCoeff%C(2,j,IVPOL) + & + TWO * AZCoeff%C(3,j,IVPOL) * iVar%w18 + & + THREE * AZCoeff%C(4,j,IVPOL) * iVar%w18**2 ) - & + TWO * AZCoeff%C(1,j,IVPOL) * AZCoeff%C(5,j,IVPOL) * iVar%w18 * & + exp_term_v * & + ( AZCoeff%C(2,j,IVPOL) * iVar%w18 + & + AZCoeff%C(3,j,IVPOL) * iVar%w18**2 + & + AZCoeff%C(4,j,IVPOL) * iVar%w18**3 ) ) * Wind_Speed_TL + END IF A2v_TL(j) = AZCoeff%C(6,j,IVPOL) * Wind_Speed_TL A1h_TL(j) = AZCoeff%C(1,j,IHPOL) * Wind_Speed_TL - A2h_TL(j) = ( AZCoeff%C(2,j,IHPOL) * ( EXP(-AZCoeff%C(6,j,IHPOL) * iVar%w18**2 ) - ONE ) * & - ( AZCoeff%C(3,j,IHPOL) + & - TWO * AZCoeff%C(4,j,IHPOL) * iVar%w18 + & - THREE * AZCoeff%C(5,j,IHPOL) * iVar%w18**2 ) - & - TWO * AZCoeff%C(2,j,IHPOL) * AZCoeff%C(6,j,IHPOL) * iVar%w18 * & - EXP(-AZCoeff%C(6,j,IHPOL) * iVar%w18**2 ) * & - ( AZCoeff%C(3,j,IHPOL) * iVar%w18 + & - AZCoeff%C(4,j,IHPOL) * iVar%w18**2 + & - AZCoeff%C(5,j,IHPOL) * iVar%w18**3 ) ) * Wind_Speed_TL + IF ( exp_clamped_h ) THEN + A2h_TL(j) = ( AZCoeff%C(2,j,IHPOL) * ( exp_term_h - ONE ) * & + ( AZCoeff%C(3,j,IHPOL) + & + TWO * AZCoeff%C(4,j,IHPOL) * iVar%w18 + & + THREE * AZCoeff%C(5,j,IHPOL) * iVar%w18**2 ) ) * Wind_Speed_TL + ELSE + A2h_TL(j) = ( AZCoeff%C(2,j,IHPOL) * ( exp_term_h - ONE ) * & + ( AZCoeff%C(3,j,IHPOL) + & + TWO * AZCoeff%C(4,j,IHPOL) * iVar%w18 + & + THREE * AZCoeff%C(5,j,IHPOL) * iVar%w18**2 ) - & + TWO * AZCoeff%C(2,j,IHPOL) * AZCoeff%C(6,j,IHPOL) * iVar%w18 * & + exp_term_h * & + ( AZCoeff%C(3,j,IHPOL) * iVar%w18 + & + AZCoeff%C(4,j,IHPOL) * iVar%w18**2 + & + AZCoeff%C(5,j,IHPOL) * iVar%w18**3 ) ) * Wind_Speed_TL + END IF END IF @@ -366,6 +417,9 @@ SUBROUTINE Azimuth_Emissivity_F6_AD( & REAL(fp), DIMENSION(N_FREQUENCIES) :: A1s1_theta_AD, A1s2_theta_AD, A2s1_theta_AD, A2s2_theta_AD REAL(fp), DIMENSION(N_FREQUENCIES) :: A1v_theta_AD , A1h_theta_AD , A2v_theta_AD , A2h_theta_AD REAL(fp) :: azimuth_component_AD(N_FREQUENCIES, N_STOKES) + REAL(fp) :: exp_arg_v, exp_arg_h + REAL(fp) :: exp_term_v, exp_term_h + LOGICAL :: exp_clamped_v, exp_clamped_h ! Initialise local adjoint variables phi_AD = ZERO @@ -482,16 +536,41 @@ SUBROUTINE Azimuth_Emissivity_F6_AD( & A1h_AD(j) = ZERO A2h_AD(j) = ZERO ELSE - Wind_Speed_AD = Wind_Speed_AD + & - ( AZCoeff%C(2,j,IHPOL) * ( EXP(-AZCoeff%C(6,j,IHPOL) * iVar%w18**2 ) - ONE ) * & - ( AZCoeff%C(3,j,IHPOL) + & - TWO * AZCoeff%C(4,j,IHPOL) * iVar%w18 + & - THREE * AZCoeff%C(5,j,IHPOL) * iVar%w18**2 ) - & - TWO * AZCoeff%C(2,j,IHPOL) * AZCoeff%C(6,j,IHPOL) * iVar%w18 * & - EXP(-AZCoeff%C(6,j,IHPOL) * iVar%w18**2 ) * & - ( AZCoeff%C(3,j,IHPOL) * iVar%w18 + & - AZCoeff%C(4,j,IHPOL) * iVar%w18**2 + & - AZCoeff%C(5,j,IHPOL) * iVar%w18**3 ) ) * A2h_AD(j) + exp_arg_v = AZCoeff%C(5,j,IVPOL) * iVar%w18**2 + IF ( exp_arg_v >= MAX_EXP_ARG ) THEN + exp_term_v = MIN_EXP + exp_clamped_v = .TRUE. + ELSE + exp_term_v = EXP(-exp_arg_v) + exp_clamped_v = .FALSE. + END IF + exp_arg_h = AZCoeff%C(6,j,IHPOL) * iVar%w18**2 + IF ( exp_arg_h >= MAX_EXP_ARG ) THEN + exp_term_h = MIN_EXP + exp_clamped_h = .TRUE. + ELSE + exp_term_h = EXP(-exp_arg_h) + exp_clamped_h = .FALSE. + END IF + + IF ( exp_clamped_h ) THEN + Wind_Speed_AD = Wind_Speed_AD + & + ( AZCoeff%C(2,j,IHPOL) * ( exp_term_h - ONE ) * & + ( AZCoeff%C(3,j,IHPOL) + & + TWO * AZCoeff%C(4,j,IHPOL) * iVar%w18 + & + THREE * AZCoeff%C(5,j,IHPOL) * iVar%w18**2 ) ) * A2h_AD(j) + ELSE + Wind_Speed_AD = Wind_Speed_AD + & + ( AZCoeff%C(2,j,IHPOL) * ( exp_term_h - ONE ) * & + ( AZCoeff%C(3,j,IHPOL) + & + TWO * AZCoeff%C(4,j,IHPOL) * iVar%w18 + & + THREE * AZCoeff%C(5,j,IHPOL) * iVar%w18**2 ) - & + TWO * AZCoeff%C(2,j,IHPOL) * AZCoeff%C(6,j,IHPOL) * iVar%w18 * & + exp_term_h * & + ( AZCoeff%C(3,j,IHPOL) * iVar%w18 + & + AZCoeff%C(4,j,IHPOL) * iVar%w18**2 + & + AZCoeff%C(5,j,IHPOL) * iVar%w18**3 ) ) * A2h_AD(j) + END IF A2h_AD(j) = ZERO Wind_Speed_AD = Wind_Speed_AD + AZCoeff%C(1,j,IHPOL)*A1h_AD(j) @@ -500,16 +579,24 @@ SUBROUTINE Azimuth_Emissivity_F6_AD( & Wind_Speed_AD = Wind_Speed_AD + AZCoeff%C(6,j,IVPOL)*A2v_AD(j) A2v_AD(j) = ZERO - Wind_Speed_AD = Wind_Speed_AD + & - ( AZCoeff%C(1,j,IVPOL) * ( EXP(-AZCoeff%C(5,j,IVPOL) * iVar%w18**2 ) - ONE ) * & - ( AZCoeff%C(2,j,IVPOL) + & - TWO * AZCoeff%C(3,j,IVPOL) * iVar%w18 + & - THREE * AZCoeff%C(4,j,IVPOL) * iVar%w18**2 ) - & - TWO * AZCoeff%C(1,j,IVPOL) * AZCoeff%C(5,j,IVPOL) * iVar%w18 * & - EXP(-AZCoeff%C(5,j,IVPOL) * iVar%w18**2 ) * & - ( AZCoeff%C(2,j,IVPOL) * iVar%w18 + & - AZCoeff%C(3,j,IVPOL) * iVar%w18**2 + & - AZCoeff%C(4,j,IVPOL) * iVar%w18**3 ) ) * A1v_AD(j) + IF ( exp_clamped_v ) THEN + Wind_Speed_AD = Wind_Speed_AD + & + ( AZCoeff%C(1,j,IVPOL) * ( exp_term_v - ONE ) * & + ( AZCoeff%C(2,j,IVPOL) + & + TWO * AZCoeff%C(3,j,IVPOL) * iVar%w18 + & + THREE * AZCoeff%C(4,j,IVPOL) * iVar%w18**2 ) ) * A1v_AD(j) + ELSE + Wind_Speed_AD = Wind_Speed_AD + & + ( AZCoeff%C(1,j,IVPOL) * ( exp_term_v - ONE ) * & + ( AZCoeff%C(2,j,IVPOL) + & + TWO * AZCoeff%C(3,j,IVPOL) * iVar%w18 + & + THREE * AZCoeff%C(4,j,IVPOL) * iVar%w18**2 ) - & + TWO * AZCoeff%C(1,j,IVPOL) * AZCoeff%C(5,j,IVPOL) * iVar%w18 * & + exp_term_v * & + ( AZCoeff%C(2,j,IVPOL) * iVar%w18 + & + AZCoeff%C(3,j,IVPOL) * iVar%w18**2 + & + AZCoeff%C(4,j,IVPOL) * iVar%w18**3 ) ) * A1v_AD(j) + END IF A1v_AD(j) = ZERO END IF diff --git a/src/SfcOptics/MW_Water/FASTEM_MWSSEM/CRTM_Fastem1.f90 b/src/SfcOptics/MW_Water/FASTEM_MWSSEM/CRTM_Fastem1.f90 index 5bb1f4d1..18fcce55 100644 --- a/src/SfcOptics/MW_Water/FASTEM_MWSSEM/CRTM_Fastem1.f90 +++ b/src/SfcOptics/MW_Water/FASTEM_MWSSEM/CRTM_Fastem1.f90 @@ -144,6 +144,9 @@ MODULE CRTM_Fastem1 0.226701E-05_fp, -.119072E-02_fp, -.263165E-04_fp, .114597E-06_fp,& 0.406300E+00_fp, .200031E-02_fp, -.781635E-05_fp/), (/59/) ) + REAL( fp ), PARAMETER :: MIN_EXP = TINY(ONE) + REAL( fp ), PARAMETER :: MAX_EXP_ARG = -LOG(MIN_EXP) + CONTAINS @@ -322,7 +325,7 @@ SUBROUTINE Fastem1(Frequency, & ! INPUT rhorzs=rhorzsr*rhorzsr+rhorzsi*rhorzsi ! calculate small scale xcorr to reflection coefficients - xcorr1=exp(emc(21)*u10mps*pc2/freqghz2) + xcorr1=exp(min(emc(21)*u10mps*pc2/freqghz2, MAX_EXP_ARG)) ! calculate large scale geometric correction ! to calculate a correction to the fresnel reflection coefficients diff --git a/src/SfcOptics/MW_Water/FASTEM_MWSSEM/CRTM_Fastem3.f90 b/src/SfcOptics/MW_Water/FASTEM_MWSSEM/CRTM_Fastem3.f90 index d135a650..f23e0191 100644 --- a/src/SfcOptics/MW_Water/FASTEM_MWSSEM/CRTM_Fastem3.f90 +++ b/src/SfcOptics/MW_Water/FASTEM_MWSSEM/CRTM_Fastem3.f90 @@ -74,6 +74,9 @@ MODULE CRTM_Fastem3 CHARACTER(*), PRIVATE, PARAMETER :: MODULE_RCS_ID = & ! Literal constants REAL(fp), PARAMETER :: HUNDRED = 100.0_fp + REAL(fp), PARAMETER :: MIN_TRANSMITTANCE = TINY(ONE) + REAL(fp), PARAMETER :: MAX_TRANSMITTANCE = ONE - TINY(ONE) + REAL(fp), PARAMETER :: MAX_OPTICAL_DEPTH = -LOG(MIN_TRANSMITTANCE) ! ------------------------- @@ -428,6 +431,8 @@ SUBROUTINE Fastem3(Frequency, & ! INPUT REAL(fp) :: opdpsfc,freqr REAL(fp) :: zrough_v,zrough_h REAL(fp) :: zreflmod_v,zreflmod_h + REAL(fp) :: transmittance_safe, denom + REAL(fp) :: transmittance_safe, denom REAL(fp) :: delta,delta2 REAL(fp) :: qdepol,emissfactor REAL(fp) :: emissfactor_v,emissfactor_h @@ -669,7 +674,9 @@ SUBROUTINE Fastem3(Frequency, & ! INPUT If ( variance < ZERO ) variance = ZERO !Compute surface to space optical depth - opdpsfc = -log(Transmittance ) / seczen + transmittance_safe = MIN(MAX(Transmittance, MIN_TRANSMITTANCE), MAX_TRANSMITTANCE) + denom = ONE - transmittance_safe + opdpsfc = -log(transmittance_safe) / seczen !Define nine predictors for the effective angle calculation zx(1) = ONE @@ -697,10 +704,8 @@ SUBROUTINE Fastem3(Frequency, & ! INPUT & + zx(9) * c(119+jcofm1*3) ) End Do - zreflmod_v = (ONE-Transmittance ** zrough_v) & - & / (ONE-Transmittance ) - zreflmod_h = (ONE-Transmittance ** zrough_h) & - & / (ONE-Transmittance ) + zreflmod_v = (ONE-transmittance_safe ** zrough_v) / denom + zreflmod_h = (ONE-transmittance_safe ** zrough_h) / denom Reflectivity(1) = zreflmod_v * (ONE-Emissivity(1)) Reflectivity(2) = zreflmod_h * (ONE-Emissivity(2)) Reflectivity(3) = ZERO @@ -859,6 +864,7 @@ SUBROUTINE Fastem3_TL(Frequency, & ! INPUT REAL(fp) :: opdpsfc_tl, freqr_tl REAL(fp) :: zrough_v_tl, zrough_h_tl REAL(fp) :: zreflmod_v_tl, zreflmod_h_tl + REAL(fp) :: transmittance_safe, transmittance_tl_safe, denom REAL(fp) :: delta_tl, delta2_tl REAL(fp) :: qdepol_tl, emissfactor_tl REAL(fp) :: emissfactor_v_tl, emissfactor_h_tl @@ -1187,8 +1193,15 @@ SUBROUTINE Fastem3_TL(Frequency, & ! INPUT Endif !Compute surface to space optical depth - opdpsfc = -log(Transmittance) / seczen - opdpsfc_tl = -Transmittance_TL / ( Transmittance * seczen ) + transmittance_safe = MIN(MAX(Transmittance, MIN_TRANSMITTANCE), MAX_TRANSMITTANCE) + IF ( transmittance_safe /= Transmittance ) THEN + transmittance_tl_safe = ZERO + ELSE + transmittance_tl_safe = Transmittance_TL + END IF + denom = ONE - transmittance_safe + opdpsfc = -log(transmittance_safe) / seczen + opdpsfc_tl = -transmittance_tl_safe / ( transmittance_safe * seczen ) !Define nine predictors for the effective angle calculation zx(1) = ONE @@ -1245,26 +1258,26 @@ SUBROUTINE Fastem3_TL(Frequency, & ! INPUT & + zx(9) * c(119+jcofm1*3) ) End Do - zreflmod_v = (ONE-Transmittance**zrough_v)/(ONE-Transmittance) - zreflmod_h = (ONE-Transmittance**zrough_h)/(ONE-Transmittance) + zreflmod_v = (ONE-transmittance_safe**zrough_v)/denom + zreflmod_h = (ONE-transmittance_safe**zrough_h)/denom - zreflmod_v_tl = Transmittance_tl *& - & (-zrough_v * Transmittance**(zrough_v-ONE) * & - & (ONE-Transmittance)+& - & ( ONE-Transmittance**zrough_v)) & - & / (ONE-Transmittance)**2 + zreflmod_v_tl = transmittance_tl_safe * & + & (-zrough_v * transmittance_safe**(zrough_v-ONE) * & + & (ONE-transmittance_safe)+& + & ( ONE-transmittance_safe**zrough_v)) & + & / denom**2 zreflmod_v_tl = zreflmod_v_tl - & - & ( Transmittance**zrough_v * Log(Transmittance) * zrough_v_tl ) / & - & (ONE-Transmittance) + & ( transmittance_safe**zrough_v * Log(transmittance_safe) * zrough_v_tl ) / & + & denom - zreflmod_h_tl = Transmittance_tl *& - & (-zrough_h * Transmittance**(zrough_h-1.0) * (1.0-Transmittance) + & - & ( 1.0-Transmittance**zrough_h) ) & - & / (1.0-Transmittance)**2 + zreflmod_h_tl = transmittance_tl_safe * & + & (-zrough_h * transmittance_safe**(zrough_h-1.0) * (1.0-transmittance_safe) + & + & ( 1.0-transmittance_safe**zrough_h) ) & + & / denom**2 zreflmod_h_tl = zreflmod_h_tl - & - & ( Transmittance**zrough_h * Log(Transmittance) * zrough_h_tl ) / & - & (1.0-Transmittance) + & ( transmittance_safe**zrough_h * Log(transmittance_safe) * zrough_h_tl ) / & + & denom Reflectivity_tl(1) = zreflmod_v_tl * (1.0-Emissivity(1)) - zreflmod_v * Emissivity_TL(1) Reflectivity_tl(2) = zreflmod_h_tl * (1.0-Emissivity(2)) - zreflmod_h * Emissivity_TL(2) @@ -1620,7 +1633,9 @@ SUBROUTINE Fastem3_AD(Frequency, & ! INPUT Endif !Compute surface to space optical depth - opdpsfc = -log(Transmittance) / seczen + transmittance_safe = MIN(MAX(Transmittance, MIN_TRANSMITTANCE), MAX_TRANSMITTANCE) + denom = ONE - transmittance_safe + opdpsfc = -log(transmittance_safe) / seczen !Define nine predictors for the effective angle calculation zx(1) = ONE @@ -1648,8 +1663,8 @@ SUBROUTINE Fastem3_AD(Frequency, & ! INPUT & + zx(8) * c(118+jcofm1*3) & & + zx(9) * c(119+jcofm1*3) ) End Do - zreflmod_v = (ONE-Transmittance**zrough_v) / (ONE-Transmittance) - zreflmod_h = (ONE-Transmittance**zrough_h) / (ONE-Transmittance) + zreflmod_v = (ONE-transmittance_safe**zrough_v) / denom + zreflmod_h = (ONE-transmittance_safe**zrough_h) / denom End If @@ -1668,24 +1683,24 @@ SUBROUTINE Fastem3_AD(Frequency, & ! INPUT emissstokes_ad(2) = emissstokes_ad(2) - reflectstokes_ad(2) * zreflmod_h emissstokes_ad(1) = emissstokes_ad(1) - reflectstokes_ad(1) * zreflmod_v zrough_h_ad = - zreflmod_h_ad * & - & ( Transmittance**zrough_h * Log(Transmittance) ) / & - & (ONE-Transmittance) + & ( transmittance_safe**zrough_h * Log(transmittance_safe) ) / & + & denom Transmittance_AD = Transmittance_AD + zreflmod_h_ad *& - & (-zrough_h * Transmittance**(zrough_h-ONE) * & - & (ONE-Transmittance) + & - & ( ONE-TRansmittance**zrough_h) ) & - & / (ONE-Transmittance)**2 + & (-zrough_h * transmittance_safe**(zrough_h-ONE) * & + & (ONE-transmittance_safe) + & + & ( ONE-transmittance_safe**zrough_h) ) & + & / denom**2 zrough_v_ad = -zreflmod_v_ad * & - & ( Transmittance**zrough_v * Log(Transmittance) ) / & - & (ONE-Transmittance) + & ( transmittance_safe**zrough_v * Log(transmittance_safe) ) / & + & denom Transmittance_AD = Transmittance_AD + zreflmod_v_ad *& - & (-zrough_v * Transmittance**(zrough_v-ONE) * & - & (ONE-Transmittance) + & - & ( ONE-Transmittance**zrough_v) ) & - & / (ONE-Transmittance)**2 + & (-zrough_v * transmittance_safe**(zrough_v-ONE) * & + & (ONE-transmittance_safe) + & + & ( ONE-transmittance_safe**zrough_v) ) & + & / denom**2 zx_ad(:) = ZERO Do jcof = 1,7 diff --git a/src/SfcOptics/MW_Water/LowFrequency_MWSSEM/CRTM_OceanEM_AMSRE.f90 b/src/SfcOptics/MW_Water/LowFrequency_MWSSEM/CRTM_OceanEM_AMSRE.f90 index 93c977df..1b95b163 100644 --- a/src/SfcOptics/MW_Water/LowFrequency_MWSSEM/CRTM_OceanEM_AMSRE.f90 +++ b/src/SfcOptics/MW_Water/LowFrequency_MWSSEM/CRTM_OceanEM_AMSRE.f90 @@ -62,6 +62,8 @@ MODULE CRTM_OCEANEM_AMSRE REAL( fp ), PARAMETER :: Salinity_Defaults = 35.0_fp REAL( fp ), PARAMETER :: GHz2Hz = 1.0E09_fp REAL( fp ), PARAMETER :: Kelvin2Celsius = 273.15_fp + REAL( fp ), PARAMETER :: MIN_EXP = TINY(ONE) + REAL( fp ), PARAMETER :: MAX_EXP_ARG = -LOG(MIN_EXP) REAL( fp ), SAVE :: MODEL_LIMIT(3) ! 1: Permittivity_limit(GHz), Foam_limit(m/s), Small scale limit(GHz) REAL( fp ), SAVE :: sdd(40,100) REAL( fp ), SAVE :: Coeff(8) @@ -1006,6 +1008,7 @@ SUBROUTINE OCEANEM_AMSRE(Frequency, & ! REAL(fp) frequency_sdd(100),wind_speed_sdd(40) REAL(fp) freq,x1,x2,y1,y2,sdd_intrp,satellite_zenith,salinity + REAL(fp) exp_arg, exp_term REAL(fp) csl,csl2 COMPLEX( fp ) :: eps_ocean,Rvv,Rhh @@ -1081,8 +1084,10 @@ SUBROUTINE OCEANEM_AMSRE(Frequency, & ! IF ( freq > MODEL_LIMIT(3) ) THEN - Rv_Small = exp( - sdd_intrp * csl2) - Rh_Small = exp( - sdd_intrp * csl2) + exp_arg = sdd_intrp * csl2 + exp_term = exp(-min(exp_arg, MAX_EXP_ARG)) + Rv_Small = exp_term + Rh_Small = exp_term ELSE Rv_Small = ONE Rh_Small = ONE @@ -1128,6 +1133,8 @@ SUBROUTINE OCEANEM_AMSRE_TL(Frequency, & REAL(fp) frequency_sdd(100),wind_speed_sdd(40) REAL(fp) freq,satellite_zenith,salinity + REAL(fp) exp_arg, exp_term + REAL(fp) exp_arg, exp_term, exp_term_tl REAL(fp) x1,x2,y1,y2,sdd_intrp REAL(fp) x1_tl,x2_tl,sdd_intrp_tl @@ -1234,10 +1241,17 @@ SUBROUTINE OCEANEM_AMSRE_TL(Frequency, & END IF IF ( freq > MODEL_LIMIT(3) ) THEN - Rv_Small = exp( - sdd_intrp * csl2) - Rh_Small = exp( - sdd_intrp * csl2) - Rv_Small_tl = -csl2 * exp( - sdd_intrp * csl2) * sdd_intrp_tl - Rh_Small_tl = -csl2 * exp( - sdd_intrp * csl2) * sdd_intrp_tl + exp_arg = sdd_intrp * csl2 + exp_term = exp(-min(exp_arg, MAX_EXP_ARG)) + Rv_Small = exp_term + Rh_Small = exp_term + IF ( exp_arg >= MAX_EXP_ARG ) THEN + exp_term_tl = ZERO + ELSE + exp_term_tl = -csl2 * exp_term * sdd_intrp_tl + END IF + Rv_Small_tl = exp_term_tl + Rh_Small_tl = exp_term_tl ELSE Rv_Small = ONE Rh_Small = ONE @@ -1414,8 +1428,10 @@ SUBROUTINE OCEANEM_AMSRE_AD(Frequency, & IF ( freq > MODEL_LIMIT(3) ) THEN - Rv_Small = exp( - sdd_intrp * csl2) - Rh_Small = exp( - sdd_intrp * csl2) + exp_arg = sdd_intrp * csl2 + exp_term = exp(-min(exp_arg, MAX_EXP_ARG)) + Rv_Small = exp_term + Rh_Small = exp_term ELSE Rv_Small = ONE Rh_Small = ONE @@ -1464,8 +1480,12 @@ SUBROUTINE OCEANEM_AMSRE_AD(Frequency, & IF ( freq > MODEL_LIMIT(3) ) THEN - sdd_intrp_ad = - csl2 * exp( - sdd_intrp * csl2) * Rh_Small_ad - sdd_intrp_ad = sdd_intrp_ad -csl2 * exp( - sdd_intrp * csl2) * Rv_Small_ad + IF ( exp_arg >= MAX_EXP_ARG ) THEN + sdd_intrp_ad = ZERO + ELSE + sdd_intrp_ad = - csl2 * exp_term * Rh_Small_ad + sdd_intrp_ad = sdd_intrp_ad -csl2 * exp_term * Rv_Small_ad + END IF ELSE sdd_intrp_ad = ZERO END IF @@ -2200,4 +2220,3 @@ END FUNCTION Compute_Fh END MODULE CRTM_OCEANEM_AMSRE - diff --git a/src/SfcOptics/MW_Water/LowFrequency_MWSSEM/RSS_Emissivity_Model.f90 b/src/SfcOptics/MW_Water/LowFrequency_MWSSEM/RSS_Emissivity_Model.f90 index 99f3aef1..f4ea9fd1 100644 --- a/src/SfcOptics/MW_Water/LowFrequency_MWSSEM/RSS_Emissivity_Model.f90 +++ b/src/SfcOptics/MW_Water/LowFrequency_MWSSEM/RSS_Emissivity_Model.f90 @@ -36,6 +36,9 @@ module RSS_Emissivity_Model ! Declare local constant Pi REAL, PARAMETER :: Pi = 3.14159265358979323846 REAL(fp), parameter :: f0=17.97510 +REAL(fp), PARAMETER :: MIN_EXP = TINY(1.0_fp) +REAL(fp), PARAMETER :: MIN_EXP_ARG = LOG(MIN_EXP) +REAL(fp), PARAMETER :: MAX_EXP_ARG = -MIN_EXP_ARG public :: fdem0_meissner_wentz, & FDEM0_MEISSNER_WENTZ_TL, & @@ -44,6 +47,24 @@ module RSS_Emissivity_Model contains + +subroutine Safe_Exp(arg, exp_val, is_clamped) + real(fp), intent(in) :: arg + real(fp), intent(out) :: exp_val + logical, intent(out) :: is_clamped + + if (arg > MAX_EXP_ARG) then + exp_val = exp(MAX_EXP_ARG) + is_clamped = .true. + elseif (arg < MIN_EXP_ARG) then + exp_val = exp(MIN_EXP_ARG) + is_clamped = .true. + else + exp_val = exp(arg) + is_clamped = .false. + end if +end subroutine Safe_Exp + subroutine fdem0_meissner_wentz(freq,tht,sst,salinity, em0) ! input: ! name parameter unit range @@ -134,6 +155,8 @@ subroutine dielectric_meissner_wentz(sst_in,s, e0s,e1s,e2s,n1s,n2s,sig) real(fp), dimension(5), parameter :: b1coef=(/0.23232E-02, -0.79208E-04, 0.36764E-05, -0.35594E-06, 0.89795E-08/) real(fp) :: e0,e1,e2,n1,n2 real(fp) :: a0,a1,a2,b1,b2 +real(fp) :: arg0,arg1 +logical :: clamp0, clamp1 real(fp) :: sig35,r15,rtr15,alpha0,alpha1 real(fp) :: sst,sst2,sst3,sst4,s2 @@ -147,9 +170,9 @@ subroutine dielectric_meissner_wentz(sst_in,s, e0s,e1s,e2s,n1s,n2s,sig) ! pure water e0 = (3.70886e4 - 8.2168e1*sst)/(4.21854e2 + sst) ! stogryn et al. - e1 = x(1) + x(2)*sst + x(3)*sst2 + e1 = x(1) + x(2)*sst + x(3)*sst2 n1 = (45.00 + sst)/(x(4) + x(5)*sst + x(6)*sst2) - e2 = x(7) + x(8)*sst + e2 = x(7) + x(8)*sst n2 = (45.00 + sst)/(x(9) + x(10)*sst + x(11)*sst2) ! saline water @@ -164,7 +187,8 @@ subroutine dielectric_meissner_wentz(sst_in,s, e0s,e1s,e2s,n1s,n2s,sig) sig = sig35*r15*rtr15 ! permittivity - a0 = exp(a0coef(1)*s + a0coef(2)*s2 + a0coef(3)*s*sst) + arg0 = a0coef(1)*s + a0coef(2)*s2 + a0coef(3)*s*sst + call Safe_Exp(arg0, a0, clamp0) e0s = a0*e0 if(sst.le.30) then @@ -174,7 +198,8 @@ subroutine dielectric_meissner_wentz(sst_in,s, e0s,e1s,e2s,n1s,n2s,sig) endif n1s = n1*b1 - a1 = exp(z(7)*s + z(8)*s2 + z(9)*s*sst) + arg1 = z(7)*s + z(8)*s2 + z(9)*s*sst + call Safe_Exp(arg1, a1, clamp1) e1s = e1*a1 ! b2 = 1.0 + s*(z(10) + z(11)*sst) @@ -359,9 +384,12 @@ SUBROUTINE DIELECTRIC_MEISSNER_WENTZ_TL(sst_in, sst_in_TL, salinity, salinity_TL sig = sig35*r15*rtr15 ! permittivity temp0 = a0coef(1)*s + a0coef(2)*s2 + a0coef(3)*s*sst - a0_TL = EXP(temp0)*(a0coef(1)*s_TL+a0coef(2)*s2_TL+a0coef(3)*(sst*s_TL+s*sst_TL)& -& ) - a0 = EXP(temp0) + call Safe_Exp(temp0, a0, clamp0) + if (clamp0) then + a0_TL = 0.0_fp + else + a0_TL = a0*(a0coef(1)*s_TL+a0coef(2)*s2_TL+a0coef(3)*(sst*s_TL+s*sst_TL)) + end if e0s_TL = e0*a0_TL + a0*e0_TL e0s = a0*e0 IF (sst .LE. 30) THEN @@ -378,8 +406,12 @@ SUBROUTINE DIELECTRIC_MEISSNER_WENTZ_TL(sst_in, sst_in_TL, salinity, salinity_TL n1s_TL = b1*n1_TL + n1*b1_TL n1s = n1*b1 temp0 = z(7)*s + z(8)*s2 + z(9)*s*sst - a1_TL = EXP(temp0)*(z(7)*s_TL+z(8)*s2_TL+z(9)*(sst*s_TL+s*sst_TL)) - a1 = EXP(temp0) + call Safe_Exp(temp0, a1, clamp1) + if (clamp1) then + a1_TL = 0.0_fp + else + a1_TL = a1*(z(7)*s_TL+z(8)*s2_TL+z(9)*(sst*s_TL+s*sst_TL)) + end if e1s_TL = a1*e1_TL + e1*a1_TL e1s = e1*a1 ! b2 = 1.0 + s*(z(10) + z(11)*sst) @@ -504,6 +536,8 @@ SUBROUTINE DIELECTRIC_MEISSNER_WENTZ_K(sst_in, salinity, d_e0s_dsst, d_e0s_ds & REAL(FP) :: d_n2_dsst, d_sig35_dsst,d_r15_ds, d_alpha0_ds REAL(FP) :: d_alpha1_ds,d_rtr15_dsst,d_rtr15_ds REAL(FP) :: d_a0_dsst,d_a0_ds,d_a1_dsst,d_a1_ds,d_a2_dsst + REAL(FP) :: arg0, arg1 + LOGICAL :: clamp0, clamp1 REAL(FP) :: d_a2_ds,d_b1_dsst,d_b1_ds,d_b2_dsst,d_b2_ds sst = sst_in s = salinity @@ -548,9 +582,15 @@ SUBROUTINE DIELECTRIC_MEISSNER_WENTZ_K(sst_in, salinity, d_e0s_dsst, d_e0s_ds & d_sig_ds = sig35 *(d_r15_ds *rtr15 + r15*d_rtr15_ds) ! permittivity - a0 = exp(a0coef(1)*s + a0coef(2)*s**2 + a0coef(3)*s*sst) - d_a0_dsst =(a0coef(3)*s) * exp(a0coef(1)*s + a0coef(2)*s**2 + a0coef(3)*s*sst) - d_a0_ds =(2*a0coef(2)*s+a0coef(1)+a0coef(3)*sst)* exp(s*(a0coef(2)*s+a0coef(1)+a0coef(3)*sst)) + arg0 = a0coef(1)*s + a0coef(2)*s**2 + a0coef(3)*s*sst + call Safe_Exp(arg0, a0, clamp0) + if (clamp0) then + d_a0_dsst = 0.0_fp + d_a0_ds = 0.0_fp + else + d_a0_dsst =(a0coef(3)*s) * a0 + d_a0_ds =(2*a0coef(2)*s+a0coef(1)+a0coef(3)*sst) * a0 + end if !e0s = a0*e0 d_e0s_dsst= d_a0_dsst * e0 +a0*d_e0_dsst @@ -569,9 +609,15 @@ SUBROUTINE DIELECTRIC_MEISSNER_WENTZ_K(sst_in, salinity, d_e0s_dsst, d_e0s_ds & d_n1s_dsst= d_n1_dsst *b1 + n1 * d_b1_dsst d_n1s_ds = n1 * d_b1_ds - a1 = exp(z(7)*s + z(8)*s**2 + z(9)*s*sst) - d_a1_dsst = exp(z(7)*s + z(8)*s**2) *z(9)*s*exp(z(9)*s*sst) - d_a1_ds = (2*z(8)*s+z(7)+z(9)*sst) * exp(s*(z(8)*s+z(7)+z(9)*sst)) + arg1 = z(7)*s + z(8)*s**2 + z(9)*s*sst + call Safe_Exp(arg1, a1, clamp1) + if (clamp1) then + d_a1_dsst = 0.0_fp + d_a1_ds = 0.0_fp + else + d_a1_dsst = (z(9)*s) * a1 + d_a1_ds = (2*z(8)*s+z(7)+z(9)*sst) * a1 + end if !e1s = e1*a1 d_e1s_dsst = d_e1_dsst*a1 +e1* d_a1_dsst diff --git a/src/SfcOptics/NESDIS_Emissivity/NESDIS_ATMS_SnowEM_Module.f90 b/src/SfcOptics/NESDIS_Emissivity/NESDIS_ATMS_SnowEM_Module.f90 index 36166174..543b128b 100644 --- a/src/SfcOptics/NESDIS_Emissivity/NESDIS_ATMS_SnowEM_Module.f90 +++ b/src/SfcOptics/NESDIS_Emissivity/NESDIS_ATMS_SnowEM_Module.f90 @@ -91,6 +91,7 @@ MODULE NESDIS_ATMS_SnowEM_Module PRIVATE PUBLIC :: NESDIS_ATMS_SNOWEM + REAL(fp), PARAMETER :: MIN_LOG_ARG = TINY(1.0_fp) CONTAINS @@ -543,7 +544,7 @@ SUBROUTINE ATMS_SNOW_ByTBTs_D(frequency,tb,ts,snow_type,em_vector) !*** adjustment from the library values emw=em(windex,snow_type) X=1.0_fp/emw - Y((/1,2,4,5/)) = log(Tb((/1,2,4,5/))/(Ts*emw((/1,2,4,5/)))) + Y((/1,2,4,5/)) = log(MAX(Tb((/1,2,4,5/))/(Ts*emw((/1,2,4,5/))), MIN_LOG_ARG)) IF(frequency > 100.0_fp) THEN XX=DOT_PRODUCT(X((/1,2,4,5/)),X((/1,2,4,5/))) XY=DOT_PRODUCT(X((/1,2,4,5/)),Y((/1,2,4,5/))) diff --git a/src/SfcOptics/NESDIS_Emissivity/NESDIS_OCEANEM_Module.f90 b/src/SfcOptics/NESDIS_Emissivity/NESDIS_OCEANEM_Module.f90 index a58b5165..b2a0a1a2 100644 --- a/src/SfcOptics/NESDIS_Emissivity/NESDIS_OCEANEM_Module.f90 +++ b/src/SfcOptics/NESDIS_Emissivity/NESDIS_OCEANEM_Module.f90 @@ -33,6 +33,8 @@ MODULE NESDIS_OCEANEM_Module REAL(fp), PARAMETER :: SST_max = 330.0_fp REAL(fp), PARAMETER :: wind_min = 0.0_fp REAL(fp), PARAMETER :: wind_max = 100.0_fp + REAL(fp), PARAMETER :: MIN_EXP = TINY(one) + REAL(fp), PARAMETER :: MAX_EXP_ARG = -LOG(MIN_EXP) CONTAINS @@ -205,6 +207,7 @@ subroutine NESDIS_OCeanEM(Frequency, & ! I real(fp) :: Emissivity_H,Emissivity_V,EH_dSST, EH_dSSW, EV_dSST, EV_dSSW real(fp) :: f,foam,g,tr,rfoam,ref,rclear + real(fp) :: exp_arg, exp_term complex mu, eps, aid1,aid2,aid3,cang,rh,rv @@ -247,7 +250,9 @@ subroutine NESDIS_OCeanEM(Frequency, & ! I else - foam=0.006_fp*(1.0_fp-exp(-f*1.0e-9/7.5_fp))*(wind-7.0_fp) + exp_arg = f*1.0e-9/7.5_fp + exp_term = exp(-min(exp_arg, MAX_EXP_ARG)) + foam=0.006_fp*(1.0_fp-exp_term)*(wind-7.0_fp) endif @@ -327,10 +332,10 @@ real function EPSP (t1,s,f) end function EPSP - real function EPSPP (t1,s,f) +real function EPSPP (t1,s,f) - real(fp) s,f,t1,t,t2,eswi,eo,eswo,a,b,d,esw,tswo,tsw,sswo,fi,ssw + real(fp) s,f,t1,t,t2,eswi,eo,eswo,a,b,d,esw,tswo,tsw,sswo,fi,ssw,exp_term t=t1-273.0_fp @@ -359,7 +364,8 @@ real function EPSPP (t1,s,f) fi = d*(2.033e-2+1.266e-4*d+2.464e-6*d**2- s*(1.849e-5-2.551e-7*d+2.551e-8*d*d)) - ssw = sswo*exp(-fi) + exp_term = exp(-min(fi, MAX_EXP_ARG)) + ssw = sswo*exp_term EPSPP = tsw*f*(esw-eswi)/(1.0_fp+(tsw*f)**2) @@ -425,9 +431,11 @@ subroutine OceanEM_TL_SSTW(degre,frequency,sst,wind,Salinity,deh_dt,deh_dw,dev_d dfoam_dw=zero else - foam=0.006_fp*(one-exp(-f*1.0e-9/7.5_fp))*(wind-7.0_fp) + exp_arg = f*1.0e-9/7.5_fp + exp_term = exp(-min(exp_arg, MAX_EXP_ARG)) + foam=0.006_fp*(one-exp_term)*(wind-7.0_fp) - dfoam_dw=0.006_fp*(one-exp(-f*1.0e-9/7.5_fp)) + dfoam_dw=0.006_fp*(one-exp_term) endif @@ -555,7 +563,7 @@ end function depsp_dt real function depspp_dt (t1,s,f) - real(fp) :: s,f,t1,t,t2,eswi,eo,eswo,a,b,d,esw,tswo,tsw,sswo,fi,ssw,epspp + real(fp) :: s,f,t1,t,t2,eswi,eo,eswo,a,b,d,esw,tswo,tsw,sswo,fi,ssw,epspp,exp_term real(fp) :: deswo_dt,da_dt,db_dt,desw_dt,dtswo_dt,dtsw_dt,dfi_dt,dssw_dt,aid @@ -602,9 +610,10 @@ real function depspp_dt (t1,s,f) - d*(1.266e-4+2.464e-6*d- s*(-2.551e-7+2.0*2.551e-8*d)) - ssw = sswo*exp(-fi) + exp_term = exp(-min(fi, MAX_EXP_ARG)) + ssw = sswo*exp_term - dssw_dt = - sswo*exp(-fi)*dfi_dt + dssw_dt = - sswo*exp_term*dfi_dt epspp = tsw*f*(esw-eswi)/(1.0_fp+(tsw*f)**2) diff --git a/src/SfcOptics/NESDIS_Emissivity/NESDIS_SEAICE_PHYEM_MODULE.f90 b/src/SfcOptics/NESDIS_Emissivity/NESDIS_SEAICE_PHYEM_MODULE.f90 index d937a122..6bd688f8 100644 --- a/src/SfcOptics/NESDIS_Emissivity/NESDIS_SEAICE_PHYEM_MODULE.f90 +++ b/src/SfcOptics/NESDIS_Emissivity/NESDIS_SEAICE_PHYEM_MODULE.f90 @@ -21,6 +21,9 @@ MODULE NESDIS_SEAICE_PHYEM_MODULE ! ----------------- ! Module parameters ! ----------------- + REAL(fp), PARAMETER :: one = 1.0_fp + REAL(fp), PARAMETER :: MIN_EXP = TINY(one) + REAL(fp), PARAMETER :: MAX_EXP_ARG = -LOG(MIN_EXP) CONTAINS @@ -178,7 +181,7 @@ subroutine permitivity(itype, temp, s, f, e_real, e_img) real(fp) :: ESWI, EW0, TPTW, FH, X, Y, TPTWFQ real(fp) :: e_real, e_img - real(fp) :: s,t4,t5,sigma + real(fp) :: s,t4,t5,sigma,exp_term real(fp) :: ebo,eb,sb,nb,eb_real,eb_imag real(fp) :: den_ice,den_sice,den_brine @@ -223,7 +226,8 @@ subroutine permitivity(itype, temp, s, f, e_real, e_img) fi = d * (2.033e-2 + 1.266e-4 * d + 2.464e-6 * d * d - & s * (1.849e-5 - 2.551e-7 * d + 2.551e-8 * d * d)) - ssw = sswo * exp(-fi) + exp_term = exp(-min(fi, MAX_EXP_ARG)) + ssw = sswo * exp_term epspp = tsw * f * (esw - ESWI) / (1.0 + (tsw * f)**2) epspp = epspp + ssw/(2.0 * pi * eo * f) diff --git a/src/Source_Functions/CRTM_Planck_Functions.f90 b/src/Source_Functions/CRTM_Planck_Functions.f90 index 43395740..b963aee4 100644 --- a/src/Source_Functions/CRTM_Planck_Functions.f90 +++ b/src/Source_Functions/CRTM_Planck_Functions.f90 @@ -16,7 +16,7 @@ MODULE CRTM_Planck_Functions ! ----------------- ! Module use statements USE Type_Kinds , ONLY: fp - USE CRTM_Parameters, ONLY: ONE + USE CRTM_Parameters, ONLY: ONE, ZERO USE CRTM_SpcCoeff , ONLY: SC ! Disable all implicit typing IMPLICIT NONE @@ -37,6 +37,7 @@ MODULE CRTM_Planck_Functions ! ---------- ! Parameters ! ---------- + REAL(fp), PARAMETER :: MIN_RADIANCE = TINY(ONE) CONTAINS @@ -477,11 +478,13 @@ SUBROUTINE CRTM_Planck_Temperature( n, l , & ! Input REAL(fp), INTENT(OUT) :: Temperature ! Local variables REAL(fp) :: Effective_Temperature + REAL(fp) :: Radiance_Safe ! Calculate the effective temperature + Radiance_Safe = MAX(Radiance, MIN_RADIANCE) Effective_Temperature = SC(n)%Planck_C2(l) / & ! ---------------------------------------------- - LOG( ( SC(n)%Planck_C1(l) / Radiance ) + ONE ) + LOG( ( SC(n)%Planck_C1(l) / Radiance_Safe ) + ONE ) ! Apply the polychromatic correction to ! obtain the true temperature @@ -588,15 +591,21 @@ SUBROUTINE CRTM_Planck_Temperature_TL( n, l , & ! Input ! Local variables REAL(fp) :: Argument REAL(fp) :: F + REAL(fp) :: Radiance_Safe ! Calculate the Planck function operator ! ! The logarithm argument - Argument = ( SC(n)%Planck_C1(l) / Radiance ) + ONE + IF ( Radiance <= MIN_RADIANCE ) THEN + Temperature_TL = ZERO + RETURN + END IF + Radiance_Safe = Radiance + Argument = ( SC(n)%Planck_C1(l) / Radiance_Safe ) + ONE ! The operator, call it F F = SC(n)%Planck_C1(l) * SC(n)%Planck_C2(l) / & ! ------------------------------------------------------------------- - ( SC(n)%Band_C2(l) * Argument * ( Radiance * LOG( Argument ) )**2 ) + ( SC(n)%Band_C2(l) * Argument * ( Radiance_Safe * LOG( Argument ) )**2 ) ! Calculate the tangent-linear temperature Temperature_TL = F * Radiance_TL @@ -707,15 +716,18 @@ SUBROUTINE CRTM_Planck_Temperature_AD( n, l , & ! Input ! Local variables REAL(fp) :: Argument REAL(fp) :: F + REAL(fp) :: Radiance_Safe ! Calculate the Planck function operator ! ! The logarithm Argument - Argument = ( SC(n)%Planck_C1(l) / Radiance ) + ONE + IF ( Radiance <= MIN_RADIANCE ) RETURN + Radiance_Safe = Radiance + Argument = ( SC(n)%Planck_C1(l) / Radiance_Safe ) + ONE ! The operator, call it F F = SC(n)%Planck_C1(l) * SC(n)%Planck_C2(l) / & ! ------------------------------------------------------------------- - ( SC(n)%Band_C2(l) * Argument * ( Radiance * LOG( Argument ) )**2 ) + ( SC(n)%Band_C2(l) * Argument * ( Radiance_Safe * LOG( Argument ) )**2 ) ! Calculate the adjoint radiance Radiance_AD = Radiance_AD + ( F * Temperature_AD ) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 692dad4c..c50f1642 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -281,12 +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) diff --git a/test/mains/regression/adjoint/test_ClearSky/Compare_Diagnostics.inc b/test/mains/regression/adjoint/test_ClearSky/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/adjoint/test_ClearSky/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/adjoint/test_ClearSky/test_ClearSky.f90 b/test/mains/regression/adjoint/test_ClearSky/test_ClearSky.f90 index 65e3b7f7..35581f17 100644 --- a/test/mains/regression/adjoint/test_ClearSky/test_ClearSky.f90 +++ b/test/mains/regression/adjoint/test_ClearSky/test_ClearSky.f90 @@ -391,6 +391,7 @@ PROGRAM test_ClearSky ELSE Message = 'Atmosphere_AD Adjoints are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats_R1( PROGRAM_NAME, 'Atmosphere_AD', Atmosphere_AD, atm_AD ) ! Write the current Atmosphere_AD results to file atmad_file = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmad_file, Atmosphere_AD, Quiet=.TRUE. ) @@ -407,6 +408,7 @@ PROGRAM test_ClearSky ELSE Message = 'Surface_AD Adjoints are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats_R1( PROGRAM_NAME, 'Surface_AD', Surface_AD, sfc_AD ) ! Write the current Surface_AD results to file sfcad_file = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfcad_file, Surface_AD, Quiet=.TRUE. ) @@ -423,6 +425,7 @@ PROGRAM test_ClearSky ELSE Message = 'RTSolution_AD results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_AD', RTSolution_AD, rts_AD ) rtsad_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_AD.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsad_File, RTSolution_AD, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -466,5 +469,6 @@ PROGRAM test_ClearSky INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_ClearSky diff --git a/test/mains/regression/adjoint/test_Simple/Compare_Diagnostics.inc b/test/mains/regression/adjoint/test_Simple/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/adjoint/test_Simple/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/adjoint/test_Simple/test_Simple.f90 b/test/mains/regression/adjoint/test_Simple/test_Simple.f90 index 4cd13f8e..a61f76af 100644 --- a/test/mains/regression/adjoint/test_Simple/test_Simple.f90 +++ b/test/mains/regression/adjoint/test_Simple/test_Simple.f90 @@ -464,6 +464,7 @@ PROGRAM test_Simple ELSE Message = 'Atmosphere_AD Adjoints are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats_R1( PROGRAM_NAME, 'Atmosphere_AD', Atmosphere_AD, atm_AD ) ! Write the current Atmosphere_AD results to file atmad_file = TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmad_file, atm_AD, Quiet=.TRUE. ) @@ -480,6 +481,7 @@ PROGRAM test_Simple ELSE Message = 'Surface_AD Adjoints are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats_R1( PROGRAM_NAME, 'Surface_AD', Surface_AD, sfc_AD ) ! Write the current Surface_AD results to file sfcad_file = TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfcad_file, Surface_AD, Quiet=.TRUE. ) @@ -496,6 +498,7 @@ PROGRAM test_Simple ELSE Message = 'RTSolution_AD results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_AD', RTSolution_AD, rts_AD ) rtsad_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_AD.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsad_File, RTSolution_AD, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -544,5 +547,6 @@ PROGRAM test_Simple INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_Simple diff --git a/test/mains/regression/forward/test_AOD/Compare_Diagnostics.inc b/test/mains/regression/forward/test_AOD/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_AOD/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_AOD/test_AOD.f90 b/test/mains/regression/forward/test_AOD/test_AOD.f90 index 950a3e22..a3cf2009 100644 --- a/test/mains/regression/forward/test_AOD/test_AOD.f90 +++ b/test/mains/regression/forward/test_AOD/test_AOD.f90 @@ -247,6 +247,7 @@ PROGRAM test_AOD ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -286,5 +287,6 @@ PROGRAM test_AOD CONTAINS INCLUDE 'Load_Atm_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_AOD diff --git a/test/mains/regression/forward/test_Aircraft/Compare_Diagnostics.inc b/test/mains/regression/forward/test_Aircraft/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_Aircraft/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_Aircraft/test_Aircraft.f90 b/test/mains/regression/forward/test_Aircraft/test_Aircraft.f90 index af70cc98..ed2aab8a 100644 --- a/test/mains/regression/forward/test_Aircraft/test_Aircraft.f90 +++ b/test/mains/regression/forward/test_Aircraft/test_Aircraft.f90 @@ -273,6 +273,7 @@ PROGRAM test_Aircraft ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -314,5 +315,6 @@ PROGRAM test_Aircraft INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_Aircraft diff --git a/test/mains/regression/forward/test_ChannelSubset/Compare_Diagnostics.inc b/test/mains/regression/forward/test_ChannelSubset/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_ChannelSubset/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_ChannelSubset/test_ChannelSubset.f90 b/test/mains/regression/forward/test_ChannelSubset/test_ChannelSubset.f90 index a1cadbb9..0d8e9d16 100644 --- a/test/mains/regression/forward/test_ChannelSubset/test_ChannelSubset.f90 +++ b/test/mains/regression/forward/test_ChannelSubset/test_ChannelSubset.f90 @@ -285,6 +285,7 @@ PROGRAM test_ChannelSubset ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -327,5 +328,6 @@ PROGRAM test_ChannelSubset INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_ChannelSubset diff --git a/test/mains/regression/forward/test_ClearSky/Compare_Diagnostics.inc b/test/mains/regression/forward/test_ClearSky/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_ClearSky/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_ClearSky/test_ClearSky.f90 b/test/mains/regression/forward/test_ClearSky/test_ClearSky.f90 index 6e93b767..9fadf149 100644 --- a/test/mains/regression/forward/test_ClearSky/test_ClearSky.f90 +++ b/test/mains/regression/forward/test_ClearSky/test_ClearSky.f90 @@ -268,6 +268,7 @@ PROGRAM test_ClearSky ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -309,5 +310,6 @@ PROGRAM test_ClearSky INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_ClearSky diff --git a/test/mains/regression/forward/test_Downwelling_Radiance/Compare_Diagnostics.inc b/test/mains/regression/forward/test_Downwelling_Radiance/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_Downwelling_Radiance/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_Downwelling_Radiance/test_Downwelling_Radiance.f90 b/test/mains/regression/forward/test_Downwelling_Radiance/test_Downwelling_Radiance.f90 index e8846d22..e0720e8d 100644 --- a/test/mains/regression/forward/test_Downwelling_Radiance/test_Downwelling_Radiance.f90 +++ b/test/mains/regression/forward/test_Downwelling_Radiance/test_Downwelling_Radiance.f90 @@ -273,6 +273,7 @@ PROGRAM test_Downwelling_Radiance ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -315,5 +316,6 @@ PROGRAM test_Downwelling_Radiance INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_Downwelling_Radiance diff --git a/test/mains/regression/forward/test_OMPoverChannels/Compare_Diagnostics.inc b/test/mains/regression/forward/test_OMPoverChannels/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_OMPoverChannels/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_OMPoverChannels/test_OMPoverChannels.F90 b/test/mains/regression/forward/test_OMPoverChannels/test_OMPoverChannels.F90 index 2f451de6..a07e6a20 100644 --- a/test/mains/regression/forward/test_OMPoverChannels/test_OMPoverChannels.F90 +++ b/test/mains/regression/forward/test_OMPoverChannels/test_OMPoverChannels.F90 @@ -284,6 +284,7 @@ PROGRAM test_OMPoverChannels ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -324,5 +325,6 @@ PROGRAM test_OMPoverChannels INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_OMPoverChannels diff --git a/test/mains/regression/forward/test_SOI/Compare_Diagnostics.inc b/test/mains/regression/forward/test_SOI/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_SOI/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_SOI/test_SOI.f90 b/test/mains/regression/forward/test_SOI/test_SOI.f90 index fc451149..de70a0cb 100644 --- a/test/mains/regression/forward/test_SOI/test_SOI.f90 +++ b/test/mains/regression/forward/test_SOI/test_SOI.f90 @@ -270,6 +270,7 @@ PROGRAM test_SOI ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -311,5 +312,6 @@ PROGRAM test_SOI INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_SOI diff --git a/test/mains/regression/forward/test_SSU/Compare_Diagnostics.inc b/test/mains/regression/forward/test_SSU/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_SSU/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_SSU/test_SSU.f90 b/test/mains/regression/forward/test_SSU/test_SSU.f90 index 2e97d18b..d78db912 100644 --- a/test/mains/regression/forward/test_SSU/test_SSU.f90 +++ b/test/mains/regression/forward/test_SSU/test_SSU.f90 @@ -281,6 +281,7 @@ PROGRAM test_SSU ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -322,5 +323,6 @@ PROGRAM test_SSU INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_SSU diff --git a/test/mains/regression/forward/test_ScatteringSwitch/Compare_Diagnostics.inc b/test/mains/regression/forward/test_ScatteringSwitch/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_ScatteringSwitch/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_ScatteringSwitch/test_ScatteringSwitch.f90 b/test/mains/regression/forward/test_ScatteringSwitch/test_ScatteringSwitch.f90 index c65c3605..f3c87d59 100644 --- a/test/mains/regression/forward/test_ScatteringSwitch/test_ScatteringSwitch.f90 +++ b/test/mains/regression/forward/test_ScatteringSwitch/test_ScatteringSwitch.f90 @@ -274,6 +274,7 @@ PROGRAM test_ScatteringSwitch ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -315,5 +316,6 @@ PROGRAM test_ScatteringSwitch INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_ScatteringSwitch diff --git a/test/mains/regression/forward/test_Simple/Compare_Diagnostics.inc b/test/mains/regression/forward/test_Simple/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_Simple/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_Simple/test_Simple.f90 b/test/mains/regression/forward/test_Simple/test_Simple.f90 index c003297e..8815dc59 100644 --- a/test/mains/regression/forward/test_Simple/test_Simple.f90 +++ b/test/mains/regression/forward/test_Simple/test_Simple.f90 @@ -283,6 +283,7 @@ PROGRAM test_Simple ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.nc' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, NetCDF=.TRUE., Quiet=.TRUE. ) @@ -326,5 +327,6 @@ PROGRAM test_Simple INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_Simple diff --git a/test/mains/regression/forward/test_User_Emissivity/Compare_Diagnostics.inc b/test/mains/regression/forward/test_User_Emissivity/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_User_Emissivity/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_User_Emissivity/test_User_Emissivity.f90 b/test/mains/regression/forward/test_User_Emissivity/test_User_Emissivity.f90 index 588ca7a1..f3a7632f 100644 --- a/test/mains/regression/forward/test_User_Emissivity/test_User_Emissivity.f90 +++ b/test/mains/regression/forward/test_User_Emissivity/test_User_Emissivity.f90 @@ -283,6 +283,7 @@ PROGRAM test_User_Emissivity ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -324,5 +325,6 @@ PROGRAM test_User_Emissivity INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_User_Emissivity diff --git a/test/mains/regression/forward/test_VerticalCoordinates/Compare_Diagnostics.inc b/test/mains/regression/forward/test_VerticalCoordinates/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_VerticalCoordinates/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_VerticalCoordinates/test_VerticalCoordinates.f90 b/test/mains/regression/forward/test_VerticalCoordinates/test_VerticalCoordinates.f90 index 135de42c..da4c2121 100644 --- a/test/mains/regression/forward/test_VerticalCoordinates/test_VerticalCoordinates.f90 +++ b/test/mains/regression/forward/test_VerticalCoordinates/test_VerticalCoordinates.f90 @@ -390,6 +390,7 @@ PROGRAM test_VerticalCoordinates ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -406,6 +407,7 @@ PROGRAM test_VerticalCoordinates ELSE Message = 'regional RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_NAM', RTSolution_NAM, rts_NAM ) ! Write the current RTSolution results to file rts_NAM_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.NAM.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_NAM_File, RTSolution_NAM, Quiet=.TRUE. ) @@ -422,6 +424,7 @@ PROGRAM test_VerticalCoordinates ELSE Message = 'global RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_GFS', RTSolution_GFS, rts_GFS ) ! Write the current RTSolution results to file rts_GFS_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.GFS.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_GFS_File, RTSolution_GFS, Quiet=.TRUE. ) @@ -469,5 +472,6 @@ PROGRAM test_VerticalCoordinates INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' INCLUDE 'Map_To_NCEP_Model_Coordinates.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_VerticalCoordinates diff --git a/test/mains/regression/forward/test_Zeeman/Compare_Diagnostics.inc b/test/mains/regression/forward/test_Zeeman/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/forward/test_Zeeman/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/forward/test_Zeeman/test_Zeeman.f90 b/test/mains/regression/forward/test_Zeeman/test_Zeeman.f90 index d8f95366..d022ad34 100644 --- a/test/mains/regression/forward/test_Zeeman/test_Zeeman.f90 +++ b/test/mains/regression/forward/test_Zeeman/test_Zeeman.f90 @@ -299,6 +299,7 @@ PROGRAM test_Zeeman ELSE Message = 'RTSolution results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution', RTSolution, rts ) ! Write the current RTSolution results to file rts_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution, Quiet=.TRUE. ) @@ -469,4 +470,6 @@ SUBROUTINE Load_AtmSfc_Data() END SUBROUTINE Load_AtmSfc_Data + INCLUDE 'Compare_Diagnostics.inc' + END PROGRAM test_Zeeman diff --git a/test/mains/regression/incsrc/Compare_Diagnostics.inc b/test/mains/regression/incsrc/Compare_Diagnostics.inc new file mode 100644 index 00000000..2fec3238 --- /dev/null +++ b/test/mains/regression/incsrc/Compare_Diagnostics.inc @@ -0,0 +1,520 @@ + ! + ! Include file containing internal subprograms to report comparison diagnostics + ! + SUBROUTINE Report_RTSolution_DiffStats_R1( Program_Name, Label, x, y ) + USE Type_Kinds, ONLY: fp + USE Message_Handler, ONLY: INFORMATION, Display_Message + USE CRTM_RTSolution_Define, ONLY: CRTM_RTSolution_type, CRTM_RTSolution_Associated + CHARACTER(*), INTENT(IN) :: Program_Name + CHARACTER(*), INTENT(IN) :: Label + TYPE(CRTM_RTSolution_type), INTENT(IN) :: x(:), y(:) + REAL(fp) :: min_diff, max_diff, max_rel + LOGICAL :: has_value + CHARACTER(256) :: Message + INTEGER :: i + + IF ( SIZE(x) /= SIZE(y) ) THEN + Message = TRIM(Label)//' diff stats: size mismatch' + CALL Display_Message( Program_Name, Message, INFORMATION ) + RETURN + END IF + + has_value = .FALSE. + + DO i = 1, SIZE(x) + IF ( .NOT. CRTM_RTSolution_Associated(x(i)) .OR. & + .NOT. CRTM_RTSolution_Associated(y(i)) ) CYCLE + + CALL Update_Diff_Scalar( x(i)%SOD, y(i)%SOD, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Surface_Emissivity, y(i)%Surface_Emissivity, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Surface_Reflectivity, y(i)%Surface_Reflectivity, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Up_Radiance, y(i)%Up_Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Down_Radiance, y(i)%Down_Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Down_Solar_Radiance, y(i)%Down_Solar_Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Surface_Planck_Radiance, y(i)%Surface_Planck_Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Total_Cloud_Cover, y(i)%Total_Cloud_Cover, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%R_clear, y(i)%R_clear, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Tb_clear, y(i)%Tb_clear, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Reflectance_clear, y(i)%Reflectance_clear, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Radiance, y(i)%Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Brightness_Temperature, y(i)%Brightness_Temperature, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Solar_Irradiance, y(i)%Solar_Irradiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Reflectance, y(i)%Reflectance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_1D( x(i)%Stokes, y(i)%Stokes, & + min_diff, max_diff, max_rel, has_value ) + + IF ( ALLOCATED(x(i)%Upwelling_Overcast_Radiance) .AND. ALLOCATED(y(i)%Upwelling_Overcast_Radiance) ) THEN + CALL Update_Diff_1D( x(i)%Upwelling_Overcast_Radiance, y(i)%Upwelling_Overcast_Radiance, & + min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Upwelling_Radiance) .AND. ALLOCATED(y(i)%Upwelling_Radiance) ) THEN + CALL Update_Diff_1D( x(i)%Upwelling_Radiance, y(i)%Upwelling_Radiance, & + min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Layer_Optical_Depth) .AND. ALLOCATED(y(i)%Layer_Optical_Depth) ) THEN + CALL Update_Diff_1D( x(i)%Layer_Optical_Depth, y(i)%Layer_Optical_Depth, & + min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Reflectivity) .AND. ALLOCATED(y(i)%Reflectivity) ) THEN + CALL Update_Diff_1D( x(i)%Reflectivity, y(i)%Reflectivity, & + min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Reflectivity_Attenuated) .AND. ALLOCATED(y(i)%Reflectivity_Attenuated) ) THEN + CALL Update_Diff_1D( x(i)%Reflectivity_Attenuated, y(i)%Reflectivity_Attenuated, & + min_diff, max_diff, max_rel, has_value ) + END IF + END DO + + IF ( .NOT. has_value ) THEN + Message = TRIM(Label)//' diff stats: no floating-point components compared' + ELSE + WRITE( Message, '(A,": min abs=",es13.6," max abs=",es13.6," max rel=",es13.6)' ) & + TRIM(Label)//' diff stats', min_diff, max_diff, max_rel + END IF + CALL Display_Message( Program_Name, Message, INFORMATION ) + END SUBROUTINE Report_RTSolution_DiffStats_R1 + + + SUBROUTINE Report_RTSolution_DiffStats( Program_Name, Label, x, y ) + USE Type_Kinds, ONLY: fp + USE Message_Handler, ONLY: INFORMATION, Display_Message + USE CRTM_RTSolution_Define, ONLY: CRTM_RTSolution_type, CRTM_RTSolution_Associated + CHARACTER(*), INTENT(IN) :: Program_Name + CHARACTER(*), INTENT(IN) :: Label + TYPE(CRTM_RTSolution_type), INTENT(IN) :: x(:,:), y(:,:) + REAL(fp) :: min_diff, max_diff, max_rel + LOGICAL :: has_value + CHARACTER(256) :: Message + INTEGER :: i, j + + IF ( ANY(SHAPE(x) /= SHAPE(y)) ) THEN + Message = TRIM(Label)//' diff stats: shape mismatch' + CALL Display_Message( Program_Name, Message, INFORMATION ) + RETURN + END IF + + has_value = .FALSE. + + DO j = 1, SIZE(x, DIM=2) + DO i = 1, SIZE(x, DIM=1) + IF ( .NOT. CRTM_RTSolution_Associated(x(i,j)) .OR. & + .NOT. CRTM_RTSolution_Associated(y(i,j)) ) CYCLE + + CALL Update_Diff_Scalar( x(i,j)%SOD, y(i,j)%SOD, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Surface_Emissivity, y(i,j)%Surface_Emissivity, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Surface_Reflectivity, y(i,j)%Surface_Reflectivity, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Up_Radiance, y(i,j)%Up_Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Down_Radiance, y(i,j)%Down_Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Down_Solar_Radiance, y(i,j)%Down_Solar_Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Surface_Planck_Radiance, y(i,j)%Surface_Planck_Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Total_Cloud_Cover, y(i,j)%Total_Cloud_Cover, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%R_clear, y(i,j)%R_clear, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Tb_clear, y(i,j)%Tb_clear, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Reflectance_clear, y(i,j)%Reflectance_clear, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Radiance, y(i,j)%Radiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Brightness_Temperature, y(i,j)%Brightness_Temperature, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Solar_Irradiance, y(i,j)%Solar_Irradiance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Reflectance, y(i,j)%Reflectance, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_1D( x(i,j)%Stokes, y(i,j)%Stokes, & + min_diff, max_diff, max_rel, has_value ) + + IF ( ALLOCATED(x(i,j)%Upwelling_Overcast_Radiance) .AND. ALLOCATED(y(i,j)%Upwelling_Overcast_Radiance) ) THEN + CALL Update_Diff_1D( x(i,j)%Upwelling_Overcast_Radiance, y(i,j)%Upwelling_Overcast_Radiance, & + min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Upwelling_Radiance) .AND. ALLOCATED(y(i,j)%Upwelling_Radiance) ) THEN + CALL Update_Diff_1D( x(i,j)%Upwelling_Radiance, y(i,j)%Upwelling_Radiance, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Layer_Optical_Depth) .AND. ALLOCATED(y(i,j)%Layer_Optical_Depth) ) THEN + CALL Update_Diff_1D( x(i,j)%Layer_Optical_Depth, y(i,j)%Layer_Optical_Depth, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Reflectivity) .AND. ALLOCATED(y(i,j)%Reflectivity) ) THEN + CALL Update_Diff_1D( x(i,j)%Reflectivity, y(i,j)%Reflectivity, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Reflectivity_Attenuated) .AND. ALLOCATED(y(i,j)%Reflectivity_Attenuated) ) THEN + CALL Update_Diff_1D( x(i,j)%Reflectivity_Attenuated, y(i,j)%Reflectivity_Attenuated, & + min_diff, max_diff, max_rel, has_value ) + END IF + END DO + END DO + + IF ( .NOT. has_value ) THEN + Message = TRIM(Label)//' diff stats: no floating-point components compared' + ELSE + WRITE( Message, '(A,": min abs=",es13.6," max abs=",es13.6," max rel=",es13.6)' ) & + TRIM(Label)//' diff stats', min_diff, max_diff, max_rel + END IF + CALL Display_Message( Program_Name, Message, INFORMATION ) + END SUBROUTINE Report_RTSolution_DiffStats + + + SUBROUTINE Report_Atmosphere_DiffStats_R1( Program_Name, Label, x, y ) + USE Type_Kinds, ONLY: fp + USE Message_Handler, ONLY: INFORMATION, Display_Message + USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type, CRTM_Atmosphere_Associated + CHARACTER(*), INTENT(IN) :: Program_Name + CHARACTER(*), INTENT(IN) :: Label + TYPE(CRTM_Atmosphere_type), INTENT(IN) :: x(:), y(:) + REAL(fp) :: min_diff, max_diff, max_rel + LOGICAL :: has_value + CHARACTER(256) :: Message + INTEGER :: i + + IF ( SIZE(x) /= SIZE(y) ) THEN + Message = TRIM(Label)//' diff stats: size mismatch' + CALL Display_Message( Program_Name, Message, INFORMATION ) + RETURN + END IF + + has_value = .FALSE. + + DO i = 1, SIZE(x) + IF ( .NOT. CRTM_Atmosphere_Associated(x(i)) .OR. & + .NOT. CRTM_Atmosphere_Associated(y(i)) ) CYCLE + + IF ( ALLOCATED(x(i)%Level_Pressure) .AND. ALLOCATED(y(i)%Level_Pressure) ) THEN + CALL Update_Diff_1D( x(i)%Level_Pressure, y(i)%Level_Pressure, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Height) .AND. ALLOCATED(y(i)%Height) ) THEN + CALL Update_Diff_1D( x(i)%Height, y(i)%Height, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Pressure) .AND. ALLOCATED(y(i)%Pressure) ) THEN + CALL Update_Diff_1D( x(i)%Pressure, y(i)%Pressure, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Temperature) .AND. ALLOCATED(y(i)%Temperature) ) THEN + CALL Update_Diff_1D( x(i)%Temperature, y(i)%Temperature, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Relative_Humidity) .AND. ALLOCATED(y(i)%Relative_Humidity) ) THEN + CALL Update_Diff_1D( x(i)%Relative_Humidity, y(i)%Relative_Humidity, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Absorber) .AND. ALLOCATED(y(i)%Absorber) ) THEN + CALL Update_Diff_2D( x(i)%Absorber, y(i)%Absorber, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i)%Cloud_Fraction) .AND. ALLOCATED(y(i)%Cloud_Fraction) ) THEN + CALL Update_Diff_1D( x(i)%Cloud_Fraction, y(i)%Cloud_Fraction, min_diff, max_diff, max_rel, has_value ) + END IF + END DO + + IF ( .NOT. has_value ) THEN + Message = TRIM(Label)//' diff stats: no floating-point components compared' + ELSE + WRITE( Message, '(A,": min abs=",es13.6," max abs=",es13.6," max rel=",es13.6)' ) & + TRIM(Label)//' diff stats', min_diff, max_diff, max_rel + END IF + CALL Display_Message( Program_Name, Message, INFORMATION ) + END SUBROUTINE Report_Atmosphere_DiffStats_R1 + + + SUBROUTINE Report_Atmosphere_DiffStats( Program_Name, Label, x, y ) + USE Type_Kinds, ONLY: fp + USE Message_Handler, ONLY: INFORMATION, Display_Message + USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type, CRTM_Atmosphere_Associated + CHARACTER(*), INTENT(IN) :: Program_Name + CHARACTER(*), INTENT(IN) :: Label + TYPE(CRTM_Atmosphere_type), INTENT(IN) :: x(:,:), y(:,:) + REAL(fp) :: min_diff, max_diff, max_rel + LOGICAL :: has_value + CHARACTER(256) :: Message + INTEGER :: i, j + + IF ( ANY(SHAPE(x) /= SHAPE(y)) ) THEN + Message = TRIM(Label)//' diff stats: shape mismatch' + CALL Display_Message( Program_Name, Message, INFORMATION ) + RETURN + END IF + + has_value = .FALSE. + + DO j = 1, SIZE(x, DIM=2) + DO i = 1, SIZE(x, DIM=1) + IF ( .NOT. CRTM_Atmosphere_Associated(x(i,j)) .OR. & + .NOT. CRTM_Atmosphere_Associated(y(i,j)) ) CYCLE + + IF ( ALLOCATED(x(i,j)%Level_Pressure) .AND. ALLOCATED(y(i,j)%Level_Pressure) ) THEN + CALL Update_Diff_1D( x(i,j)%Level_Pressure, y(i,j)%Level_Pressure, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Height) .AND. ALLOCATED(y(i,j)%Height) ) THEN + CALL Update_Diff_1D( x(i,j)%Height, y(i,j)%Height, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Pressure) .AND. ALLOCATED(y(i,j)%Pressure) ) THEN + CALL Update_Diff_1D( x(i,j)%Pressure, y(i,j)%Pressure, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Temperature) .AND. ALLOCATED(y(i,j)%Temperature) ) THEN + CALL Update_Diff_1D( x(i,j)%Temperature, y(i,j)%Temperature, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Relative_Humidity) .AND. ALLOCATED(y(i,j)%Relative_Humidity) ) THEN + CALL Update_Diff_1D( x(i,j)%Relative_Humidity, y(i,j)%Relative_Humidity, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Absorber) .AND. ALLOCATED(y(i,j)%Absorber) ) THEN + CALL Update_Diff_2D( x(i,j)%Absorber, y(i,j)%Absorber, min_diff, max_diff, max_rel, has_value ) + END IF + IF ( ALLOCATED(x(i,j)%Cloud_Fraction) .AND. ALLOCATED(y(i,j)%Cloud_Fraction) ) THEN + CALL Update_Diff_1D( x(i,j)%Cloud_Fraction, y(i,j)%Cloud_Fraction, min_diff, max_diff, max_rel, has_value ) + END IF + END DO + END DO + + IF ( .NOT. has_value ) THEN + Message = TRIM(Label)//' diff stats: no floating-point components compared' + ELSE + WRITE( Message, '(A,": min abs=",es13.6," max abs=",es13.6," max rel=",es13.6)' ) & + TRIM(Label)//' diff stats', min_diff, max_diff, max_rel + END IF + CALL Display_Message( Program_Name, Message, INFORMATION ) + END SUBROUTINE Report_Atmosphere_DiffStats + + + SUBROUTINE Report_Surface_DiffStats_R1( Program_Name, Label, x, y ) + USE Type_Kinds, ONLY: fp + USE Message_Handler, ONLY: INFORMATION, Display_Message + USE CRTM_Surface_Define, ONLY: CRTM_Surface_type + CHARACTER(*), INTENT(IN) :: Program_Name + CHARACTER(*), INTENT(IN) :: Label + TYPE(CRTM_Surface_type), INTENT(IN) :: x(:), y(:) + REAL(fp) :: min_diff, max_diff, max_rel + LOGICAL :: has_value + CHARACTER(256) :: Message + INTEGER :: i + + IF ( SIZE(x) /= SIZE(y) ) THEN + Message = TRIM(Label)//' diff stats: size mismatch' + CALL Display_Message( Program_Name, Message, INFORMATION ) + RETURN + END IF + + has_value = .FALSE. + + DO i = 1, SIZE(x) + CALL Update_Diff_Scalar( x(i)%Land_Coverage, y(i)%Land_Coverage, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Water_Coverage, y(i)%Water_Coverage, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Snow_Coverage, y(i)%Snow_Coverage, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Ice_Coverage, y(i)%Ice_Coverage, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Land_Temperature, y(i)%Land_Temperature, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Soil_Moisture_Content, y(i)%Soil_Moisture_Content, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Canopy_Water_Content, y(i)%Canopy_Water_Content, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Vegetation_Fraction, y(i)%Vegetation_Fraction, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Soil_Temperature, y(i)%Soil_Temperature, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%LAI, y(i)%LAI, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Water_Temperature, y(i)%Water_Temperature, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Wind_Speed, y(i)%Wind_Speed, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Wind_Direction, y(i)%Wind_Direction, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Salinity, y(i)%Salinity, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Snow_Temperature, y(i)%Snow_Temperature, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Snow_Depth, y(i)%Snow_Depth, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Snow_Density, y(i)%Snow_Density, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Snow_Grain_Size, y(i)%Snow_Grain_Size, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Ice_Temperature, y(i)%Ice_Temperature, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Ice_Thickness, y(i)%Ice_Thickness, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Ice_Density, y(i)%Ice_Density, & + min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i)%Ice_Roughness, y(i)%Ice_Roughness, & + min_diff, max_diff, max_rel, has_value ) + END DO + + IF ( .NOT. has_value ) THEN + Message = TRIM(Label)//' diff stats: no floating-point components compared' + ELSE + WRITE( Message, '(A,": min abs=",es13.6," max abs=",es13.6," max rel=",es13.6)' ) & + TRIM(Label)//' diff stats', min_diff, max_diff, max_rel + END IF + CALL Display_Message( Program_Name, Message, INFORMATION ) + END SUBROUTINE Report_Surface_DiffStats_R1 + + + SUBROUTINE Report_Surface_DiffStats( Program_Name, Label, x, y ) + USE Type_Kinds, ONLY: fp + USE Message_Handler, ONLY: INFORMATION, Display_Message + USE CRTM_Surface_Define, ONLY: CRTM_Surface_type + CHARACTER(*), INTENT(IN) :: Program_Name + CHARACTER(*), INTENT(IN) :: Label + TYPE(CRTM_Surface_type), INTENT(IN) :: x(:,:), y(:,:) + REAL(fp) :: min_diff, max_diff, max_rel + LOGICAL :: has_value + CHARACTER(256) :: Message + INTEGER :: i, j + + IF ( ANY(SHAPE(x) /= SHAPE(y)) ) THEN + Message = TRIM(Label)//' diff stats: shape mismatch' + CALL Display_Message( Program_Name, Message, INFORMATION ) + RETURN + END IF + + has_value = .FALSE. + + DO j = 1, SIZE(x, DIM=2) + DO i = 1, SIZE(x, DIM=1) + CALL Update_Diff_Scalar( x(i,j)%Land_Coverage, y(i,j)%Land_Coverage, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Water_Coverage, y(i,j)%Water_Coverage, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Snow_Coverage, y(i,j)%Snow_Coverage, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Ice_Coverage, y(i,j)%Ice_Coverage, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Land_Temperature, y(i,j)%Land_Temperature, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Soil_Moisture_Content,y(i,j)%Soil_Moisture_Content,min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Canopy_Water_Content, y(i,j)%Canopy_Water_Content, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Vegetation_Fraction, y(i,j)%Vegetation_Fraction, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Soil_Temperature, y(i,j)%Soil_Temperature, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%LAI, y(i,j)%LAI, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Water_Temperature, y(i,j)%Water_Temperature, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Wind_Speed, y(i,j)%Wind_Speed, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Wind_Direction, y(i,j)%Wind_Direction, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Salinity, y(i,j)%Salinity, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Snow_Temperature, y(i,j)%Snow_Temperature, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Snow_Depth, y(i,j)%Snow_Depth, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Snow_Density, y(i,j)%Snow_Density, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Snow_Grain_Size, y(i,j)%Snow_Grain_Size, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Ice_Temperature, y(i,j)%Ice_Temperature, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Ice_Thickness, y(i,j)%Ice_Thickness, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Ice_Density, y(i,j)%Ice_Density, min_diff, max_diff, max_rel, has_value ) + CALL Update_Diff_Scalar( x(i,j)%Ice_Roughness, y(i,j)%Ice_Roughness, min_diff, max_diff, max_rel, has_value ) + END DO + END DO + + IF ( .NOT. has_value ) THEN + Message = TRIM(Label)//' diff stats: no floating-point components compared' + ELSE + WRITE( Message, '(A,": min abs=",es13.6," max abs=",es13.6," max rel=",es13.6)' ) & + TRIM(Label)//' diff stats', min_diff, max_diff, max_rel + END IF + CALL Display_Message( Program_Name, Message, INFORMATION ) + END SUBROUTINE Report_Surface_DiffStats + + + SUBROUTINE Update_Diff_Scalar( x, y, min_diff, max_diff, max_rel, has_value ) + USE Type_Kinds, ONLY: fp + REAL(fp), INTENT(IN) :: x + REAL(fp), INTENT(IN) :: y + REAL(fp), INTENT(IN OUT) :: min_diff + REAL(fp), INTENT(IN OUT) :: max_diff + REAL(fp), INTENT(IN OUT) :: max_rel + LOGICAL, INTENT(IN OUT) :: has_value + REAL(fp) :: diff + REAL(fp) :: rel + REAL(fp) :: denom + + diff = ABS(x - y) + denom = MAX(ABS(y), TINY(1.0_fp)) + rel = diff / denom + IF ( .NOT. has_value ) THEN + min_diff = diff + max_diff = diff + max_rel = rel + has_value = .TRUE. + ELSE + min_diff = MIN(min_diff, diff) + max_diff = MAX(max_diff, diff) + max_rel = MAX(max_rel, rel) + END IF + END SUBROUTINE Update_Diff_Scalar + + + SUBROUTINE Update_Diff_1D( x, y, min_diff, max_diff, max_rel, has_value ) + USE Type_Kinds, ONLY: fp + REAL(fp), INTENT(IN) :: x(:) + REAL(fp), INTENT(IN) :: y(:) + REAL(fp), INTENT(IN OUT) :: min_diff + REAL(fp), INTENT(IN OUT) :: max_diff + REAL(fp), INTENT(IN OUT) :: max_rel + LOGICAL, INTENT(IN OUT) :: has_value + REAL(fp) :: min_local + REAL(fp) :: max_local + REAL(fp) :: max_rel_local + REAL(fp) :: denom_floor + + IF ( SIZE(x) == 0 ) RETURN + min_local = MINVAL( ABS(x - y) ) + max_local = MAXVAL( ABS(x - y) ) + denom_floor = TINY(1.0_fp) + max_rel_local = MAXVAL( ABS(x - y) / MAX(ABS(y), denom_floor) ) + + IF ( .NOT. has_value ) THEN + min_diff = min_local + max_diff = max_local + max_rel = max_rel_local + has_value = .TRUE. + ELSE + min_diff = MIN(min_diff, min_local) + max_diff = MAX(max_diff, max_local) + max_rel = MAX(max_rel, max_rel_local) + END IF + END SUBROUTINE Update_Diff_1D + + + SUBROUTINE Update_Diff_2D( x, y, min_diff, max_diff, max_rel, has_value ) + USE Type_Kinds, ONLY: fp + REAL(fp), INTENT(IN) :: x(:,:) + REAL(fp), INTENT(IN) :: y(:,:) + REAL(fp), INTENT(IN OUT) :: min_diff + REAL(fp), INTENT(IN OUT) :: max_diff + REAL(fp), INTENT(IN OUT) :: max_rel + LOGICAL, INTENT(IN OUT) :: has_value + REAL(fp) :: min_local + REAL(fp) :: max_local + REAL(fp) :: max_rel_local + REAL(fp) :: denom_floor + + IF ( SIZE(x) == 0 ) RETURN + min_local = MINVAL( ABS(x - y) ) + max_local = MAXVAL( ABS(x - y) ) + denom_floor = TINY(1.0_fp) + max_rel_local = MAXVAL( ABS(x - y) / MAX(ABS(y), denom_floor) ) + + IF ( .NOT. has_value ) THEN + min_diff = min_local + max_diff = max_local + max_rel = max_rel_local + has_value = .TRUE. + ELSE + min_diff = MIN(min_diff, min_local) + max_diff = MAX(max_diff, max_local) + max_rel = MAX(max_rel, max_rel_local) + END IF + END SUBROUTINE Update_Diff_2D diff --git a/test/mains/regression/k_matrix/test_AOD/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_AOD/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_AOD/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 b/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 index 0de8e78c..87bbe17d 100644 --- a/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 +++ b/test/mains/regression/k_matrix/test_AOD/test_AOD.f90 @@ -286,6 +286,7 @@ PROGRAM test_AOD ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -330,5 +331,6 @@ PROGRAM test_AOD CONTAINS INCLUDE 'Load_Atm_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_AOD diff --git a/test/mains/regression/k_matrix/test_ChannelSubset/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_ChannelSubset/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_ChannelSubset/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/k_matrix/test_ChannelSubset/test_ChannelSubset.f90 b/test/mains/regression/k_matrix/test_ChannelSubset/test_ChannelSubset.f90 index d8f34991..6ced9fb5 100644 --- a/test/mains/regression/k_matrix/test_ChannelSubset/test_ChannelSubset.f90 +++ b/test/mains/regression/k_matrix/test_ChannelSubset/test_ChannelSubset.f90 @@ -430,6 +430,7 @@ PROGRAM Example6_ChannelSubset ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -446,6 +447,7 @@ PROGRAM Example6_ChannelSubset ELSE Message = 'Surface_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_K', Surface_K, sfc_k ) ! Write the current Surface_K results to file sfck_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_file, Surface_K, Quiet=.TRUE. ) @@ -462,6 +464,7 @@ PROGRAM Example6_ChannelSubset ELSE Message = 'RTSolution_K results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_K', RTSolution_K, rts_k ) rtsk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_K.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsk_File, RTSolution_K, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -515,5 +518,6 @@ PROGRAM Example6_ChannelSubset INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM Example6_ChannelSubset diff --git a/test/mains/regression/k_matrix/test_ClearSky/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_ClearSky/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_ClearSky/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/k_matrix/test_ClearSky/test_ClearSky.f90 b/test/mains/regression/k_matrix/test_ClearSky/test_ClearSky.f90 index 019f347a..6eed6449 100644 --- a/test/mains/regression/k_matrix/test_ClearSky/test_ClearSky.f90 +++ b/test/mains/regression/k_matrix/test_ClearSky/test_ClearSky.f90 @@ -392,6 +392,7 @@ PROGRAM test_ClearSky ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -408,6 +409,7 @@ PROGRAM test_ClearSky ELSE Message = 'Surface_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_K', Surface_K, sfc_k ) ! Write the current Surface_K results to file sfck_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_file, Surface_K, Quiet=.TRUE. ) @@ -424,6 +426,7 @@ PROGRAM test_ClearSky ELSE Message = 'RTSolution_K results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_K', RTSolution_K, rts_k ) rtsk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_K.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsk_File, RTSolution_K, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -474,5 +477,6 @@ PROGRAM test_ClearSky INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_ClearSky diff --git a/test/mains/regression/k_matrix/test_SOI/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_SOI/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_SOI/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/k_matrix/test_SOI/test_SOI.f90 b/test/mains/regression/k_matrix/test_SOI/test_SOI.f90 index f912a3d9..f08e74d7 100644 --- a/test/mains/regression/k_matrix/test_SOI/test_SOI.f90 +++ b/test/mains/regression/k_matrix/test_SOI/test_SOI.f90 @@ -400,6 +400,7 @@ PROGRAM test_SOI ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -416,6 +417,7 @@ PROGRAM test_SOI ELSE Message = 'Surface_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_K', Surface_K, sfc_k ) ! Write the current Surface_K results to file sfck_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_file, Surface_K, Quiet=.TRUE. ) @@ -432,6 +434,7 @@ PROGRAM test_SOI ELSE Message = 'RTSolution_K results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_K', RTSolution_K, rts_k ) rtsk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_K.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsk_File, RTSolution_K, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -477,5 +480,6 @@ PROGRAM test_SOI INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_SOI diff --git a/test/mains/regression/k_matrix/test_SSU/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_SSU/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_SSU/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/k_matrix/test_SSU/test_SSU.f90 b/test/mains/regression/k_matrix/test_SSU/test_SSU.f90 index 9cebcd20..ca357730 100644 --- a/test/mains/regression/k_matrix/test_SSU/test_SSU.f90 +++ b/test/mains/regression/k_matrix/test_SSU/test_SSU.f90 @@ -406,6 +406,7 @@ PROGRAM test_SSU ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -422,6 +423,7 @@ PROGRAM test_SSU ELSE Message = 'Surface_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_K', Surface_K, sfc_k ) ! Write the current Surface_K results to file sfck_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_file, Surface_K, Quiet=.TRUE. ) @@ -438,6 +440,7 @@ PROGRAM test_SSU ELSE Message = 'RTSolution_K results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_K', RTSolution_K, rts_k ) rtsk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_K.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsk_File, RTSolution_K, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -485,5 +488,6 @@ PROGRAM test_SSU INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_SSU diff --git a/test/mains/regression/k_matrix/test_ScatteringSwitch/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_ScatteringSwitch/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_ScatteringSwitch/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/k_matrix/test_ScatteringSwitch/test_ScatteringSwitch.f90 b/test/mains/regression/k_matrix/test_ScatteringSwitch/test_ScatteringSwitch.f90 index 1fccdc0d..cd9983c0 100644 --- a/test/mains/regression/k_matrix/test_ScatteringSwitch/test_ScatteringSwitch.f90 +++ b/test/mains/regression/k_matrix/test_ScatteringSwitch/test_ScatteringSwitch.f90 @@ -401,6 +401,7 @@ PROGRAM test_ScatteringSwitch ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -417,6 +418,7 @@ PROGRAM test_ScatteringSwitch ELSE Message = 'Surface_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_K', Surface_K, sfc_k ) ! Write the current Surface_K results to file sfck_File = TRIM(PROGRAM_NAME)//'_'//TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_file, Surface_K, Quiet=.TRUE. ) @@ -433,6 +435,7 @@ PROGRAM test_ScatteringSwitch ELSE Message = 'RTSolution_K results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_K', RTSolution_K, rts_k ) rtsk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_K.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsk_File, RTSolution_K, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -479,6 +482,7 @@ PROGRAM test_ScatteringSwitch INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_ScatteringSwitch diff --git a/test/mains/regression/k_matrix/test_Simple/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_Simple/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_Simple/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file 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..1b49b731 100644 --- a/test/mains/regression/k_matrix/test_Simple/test_Simple.f90 +++ b/test/mains/regression/k_matrix/test_Simple/test_Simple.f90 @@ -404,6 +404,7 @@ PROGRAM test_Simple ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -420,6 +421,7 @@ PROGRAM test_Simple ELSE Message = 'Surface_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_K', Surface_K, sfc_k ) ! Write the current Surface_K results to file sfck_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_file, Surface_K, Quiet=.TRUE. ) @@ -436,6 +438,7 @@ PROGRAM test_Simple ELSE Message = 'RTSolution_K results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_K', RTSolution_K, rts_k ) rtsk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_K.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsk_File, RTSolution_K, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -486,5 +489,6 @@ PROGRAM test_Simple INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_Simple diff --git a/test/mains/regression/k_matrix/test_User_Emissivity/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_User_Emissivity/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_User_Emissivity/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/k_matrix/test_User_Emissivity/test_User_Emissivity.f90 b/test/mains/regression/k_matrix/test_User_Emissivity/test_User_Emissivity.f90 index 4e3162cb..ab1a6a05 100644 --- a/test/mains/regression/k_matrix/test_User_Emissivity/test_User_Emissivity.f90 +++ b/test/mains/regression/k_matrix/test_User_Emissivity/test_User_Emissivity.f90 @@ -412,6 +412,7 @@ PROGRAM test_User_Emissivity ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -428,6 +429,7 @@ PROGRAM test_User_Emissivity ELSE Message = 'Surface_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_K', Surface_K, sfc_k ) ! Write the current Surface_K results to file sfck_File = TRIM(PROGRAM_NAME)//'_'//TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_file, Surface_K, Quiet=.TRUE. ) @@ -444,6 +446,7 @@ PROGRAM test_User_Emissivity ELSE Message = 'RTSolution_K results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_K', RTSolution_K, rts_k ) rtsk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_K.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsk_File, RTSolution_K, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -490,6 +493,7 @@ PROGRAM test_User_Emissivity INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_User_Emissivity diff --git a/test/mains/regression/k_matrix/test_VerticalCoordinates/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_VerticalCoordinates/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_VerticalCoordinates/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/k_matrix/test_VerticalCoordinates/test_VerticalCoordinates.f90 b/test/mains/regression/k_matrix/test_VerticalCoordinates/test_VerticalCoordinates.f90 index 27b2e7fc..f9dfe2ca 100644 --- a/test/mains/regression/k_matrix/test_VerticalCoordinates/test_VerticalCoordinates.f90 +++ b/test/mains/regression/k_matrix/test_VerticalCoordinates/test_VerticalCoordinates.f90 @@ -611,6 +611,7 @@ PROGRAM test_VerticalCoordinates ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -627,6 +628,7 @@ PROGRAM test_VerticalCoordinates ELSE Message = 'Atmosphere_NAM_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_NAM_K', Atmosphere_NAM_K, atm_NAM_k ) ! Write the current Atmosphere_NAM_K results to file atmk_NAM_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.NAM.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_NAM_file, Atmosphere_NAM_K, Quiet=.TRUE. ) @@ -643,6 +645,7 @@ PROGRAM test_VerticalCoordinates ELSE Message = 'Atmosphere_GFS_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_GFS_K', Atmosphere_GFS_K, atm_GFS_k ) ! Write the current Atmosphere_GFS_K results to file atmk_GFS_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.GFS.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_GFS_file, Atmosphere_GFS_K, Quiet=.TRUE. ) @@ -661,6 +664,7 @@ PROGRAM test_VerticalCoordinates ELSE Message = 'Surface_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_K', Surface_K, sfc_k ) ! Write the current Surface_K results to file sfck_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_file, Surface_K, Quiet=.TRUE. ) @@ -677,6 +681,7 @@ PROGRAM test_VerticalCoordinates ELSE Message = 'Surface_NAM_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_NAM_K', Surface_NAM_K, sfc_NAM_k ) ! Write the current Surface_NAM_K results to file sfck_NAM_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.NAM.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_NAM_file, Surface_NAM_K, Quiet=.TRUE. ) @@ -693,6 +698,7 @@ PROGRAM test_VerticalCoordinates ELSE Message = 'Surface_GFS_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_GFS_K', Surface_GFS_K, sfc_GFS_k ) ! Write the current Surface_GFS_K results to file sfck_GFS_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.GFS.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_GFS_file, Surface_GFS_K, Quiet=.TRUE. ) @@ -752,5 +758,6 @@ PROGRAM test_VerticalCoordinates INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' INCLUDE 'Map_To_NCEP_Model_Coordinates.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_VerticalCoordinates diff --git a/test/mains/regression/k_matrix/test_Zeeman/Compare_Diagnostics.inc b/test/mains/regression/k_matrix/test_Zeeman/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/k_matrix/test_Zeeman/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/k_matrix/test_Zeeman/test_Zeeman.f90 b/test/mains/regression/k_matrix/test_Zeeman/test_Zeeman.f90 index 2320f92f..f359328d 100644 --- a/test/mains/regression/k_matrix/test_Zeeman/test_Zeeman.f90 +++ b/test/mains/regression/k_matrix/test_Zeeman/test_Zeeman.f90 @@ -422,6 +422,7 @@ PROGRAM test_Zeeman ELSE Message = 'Atmosphere_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Atmosphere_DiffStats( PROGRAM_NAME, 'Atmosphere_K', Atmosphere_K, atm_k ) ! Write the current Atmosphere_K results to file atmk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin' Error_Status = CRTM_Atmosphere_WriteFile( atmk_file, Atmosphere_K, Quiet=.TRUE. ) @@ -438,6 +439,7 @@ PROGRAM test_Zeeman ELSE Message = 'Surface_K Jacobians are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_Surface_DiffStats( PROGRAM_NAME, 'Surface_K', Surface_K, sfc_k ) ! Write the current Surface_K results to file sfck_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin' Error_Status = CRTM_Surface_WriteFile( sfck_file, Surface_K, Quiet=.TRUE. ) @@ -454,6 +456,7 @@ PROGRAM test_Zeeman ELSE Message = 'RTSolution_K results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_K', RTSolution_K, rts_k ) rtsk_File = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.RTSolution_K.bin' Error_Status = CRTM_RTSolution_WriteFile( rtsk_File, RTSolution_K, Quiet=.TRUE. ) IF ( Error_Status /= SUCCESS ) THEN @@ -627,4 +630,6 @@ SUBROUTINE Load_AtmSfc_Data() END SUBROUTINE Load_AtmSfc_Data + INCLUDE 'Compare_Diagnostics.inc' + END PROGRAM test_Zeeman diff --git a/test/mains/regression/tangent_linear/test_ClearSky/Compare_Diagnostics.inc b/test/mains/regression/tangent_linear/test_ClearSky/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/tangent_linear/test_ClearSky/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file diff --git a/test/mains/regression/tangent_linear/test_ClearSky/test_ClearSky.f90 b/test/mains/regression/tangent_linear/test_ClearSky/test_ClearSky.f90 index a063dc81..1bd35e03 100644 --- a/test/mains/regression/tangent_linear/test_ClearSky/test_ClearSky.f90 +++ b/test/mains/regression/tangent_linear/test_ClearSky/test_ClearSky.f90 @@ -310,6 +310,7 @@ PROGRAM test_ClearSky ELSE Message = 'RTSolution_TL results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_TL', RTSolution_TL, rts_TL ) ! Write the current RTSolution results to file rts_File = TRIM(Sensor_Id)//'.RTSolution_TL.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution_TL, Quiet=.TRUE. ) @@ -354,5 +355,6 @@ PROGRAM test_ClearSky INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_ClearSky diff --git a/test/mains/regression/tangent_linear/test_Simple/Compare_Diagnostics.inc b/test/mains/regression/tangent_linear/test_Simple/Compare_Diagnostics.inc new file mode 120000 index 00000000..7cb1bb33 --- /dev/null +++ b/test/mains/regression/tangent_linear/test_Simple/Compare_Diagnostics.inc @@ -0,0 +1 @@ +../../incsrc/Compare_Diagnostics.inc \ No newline at end of file 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..68325874 100644 --- a/test/mains/regression/tangent_linear/test_Simple/test_Simple.f90 +++ b/test/mains/regression/tangent_linear/test_Simple/test_Simple.f90 @@ -308,6 +308,7 @@ PROGRAM test_Simple ELSE Message = 'RTSolution_TL results are different!' CALL Display_Message( PROGRAM_NAME, Message, FAILURE ) + CALL Report_RTSolution_DiffStats( PROGRAM_NAME, 'RTSolution_TL', RTSolution_TL, rts_TL ) ! Write the current RTSolution results to file rts_File = TRIM(Sensor_Id)//'.RTSolution_TL.bin' Error_Status = CRTM_RTSolution_WriteFile( rts_File, RTSolution_TL, Quiet=.TRUE. ) @@ -350,5 +351,6 @@ PROGRAM test_Simple INCLUDE 'Load_Atm_Data.inc' INCLUDE 'Load_Sfc_Data.inc' + INCLUDE 'Compare_Diagnostics.inc' END PROGRAM test_Simple