[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]

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