[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