[Dart-dev] [4860] DART/trunk/assim_tools/assim_tools_mod.f90: proposed update to assim_tools_mod. differen
nancy at ucar.edu
nancy at ucar.edu
Fri Apr 15 13:19:53 MDT 2011
Revision: 4860
Author: nancy
Date: 2011-04-15 13:19:53 -0600 (Fri, 15 Apr 2011)
Log Message:
-----------
proposed update to assim_tools_mod. different method to compute variance & covariance; _r8 on constants, additional limit tests for values
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-15 19:16:58 UTC (rev 4859)
+++ DART/trunk/assim_tools/assim_tools_mod.f90 2011-04-15 19:19:53 UTC (rev 4860)
@@ -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,42 @@
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
- low_indx = floor((scorrel + 0.995_r8) / 0.01_r8 + 1.0_r8)
- low_correl = -0.995_r8 + (low_indx - 1) * 0.01
+ !low_indx = floor((scorrel + 0.995_r8) / 0.01_r8 + 1.0_r8)
+ !if(low_indx < 1) low_indx = 1
+ !if(low_indx > 199) low_indx = 199
+ low_indx = min(floor((scorrel + 0.995_r8) / 0.01_r8 + 1.0_r8), 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