[Dart-dev] DART/branches Revision: 11857

[dart at ucar.edu]

nancy at ucar.edu
2017-08-02 09:26:45 -0600 (Wed, 02 Aug 2017)
280
first stab at adding the GIGG code to filter_assim().
does not work - i don't understand the T1Rr input
and it's generating negative values for the scale
parameter into the random gamma/inv gamma generator.

but it does compile!  craig b. will have to help
with test cases here.




Modified: DART/branches/rma_gigg/assimilation_code/modules/assimilation/assim_tools_mod.f90
===================================================================
--- DART/branches/rma_gigg/assimilation_code/modules/assimilation/assim_tools_mod.f90	2017-08-01 23:03:11 UTC (rev 11856)
+++ DART/branches/rma_gigg/assimilation_code/modules/assimilation/assim_tools_mod.f90	2017-08-02 15:26:45 UTC (rev 11857)
@@ -14,10 +14,10 @@
 use  utilities_mod,       only : file_exist, get_unit, check_namelist_read, do_output,    &
                                  find_namelist_in_file, register_module, error_handler,   &
                                  E_ERR, E_MSG, nmlfileunit, do_nml_file, do_nml_term,     &
-                                 open_file, close_file, timestamp
+                                 open_file, close_file, timestamp, to_upper
 use       sort_mod,       only : index_sort 
 use random_seq_mod,       only : random_seq_type, random_gaussian, init_random_seq,       &
-                                 random_uniform
+                                 random_uniform, random_gamma, random_inverse_gamma
 
 use obs_sequence_mod,     only : obs_sequence_type, obs_type, get_num_copies, get_num_qc, &
                                  init_obs, get_obs_from_key, get_obs_def, get_obs_values, &
@@ -88,6 +88,10 @@
 logical                :: first_inc_ran_call = .true.
 type (random_seq_type) :: inc_ran_seq
 
+! sequence for the GIGG filter
+logical                :: first_ran_gigg = .true.
+type (random_seq_type) :: gigg_ran_seq
+
 integer                :: num_types = 0
 real(r8), allocatable  :: cutoff_list(:)
 logical                :: has_special_cutoffs
@@ -108,6 +112,12 @@
 ! and fill this 2d impact table. 
 real(r8), allocatable  :: obs_impact_table(:,:)
 
+! convert string to an integer in the init routine
+integer, parameter :: GAUSS_GAUSS    = 1
+integer, parameter :: GAMMA_INVGAMMA = 2
+integer, parameter :: INVGAMMA_GAMMA = 3
+integer :: GIGG_type = GAUSS_GAUSS   ! should it be something else?
+
 ! version controlled file description for error handling, do not edit
 character(len=256), parameter :: source   = &
    "$URL$"
@@ -126,6 +136,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 (see Poterjoy 2016), in assim_tools_mod.pf.f90)
+!     10 = GIGG, Craig Bishop. Needs additional GIGG_distribution_type flag as follows:
+!          "Gauss_Gauss"     gaussian prior, gaussian observation likelihood (default?)
+!          "Gamma_InvGamma"  gamma prior, inverse-gamma observation likelihood
+!          "InvGamma_Gamma"  inverse-gamma prior, gamma observation likelihood
+!      
 !
 !  special_localization_obs_types -> Special treatment for the specified observation types
 !  special_localization_cutoffs   -> Different cutoff value for each specified obs type
@@ -138,6 +154,7 @@
 integer  :: adaptive_localization_threshold = -1
 real(r8) :: adaptive_cutoff_floor           = 0.0_r8
 integer  :: print_every_nth_obs             = 0
+character(len=64) :: GIGG_distribution_type = "Gauss_Gauss"   ! DOES THIS MAKE SENSE AS DEFAULT?
 
 ! since this is in the namelist, it has to have a fixed size.
 integer, parameter   :: MAX_ITEMS = 300
@@ -211,6 +228,7 @@
    allow_missing_in_clm, distribute_mean, close_obs_caching,               &
    adjust_obs_impact, obs_impact_filename, allow_any_impact_values,        &
    convert_all_state_verticals_first, convert_all_obs_verticals_first,     &
+   GIGG_distribution_type,                                                 &
    lanai_bitwise ! don't document this one -- only used for regression tests
 
 !============================================================================
@@ -224,6 +242,7 @@
 integer :: iunit, io, i, j
 integer :: num_special_cutoff, type_index
 logical :: cache_override = .false.
+character(len=64) :: GIGG_upper
 
 call register_module(source, revision, revdate)
 
@@ -312,6 +331,28 @@
 is_doing_vertical_conversion = (has_vertical_choice() .and. vertical_localization_on() .and. &
                                 .not. lanai_bitwise)
 
+! set GIGG_type 1, 2 or 3 here via select
+!          "Gauss_Gauss"     gaussian prior, gaussian observation likelihood (default)
+!          "Gamma_InvGamma"  gamma prior, inverse-gamma observation likelihood
+!          "InvGamma_Gamma"  inverse-gamma prior, gamma observation likelihood
+
+GIGG_upper = GIGG_distribution_type
+call to_upper(GIGG_upper)   ! uppercases in place
+
+select case (GIGG_upper)
+ case("GAUSS_GAUSS")
+   GIGG_type = GAUSS_GAUSS
+ case("GAMMA_INVGAMMA")
+   GIGG_type = GAMMA_INVGAMMA 
+ case("INVGAMMA_GAMMA")
+   GIGG_type = INVGAMMA_GAMMA 
+ case default
+   write(msgstring, *) 'Unrecognized GIGG filter type string: "'//trim(GIGG_distribution_type)//'"'
+   call error_handler(E_ERR, 'assim_tools_init: ', &
+                      'Valid types are "Gauss_Gauss", "Gamma_InvGamma", and "InvGamma_Gamma"', &