[Dart-dev] DART/branches Revision: 12616

[dart at ucar.edu]

jla at ucar.edu
2018-06-01 13:42:35 -0600 (Fri, 01 Jun 2018)
104
Gigg filter modifications. Tested with Craig and confirmed correct for single prior observation tests.




Modified: DART/branches/gigg/assimilation_code/modules/assimilation/assim_tools_mod.f90
===================================================================
--- DART/branches/gigg/assimilation_code/modules/assimilation/assim_tools_mod.f90	2018-06-01 19:07:17 UTC (rev 12615)
+++ DART/branches/gigg/assimilation_code/modules/assimilation/assim_tools_mod.f90	2018-06-01 19:42:35 UTC (rev 12616)
@@ -133,10 +133,11 @@
 !      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.  Can do the following subtypes:
-!          "Gauss_Gauss"     gaussian prior, gaussian observation likelihood (default?)
-!          "Gamma_InvGamma"  gamma prior, inverse-gamma observation likelihood
-!          "InvGamma_Gamma"  inverse-gamma prior, gamma observation likelihood
+!
+!     100 = GIGG, Craig Bishop. Detect interior type from observation quantity (NOT YET IMPLEMENTED) 
+!     101 ="Gauss_Gauss"     gaussian prior, gaussian observation likelihood (default?)
+!     102 ="Gamma_InvGamma"  gamma prior, inverse-gamma observation likelihood
+!     103 ="InvGamma_Gamma"  inverse-gamma prior, gamma observation likelihood
 !      
 !
 !  special_localization_obs_types -> Special treatment for the specified observation types
@@ -1341,11 +1342,11 @@
       call obs_increment_boxcar(ens, ens_size, obs, obs_var, obs_inc, rel_weights)
     case (8) 
       call obs_increment_rank_histogram(ens, ens_size, prior_var, obs, obs_var, obs_inc)
-    case (10)
+    case (100, 101, 102, 103)
       call obs_increment_GIGG(ens, ens_size, prior_mean, prior_var, obs, obs_var, quantity, obs_inc)
     case default
       call error_handler(E_ERR,'obs_increment', &
-                 'Unexpected error: Illegal value of filter_kind in assim_tools namelist [1-8,10 OK]', &
+                 'Unexpected error: Illegal value of filter_kind in assim_tools namelist [1-8, 100-103 OK]', &
                  source, revision, revdate)
    end select
 endif
@@ -2426,6 +2427,8 @@
 real(r8) :: alpha,theta,lambda,beta,k_gamma
 real(r8) :: nanal_perts(ens_size),nfcst_perts(ens_size),ny_PO(ens_size),nfcst
 real(r8) :: y_GIG_PO(ens_size),y_IGG_PO(ens_size)
+real(r8) :: k_gig, k_post, k_prior, rescale, theta_gig, wk_samp_mean, correction_ratio
+real(r8) :: y_anal_wk(ens_size), y_anal(ens_size), prior_mean_temp
 
 
 
@@ -2432,6 +2435,17 @@
 ! See comment by other init_random_seq() call.  Will preproduce results
 ! exactly for the same number of mpi tasks; will not if # tasks changes.
 
+! Have three namelist variants for filter_kind that come here
+! 100 means that the choice of filter kind within gigg is a function of quantity (not yet implemented)
+! 101 means just do EAKF
+! 102 means GIG (gamma prior, inverse gamma likelihood)
+! 103 means IGG (inverse gamma prior, gamma likelihood)
+
+if(filter_kind == 100) then
+   call error_handler(E_ERR, 'obs_increment_GIGG:', 'filter_kind 100 not yet implemented', &
+                      source, revision, revdate)
+endif
+
 ! seed the random number generator the first time through this code.
 if(first_ran_gigg) then
    call init_random_seq(gigg_ran_seq, my_task_id() + 1)
@@ -2443,8 +2457,8 @@
 ! prior_var), and the quantity of interest (quantity, e.g. QTY_TEMPERATURE, QTY_U_WIND_COMPONENT, etc)
 ! and set GIGG_type here.  also set T1Rr
 
-GIGG_type = 1 ! or 2 or 3
-T1Rr = 1.0_r8
+GIGG_type = filter_kind - 100 ! or 2 or 3
+T1Rr = obs_var
 
 select case (GIGG_type)
  case (1)  ! gaussian ens prior, gaussian obs likelihood
@@ -2458,6 +2472,7 @@
    obs_inc = a * (ens - prior_mean) + new_mean - ens
 
  case (2)   ! gamma ens prior, inverse gamma obs likelihood
+
    !Define type 2 prior and ob relative variance
    T2Pr=prior_var/(prior_mean**2+prior_var)
    T2Rr=1/((1/T1Rr)+1)
@@ -2480,12 +2495,14 @@
    theta_gig=prior_mean/k_prior  !note that theta_gig=theta_mean as required by summation theorem
    wk_samp_mean=0                !summation variable to compute sample mean
    do i=1,ens_size
-if (theta_gig < 0) then
-   print *, 'correcting negative theta: ', theta
-   theta_gig = 10**(-10)/k_gig
-endif
-       y_anal_wk(i)=rescale*(ens(i)+random_gamma(gigg_ran_seq,k_gig,theta_gig))
-       wk_samp_mean=wk_samp_mean+y_anal_wk(i)
+      if (prior_mean < 0) then
+         print *, 'correcting negative prior_mean: ', prior_mean
+         !!!theta_gig = 10**(-10)/k_gig
+         prior_mean_temp = 10e-4_r8
+         theta_gig=prior_mean_temp/k_prior  !note that theta_gig=theta_mean as required by summation theorem
+      endif
+      y_anal_wk(i)=rescale*(ens(i)+random_gamma(gigg_ran_seq, k_gig, theta_gig))
+      wk_samp_mean=wk_samp_mean+y_anal_wk(i)
    enddo
    !correct degradation in posterior mean associated with sampling
    wk_samp_mean=wk_samp_mean/ens_size