diff --git a/common/dealloc_output_data.F90 b/common/dealloc_output_data.F90 index 20a03429..3afdf9fa 100644 --- a/common/dealloc_output_data.F90 +++ b/common/dealloc_output_data.F90 @@ -272,6 +272,15 @@ subroutine dealloc_output_data_secondary(data) deallocate(data%channels_vmax) deallocate(data%channels) + if (associated(data%Sy)) then + deallocate(data%Sy) + deallocate(data%vid_Sy) + deallocate(data%Sy_scale) + deallocate(data%Sy_offset) + deallocate(data%Sy_vmin) + deallocate(data%Sy_vmax) + end if + deallocate(data%vid_y0) deallocate(data%y0_scale) deallocate(data%y0_offset) diff --git a/common/def_output_primary.F90 b/common/def_output_primary.F90 index 2ab5b5dd..a4a6a44e 100644 --- a/common/def_output_primary.F90 +++ b/common/def_output_primary.F90 @@ -1582,8 +1582,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'number of retrieval iterations', & standard_name = '', & fill_value = byte_fill_value, & - scale_factor = output_data%niter_scale, & - add_offset = output_data%niter_offset, & valid_min = output_data%niter_vmin, & valid_max = output_data%niter_vmax, & units = '1', & @@ -1710,8 +1708,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'land/sea flag', & standard_name = 'land_binary_mask', & fill_value = byte_fill_value, & - scale_factor = output_data%lsflag_scale, & - add_offset = output_data%lsflag_offset, & valid_min = output_data%lsflag_vmin, & valid_max = output_data%lsflag_vmax, & units = '1', & @@ -1732,8 +1728,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'land use flag', & standard_name = 'land_use_mask', & fill_value = byte_fill_value, & - scale_factor = output_data%lusflag_scale, & - add_offset = output_data%lusflag_offset, & valid_min = output_data%lusflag_vmin, & valid_max = output_data%lusflag_vmax, & units = '1', & @@ -1801,8 +1795,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'illumination flag', & standard_name = '', & fill_value = byte_fill_value, & - scale_factor = output_data%illum_scale, & - add_offset = output_data%illum_offset, & valid_min = output_data%illum_vmin, & valid_max = output_data%illum_vmax, & units = '1', & @@ -1834,8 +1826,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'Pavolonis cloud type', & standard_name = '', & fill_value = byte_fill_value, & - scale_factor = output_data%cldtype_scale, & - add_offset = output_data%cldtype_offset, & valid_min = output_data%cldtype_vmin, & valid_max = output_data%cldtype_vmax, & units = '1', & @@ -1857,8 +1847,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'Neural net cloud mask (radiance based)', & standard_name = '', & fill_value = byte_fill_value, & - scale_factor = output_data%cldmask_scale, & - add_offset = output_data%cldmask_offset, & valid_min = output_data%cldmask_vmin, & valid_max = output_data%cldmask_vmax, & units = '1', & @@ -1903,8 +1891,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'Neural net cloud phase mask (radiance based)', & standard_name = '', & fill_value = byte_fill_value, & - scale_factor = output_data%ann_phase_scale, & - add_offset = output_data%ann_phase_offset, & valid_min = output_data%ann_phase_vmin, & valid_max = output_data%ann_phase_vmax, & units = '1', & @@ -1986,8 +1972,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & standard_name = 'thermodynamic_phase_of_cloud_water_particles_' // & 'at_cloud_top', & fill_value = byte_fill_value, & - scale_factor = output_data%phase_scale, & - add_offset = output_data%phase_offset, & valid_min = output_data%phase_vmin, & valid_max = output_data%phase_vmax, & flag_values = input_dummy2, & @@ -2011,8 +1995,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & standard_name = 'thermodynamic_phase_of_cloud_water_particles_' // & 'at_cloud_top', & fill_value = byte_fill_value, & - scale_factor = output_data%phase_pavolonis_scale, & - add_offset = output_data%phase_pavolonis_offset, & valid_min = output_data%phase_pavolonis_vmin, & valid_max = output_data%phase_pavolonis_vmax, & units = '1', & @@ -2035,8 +2017,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'instrument channel index', & standard_name = '', & fill_value = byte_fill_value, & - scale_factor = output_data%y_id_scale, & - add_offset = output_data%y_id_offset, & valid_min = output_data%y_id_vmin, & valid_max = output_data%y_id_vmax, & deflate_level = deflate_level, & @@ -2054,8 +2034,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'instrument channel flags', & standard_name = '', & fill_value = byte_fill_value, & - scale_factor = output_data%ch_is_scale, & - add_offset = output_data%ch_is_offset, & valid_min = output_data%ch_is_vmin, & valid_max = output_data%ch_is_vmax, & deflate_level = deflate_level, & @@ -2073,8 +2051,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'surface reflectance output flags', & standard_name = '', & fill_value = byte_fill_value, & - scale_factor = output_data%rho_flags_scale, & - add_offset = output_data%rho_flags_offset, & valid_min = output_data%rho_flags_vmin, & valid_max = output_data%rho_flags_vmax, & deflate_level = deflate_level, & @@ -2094,8 +2070,6 @@ subroutine def_output_primary(ncid, dim3d_var, output_data, indexing, & long_name = 'instrument channel view index', & standard_name = '', & fill_value = byte_fill_value, & - scale_factor = output_data%view_id_scale, & - add_offset = output_data%view_id_offset, & valid_min = output_data%view_id_vmin, & valid_max = output_data%view_id_vmax, & deflate_level = deflate_level, & diff --git a/common/def_output_secondary.F90 b/common/def_output_secondary.F90 index 3d27a5c5..67b69eef 100644 --- a/common/def_output_secondary.F90 +++ b/common/def_output_secondary.F90 @@ -750,9 +750,17 @@ subroutine def_output_secondary(ncid, dim3d_var, output_data, indexing, & input_dummy2 = 'measurement uncertainty in channel no '// & trim(adjustl(input_num)) if (btest(indexing%Ch_Is(i), ThermalBit)) then - input_dummy3 = 'kelvin' + input_dummy='measurement_uncertainty_in_channel_no_'// & + trim(adjustl(input_num)) + input_dummy2='measurement uncertainty in channel no '// & + trim(adjustl(input_num)) + input_dummy3='kelvin' else - input_dummy3 = '1' + input_dummy='measurement_uncertainty_in_channel_no_'// & + trim(adjustl(input_num)) + input_dummy2='measurement uncertainty in channel no '// & + trim(adjustl(input_num)) + input_dummy3='1' end if call ncdf_def_var_short_packed_float( & ncid, & diff --git a/src/ctrl.F90 b/src/ctrl.F90 index 981dad37..da321522 100644 --- a/src/ctrl.F90 +++ b/src/ctrl.F90 @@ -285,6 +285,8 @@ module Ctrl_m logical :: process_aerosol_only logical :: all_channels_same_view logical :: use_ann_phase + logical :: use_new_meas_error ! Enable or disable new LUT + ! uncertainty characterisation integer :: NTypes_to_process ! # of valid values in above integer(byte) :: Types_to_process(MaxTypes) ! Pavolonis (or other) ! type codes for pixels to diff --git a/src/get_measurements.F90 b/src/get_measurements.F90 index ec35b2ed..c5249ed4 100644 --- a/src/get_measurements.F90 +++ b/src/get_measurements.F90 @@ -98,7 +98,11 @@ ! that uses the observed radiance rather than the solar constant to convert ! Solar%NeHomog/Coreg from reflectance into radiance. Additionally, add ! Thermal%NeHomog/Coreg to IDay pixels as well. -! +! 2023/10/10, GT: Introduced Ctrl%use_new_meas_error as the method of switching +! between the old measurement uncertainty characterisation (associated with +! the old text-based LUTs) and the new one. Note that this doesn't affect +! the calculations involving NeHomog or Coreg, as the new LUTs simply don't +! support the SAD-file method of setting these values. ! Bugs: ! None known. !------------------------------------------------------------------------------- @@ -185,7 +189,7 @@ subroutine Get_Measurements(Ctrl, SAD_Chan, SPixel, MSI_Data, status) do i = 1, SPixel%Ind%Ny ii = SPixel%spixel_y_to_ctrl_y_index(i) if (SAD_Chan(ii)%Thermal%Flag > 0) then - if (len_trim(Ctrl%LUTClass) == 3) then + if (.not. Ctrl%use_new_meas_error) then ! Old LUTs approach (fixed uncertainty) ! Note: NEBT is already squared when old LUT is read in. SPixel%Sy(i,i) = SAD_Chan(ii)%Thermal%NEBT @@ -199,7 +203,7 @@ subroutine Get_Measurements(Ctrl, SAD_Chan, SPixel, MSI_Data, status) (dR_dT0(1) / dR_dTm(1))) ** 2 end if else - if (len_trim(Ctrl%LUTClass) == 3) then + if (.not. Ctrl%use_new_meas_error) then ! Old LUTs approach (fixed uncertainty) ! Note: NEDR is squared when old LUT is read in. SPixel%Sy(i,i) = SAD_Chan(ii)%Solar%NEDR @@ -207,8 +211,9 @@ subroutine Get_Measurements(Ctrl, SAD_Chan, SPixel, MSI_Data, status) ! New LUTs approach (varying uncertainty) ! Convert reflectance to radiance j = SPixel%ViewIdx(i) - Lx = (SPixel%Ym(i) * cos(SPixel%Geom%SolZen(j) * d2r) * & - SAD_Chan(ii)%Solar%F0) / Pi + !Lx = (SPixel%Ym(i) * cos(SPixel%Geom%SolZen(j) * d2r) * & + ! SAD_Chan(ii)%Solar%F0) / Pi + Lx = (100.0 * SPixel%Ym(i) * SAD_Chan(ii)%Solar%F0) / Pi ! Compute uncertainty in terms of radiance if (SAD_Chan(ii)%Solar%SNR > 0.0) then if (MSI_Data%cal_gain(ii) > 0.0) then @@ -227,9 +232,10 @@ subroutine Get_Measurements(Ctrl, SAD_Chan, SPixel, MSI_Data, status) end if ! Convert radiance uncertainty to reflectance uncertainty ! then square it to convert to reflectance variance - SPixel%Sy(i,i) = ((Pi * dLx) / & - (cos(SPixel%Geom%SolZen(j) * d2r) * & - SAD_Chan(ii)%Solar%F0)) ** 2 + !SPixel%Sy(i,i) = ((Pi * dLx) / & + ! (cos(SPixel%Geom%SolZen(j) * d2r) * & + ! SAD_Chan(ii)%Solar%F0)) ** 2 + SPixel%Sy(i,i) = ((Pi * dLx) / (100.0 * SAD_Chan(ii)%Solar%F0)) ** 2 end if end if end do diff --git a/src/read_driver.F90 b/src/read_driver.F90 index 392fc177..f65a9616 100644 --- a/src/read_driver.F90 +++ b/src/read_driver.F90 @@ -148,6 +148,8 @@ ! 2019/08/14, SP: Add Fengyun4A support. ! 2021/04/06, AP: New LUT names. ! 2022/01/27, GT: Added CTP input file to Ctrl%FID structure. +! 2023/10/10, GT: Added a check to override an attempt to use the new-LUT +! uncertainty calculation with old text-based LUTs. ! ! Bugs: ! None known. @@ -239,6 +241,7 @@ subroutine Read_Driver(Ctrl, global_atts, source_atts) Ctrl%process_cloudy_only = .true. Ctrl%process_aerosol_only = .false. Ctrl%force_nighttime_retrieval = .false. + Ctrl%use_new_meas_error = .true. !---------------------------------------------------------------------------- @@ -418,10 +421,19 @@ subroutine Read_Driver(Ctrl, global_atts, source_atts) end if end if + ! Check if we're using the old text-based LUTs, which do not include + ! the new measurement uncertainty characterisation + if (len_trim(Ctrl%LUTClass) == 3 .and. Ctrl%use_new_meas_error) then + Ctrl%use_new_meas_error = .false. + write(*,*) 'Warning: Using old text-LUTs. Old uncertainty characterisation forced.' + end if + if (Ctrl%verbose) & + write(*,*) 'use_new_meas_error: ', Ctrl%use_new_meas_error + if (Ctrl%Class == -1) then ! Class not set, so deduce it from the LUTClass and/or Approach if (Ctrl%Approach == AppCld1L .and. (Ctrl%LUTClass(1:3) == 'WAT' .or. & - Ctrl%LUTClass(1:12) == 'liquid-water')) then + index(Ctrl%LUTClass, 'liquid-water') .gt. 0)) then Ctrl%Class = ClsCldWat else if (Ctrl%Approach == AppCld1L .and. (Ctrl%LUTClass(1:3) == 'ICE' .or. & Ctrl%LUTClass(1:9) == 'water-ice')) then @@ -519,8 +531,8 @@ subroutine Read_Driver(Ctrl, global_atts, source_atts) !----------------------- Ctrl%EqMPN -------------------- Ctrl%EqMPN%SySelm = switch_app(a, Default=SelmAux) - Ctrl%EqMPN%Homog = switch_app(a, Default=.true., Aer=.false.) - Ctrl%EqMPN%Coreg = switch_app(a, Default=.true., Aer=.false.) + Ctrl%EqMPN%Homog = switch_app(a, Default=.true., Aer=.false.) + Ctrl%EqMPN%Coreg = switch_app(a, Default=.true., Aer=.false.) !----------------------- Ctrl%Invpar ------------------- Ctrl%Invpar%ConvTest = switch_app(a, Default=.false., Aer=.true.) diff --git a/src/read_sad_lut.F90 b/src/read_sad_lut.F90 index 4798e508..08f57e14 100644 --- a/src/read_sad_lut.F90 +++ b/src/read_sad_lut.F90 @@ -15,6 +15,9 @@ ! grids consistent in dimension size, spacing and vertex values. ! 2017/03/16, GT: Changes for single-view aerosol retrieval mode. ! 2021/03/08, AP: Gather grid dimensions into LUT_Grid_t +! 2023/10/10, GT: Added calculation of old NedR (solar channel noise +! equivalent radiance) Read_NCDF_SAD_LUT(), to allow use of old +! uncertainty calculation with new LUTs. ! ! Bugs: ! - Arrays are allocated excessively large as dimensions of the text tables @@ -1225,8 +1228,12 @@ subroutine Read_NCDF_SAD_LUT(Ctrl, LUTFilename, SAD_Chan, SAD_LUT) ! Uncertainty = a**2 L**2 + b**2 L + c**2 call ncdf_read_array(fid, "ru"//char(j+96), chan_tmp_real) do i = 1, Ctrl%Ind%NSolar + ! ru values stored in LUT already squared?! SAD_Chan(Ctrl%Ind%YSolar(i))%Solar%ru2(j) = & - chan_tmp_real(solar_indices(i)) ** 2 + chan_tmp_real(solar_indices(i)) !** 2 + ! Ensure SNR is explicitly set to 0 (used to determine if ru2 + ! should be used in get_measurments) + SAD_Chan(Ctrl%Ind%YSolar(i))%Solar%SNR = 0.0 end do end do else @@ -1237,6 +1244,13 @@ subroutine Read_NCDF_SAD_LUT(Ctrl, LUTFilename, SAD_Chan, SAD_LUT) chan_tmp_real(solar_indices(i)) end do end if + call ncdf_read_array(fid, "oldnefr", chan_tmp_real) + do i = 1, Ctrl%Ind%NSolar + SAD_Chan(Ctrl%Ind%YSolar(i))%Solar%NedR = & + ( chan_tmp_real(solar_indices(i)) / & + SAD_Chan(Ctrl%Ind%YSolar(i))%Solar%F0 )**2 + end do + deallocate(chan_tmp_real) end if