Skip to content

Commit

Permalink
Removed an amplitude adjustment factor from inv_bnrh_cdf. This was a …
Browse files Browse the repository at this point in the history
…legacy bit of code

that was required to bitwise reproduce earlier versions. About half of ensemble sizes
continue to bitwise reproduce with this change using Jeff's lorenz96_tracer_mod tests.
For the others, the amp_adj was different from 1 at the smallest writeable decimal
point and the tests ran producing reasonable, but no longer backward compatible,
results.
  • Loading branch information
jlaucar committed Dec 18, 2023
1 parent 4fd58ef commit 4e27dcf
Showing 1 changed file with 7 additions and 13 deletions.
20 changes: 7 additions & 13 deletions assimilation_code/modules/assimilation/bnrh_distribution_mod.f90
Original file line number Diff line number Diff line change
Expand Up @@ -409,24 +409,19 @@ subroutine inv_bnrh_cdf(quantiles, ens_size, sort_ens, &
x(i) = lower_bound + (curr_q / q(1)) * (upper_state - lower_bound)
else
! Find the mass at the lower bound (which could be unbounded)
! NOTE: The amplitude here should be one since there is no likelihood. However, there is
! round-off error that occurs in the statement below. In the long term, amp_adj should be
! removed from this code block and the onn for the upper region. However, removing it now
! would require resetting the baseline for the large number of baseline archived experiments.
amp_adj = q(1) / del_q
if(bounded_below) then
lower_mass = amp_adj * tail_amp_left * &
lower_mass = tail_amp_left * &
normal_cdf(lower_bound, tail_mean_left, tail_sd_left)
else
lower_mass = 0.0_r8
endif
! Find the mass at the upper bound (ensemble member 1)
upper_mass = amp_adj * tail_amp_left * &
upper_mass = tail_amp_left * &
normal_cdf(sort_ens(1), tail_mean_left, tail_sd_left)
! What fraction of this mass difference should we go?
fract = curr_q / q(1)
target_mass = lower_mass + fract * (upper_mass - lower_mass)
x(i) = inv_weighted_normal_cdf(amp_adj*tail_amp_left, tail_mean_left, &
x(i) = inv_weighted_normal_cdf(tail_amp_left, tail_mean_left, &
tail_sd_left, target_mass)
endif

Expand All @@ -441,20 +436,19 @@ subroutine inv_bnrh_cdf(quantiles, ens_size, sort_ens, &
else
! Upper tail is (bounded) normal
! Find the mass at the upper bound (which could be unbounded)
amp_adj = (1.0_r8 - q(ens_size)) / del_q
if(bounded_above) then
upper_mass = amp_adj * tail_amp_right * &
upper_mass = tail_amp_right * &
normal_cdf(upper_bound, tail_mean_right, tail_sd_right)
else
upper_mass = amp_adj * 1.0_r8
upper_mass = 1.0_r8
endif
! Find the mass at the lower edge of the region (ensemble member n)
lower_mass = amp_adj * tail_amp_right * &
lower_mass = tail_amp_right * &
normal_cdf(sort_ens(ens_size), tail_mean_right, tail_sd_right)
! What fraction of the last interval do we need to move
fract = (curr_q - q(ens_size)) / (1.0_r8 - q(ens_size))
target_mass = lower_mass + fract * (upper_mass - lower_mass)
x(i) = inv_weighted_normal_cdf(amp_adj * tail_amp_right, tail_mean_right, &
x(i) = inv_weighted_normal_cdf(tail_amp_right, tail_mean_right, &
tail_sd_right, target_mass)
endif

Expand Down

0 comments on commit 4e27dcf

Please sign in to comment.