From 6731e2843196e4a5a067e5d4259679631c134a89 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Mon, 9 Dec 2024 11:51:18 -0700 Subject: [PATCH 1/7] Existing code did not take appropriate action when the bnrh transform failed because the sample ensemble covariance was not positive. The existing code simply returned the input ensemble but did not indicate that an error had occurred. This resulted in the untransformed ensemble being used as if it had been transformed in filter_mod.f90 for inflation and in assim_tools_mod.f90 for assimilation. In general, this was not a problem because zero covariance in the ensemble led to both inflation and assimilation doing nothing. However, due to numerical round-off, it was possible to get cases where the inflation could result in violating the bounds of the bnrh distribution. This commit attempts to correct this behavior by updates to probit_transform_mod.f90, filter_mod.f90, and assim_tools_mod.f90. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit probit_transform_mod.f90: An error return argument was added to subroutine to_probit_bounded_normal_rh. It returns a 1 if the ensemble standard deviation was not positive (failure) and a 0 otherwise. This returned argument is passed up through transform_to_probit and transform_all_to_probit, where it becomes an array with one entry for each state variable being transformed. In the case of an error, the returned ‘transformed’ ensemble is simply the input original ensemble. The check for a non-positive standard deviation was eliminated in the case where ‘use_input_p’ is true since this should never be called if the original transformation failed. Subroutine transform_all_from_probit was modified to include an additional logical array that indicates whether each of the state variables to be transformed was successfully transformed to probit space. The inverse transform is only computed for the state variables for which the original transform was successful. filter_mod.f90: All changes are in subroutine filter_ensemble_inflate. The inflation is only performed if the error return from the transform_to_probit is 0. Note that this is called for each group if a group filter is being run and it would be possible for some groups to transform successfully while others failed. This case is handled properly. assim_tools_mod.f90: Two logical arrays were added to keep track of which transforms fail for the observed variables (obs_probit_trans_ok) and state variables (state_probit_trans_ok); these are allocated on entry into subroutine filter_assim and deallocated on exit. The result of the transforms is only copied into the appropriate ensemble_handle storage if the transform was successful. When observations are sequentially assimilated, the observation ensemble is only converted back to physical space if its original transform was successful. The observation prior and posterior are transformed back to probit space in the SEQUENTIAL_OBS loop. If this transform fails for any group, the loop is cycled for this observation so that it has no impact. In the STATE_UPDATE loop, if the transform to probit space for the current state variable failed, the loop is cycled. The same is done in the OBS_UPDATE loop. Note that groups are not relevant to this step. Finally, when transform_all_from_probit is called for state variables, only those that were successfully transformed originally are transformed back. --- .../modules/assimilation/assim_tools_mod.f90 | 49 ++++++--- .../modules/assimilation/filter_mod.f90 | 26 +++-- .../assimilation/probit_transform_mod.f90 | 102 +++++++++++------- 3 files changed, 111 insertions(+), 66 deletions(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index 5787c0beb..f5c7147ba 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -346,7 +346,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & integer(i8), allocatable :: my_state_indx(:) integer(i8), allocatable :: my_obs_indx(:) -integer :: my_num_obs, i, j, owner, owners_index, my_num_state +integer :: my_num_obs, i, j, owner, owners_index, my_num_state, ierr integer :: obs_mean_index, obs_var_index integer :: grp_beg(num_groups), grp_end(num_groups), grp_size, grp_bot, grp_top, group integer :: num_close_obs, obs_index, num_close_states @@ -373,6 +373,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & type(obs_def_type) :: obs_def type(time_type) :: obs_time +logical, allocatable :: obs_probit_trans_ok(:), state_probit_trans_ok(:) logical :: allow_missing_in_state logical :: local_single_ss_inflate logical :: local_varying_ss_inflate @@ -395,13 +396,15 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & my_obs_indx( obs_ens_handle%my_num_vars), & my_obs_kind( obs_ens_handle%my_num_vars), & my_obs_type( obs_ens_handle%my_num_vars), & - my_obs_loc( obs_ens_handle%my_num_vars)) + my_obs_loc( obs_ens_handle%my_num_vars), & + obs_probit_trans_ok(obs_ens_handle%my_num_vars)) allocate(close_state_dist( ens_handle%my_num_vars), & close_state_ind( ens_handle%my_num_vars), & my_state_indx( ens_handle%my_num_vars), & my_state_kind( ens_handle%my_num_vars), & - my_state_loc( ens_handle%my_num_vars)) + my_state_loc( ens_handle%my_num_vars), & + state_probit_trans_ok(ens_handle%my_num_vars)) ! end alloc ! Initialize assim_tools_module if needed @@ -503,8 +506,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Convert all my state variables to appropriate probit space call transform_to_probit(ens_size, ens_handle%copies(1:ens_size, i), dist_for_state, & state_dist_params(i), probit_ens, .false., & - bounded_below, bounded_above, lower_bound, upper_bound) - ens_handle%copies(1:ens_size, i) = probit_ens + bounded_below, bounded_above, lower_bound, upper_bound, ierr) + state_probit_trans_ok(i) = (ierr == 0) + if(state_probit_trans_ok(1)) ens_handle%copies(1:ens_size, i) = probit_ens end do !> optionally convert all state location verticals @@ -541,8 +545,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Convert all my obs (extended state) variables to appropriate probit space call transform_to_probit(ens_size, obs_ens_handle%copies(1:ens_size, i), dist_for_obs, & obs_dist_params(i), probit_ens, .false., & - bounded_below, bounded_above, lower_bound, upper_bound) - obs_ens_handle%copies(1:ens_size, i) = probit_ens + bounded_below, bounded_above, lower_bound, upper_bound, ierr) + obs_probit_trans_ok(i) = (ierr == 0) + if(obs_probit_trans_ok(i)) obs_ens_handle%copies(1:ens_size, i) = probit_ens endif end do @@ -635,9 +640,11 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & orig_obs_prior_var = obs_ens_handle%copies(OBS_PRIOR_VAR_START: & OBS_PRIOR_VAR_END, owners_index) - ! If QC is okay, convert this observation ensemble from probit to regular space - call transform_from_probit(ens_size, obs_ens_handle%copies(1:ens_size, owners_index) , & - obs_dist_params(owners_index), obs_ens_handle%copies(1:ens_size, owners_index)) + ! Convert this observation ensemble from probit back to regular space if to_probit was successful + if(obs_probit_trans_ok(owners_index)) then + call transform_from_probit(ens_size, obs_ens_handle%copies(1:ens_size, owners_index) , & + obs_dist_params(owners_index), obs_ens_handle%copies(1:ens_size, owners_index)) + endif obs_prior = obs_ens_handle%copies(1:ens_size, owners_index) endif IF_QC_IS_OKAY @@ -698,10 +705,14 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Convert the prior and posterior for this observation to probit space call transform_to_probit(grp_size, obs_prior(grp_bot:grp_top), dist_for_obs, & temp_dist_params, probit_obs_prior(grp_bot:grp_top), .false., & - bounded_below, bounded_above, lower_bound, upper_bound) + bounded_below, bounded_above, lower_bound, upper_bound, ierr) + ! If probit transform fails for any group, this observation cannot be assimimlated + if(ierr /= 0) cycle SEQUENTIAL_OBS call transform_to_probit(grp_size, obs_post(grp_bot:grp_top), dist_for_obs, & temp_dist_params, probit_obs_post(grp_bot:grp_top), .true., & - bounded_below, bounded_above, lower_bound, upper_bound) + bounded_below, bounded_above, lower_bound, upper_bound, ierr) + ! If probit transform fails for any group, this observation cannot be assimimlated + if(ierr /= 0) cycle SEQUENTIAL_OBS ! Free up the storage used for this transform call deallocate_distribution_params(temp_dist_params) @@ -761,6 +772,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & STATE_UPDATE: do j = 1, num_close_states state_index = close_state_ind(j) + ! If transform to probit failed for this state variable, don't update it + if(.not. state_probit_trans_ok(state_index)) cycle STATE_UPDATE + if ( allow_missing_in_state ) then ! Don't allow update of state ensemble with any missing values if (any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)) cycle STATE_UPDATE @@ -796,6 +810,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & OBS_UPDATE: do j = 1, num_close_obs obs_index = close_obs_ind(j) + ! If transform to probit failed for this observed variable, don't update it + if(.not. obs_probit_trans_ok(obs_index)) cycle OBS_UPDATE + ! Only have to update obs that have not yet been used if(my_obs_indx(obs_index) > i) then @@ -820,7 +837,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & ! Do the inverse probit transform for state variables call transform_all_from_probit(ens_size, ens_handle%my_num_vars, ens_handle%copies, & - state_dist_params, ens_handle%copies) + state_dist_params, ens_handle%copies, state_probit_trans_ok) ! Every pe needs to get the current my_inflate and my_inflate_sd back if(local_single_ss_inflate) then @@ -876,13 +893,15 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & my_obs_type, & close_obs_ind, & vstatus, & - my_obs_loc) + my_obs_loc, & + state_probit_trans_ok) deallocate(close_state_dist, & my_state_indx, & close_state_ind, & my_state_kind, & - my_state_loc) + my_state_loc, & + obs_probit_trans_ok) ! end dealloc diff --git a/assimilation_code/modules/assimilation/filter_mod.f90 b/assimilation_code/modules/assimilation/filter_mod.f90 index 710091f33..57e137942 100644 --- a/assimilation_code/modules/assimilation/filter_mod.f90 +++ b/assimilation_code/modules/assimilation/filter_mod.f90 @@ -1561,7 +1561,7 @@ subroutine filter_ensemble_inflate(ens_handle, inflate_copy, inflate, ENS_MEAN_C real(r8) :: probit_ens(ens_size), probit_ens_mean logical :: bounded_below, bounded_above real(r8) :: lower_bound, upper_bound -integer :: dist_type +integer :: dist_type, ierr ! Inflate each group separately; Divide ensemble into num_groups groups grp_size = ens_size / num_groups @@ -1603,16 +1603,20 @@ subroutine filter_ensemble_inflate(ens_handle, inflate_copy, inflate, ENS_MEAN_C call transform_to_probit(grp_size, ens_handle%copies(grp_bot:grp_top, j), & dist_type, dist_params, probit_ens(1:grp_size), .false., & - bounded_below, bounded_above, lower_bound, upper_bound) - - ! Compute the ensemble mean in transformed space - probit_ens_mean = sum(probit_ens(1:grp_size)) / grp_size - ! Inflate in probit space - call inflate_ens(inflate, probit_ens(1:grp_size), probit_ens_mean, & - ens_handle%copies(inflate_copy, j)) - ! Transform back from probit space - call transform_from_probit(grp_size, probit_ens(1:grp_size), & - dist_params, ens_handle%copies(grp_bot:grp_top, j)) + bounded_below, bounded_above, lower_bound, upper_bound, ierr) + + ! Only inflate if successful return from the probit transform + if(ierr == 0) then + ! Compute the ensemble mean in transformed space + probit_ens_mean = sum(probit_ens(1:grp_size)) / grp_size + ! Inflate in probit space + call inflate_ens(inflate, probit_ens(1:grp_size), probit_ens_mean, & + ens_handle%copies(inflate_copy, j)) + ! Transform back from probit space + call transform_from_probit(grp_size, probit_ens(1:grp_size), & + dist_params, ens_handle%copies(grp_bot:grp_top, j)) + endif + end do endif end do diff --git a/assimilation_code/modules/assimilation/probit_transform_mod.f90 b/assimilation_code/modules/assimilation/probit_transform_mod.f90 index c5645c15a..984ca8bf2 100644 --- a/assimilation_code/modules/assimilation/probit_transform_mod.f90 +++ b/assimilation_code/modules/assimilation/probit_transform_mod.f90 @@ -62,7 +62,7 @@ module probit_transform_mod !------------------------------------------------------------------------ subroutine transform_all_to_probit(ens_size, num_vars, state_ens, distribution_type, & - p, probit_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound) + p, probit_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr) integer, intent(in) :: ens_size integer, intent(in) :: num_vars @@ -73,6 +73,7 @@ subroutine transform_all_to_probit(ens_size, num_vars, state_ens, distribution_t logical, intent(in) :: use_input_p logical, intent(in) :: bounded_below, bounded_above real(r8), intent(in) :: lower_bound, upper_bound +integer, intent(out) :: ierr(num_vars) ! NOTE THAT WILL MAKE HELEN CRAZY: THIS WORKS WITH THE INPUT CALLING ARGUMENTS FOR STATE_ENS AND @@ -89,8 +90,8 @@ subroutine transform_all_to_probit(ens_size, num_vars, state_ens, distribution_t do i = 1, num_vars call transform_to_probit(ens_size, state_ens(1:ens_size, i), distribution_type(i), & - p(i), temp_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound) - probit_ens(1:ens_size, i) = temp_ens + p(i), temp_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr(i)) + if(ierr(i) == 0) probit_ens(1:ens_size, i) = temp_ens end do end subroutine transform_all_to_probit @@ -98,7 +99,8 @@ end subroutine transform_all_to_probit !------------------------------------------------------------------------ subroutine transform_to_probit(ens_size, state_ens_in, distribution_type, p, & - probit_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound) + probit_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, & + ierr) integer, intent(in) :: ens_size real(r8), intent(in) :: state_ens_in(ens_size) @@ -108,6 +110,7 @@ subroutine transform_to_probit(ens_size, state_ens_in, distribution_type, p, & logical, intent(in) :: use_input_p logical, intent(in) :: bounded_below, bounded_above real(r8), intent(in) :: lower_bound, upper_bound +integer, intent(out) :: ierr real(r8) :: state_ens(ens_size) real(r8) :: probit_ens_temp(ens_size), state_ens_temp(ens_size), diff(ens_size) @@ -117,6 +120,9 @@ subroutine transform_to_probit(ens_size, state_ens_in, distribution_type, p, & ! If not initialized, read in the namelist if(.not. module_initialized) call initialize_probit_transform +! Default is no error, ierr is 0 +ierr = 0 + ! Fix bounds violations if requested if(fix_bound_violations) then do i = 1, ens_size @@ -145,7 +151,7 @@ subroutine transform_to_probit(ens_size, state_ens_in, distribution_type, p, & lower_bound, upper_bound) elseif(p%distribution_type == BOUNDED_NORMAL_RH_DISTRIBUTION) then call to_probit_bounded_normal_rh(ens_size, state_ens, p, probit_ens, & - use_input_p, bounded_below, bounded_above, lower_bound, upper_bound) + use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr) !---------------------------------------------------------------------------------- ! The following code block tests that the to/from probit calls are nearly inverse @@ -153,33 +159,37 @@ subroutine transform_to_probit(ens_size, state_ens_in, distribution_type, p, & if(do_inverse_check) then if(.not. use_input_p) then call to_probit_bounded_normal_rh(ens_size, state_ens, p_temp, probit_ens_temp, & - use_input_p, bounded_below, bounded_above, lower_bound, upper_bound) - call from_probit_bounded_normal_rh(ens_size, probit_ens_temp, p_temp, state_ens_temp) - diff = state_ens - state_ens_temp - if(abs(maxval(diff)) > 1.0e-8_r8) then - write(*, *) 'Maximum allowed value of probit to/from difference exceeded' - write(*, *) 'Location of minimum ensemble member ', minloc(state_ens) - write(*, *) 'Location of maximum ensemble member ', maxloc(state_ens) - do i = 1, ens_size - write(*, *) i, state_ens(i), state_ens_temp(i), diff(i) - enddo - stop + use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr) + if(ierr == 0) then + call from_probit_bounded_normal_rh(ens_size, probit_ens_temp, p_temp, state_ens_temp) + diff = state_ens - state_ens_temp + if(abs(maxval(diff)) > 1.0e-8_r8) then + write(*, *) 'Maximum allowed value of probit to/from difference exceeded' + write(*, *) 'Location of minimum ensemble member ', minloc(state_ens) + write(*, *) 'Location of maximum ensemble member ', maxloc(state_ens) + do i = 1, ens_size + write(*, *) i, state_ens(i), state_ens_temp(i), diff(i) + enddo + stop + endif endif endif if(use_input_p) then call to_probit_bounded_normal_rh(ens_size, state_ens, p, probit_ens_temp, & - use_input_p, bounded_below, bounded_above, lower_bound, upper_bound) - call from_probit_bounded_normal_rh(ens_size, probit_ens_temp, p, state_ens_temp) - diff = state_ens - state_ens_temp - if(abs(maxval(diff)) > 1.0e-8_r8) then - write(*, *) 'Maximum allowed value of probit to/from difference for input p exceeded' - write(*, *) 'Location of minimum ensemble member ', minloc(state_ens) - write(*, *) 'Location of maximum ensemble member ', maxloc(state_ens) - do i = 1, ens_size - write(*, *) i, state_ens(i), state_ens_temp(i), diff(i) - enddo - stop + use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr) + if(ierr == 0) then + call from_probit_bounded_normal_rh(ens_size, probit_ens_temp, p, state_ens_temp) + diff = state_ens - state_ens_temp + if(abs(maxval(diff)) > 1.0e-8_r8) then + write(*, *) 'Maximum allowed value of probit to/from difference for input p exceeded' + write(*, *) 'Location of minimum ensemble member ', minloc(state_ens) + write(*, *) 'Location of maximum ensemble member ', maxloc(state_ens) + do i = 1, ens_size + write(*, *) i, state_ens(i), state_ens_temp(i), diff(i) + enddo + stop + endif endif endif @@ -325,7 +335,7 @@ end subroutine to_probit_beta !------------------------------------------------------------------------ subroutine to_probit_bounded_normal_rh(ens_size, state_ens, p, probit_ens, & - use_input_p, bounded_below, bounded_above, lower_bound, upper_bound) + use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr) ! Note that this is just for transforming back and forth, not for doing the RHF observation update ! This means that we know a prior that the quantiles associated with the initial ensemble are @@ -334,6 +344,11 @@ subroutine to_probit_bounded_normal_rh(ens_size, state_ens, p, probit_ens, & ! How to handle identical ensemble members is an open question for now. This is also a problem ! for ensemble members that are identical to one of the bounds. +! If the stanard deviation computed for the sample ensemble is not positive, don't know how +! to do a transform? This can happen when the standard deviation computation is 0 or negative +! due to computational precision errors. Could try to work around this, but challenge with knowing +! what to do on the tails where a normal distribution requires a standard deviation. + integer, intent(in) :: ens_size real(r8), intent(in) :: state_ens(ens_size) type(distribution_params_type), intent(inout) :: p @@ -341,31 +356,31 @@ subroutine to_probit_bounded_normal_rh(ens_size, state_ens, p, probit_ens, & logical, intent(in) :: use_input_p logical, intent(in) :: bounded_below, bounded_above real(r8), intent(in) :: lower_bound, upper_bound +integer, intent(out) :: ierr ! Probit transform for bounded normal rh. integer :: i real(r8) :: quantile(ens_size) -if(use_input_p) then - ! Do not know what to do if sd of original ensemble is 0 (or small, work on this later) - if(get_bnrh_sd(p) <= 0.0_r8) then - ! Just return the original ensemble - probit_ens = state_ens - return - endif +! Successful return has ierr 0 +ierr = 0 +if(use_input_p) then + ! Need tthis to fail if p isn't available. Is that possible? ! Get the quantiles for each of the ensemble members in a BNRH distribution call bnrh_cdf_initialized_vector(state_ens, ens_size, p, quantile) - else ! Get all the info about the rank histogram cdf call bnrh_cdf_params(state_ens, ens_size, bounded_below, bounded_above, & lower_bound, upper_bound, p, quantile) - ! Do not know what to do if sd is 0 (or small, work on this later) + ! Fail if sd is not positive (or small, work on this later) if(get_bnrh_sd(p) <= 0.0_r8) then + ierr = 1 ! Just return the original ensemble probit_ens = state_ens + ! Free up the storage that would have been used for the transformed variables + call deallocate_distribution_params(p) return endif @@ -545,21 +560,28 @@ subroutine to_probit_kde(ens_size, state_ens, p, probit_ens, use_input_p, & end subroutine to_probit_kde -subroutine transform_all_from_probit(ens_size, num_vars, probit_ens, p, state_ens) +subroutine transform_all_from_probit(ens_size, num_vars, probit_ens, p, state_ens, & + transform_ok) integer, intent(in) :: ens_size integer, intent(in) :: num_vars real(r8), intent(in) :: probit_ens(:, :) type(distribution_params_type), intent(inout) :: p(num_vars) real(r8), intent(out) :: state_ens(:, :) +logical, intent(in) :: transform_ok(num_vars) ! Transform back to the original space integer :: i real(r8) :: temp_ens(ens_size) do i = 1, num_vars - call transform_from_probit(ens_size, probit_ens(1:ens_size, i), p(i), temp_ens) - state_ens(1:ens_size, i) = temp_ens + if(transform_ok(i)) then + call transform_from_probit(ens_size, probit_ens(1:ens_size, i), p(i), temp_ens) + state_ens(1:ens_size, i) = temp_ens + else + ! Transform can't be done, so return the input + state_ens(1:ens_size, i) = probit_ens(1:ens_size, i) + endif end do end subroutine transform_all_from_probit From 844163f13baba3b577e1c4b7b03357ff23216376 Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Mon, 9 Dec 2024 14:15:55 -0700 Subject: [PATCH 2/7] Corrected use of 1 in a loop over i that was checking to see if each variable had been successfully tranformed. --- assimilation_code/modules/assimilation/assim_tools_mod.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/assimilation_code/modules/assimilation/assim_tools_mod.f90 b/assimilation_code/modules/assimilation/assim_tools_mod.f90 index f5c7147ba..b489f738b 100644 --- a/assimilation_code/modules/assimilation/assim_tools_mod.f90 +++ b/assimilation_code/modules/assimilation/assim_tools_mod.f90 @@ -508,7 +508,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, & state_dist_params(i), probit_ens, .false., & bounded_below, bounded_above, lower_bound, upper_bound, ierr) state_probit_trans_ok(i) = (ierr == 0) - if(state_probit_trans_ok(1)) ens_handle%copies(1:ens_size, i) = probit_ens + if(state_probit_trans_ok(i)) ens_handle%copies(1:ens_size, i) = probit_ens end do !> optionally convert all state location verticals From 1f2b7e1f97945693dff69b30ba532891185e9fed Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Mon, 23 Dec 2024 07:32:37 -0700 Subject: [PATCH 3/7] Explicitly copy the input file into the output if the transform fails for the transform_all case. This version has been tested with the lorenz_96_tracer tests. All tests run to completion. All of the RMSE tests reproduce the baseline values except for: Case 1 ensemble size 6 and 15 Case 2 ensemble size 6 and 15 Case 3 ensemble size 6 and 10 --- .../modules/assimilation/probit_transform_mod.f90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/assimilation_code/modules/assimilation/probit_transform_mod.f90 b/assimilation_code/modules/assimilation/probit_transform_mod.f90 index 984ca8bf2..2b66b84fe 100644 --- a/assimilation_code/modules/assimilation/probit_transform_mod.f90 +++ b/assimilation_code/modules/assimilation/probit_transform_mod.f90 @@ -91,7 +91,12 @@ subroutine transform_all_to_probit(ens_size, num_vars, state_ens, distribution_t do i = 1, num_vars call transform_to_probit(ens_size, state_ens(1:ens_size, i), distribution_type(i), & p(i), temp_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr(i)) - if(ierr(i) == 0) probit_ens(1:ens_size, i) = temp_ens + if(ierr(i) == 0) + probit_ens(1:ens_size, i) = temp_ens + else + ! If transform failed, return the input state + probit_ens(1:ens_size, i) = state_ens(1:ens_size, i) + endif end do end subroutine transform_all_to_probit @@ -366,7 +371,7 @@ subroutine to_probit_bounded_normal_rh(ens_size, state_ens, p, probit_ens, & ierr = 0 if(use_input_p) then - ! Need tthis to fail if p isn't available. Is that possible? + ! Need this to fail if p isn't available. Is that possible? ! Get the quantiles for each of the ensemble members in a BNRH distribution call bnrh_cdf_initialized_vector(state_ens, ens_size, p, quantile) else From f0e7f99138df52600ef36caf0af1187ba81da83c Mon Sep 17 00:00:00 2001 From: Jeff Anderson Date: Mon, 23 Dec 2024 12:13:55 -0700 Subject: [PATCH 4/7] Corrected syntax error in probit_tranform_mod.f90. Reran tests with results the same as reported in previous commit. --- assimilation_code/modules/assimilation/probit_transform_mod.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/assimilation_code/modules/assimilation/probit_transform_mod.f90 b/assimilation_code/modules/assimilation/probit_transform_mod.f90 index 2b66b84fe..051b89f69 100644 --- a/assimilation_code/modules/assimilation/probit_transform_mod.f90 +++ b/assimilation_code/modules/assimilation/probit_transform_mod.f90 @@ -91,7 +91,7 @@ subroutine transform_all_to_probit(ens_size, num_vars, state_ens, distribution_t do i = 1, num_vars call transform_to_probit(ens_size, state_ens(1:ens_size, i), distribution_type(i), & p(i), temp_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr(i)) - if(ierr(i) == 0) + if(ierr(i) == 0) then probit_ens(1:ens_size, i) = temp_ens else ! If transform failed, return the input state From 894c9f2d3ad9e217b8d2ac355f8823f824b756d7 Mon Sep 17 00:00:00 2001 From: Helen Kershaw Date: Mon, 27 Jan 2025 10:32:33 -0500 Subject: [PATCH 5/7] chore: remove unused routine --- .../assimilation/probit_transform_mod.f90 | 44 +------------------ 1 file changed, 1 insertion(+), 43 deletions(-) diff --git a/assimilation_code/modules/assimilation/probit_transform_mod.f90 b/assimilation_code/modules/assimilation/probit_transform_mod.f90 index 9e65dcb4f..215362fa7 100644 --- a/assimilation_code/modules/assimilation/probit_transform_mod.f90 +++ b/assimilation_code/modules/assimilation/probit_transform_mod.f90 @@ -37,7 +37,7 @@ module probit_transform_mod implicit none private -public :: transform_to_probit, transform_from_probit, transform_all_to_probit, & +public :: transform_to_probit, transform_from_probit, & transform_all_from_probit character(len=512) :: errstring @@ -61,48 +61,6 @@ module probit_transform_mod !------------------------------------------------------------------------ -subroutine transform_all_to_probit(ens_size, num_vars, state_ens, distribution_type, & - p, probit_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr) - -integer, intent(in) :: ens_size -integer, intent(in) :: num_vars -real(r8), intent(in) :: state_ens(:, :) -integer, intent(in) :: distribution_type(num_vars) -type(distribution_params_type), intent(inout) :: p(num_vars) -real(r8), intent(out) :: probit_ens(:, :) -logical, intent(in) :: use_input_p -logical, intent(in) :: bounded_below, bounded_above -real(r8), intent(in) :: lower_bound, upper_bound -integer, intent(out) :: ierr(num_vars) - - -! NOTE THAT WILL MAKE HELEN CRAZY: THIS WORKS WITH THE INPUT CALLING ARGUMENTS FOR STATE_ENS AND -! PROBIT_ENS BEING THE SAME. A TEMP IS USED TO AVOID OVERWRITING ISSUES. IS THIS YUCKY? - -! Note that the input and output arrays may have extra copies (first subscript). Passing sections of a -! leading index could be inefficient for time and storage, so avoiding that for now. - -! Assumes that the bounds are the same for any variables that are BNRH for now -! The bounds variables are not used for the normal case or the case where the input p is used - -integer :: i -real(r8) :: temp_ens(ens_size) - -do i = 1, num_vars - call transform_to_probit(ens_size, state_ens(1:ens_size, i), distribution_type(i), & - p(i), temp_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, ierr(i)) - if(ierr(i) == 0) then - probit_ens(1:ens_size, i) = temp_ens - else - ! If transform failed, return the input state - probit_ens(1:ens_size, i) = state_ens(1:ens_size, i) - endif -end do - -end subroutine transform_all_to_probit - -!------------------------------------------------------------------------ - subroutine transform_to_probit(ens_size, state_ens_in, distribution_type, p, & probit_ens, use_input_p, bounded_below, bounded_above, lower_bound, upper_bound, & ierr) From 2273740ae63bb7f441aa080b48af947c12279f02 Mon Sep 17 00:00:00 2001 From: Helen Kershaw Date: Mon, 27 Jan 2025 12:12:38 -0500 Subject: [PATCH 6/7] fix: use error_handler to terminate program rather than stop ALLMSG so tasks with the greater than the max difference print the the information. --- .../assimilation/probit_transform_mod.f90 | 27 +++++++++++-------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/assimilation_code/modules/assimilation/probit_transform_mod.f90 b/assimilation_code/modules/assimilation/probit_transform_mod.f90 index 215362fa7..a1bf8f456 100644 --- a/assimilation_code/modules/assimilation/probit_transform_mod.f90 +++ b/assimilation_code/modules/assimilation/probit_transform_mod.f90 @@ -11,7 +11,7 @@ module probit_transform_mod use sort_mod, only : index_sort -use utilities_mod, only : E_ERR, error_handler, do_nml_file, do_nml_term, nmlfileunit, & +use utilities_mod, only : E_ERR, E_ALLMSG, error_handler, do_nml_file, do_nml_term, nmlfileunit, & find_namelist_in_file, check_namelist_read use distribution_params_mod, only : distribution_params_type, deallocate_distribution_params, & @@ -79,6 +79,7 @@ subroutine transform_to_probit(ens_size, state_ens_in, distribution_type, p, & real(r8) :: probit_ens_temp(ens_size), state_ens_temp(ens_size), diff(ens_size) type(distribution_params_type) :: p_temp integer :: i +character(len=32), parameter :: routine = 'transform_to_probit' ! If not initialized, read in the namelist if(.not. module_initialized) call initialize_probit_transform @@ -125,13 +126,15 @@ subroutine transform_to_probit(ens_size, state_ens_in, distribution_type, p, & call from_probit_bounded_normal_rh(ens_size, probit_ens_temp, p_temp, state_ens_temp) diff = state_ens - state_ens_temp if(abs(maxval(diff)) > 1.0e-8_r8) then - write(*, *) 'Maximum allowed value of probit to/from difference exceeded' - write(*, *) 'Location of minimum ensemble member ', minloc(state_ens) - write(*, *) 'Location of maximum ensemble member ', maxloc(state_ens) + write(errstring, *) 'Location of minimum ensemble member ', minloc(state_ens) + call error_handler(E_ALLMSG, routine, errstring) + write(errstring, *) 'Location of maximum ensemble member ', maxloc(state_ens) + call error_handler(E_ALLMSG, routine, errstring) do i = 1, ens_size - write(*, *) i, state_ens(i), state_ens_temp(i), diff(i) + write(errstring, *) i, state_ens(i), state_ens_temp(i), diff(i) + call error_handler(E_ALLMSG, routine, errstring) enddo - stop + call error_handler(E_ERR, routine, 'Maximum allowed value of probit to/from difference exceeded') endif endif endif @@ -143,13 +146,15 @@ subroutine transform_to_probit(ens_size, state_ens_in, distribution_type, p, & call from_probit_bounded_normal_rh(ens_size, probit_ens_temp, p, state_ens_temp) diff = state_ens - state_ens_temp if(abs(maxval(diff)) > 1.0e-8_r8) then - write(*, *) 'Maximum allowed value of probit to/from difference for input p exceeded' - write(*, *) 'Location of minimum ensemble member ', minloc(state_ens) - write(*, *) 'Location of maximum ensemble member ', maxloc(state_ens) + write(errstring, *) 'Location of minimum ensemble member ', minloc(state_ens) + call error_handler(E_ALLMSG, routine, errstring) + write(errstring, *) 'Location of maximum ensemble member ', maxloc(state_ens) + call error_handler(E_ALLMSG, routine, errstring) do i = 1, ens_size - write(*, *) i, state_ens(i), state_ens_temp(i), diff(i) + write(errstring, *) i, state_ens(i), state_ens_temp(i), diff(i) + call error_handler(E_ALLMSG, routine, errstring) enddo - stop + call error_handler(E_ERR, routine, 'Maximum allowed value of probit to/from difference for input p exceeded') endif endif From 437677f4d10c79165050364f506fdd540901c82f Mon Sep 17 00:00:00 2001 From: Helen Kershaw Date: Thu, 30 Jan 2025 15:32:45 -0500 Subject: [PATCH 7/7] bump vesion and changelog for release --- CHANGELOG.rst | 5 +++++ conf.py | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 64b31d3ec..537531f11 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,11 @@ individual files. The changes are now listed with the most recent at the top. +**January 30 2025 :: Bug-fix: Explicitly handle BNRHF transform failures. Tag v11.10.1** + +- Probit transform failure is caught and an error code is returned +- filter_mod and assim_tools_mod skip variables that fail the transform + **January 23 2025 :: DART_LAB QCEFF. Tag v11.10.0** - Updated DART_LAB to include QCEFF diff --git a/conf.py b/conf.py index f9724da98..87389e14e 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '11.10.0' +release = '11.10.1' root_doc = 'index' # -- General configuration ---------------------------------------------------