[Dart-dev] [5576] DART/trunk/assim_tools/assim_tools_mod.f90: Patch for the computation of correlation for inflation/sampling error
nancy at ucar.edu
nancy at ucar.edu
Wed Feb 29 11:49:12 MST 2012
Revision: 5576
Author: nancy
Date: 2012-02-29 11:49:12 -0700 (Wed, 29 Feb 2012)
Log Message:
-----------
Patch for the computation of correlation for inflation/sampling error
when running with R4 precision. With very small values it was possible
for both A and B to test larger than 0, but A*B was small enough to
underflow and be 0. The original code was doing the tests separately
and not testing the product before dividing by it. With very small vals
it was possible to divide by 0 and get NaN or Inf (depending on the numerator).
Inf would get reset to 1.0 in the lines below but NaN would remain as
NaN and percolate through the code to the output state vector unless you
were running with floating point exceptions enabled (which usually carries
a performance penalty).
This revised code checks the values in a slightly different order to
avoid computation if we are going to be setting the correlation to 0
anyway, and it is more careful to check for 0s in the face of underflow
with R4 precision. There is still a risk of overflow in this code with
very large data values; we do not recommend running with reals compiled
R4 if your data includes values with very large (eg 1e15+) 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 2012-02-28 23:11:33 UTC (rev 5575)
+++ DART/trunk/assim_tools/assim_tools_mod.f90 2012-02-29 18:49:12 UTC (rev 5576)
@@ -1429,7 +1429,7 @@
real(r8), intent(inout) :: net_a
real(r8), optional, intent(inout) :: correl_out
-real(r8) :: obs_state_cov
+real(r8) :: obs_state_cov, intermed
real(r8) :: restoration_inc(ens_size), state_mean, state_var, correl
real(r8) :: factor, exp_true_correl, mean_factor
@@ -1444,14 +1444,28 @@
endif
! If correl_out is present, need correl for adaptive inflation
-! Also needed for file correction below
+! Also needed for file correction below.
+
+! WARNING: we have had several different numerical problems in this
+! section, especially with users running in single precision floating point.
+! Be very cautious if changing any code in this section, taking into
+! account underflow and overflow for 32 bit floats.
+
if(present(correl_out) .or. sampling_error_correction) 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
+ if (obs_state_cov == 0.0_r8 .or. obs_prior_var <= 0.0_r8) then
correl = 0.0_r8
else
- correl = obs_state_cov / sqrt(obs_prior_var * state_var)
+ state_var = sum((state - state_mean)**2) / (ens_size - 1)
+ if (state_var <= 0.0_r8) then
+ correl = 0.0_r8
+ else
+ intermed = sqrt(obs_prior_var) * sqrt(state_var)
+ if (intermed <= 0.0_r8) then
+ correl = 0.0_r8
+ else
+ correl = obs_state_cov / intermed
+ endif
+ endif
endif
if(correl > 1.0_r8) correl = 1.0_r8
if(correl < -1.0_r8) correl = -1.0_r8
More information about the Dart-dev
mailing list