[Dart-dev] [3654] DART/trunk/assim_tools:
Added namelist control to allow rectangular or trapezoidal
nancy at ucar.edu
nancy at ucar.edu
Fri Nov 14 10:55:23 MST 2008
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20081114/c59f427a/attachment.html
-------------- next part --------------
Modified: DART/trunk/assim_tools/assim_tools_mod.f90
===================================================================
--- DART/trunk/assim_tools/assim_tools_mod.f90 2008-11-12 20:17:52 UTC (rev 3653)
+++ DART/trunk/assim_tools/assim_tools_mod.f90 2008-11-14 17:55:23 UTC (rev 3654)
@@ -93,10 +93,13 @@
logical :: sampling_error_correction = .false.
integer :: adaptive_localization_threshold = -1
integer :: print_every_nth_obs = 0
+! Following only relevant for filter_kind = 8
+logical :: rectangular_quadrature = .true.
+logical :: gaussian_likelihood_tails = .false.
namelist / assim_tools_nml / filter_kind, cutoff, sort_obs_inc, &
spread_restoration, sampling_error_correction, adaptive_localization_threshold, &
- print_every_nth_obs
+ print_every_nth_obs, rectangular_quadrature, gaussian_likelihood_tails
!============================================================================
@@ -1547,21 +1550,35 @@
obs, obs_var, obs_inc)
!------------------------------------------------------------------------
!
-! Initiated 4 June, 2008
-! Attempt to clean up, improve and simplify original boxcar filter idea.
-! This will be less like a particle filter. Idea is that prior is represented
-! in the same way. In the interior boxes it is locally uniform. In the tail
-! boxes it is a gaussian although this will be defined slightly differently.
-! The likelihood is assumed gaussian here although it would be possible to
-! extend to something non-Guassian if desired for certain types of obs
-! error.
-! The product is taken over each of the interior boxes by computing the mass
-! of the likelihood funtion that fits in that box.
-! The product for the tails is trickier...
-!
-! This code is under development. Please contact Jeff Anderson (jla at ucar.edu)
-! if you want to try it out.
+! Revised 14 November 2008
+!
+! Does observation space update by approximating the prior distribution by
+! a rank histogram. Prior and posterior are assumed to have 1/(n+1) probability
+! mass between each ensemble member. The tails are assumed to be gaussian with
+! a variance equal to sample variance of the entire ensemble and a mean
+! selected so that 1/(n+1) of the mass is in each tail.
+!
+! The likelihood between the extreme ensemble members is approximated by
+! quadrature. Two options are available and controlled by the namelist entry
+! rectangular_quadrature. If this namelist is true than the likelihood between
+! a pair of ensemble members is assumed to be uniform with the average of
+! the likelihood computed at the two ensemble members. If it is false then
+! the likelihood between two ensemble members is approximated by a line
+! connecting the values of the likelihood computed at each of the ensemble
+! members (trapezoidal quadrature).
+!
+! Two options are available for approximating the likelihood on the tails.
+! If gaussian_likelihood_tails is true that the likelihood is assumed to
+! be N(obs, obs_var) on the tails. If this is false, then the likelihood
+! on the tails is taken to be uniform (to infinity) with the value at the
+! outermost ensemble members.
+!
+! A product of the approximate prior and approximate posterior is taken
+! and new ensemble members are located so that 1/(n+1) of the mass is between
+! each member and on the tails.
+! This code is still under development. Please contact Jeff Anderson at
+! jla at ucar.edu if you are interested in trying it.
integer, intent(in) :: ens_size
real(r8), intent(in) :: ens(ens_size), prior_mean, prior_var, obs, obs_var
@@ -1609,91 +1626,66 @@
1.0_r8 / (ens_size + 1.0_r8), dist_for_unit_sd)
dist_for_unit_sd = -1.0_r8 * dist_for_unit_sd
-!*********** Following block is original from early June 2008 **************
-! Compute standard deviation of lower and upper tail gaussians
-! Mean of tail distributions is ensemble mean
-! Standard deviation is selected so that exactly 1/(n+1) of mass lies
-! outside the outermost ensemble members.
-! Computing the sd for the lower tail distribution
-! This one uses the distance from the prior mean to the 1st ensemble member
-! to design a tail that is N(prior_mean, left_var) so that first ensemble
-! member is where cdf is 1/(n + 1). For small ensembles this has a tendency
-! to lose variance on the tails through sampling error.
-!!!dist_to_first = prior_mean - x(1)
-!!!left_mean = prior_mean
-!!!left_sd = dist_to_first / dist_for_unit_sd
-!!!left_var = left_sd**2
-!*********** End block *****************************************************
-
-!*********** CANDIDATE REPLACEMENT BLOCK July 2008 *************************
! Have variance of tails just be sample prior variance
! Mean is adjusted so that 1/(n+1) is outside
left_mean = x(1) + dist_for_unit_sd * prior_sd
left_var = prior_var
left_sd = prior_sd
-!*********** END BLOCK ***********************************
-
-!*********** Following block is original from early June 2008 **************
-! Get sd for the upper tail distribution
-!!!dist_to_last = x(ens_size) - prior_mean
-! Proportion to find sd so that first ensemble member is where cdf is 1/(n + 1)
-!!!right_mean = prior_mean
-!!!right_sd = dist_to_last / dist_for_unit_sd
-!!!right_var = right_sd**2
-!*********** End block ****************************************************
-
-!*********** CANDIDATE REPLACEMENT BLOCK ***********************************
-! Have variance of tails just be sample prior variance
-! Mean is adjusted so that 1/(n+1) is outside
+! Same for right tail
right_mean = x(ens_size) - dist_for_unit_sd * prior_sd
right_var = prior_var
right_sd = prior_sd
-!*********** END BLOCK ***********************************
-!*************** Block to do Gaussian-Gaussian on tail **************
-! Compute the product of the obs likelihood gaussian with the priors
-! Left tail gaussian first
-!!!var_ratio = obs_var / (left_var + obs_var)
-!!!new_var_left = var_ratio * left_var
-!!!new_sd_left = sqrt(new_var_left)
-!!!new_mean_left = var_ratio * (left_mean + left_var*obs / obs_var)
-! REMEMBER, this product has an associated weight which must be taken into account
-! See Anderson and Anderson for this weight term (or tutorial kernel filter)
-!!!prod_weight_left = exp(-0.5_r8 * (left_mean**2 / left_var + &
- !!!obs**2 / obs_var - new_mean_left**2 / new_var_left)) / &
- !!!sqrt(2.0_r8 * PI * (left_var + obs_var))
-! Determine how much mass is in the updated tails by computing gaussian cdf
-!!!mass(1) = norm_cdf(x(1), new_mean_left, new_sd_left) * prod_weight_left
-!************End Block to do Gaussian-Gaussian on tail **************
+if(gaussian_likelihood_tails) then
+ !*************** Block to do Gaussian-Gaussian on tail **************
+ ! Compute the product of the obs likelihood gaussian with the priors
+ ! Left tail gaussian first
+ var_ratio = obs_var / (left_var + obs_var)
+ new_var_left = var_ratio * left_var
+ new_sd_left = sqrt(new_var_left)
+ new_mean_left = var_ratio * (left_mean + left_var*obs / obs_var)
+ ! REMEMBER, this product has an associated weight which must be taken into account
+ ! See Anderson and Anderson for this weight term (or tutorial kernel filter)
+ ! NOTE: The constant term has been left off the likelihood so we don't have
+ ! to divide by sqrt(2 PI) in this expression
+ prod_weight_left = exp(-0.5_r8 * (left_mean**2 / left_var + &
+ obs**2 / obs_var - new_mean_left**2 / new_var_left)) / &
+ sqrt(left_var + obs_var)
+ ! Determine how much mass is in the updated tails by computing gaussian cdf
+ mass(1) = norm_cdf(x(1), new_mean_left, new_sd_left) * prod_weight_left
-! Flat tails: THIS REMOVES ALL ASSUMPTIONS ABOUT LIKELIHOOD AND CUTS COST
-new_var_left = left_var
-new_sd_left = left_sd
-new_mean_left = left_mean
-prod_weight_left = like(1)
-mass(1) = like(1) / (ens_size + 1.0_r8)
+ ! Same for the right tail
+ var_ratio = obs_var / (right_var + obs_var)
+ new_var_right = var_ratio * right_var
+ new_sd_right = sqrt(new_var_right)
+ new_mean_right = var_ratio * (right_mean + right_var*obs / obs_var)
+ ! NOTE: The constant term has been left off the likelihood so we don't have
+ ! to divide by sqrt(2 PI) in this expression
+ prod_weight_right = exp(-0.5_r8 * (right_mean**2 / right_var + &
+ obs**2 / obs_var - new_mean_right**2 / new_var_right)) / &
+ sqrt(right_var + obs_var)
+ ! Determine how much mass is in the updated tails by computing gaussian cdf
+ mass(ens_size + 1) = (1.0_r8 - norm_cdf(x(ens_size), new_mean_right, &
+ new_sd_right)) * prod_weight_right
+ !************ End Block to do Gaussian-Gaussian on tail **************
+else
+ !*************** Block to do flat tail for likelihood ****************
+ ! Flat tails: THIS REMOVES ASSUMPTIONS ABOUT LIKELIHOOD AND CUTS COST
+ new_var_left = left_var
+ new_sd_left = left_sd
+ new_mean_left = left_mean
+ prod_weight_left = like(1)
+ mass(1) = like(1) / (ens_size + 1.0_r8)
-!*************** Block to do Gaussian-Gaussian on tail **************
-! Now the right tail
-!!!var_ratio = obs_var / (right_var + obs_var)
-!!!new_var_right = var_ratio * right_var
-!!!new_sd_right = sqrt(new_var_right)
-!!!new_mean_right = var_ratio * (right_mean + right_var*obs / obs_var)
-!!!prod_weight_right = exp(-0.5_r8 * (right_mean**2 / right_var + &
- !!!obs**2 / obs_var - new_mean_right**2 / new_var_right)) / &
- !!!sqrt(2.0_r8 * PI * (right_var + obs_var))
-! Determine how much mass is in the updated tails by computing gaussian cdf
-!!!mass(ens_size + 1) = (1.0_r8 - norm_cdf(x(ens_size), new_mean_right, &
- !!!new_sd_right)) * prod_weight_right
-!************End Block to do Gaussian-Gaussian on tail **************
+ ! Same for right tail
+ new_var_right = right_var
+ new_sd_right = right_sd
+ new_mean_right = right_mean
+ prod_weight_right = like(ens_size)
+ mass(ens_size + 1) = like(ens_size) / (ens_size + 1.0_r8)
+ !*************** End block to do flat tail for likelihood ****************
+endif
-! Flat tails: THIS REMOVES ALL ASSUMPTIONS ABOUT LIKELIHOOD
-new_var_right = right_var
-new_sd_right = right_sd
-new_mean_right = right_mean
-prod_weight_right = like(ens_size)
-mass(ens_size + 1) = like(ens_size) / (ens_size + 1.0_r8)
-
! The mass in each interior box is the height times the width
! The height of the likelihood is like_dense
! For the prior, mass is 1/(n+1), and mass = height x width so...
@@ -1755,55 +1747,48 @@
! Find the box that this mass is in
if(umass >= cumul_mass(j) .and. umass <= cumul_mass(j + 1)) then
- !********* Block for rectangular quadrature *******************
- ! Linearly interpolate in mass
- new_ens(i) = x(j) + ((umass - cumul_mass(j)) / &
- (cumul_mass(j+1) - cumul_mass(j))) * (x(j + 1) - x(j))
- if(1 == 1) goto 1213
- !********* End block for rectangular quadrature *******************
+ if(rectangular_quadrature) then
+ !********* Block for rectangular quadrature *******************
+ ! Linearly interpolate in mass
+ new_ens(i) = x(j) + ((umass - cumul_mass(j)) / &
+ (cumul_mass(j+1) - cumul_mass(j))) * (x(j + 1) - x(j))
+ !********* End block for rectangular quadrature *******************
- !********* Block for trapezoidal interpolation *******************
- ! Assume that mass has linear profile, quadratic interpolation
- ! If two ensemble members are the same, just keep that value
- if(height(j + 1) < 0) then
- new_ens(i) = x(j)
else
- ! Height on left side and right side
- hleft = height(j + 1) * like(j) / mass_sum
- hright = height(j + 1) * like(j + 1) / mass_sum
- ! Will solve a quadratic for desired x-x(j)
- ! a is 0.5(hright - hleft) / (x(j+1) - x(j))
- a = 0.5_r8 * (hright - hleft) / (x(j+1) - x(j))
- ! b is hleft
- b = hleft
- ! c is cumul_mass(j) - umass
- c = cumul_mass(j) - umass
- ! Use stable quadratic solver
- call solve_quadratic(a, b, c, r1, r2)
- !!!write(*, *) 'two quadratic roots ', r1, r2
- adj_r1 = r1 + x(j)
- adj_r2 = r2 + x(j)
- if(adj_r1 >= x(j) .and. adj_r1 <= x(j+1)) then
- new_ens(i) = adj_r1
- elseif (adj_r2 >= x(j) .and. adj_r2 <= x(j+1)) then
- new_ens(i) = adj_r2
+
+ !********* Block for trapezoidal interpolation *******************
+ ! Assume that mass has linear profile, quadratic interpolation
+ ! If two ensemble members are the same, just keep that value
+ if(height(j + 1) < 0) then
+ new_ens(i) = x(j)
else
- write(*, *) 'Did not get a satisfactory quadratic root'
- write(*, *) 'j is ', j
- write(*, *) 'Adjusted roots are ', adj_r1, adj_r2
- write(*, *) 'bounds on update ', x(j), x(j+1)
- write(*, *) 'hleft, hright ', hleft, hright
- write(*, *) 'cumul mass ', cumul_mass(j), cumul_mass(j + 1)
- write(*, *) 'umass ', umass
- write(*, *) 'mass in box ', cumul_mass(j+1) - cumul_mass(j)
- write(*, *) 'comp mass ', (x(j+1) - x(j)) * (hleft + hright) / 2.0_r8
- write(*, *) 'mass_sum', mass_sum
- stop
+ ! Height on left side and right side
+ hleft = height(j + 1) * like(j) / mass_sum
+ hright = height(j + 1) * like(j + 1) / mass_sum
+ ! Will solve a quadratic for desired x-x(j)
+ ! a is 0.5(hright - hleft) / (x(j+1) - x(j))
+ a = 0.5_r8 * (hright - hleft) / (x(j+1) - x(j))
+ ! b is hleft
+ b = hleft
+ ! c is cumul_mass(j) - umass
+ c = cumul_mass(j) - umass
+ ! Use stable quadratic solver
+ call solve_quadratic(a, b, c, r1, r2)
+ adj_r1 = r1 + x(j)
+ adj_r2 = r2 + x(j)
+ if(adj_r1 >= x(j) .and. adj_r1 <= x(j+1)) then
+ new_ens(i) = adj_r1
+ elseif (adj_r2 >= x(j) .and. adj_r2 <= x(j+1)) then
+ new_ens(i) = adj_r2
+ else
+ errstring = 'Did not get a satisfactory quadratic root'
+ call error_handler(E_ERR, 'obs_increment_boxcar2', errstring, &
+ source, revision, revdate)
+ endif
endif
- endif
- !********* End block for quadratic interpolation *******************
+ !********* End block for quadratic interpolation *******************
- 1213 continue
+ endif
! Don't need to search lower boxes again
lowest_box = j
@@ -1813,7 +1798,7 @@
endif
end do
-! Now, need to convert to increments for unsorted
+! Convert to increments for unsorted
do i = 1, ens_size
obs_inc(e_ind(i)) = new_ens(i) - x(i)
end do
Modified: DART/trunk/assim_tools/assim_tools_mod.nml
===================================================================
--- DART/trunk/assim_tools/assim_tools_mod.nml 2008-11-12 20:17:52 UTC (rev 3653)
+++ DART/trunk/assim_tools/assim_tools_mod.nml 2008-11-14 17:55:23 UTC (rev 3654)
@@ -5,5 +5,7 @@
spread_restoration = .false.,
sampling_error_correction = .false.,
adaptive_localization_threshold = -1,
- print_every_nth_obs = 0 /
+ print_every_nth_obs = 0,
+ rectangular_quadrature = .true.,
+ gaussian_likelihood_tails = .false. /
More information about the Dart-dev
mailing list