[Dart-dev] DART/branches Revision: 10757

[dart at ucar.edu]

jamesmcc at ucar.edu
2016-11-16 15:33:46 -0700 (Wed, 16 Nov 2016)
37
Moving to Poterjoy particle filter. 



Modified: DART/branches/wrfHydro/assim_tools/assim_tools_mod.f90
===================================================================
--- DART/branches/wrfHydro/assim_tools/assim_tools_mod.f90	2016-11-16 22:09:25 UTC (rev 10756)
+++ DART/branches/wrfHydro/assim_tools/assim_tools_mod.f90	2016-11-16 22:33:46 UTC (rev 10757)
@@ -8,7 +8,6 @@
 
  
 ! A variety of operations required by assimilation.
-
 use      types_mod,       only : r8, digits12, PI, missing_r8
 use  utilities_mod,       only : file_exist, get_unit, check_namelist_read, do_output,    &
                                  find_namelist_in_file, register_module, error_handler,   &
@@ -41,7 +40,7 @@
                                  prepare_to_update_copies, map_pe_to_task
 
 use mpi_utilities_mod,    only : my_task_id, broadcast_send, broadcast_recv,              & 
-                                 sum_across_tasks
+                                 sum_across_tasks, task_sync
 
 use adaptive_inflate_mod, only : do_obs_inflate,  do_single_ss_inflate,                   &
                                  do_varying_ss_inflate, get_inflate, set_inflate,         &
@@ -73,7 +72,6 @@
 real(r8), allocatable  :: cutoff_list(:)
 logical                :: has_special_cutoffs
 logical                :: close_obs_caching = .true.
-real(r8), parameter    :: small = epsilon(1.0_r8)   ! threshold for avoiding NaNs/Inf
 
 character(len = 255)   :: msgstring, msgstring2, msgstring3
 
@@ -99,10 +97,12 @@
 !      5 = random draw from posterior
 !      6 = deterministic draw from posterior with fixed kurtosis
 !      8 = Rank Histogram Filter (see Anderson 2011)
+!      9 = Localized particle filter (Poterjoy Nov. 2014)
 !
 !  special_localization_obs_types -> Special treatment for the specified observation types
 !  special_localization_cutoffs   -> Different cutoff value for each specified obs type
 !
+! JPOTERJOY: added namelist variables frac_neff, pf_kddm, kddm_bins
 integer  :: filter_kind                     = 1
 real(r8) :: cutoff                          = 0.2_r8
 logical  :: sort_obs_inc                    = .false.
@@ -111,6 +111,10 @@
 integer  :: adaptive_localization_threshold = -1
 real(r8) :: adaptive_cutoff_floor           = 0.0_r8
 integer  :: print_every_nth_obs             = 0
+real(r8) :: frac_neff                       = 0.30_r8
+integer  :: pf_kddm                         = 1
+integer  :: kddm_bins                       = 5
+real(r8), parameter ::             max_infl = 99999.0_r8
 
 ! since this is in the namelist, it has to have a fixed size.
 integer, parameter   :: MAX_ITEMS = 300
@@ -135,10 +139,12 @@
 ! sections.  to try out the alternatives, set this to .false.
 logical  :: only_area_adapt  = .true.
 
+! JPOTERJOY: new namelist variables
 namelist / assim_tools_nml / filter_kind, cutoff, sort_obs_inc, &
    spread_restoration, sampling_error_correction,                          & 
    adaptive_localization_threshold, adaptive_cutoff_floor,                 &
    print_every_nth_obs, rectangular_quadrature, gaussian_likelihood_tails, &
+   pf_kddm, frac_neff, kddm_bins,                                          &
    output_localization_diagnostics, localization_diagnostics_file,         &
    special_localization_obs_types, special_localization_cutoffs,           &
    allow_missing_in_clm
@@ -207,7 +213,7 @@
                          text3=trim(msgstring))
    endif
    j = j + 1
-enddo
+end do
 num_special_cutoff = j
 
 if (num_special_cutoff > 0) has_special_cutoffs = .true.
@@ -314,8 +320,17 @@
 real(r8) :: close_state_dist(ens_handle%my_num_vars)
 real(r8) :: last_close_obs_dist(obs_ens_handle%my_num_vars)
 real(r8) :: last_close_state_dist(ens_handle%my_num_vars)
-real(r8) :: diff_sd, outlier_ratio
 
+! JPOTERJOY: added variables
+character(8)  :: date
+character(10) :: time
+real(r8) :: prior_var, prior_mean, orig_obs_prior(ens_size)
+real(r8) :: ens_init(ens_size,ens_handle%my_num_vars)
+real(r8) :: obs_ens_init(ens_size,obs_ens_handle%my_num_vars)
+real(r8) :: wo(ens_size,ens_handle%my_num_vars), hwo(ens_size,obs_ens_handle%my_num_vars)
+real(r8) :: w(ens_size), ens_mean, ens_var, ws, obs_err_infl
+integer  :: indx(ens_size)
+
 integer  :: my_num_obs, i, j, owner, owners_index, my_num_state
 integer  :: my_obs_indx(obs_ens_handle%my_num_vars), my_state_indx(ens_handle%my_num_vars)
 integer  :: this_obs_key, obs_mean_index, obs_var_index
@@ -343,12 +358,11 @@
 type(obs_def_type)   :: obs_def
 type(time_type)      :: obs_time, this_obs_time
 
-logical :: do_adapt_inf_update
-logical :: missing_in_state
 ! for performance, local copies