[Dart-dev] DART/branches Revision: 10959

[dart at ucar.edu]

thoar at ucar.edu
2017-01-26 11:17:32 -0700 (Thu, 26 Jan 2017)
250
The sampling error correction code has been vetted
and an obsolete test program removed.
Instead of using 
system_simulation/final_full_precomputed_tables/final_full.xx
all the ensemble sizes will be available in 
sampling_error_correction_table.nc




Modified: DART/branches/rma_trunk/system_simulation/gen_sampling_err_table.f90
===================================================================
--- DART/branches/rma_trunk/system_simulation/gen_sampling_err_table.f90	2017-01-26 18:09:52 UTC (rev 10958)
+++ DART/branches/rma_trunk/system_simulation/gen_sampling_err_table.f90	2017-01-26 18:17:32 UTC (rev 10959)
@@ -5,9 +5,9 @@
 ! $Id$
 
 !> Correct covariances for fixed ensemble sizes.
-!> See Anderson, J. L., 2011: Localization and Sampling Error Correction
-!>   in Ensemble Kalman Filter Data Assimilation. 
-!> Submitted for publication, Jan 2011.  Contact author.
+!> Ref: Anderson, J., 2012: 
+!> Localization and Sampling Error Correction in Ensemble Kalman Filter Data Assimilation.
+!> Mon. Wea. Rev., 140, 2359-2371, doi: 10.1175/MWR-D-11-00013.1. 
 
 !> this version is a sparse array - the ens_index array lists the ensemble
 !> sizes that have been computed, using the unlimited dimension.  to generate
@@ -28,8 +28,7 @@
                            initialize_utilities, finalize_utilities, &
                            find_namelist_in_file, check_namelist_read, &
                            do_nml_file, do_nml_term, nmlfileunit
-use random_seq_mod, only : random_seq_type, init_random_seq, random_gaussian, &          
-                           twod_gaussians, random_uniform
+use random_seq_mod, only : random_seq_type, init_random_seq, twod_gaussians
 
 use netcdf
 
@@ -43,14 +42,9 @@
 
 
 integer, parameter :: num_times   = 1
-! FIXME BEFORE COMMITTING, 1 of 2:
-integer, parameter :: num_samples = 100000  ! for developement/debugging
-!integer, parameter :: num_samples = 100000000   ! real value which should be used
+integer, parameter :: num_samples = 100000000 ! large number for statistical rigor
 
-! FIXME BEFORE COMMITTING, 2 of 2:
-!character(len=128) :: output_filename = 'sampling_error_correction_table.nc'
-! for debugging, make this short:
-character(len=128) :: output_filename = 'sec.nc'
+character(len=128) :: output_filename = 'sampling_error_correction_table.nc'
 integer, parameter :: nentries = 200
 
 integer, parameter :: MAX_LIST_LEN = 500
@@ -130,10 +124,10 @@
 contains
 
 !----------------------------------------------------------------
-! main computational routine
-!----------------------------------------------------------------
+!> main computational routine
 
 subroutine compute_table(this_size, nentries, bin_count, true_correl_mean, alpha)
+
 integer,  intent(in)  :: this_size
 integer,  intent(in)  :: nentries
 integer,  intent(out) :: bin_count(0:nentries+1)
@@ -140,10 +134,9 @@
 real(r8), intent(out) :: true_correl_mean(nentries)
 real(r8), intent(out) :: alpha(nentries)
 
-
 real(r8), allocatable  :: pairs(:,:), temp(:)
 
-real(r8) :: zero_2(2) = 0.0, cov(2, 2)
+real(r8) :: zero_2(2) = 0.0_r8, cov(2, 2)
 real(r8) :: t_correl, sample_correl, beta
 real(r8) :: s_mean(2), s_var(2), reg_mean, reg_sd, t_sd1, t_sd2
 real(r8) :: tcorrel_sum(nentries), reg_sum(nentries), reg_2_sum(nentries)
@@ -199,6 +192,7 @@
 
       ! Keep a sum of the sample_correl and squared to compute standard deviation
       tcorrel_sum(bin_num) = tcorrel_sum(bin_num) + t_correl
+
       !-----------------
       ! Also interested in finding out what the spurious variance reduction factor is
       ! First need to compute sample s.d. for the obs and unobs variable
@@ -258,8 +252,7 @@
 end subroutine
 
 !----------------------------------------------------------------
-! stat routines
-!----------------------------------------------------------------
+!> stat routines
 
 subroutine comp_correl(ens, n, correl)
 
@@ -269,7 +262,6 @@
 
 real(r8) :: sum_x, sum_y, sum_xy, sum_x2, sum_y2
 
-
 sum_x  = sum(ens(2, :))
 sum_y  = sum(ens(1, :))
 sum_xy = sum(ens(2, :) * ens(1, :))
@@ -284,6 +276,7 @@
 end subroutine comp_correl
 
 !----------------------------------------------------------------
+!>