[Dart-dev] [4869] DART/trunk/assim_tools/assim_tools_mod.f90: Several critical fixes for those users who redefine R8 to be R4.
nancy at ucar.edu
nancy at ucar.edu
Tue Apr 19 16:57:13 MDT 2011
Revision: 4869
Author: nancy
Date: 2011-04-19 16:57:13 -0600 (Tue, 19 Apr 2011)
Log Message:
-----------
Several critical fixes for those users who redefine R8 to be R4.
1. Compute variance and covariance in more numerically stable ways
if running with single precision reals. The previous code was fine
when running with real*8, but the updated code retains more accuracy
when running real*4.
2. Fix minor issues with the Sampling Error Correction; a few of
the
Modified Paths:
--------------
DART/trunk/assim_tools/assim_tools_mod.f90
-------------- next part --------------
Modified: DART/trunk/assim_tools/assim_tools_mod.f90
===================================================================
--- DART/trunk/assim_tools/assim_tools_mod.f90 2011-04-19 22:40:49 UTC (rev 4868)
+++ DART/trunk/assim_tools/assim_tools_mod.f90 2011-04-19 22:57:13 UTC (rev 4869)
@@ -236,7 +236,7 @@
write(errstring, '(A,I10,A)') 'Using adaptive localization, threshold ', &
adaptive_localization_threshold, ' obs'
call error_handler(E_MSG,'assim_tools_init:', errstring)
- if(adaptive_cutoff_floor > 0) then
+ if(adaptive_cutoff_floor > 0.0_r8) then
write(errstring, '(A,F18.6)') 'Minimum cutoff will not go below ', &
adaptive_cutoff_floor
call error_handler(E_MSG,'assim_tools_init:', 'Using adaptive localization cutoff floor.', &
@@ -248,6 +248,10 @@
call error_handler(E_MSG,'assim_tools_init:', 'Writing localization diagnostics to file:')
call error_handler(E_MSG,'assim_tools_init:', trim(localization_diagnostics_file))
endif
+
+ if(sampling_error_correction) then
+ call error_handler(E_MSG,'assim_tools_init:', 'Using Sampling Error Correction')
+ endif
endif
end subroutine assim_tools_init
@@ -580,8 +584,9 @@
grp_bot = grp_beg(group)
grp_top = grp_end(group)
obs_prior_mean(group) = sum(obs_prior(grp_bot:grp_top)) / grp_size
- obs_prior_var(group) = sum(obs_prior(grp_bot:grp_top) * obs_prior(grp_bot:grp_top)) - &
- grp_size * obs_prior_mean(group)**2
+ obs_prior_var(group) = sum((obs_prior(grp_bot:grp_top) - obs_prior_mean(group))**2) / &
+ (grp_size - 1)
+ if (obs_prior_var(group) < 0.0_r8) obs_prior_var(group) = 0.0_r8
end do
! Need to get obs density first in case of adaptive localization
@@ -738,7 +743,7 @@
obs_prior_var(group), obs_inc(grp_bot:grp_top), &
ens_handle%copies(grp_bot:grp_top, state_index), grp_size, &
increment(grp_bot:grp_top), reg_coef(group), net_a(group), correl(group))
- else
+ else
call update_from_obs_inc(obs_prior(grp_bot:grp_top), obs_prior_mean(group), &
obs_prior_var(group), obs_inc(grp_bot:grp_top), &
ens_handle%copies(grp_bot:grp_top, state_index), grp_size, &
@@ -888,7 +893,7 @@
num_close_obs_buffered + num_close_states_buffered
if (num_close_obs_buffered+num_close_obs_calls_made+ &
num_close_states_buffered+num_close_states_calls_made > 0) then
- print *, "Percent saved: ", 100. * &
+ print *, "Percent saved: ", 100.0_r8 * &
(real(num_close_obs_buffered+num_close_states_buffered, r8) / &
(num_close_obs_calls_made+num_close_obs_buffered+ &
num_close_states_calls_made+num_close_states_buffered))
@@ -1422,15 +1427,15 @@
real(r8), intent(inout) :: net_a
real(r8), optional, intent(inout) :: correl_out
-real(r8) :: t(ens_size), obs_state_cov
+real(r8) :: obs_state_cov
real(r8) :: restoration_inc(ens_size), state_mean, state_var, correl
real(r8) :: factor, exp_true_correl, mean_factor
! For efficiency, just compute regression coefficient here unless correl is needed
-t = obs - obs_prior_mean
-obs_state_cov = sum(t * state)
+state_mean = sum(state) / ens_size
+obs_state_cov = sum( (state - state_mean) * (obs - obs_prior_mean) ) / (ens_size - 1)
-if (obs_prior_var /= 0.0_r8) then
+if (obs_prior_var > 0.0_r8) then
reg_coef = obs_state_cov/obs_prior_var
else
reg_coef = 0.0_r8
@@ -1439,8 +1444,9 @@
! If correl_out is present, need correl for adaptive inflation
! Also needed for file correction below
if(present(correl_out) .or. sampling_error_correction) then
- state_var = sum(state * state) - sum(state)**2 / ens_size
- if(obs_prior_var * state_var <= 0.0_r8) then
+ state_var = sum((state - state_mean)**2) / (ens_size - 1)
+ if (state_var < 0.0_r8) state_var = 0.0_r8
+ if ((obs_prior_var <= 0.0_r8) .or. (state_var <= 0.0_r8)) then
correl = 0.0_r8
else
correl = obs_state_cov / sqrt(obs_prior_var * state_var)
@@ -1455,7 +1461,6 @@
! Get the expected actual correlation and the regression weight reduction factor
if(sampling_error_correction) then
call get_correction_from_file(ens_size, correl, mean_factor, exp_true_correl)
- !write(*, *) correl, exp_true_correl, mean_factor
! Watch out for division by zero; if correl is really small regression is safely 0
if(abs(correl) > 0.001) then
reg_coef = reg_coef * (exp_true_correl / correl) * mean_factor
@@ -1488,7 +1493,7 @@
! the reciprocal of the ensemble_size (slope = 0.80 / ens_size). These are empirical
! for now. See also README in spread_restoration_paper documentation.
!!!factor = 1.0_r8 / (1.0_r8 + (net_a - 1.0_r8) * (0.8_r8 / ens_size)) - 1.0_r8
- factor = 1.0_r8 / (1.0_r8 + (net_a - 1.0_r8) / (-2.4711 + 1.6386 * ens_size)) - 1.0_r8
+ factor = 1.0_r8 / (1.0_r8 + (net_a - 1.0_r8) / (-2.4711_r8 + 1.6386_r8 * ens_size)) - 1.0_r8
!!!factor = 1.0_r8 / (1.0_r8 + (net_a**2 - 1.0_r8) * (-0.0111_r8 + .8585_r8 / ens_size)) - 1.0_r8
! Variance restoration
@@ -1497,6 +1502,15 @@
state_inc = state_inc + restoration_inc
endif
+!! NOTE: if requested to be returned, correl_out is set further up in the
+!! code, before the sampling error correction, if enabled, is applied.
+!! this means it's returning a different larger value than the correl
+!! being returned here. it's used by the adaptive inflation and so the
+!! inflation will see a slightly different correlation value. it isn't
+!! clear that this is a bad thing; it means the inflation might be a bit
+!! larger than it would otherwise. before we move any code this would
+!! need to be studied to see what the real impact would be.
+
end subroutine update_from_obs_inc
@@ -1554,8 +1568,7 @@
first_get_correction = .false.
endif
-
-! First quick modification of correl to expected true correl for test (should interp)
+! Interpolate to get values of expected correlation and mean_factor
if(scorrel < -1.0_r8) then
correl = -1.0_r8
mean_factor = 1.0_r8
@@ -1563,29 +1576,47 @@
correl = 1.0_r8
mean_factor = 1.0_r8
else if(scorrel <= -0.995_r8) then
- fract = (scorrel + 1.0_r8) / 0.05_r8
+ fract = (scorrel + 1.0_r8) / 0.005_r8
correl = (exp_true_correl(1) + 1.0_r8) * fract - 1.0_r8
mean_factor = (alpha(1) - 1.0_r8) * fract + 1.0_r8
else if(scorrel >= 0.995_r8) then
- fract = (scorrel - 0.995_r8) / 0.05_r8
+ fract = (scorrel - 0.995_r8) / 0.005_r8
correl = (1.0_r8 - exp_true_correl(200)) * fract + exp_true_correl(200)
mean_factor = (1.0_r8 - alpha(200)) * fract + alpha(200)
else
+ ! given the ifs above, the floor() computation below for low_indx
+ ! should always result in a value in the range 1 to 199. but if this
+ ! code is compiled with r8=r4 (single precision reals) it turns out
+ ! to be possible to get values a few bits below 0 which results in
+ ! a very large negative integer. the limit tests below ensure the
+ ! index stays in a legal range.
low_indx = floor((scorrel + 0.995_r8) / 0.01_r8 + 1.0_r8)
- low_correl = -0.995_r8 + (low_indx - 1) * 0.01
+ if (low_indx < 1) low_indx = 1
+ if (low_indx > 199) low_indx = 199
+ low_correl = -0.995_r8 + (low_indx - 1) * 0.01_r8
low_exp_correl = exp_true_correl(low_indx)
low_alpha = alpha(low_indx)
high_indx = low_indx + 1
- high_correl = low_correl + 0.01
+ high_correl = low_correl + 0.01_r8
high_exp_correl = exp_true_correl(high_indx)
- high_alpha = alpha(low_indx)
+ high_alpha = alpha(high_indx)
fract = (scorrel - low_correl) / (high_correl - low_correl)
correl = (high_exp_correl - low_exp_correl) * fract + low_exp_correl
mean_factor = (high_alpha - low_alpha) * fract + low_alpha
endif
-expected_true_correl = correl
+expected_true_correl = correl
+! Don't want Monte Carlo interpolation problems to put us outside of a
+! ratio between 0 and 1 for expected_true_correl / sample_correl
+! If they have different signs, expected should just be 0
+if(expected_true_correl * scorrel <= 0.0_r8) then
+ expected_true_correl = 0.0_r8
+else if(abs(expected_true_correl) > abs(scorrel)) then
+ ! If same sign, expected should not be bigger in absolute value
+ expected_true_correl = scorrel
+endif
+
end subroutine get_correction_from_file
More information about the Dart-dev
mailing list