Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New uncertainty characterisation updates #100

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

New uncertainty characterisation updates #100

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions common/dealloc_output_data.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
26 changes: 0 additions & 26 deletions common/def_output_primary.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good idea. Should also delete them from the output_data structure.

valid_min = output_data%niter_vmin, &
valid_max = output_data%niter_vmax, &
units = '1', &
Expand Down Expand Up @@ -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', &
Expand All @@ -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', &
Expand Down Expand Up @@ -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', &
Expand Down Expand Up @@ -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', &
Expand All @@ -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', &
Expand Down Expand Up @@ -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', &
Expand Down Expand Up @@ -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, &
Expand All @@ -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', &
Expand All @@ -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, &
Expand All @@ -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, &
Expand All @@ -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, &
Expand All @@ -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, &
Expand Down
12 changes: 10 additions & 2 deletions common/def_output_secondary.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't these just repeat the block immediately before the if?

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, &
Expand Down
2 changes: 2 additions & 0 deletions src/ctrl.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
22 changes: 14 additions & 8 deletions src/get_measurements.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
!-------------------------------------------------------------------------------
Expand Down Expand Up @@ -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
Expand All @@ -199,16 +203,17 @@ 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
else
! 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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment explaining why we're removing the \cos\theta. We won't remember next time we go to war with The Angles.

! Compute uncertainty in terms of radiance
if (SAD_Chan(ii)%Solar%SNR > 0.0) then
if (MSI_Data%cal_gain(ii) > 0.0) then
Expand All @@ -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
Expand Down
18 changes: 15 additions & 3 deletions src/read_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.


!----------------------------------------------------------------------------
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.)
Expand Down
16 changes: 15 additions & 1 deletion src/read_sad_lut.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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?!
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The units in the files agree with you here. I shall check what Don meant.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don says the values are not squared and so they need to be here.

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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't this done on line 1198?

end do
end do
else
Expand All @@ -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

Expand Down