[Dart-dev] DART/branches Revision: 12692
dart at ucar.edu
dart at ucar.edu
Wed Jun 27 16:31:30 MDT 2018
jla at ucar.edu
2018-06-27 16:31:30 -0600 (Wed, 27 Jun 2018)
24
Cleanup of gigg filter.
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-27 22:29:37 UTC (rev 12691)
+++ DART/branches/gigg/assimilation_code/modules/assimilation/assim_tools_mod.f90 2018-06-27 22:31:30 UTC (rev 12692)
@@ -135,7 +135,7 @@
! (9 = Localized Particle Filter (see Poterjoy 2016), in assim_tools_mod.pf.f90)
!
! 100 = GIGG, Craig Bishop. Detect interior type from observation quantity (NOT YET IMPLEMENTED)
-! 101 ="Gauss_Gauss" gaussian prior, gaussian observation likelihood (default?)
+! 101 ="Gauss_Gauss" gaussian prior, gaussian observation likelihood
! 102 ="Gamma_InvGamma" gamma prior, inverse-gamma observation likelihood
! 103 ="InvGamma_Gamma" inverse-gamma prior, gamma observation likelihood
!
@@ -2420,13 +2420,13 @@
real(r8) :: new_mean, var_ratio, a
real(r8) :: temp_mean, temp_var, new_ens(ens_size), new_var
real(r8) :: T1Rr, T2Rr
-real(r8) :: wk1,wk2,T2Pr,inv_prior_mean
-real(r8) :: GIG_gain,T1Rr_GIG_PO,y_GIG_mean_PO,y_GIG_sample_mean_PO,nobs_GIG
-real(r8) :: IGG_gain,T1Rr_IGG_PO,y_IGG_mean_PO,y_IGG_sample_mean_PO,nobs_IGG,nanal_var_IGG,anal_var_IGG
-real(r8) :: y_GIG_var_PO,y_IGG_var_PO,T2Rr_IGG_PO
-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) :: wk1, wk2, T2Pr, inv_prior_mean
+real(r8) :: GIG_gain, T1Rr_GIG_PO, y_GIG_mean_PO, y_GIG_sample_mean_PO, nobs_GIG
+real(r8) :: IGG_gain, T1Rr_IGG_PO, y_IGG_mean_PO, y_IGG_sample_mean_PO, nobs_IGG, nanal_var_IGG, anal_var_IGG
+real(r8) :: y_GIG_var_PO, y_IGG_var_PO, T2Rr_IGG_PO
+real(r8) :: alpha, theta, lambda, beta, k_gamma, scale
+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
@@ -2455,7 +2455,8 @@
! FIXME: compute something based on the obs value (obs), the obs variance (obs_var),
! the prior ensemble distribution (ens(ens_size), the prior mean and variance (prior_mean,
! 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
+! and set GIGG_type here.
+! The value of the type 1 relative error variance is assumed to be in the obs_seq file at this point
GIGG_type = filter_kind - 100 ! or 2 or 3
T1Rr = obs_var
@@ -2474,92 +2475,94 @@
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)
+ T2Pr = prior_var / (prior_mean**2 + prior_var)
+ T2Rr = 1.0_r8 / ((1.0_r8 / T1Rr) + 1.0_r8)
!Define the GIG gain
- GIG_gain=T2Pr/(T2Pr+T2Rr)
- inv_prior_mean=1/prior_mean
+ GIG_gain = T2Pr / (T2Pr + T2Rr)
+ inv_prior_mean = 1.0_r8 / prior_mean
!Compute new_mean - the posterior mean. Equation (6) from Bishop (2016)
- wk1=inv_prior_mean+GIG_gain*((1/obs)-(T2Rr+1)*inv_prior_mean)
- new_mean=1/wk1;
+ wk1 = inv_prior_mean + GIG_gain * ((1.0_r8 / obs) - (T2Rr + 1) * inv_prior_mean)
+ new_mean = 1.0_r8 /wk1
!Compute obs_inc - the correction to the existing ensemble
- !first compute required shape parameters
- k_prior=(1/T2Pr)-1 !shape parameter for gamma approximation to prior
- k_gig=(1/T2Rr)+1 !shape parameter for perturbations to be added to forecast
- k_post=k_prior+k_gig !shape parameter for posterior distribution of gig assumptions satisfied
- !now compute rescaling factor to ensure perturbed forecasts have posterior mean
- rescale=(k_prior*new_mean)/(k_post*prior_mean)
+ !first compute required shape parameters
+ k_prior = (1.0_r8 / T2Pr) - 1.0_r8 ! shape parameter for gamma approximation to prior
+ k_gig = (1.0_r8 / T2Rr) + 1.0_r8 ! shape parameter for perturbations to be added to forecast
+ k_post = k_prior + k_gig ! shape parameter for posterior distribution of gig assumptions satisfied
+ ! Compute rescaling factor to ensure perturbed forecasts have posterior mean
+ rescale = (k_prior*new_mean) / (k_post*prior_mean)
!Now perturb the forecast perturbations so that their shape parameter would be the same as the posterior
- 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 (prior_mean < 0) then
- print *, 'correcting negative prior_mean: ', prior_mean
- !!!theta_gig = 10**(-10)/k_gig
+ theta_gig = prior_mean / k_prior ! note that theta_gig=theta_mean as required by summation theorem
+ wk_samp_mean = 0.0_r8 ! summation variable to compute sample mean
+ do i = 1, ens_size
+ if (prior_mean < 0.0_r8) then
prior_mean_temp = 10e-4_r8
- theta_gig=prior_mean_temp/k_prior !note that theta_gig=theta_mean as required by summation theorem
+ msgstring = 'Correcting negative prior_mean for GIG filter '
+ call error_handler(E_MSG, 'obs_increment_gigg:', msgstring)
+ 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)
+ 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
More information about the Dart-dev
mailing list