[Dart-dev] [3438] DART/trunk/assim_tools/assim_tools_mod.f90: Modified the computation of the tails for the priors in the

nancy at ucar.edu nancy at ucar.edu
Tue Jul 1 12:16:21 MDT 2008


An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080701/97cca3e5/attachment.html
-------------- next part --------------
Modified: DART/trunk/assim_tools/assim_tools_mod.f90
===================================================================
--- DART/trunk/assim_tools/assim_tools_mod.f90	2008-06-30 18:00:04 UTC (rev 3437)
+++ DART/trunk/assim_tools/assim_tools_mod.f90	2008-07-01 18:16:21 UTC (rev 3438)
@@ -1401,8 +1401,9 @@
 
 integer  :: i, e_ind(ens_size), lowest_box, j
 real(r8) :: a
-real(r8) :: sx, prior_mean, prior_var
+real(r8) :: sx, prior_mean, prior_var, prior_var_d2
 real(r8) :: var_ratio, new_var, new_sd, umass, left_weight, right_weight
+real(r8) :: new_mean, temp_mean, temp_var
 real(r8) :: mass(2*ens_size), weight(ens_size), cumul_mass(0:2*ens_size)
 real(r8) :: new_mean_left, new_mean_right, prod_weight_left, prod_weight_right
 real(r8) :: new_ens(ens_size), mass_sum, const_term
@@ -1449,21 +1450,21 @@
 ! Since 1/2 of a normal is outside, need to multiply by 2 / (ens_size + 1)
 
 ! Need some sort of width for the boundary kernel, try 1/2 the VAR for now
-prior_var = prior_var / 2.0_r8
+prior_var_d2 = prior_var / 2.0_r8
 
 ! Compute the product of the obs error gaussian with the prior gaussian (EAKF)
 ! Left wing first
-var_ratio = obs_var / (prior_var + obs_var)
-new_var = var_ratio * prior_var
+var_ratio = obs_var / (prior_var_d2 + obs_var)
+new_var = var_ratio * prior_var_d2
 new_sd = sqrt(new_var)
-new_mean_left  = var_ratio * (ens(e_ind(1))  + prior_var*obs / obs_var)
-new_mean_right  = var_ratio * (ens(e_ind(ens_size))  + prior_var*obs / obs_var)
+new_mean_left  = var_ratio * (ens(e_ind(1))  + prior_var_d2*obs / obs_var)
+new_mean_right  = var_ratio * (ens(e_ind(ens_size))  + prior_var_d2*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 =  2.71828_r8 ** (-0.5_r8 * (ens(e_ind(1))**2 / prior_var + &
+prod_weight_left =  2.71828_r8 ** (-0.5_r8 * (ens(e_ind(1))**2 / prior_var_d2 + &
       obs**2 / obs_var - new_mean_left**2 / new_var)) / sqrt(2.0_r8 * PI)
 
-prod_weight_right =  2.71828_r8 ** (-0.5_r8 * (ens(e_ind(ens_size))**2 / prior_var + &
+prod_weight_right =  2.71828_r8 ** (-0.5_r8 * (ens(e_ind(ens_size))**2 / prior_var_d2 + &
       obs**2 / obs_var - new_mean_right**2 / new_var)) / sqrt(2.0_r8 * PI)
 
 ! Split into 2*ens_size domains; mass in each is computed
@@ -1544,7 +1545,7 @@
 
 subroutine obs_increment_boxcar2(ens, ens_size, 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
@@ -1552,7 +1553,7 @@
 ! 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.
+! 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...
@@ -1566,9 +1567,9 @@
 real(r8), intent(out) :: obs_inc(ens_size)
 
 integer  :: i, e_ind(ens_size), lowest_box, j
-real(r8) :: prior_mean, prior_var
+real(r8) :: prior_mean, prior_var, prior_sd
 real(r8) :: var_ratio, new_var, umass, left_mass, right_mass
-real(r8) :: left_sd, left_var, right_sd, right_var
+real(r8) :: left_sd, left_var, right_sd, right_var, left_mean, right_mean
 real(r8) :: new_mean, temp_mean, temp_var
 real(r8) :: mass(ens_size + 1), like(ens_size), cumul_mass(0:ens_size + 1)
 real(r8) :: nmass(ens_size + 1)
@@ -1606,38 +1607,67 @@
 prior_mean = sum(ens) / ens_size
 ! prior_var only used for diagnostic evaluation during development
 prior_var  = sum((ens - prior_mean)**2) / (ens_size - 1)
+prior_sd = sqrt(prior_var)
 
 
 ! 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
+! 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
-dist_to_first = prior_mean - x(1)
+! Computing the sd for the lower tail distribution 
 ! For unit normal, find distance from mean to where cdf is 1/(n+1)
 ! Lots of this can be done once in first call and then saved
 call weighted_norm_inv(1.0_r8, 0.0_r8, 1.0_r8, &
    1.0_r8 / (ens_size + 1.0_r8), dist_for_unit_sd)
 dist_for_unit_sd = -1.0_r8 * dist_for_unit_sd
-! Proportion to find sd so that first ensemble member is where cdf is 1/(n + 1)
-left_sd = dist_to_first / dist_for_unit_sd
-left_var = left_sd**2
 
+
+!*********** Following block is original from early June 2008 **************
+! 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
+!!!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_sd = dist_to_last / dist_for_unit_sd
-right_var = right_sd**2
+!!!right_mean = prior_mean
+!!!right_sd = dist_to_last / dist_for_unit_sd
+!!!right_var = right_sd**2
+!*********** End block ****************************************************
 
-! Compute the product of the obs likelihood gaussian with the priors
+!*********** CANDIDATE REPLACEMENT BLOCK ***********************************
+! Have variance of tails just be sample prior variance
+! Mean is adjusted so that 1/(n+1) is outside
+right_mean = x(ens_size) - dist_for_unit_sd * prior_sd
+right_var = prior_var
+right_sd = prior_sd
+!*********** END BLOCK ***********************************
+
+
+! 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 * (prior_mean  + left_var*obs / obs_var)
+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 * (prior_mean**2 / left_var + &
+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))
 
@@ -1645,8 +1675,8 @@
 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 * (prior_mean  + right_var*obs / obs_var)
-prod_weight_right =  exp(-0.5_r8 * (prior_mean**2 / right_var + &
+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))
 
@@ -1662,7 +1692,7 @@
 ! The height of the prior is 1 / ((n+1) width);   multiplying by width leaves 1/(n+1)
 
 ! In prior, have 1/(n+1) mass in each bin, multiply by mean likelihood density
-! to get approximate mass in updated bin
+! to get approximate mass in updated bin 
 do i = 2, ens_size
    mass(i) = like_dense(i) / (ens_size + 1.0_r8)
    ! Height of prior in this bin is mass/width
@@ -1712,7 +1742,7 @@
          if(umass >= cumul_mass(j) .and. umass <= cumul_mass(j + 1)) then
 
             !!! Commented block does rectangular quadrature
-            !!! Without this comment does trapezoidal quadrature
+            !!! Block below does trapezoidal quadrature that is more accurate
             !!! 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))
@@ -1733,7 +1763,7 @@
             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
+            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
@@ -1750,9 +1780,9 @@
                write(*, *) 'mass_sum', mass_sum
                stop
             endif
+            
+           1213 continue 
 
-            1213 continue
-
             ! Don't need to search lower boxes again
             lowest_box = j
             exit FIND_BOX


More information about the Dart-dev mailing list