gamma_distribution_mod.f90: remove partial support for generalized ga…

…mma distributions.

jla commit message:
A comment now makes it clear that this module only supports the standard gamma
distribution that is bounded below by 0 and unbounded above. Subroutines gamma_cdf, inverse_gamma_cdf,
set_gamma_params_from_ens, and inv_gamma_first_guess no longer have bounded_below, bounded_above,
lower_bound, and upper_bound as arguments. inv_gamma_cdf now sets the bounded_below, bounded_above,
lower_bound and upper_bound parameters to the correct values for the gamma.

Functions inv_gamma_cdf_params and gamma_cdf_params now include an error check to make sure that the
lower and upper bounds in the distribution_params_type have been set to 0 and missing_r8 as other
values are not supported.

Subroutine set_gamma_params_from_ens has changed the distribution_params_type to an intent out
argument and defines all six parameters correctly. It also sets the parameter type to
GAMMA_DISTRIBUTION.

assimilation_code/modules/assimilation/assim_tools_mod.f90

@@ -1065,7 +1065,7 @@ subroutine obs_increment_gamma(ens, ens_size, prior_mean, prior_var, obs, obs_va
 ! Compute the prior quantiles of each ensemble member in the prior gamma distribution
 call gamma_mn_var_to_shape_scale(prior_mean, prior_var, prior_shape, prior_scale)
 do i = 1, ens_size
-   q(i) = gamma_cdf(ens(i), prior_shape, prior_scale, .true., .false., 0.0_r8, missing_r8) 
+   q(i) = gamma_cdf(ens(i), prior_shape, prior_scale) 
 end do
 
 ! Compute the statistics of the continous posterior distribution
@@ -1082,7 +1082,7 @@ subroutine obs_increment_gamma(ens, ens_size, prior_mean, prior_var, obs, obs_va
 
 ! Now invert the quantiles with the posterior distribution
 do i = 1, ens_size
-   post(i) = inv_gamma_cdf(q(i), post_shape, post_scale, .true., .false., 0.0_r8, missing_r8)
+   post(i) = inv_gamma_cdf(q(i), post_shape, post_scale)
 end do
 
 obs_inc = post - ens

assimilation_code/modules/assimilation/gamma_distribution_mod.f90

@@ -2,6 +2,9 @@
 ! by UCAR, "as is", without charge, subject to all terms of use at
 ! http://www.image.ucar.edu/DAReS/DART/DART_download
 
+! This module only supports standard gamma distributions for which the 
+! lower bound is 0.
+
 module gamma_distribution_mod
 
 use types_mod,               only : r8, PI, missing_r8
@@ -10,7 +13,7 @@ module gamma_distribution_mod
 
 use normal_distribution_mod, only : normal_cdf, inv_cdf
 
-use distribution_params_mod, only : distribution_params_type
+use distribution_params_mod, only : distribution_params_type, GAMMA_DISTRIBUTION
 
 use random_seq_mod,          only : random_seq_type, random_uniform
 
@@ -60,18 +63,17 @@ subroutine test_gamma
 ! Compare to matlab
 write(*, *) 'Absolute value of differences should be less than 1e-15'
 do i = 1, 7
-   pdf_diff(i) = gamma_pdf(mx(i), mshape(i), mscale(i)) - mpdf(i)
-   cdf_diff(i) = gamma_cdf(mx(i), mshape(i), mscale(i), .true., .false., 0.0_r8, missing_r8) - mcdf(i)
+   pdf_diff(i) = abs(gamma_pdf(mx(i), mshape(i), mscale(i)) - mpdf(i))
+   cdf_diff(i) = abs(gamma_cdf(mx(i), mshape(i), mscale(i)) - mcdf(i))
    write(*, *) i, pdf_diff(i), cdf_diff(i)
 end do
-
-if(abs(maxval(pdf_diff)) < 1e-15_r8 .and. abs(maxval(cdf_diff)) < 1e-15_r8) then
+if(maxval(pdf_diff) < 1e-15_r8 .and. maxval(cdf_diff) < 1e-15_r8) then
    write(*, *) 'Matlab Comparison Tests: PASS'
 else
    write(*, *) 'Matlab Compariosn Tests: FAIL'
 endif
 
-! Input a mean and variance
+! Test many x values for cdf and inverse cdf for a single mean and sd (shape and scale)
 mean = 10.0_r8
 sd = 1.0_r8
 variance = sd**2
@@ -84,22 +86,20 @@ subroutine test_gamma
 max_diff = -1.0_r8
 do i = 0, 1000
    x = mean + ((i - 500.0_r8) / 500.0_r8) * 5.0_r8 * sd
-   y = gamma_cdf(x, gamma_shape, gamma_scale, .true., .false., 0.0_r8, missing_r8)
-   inv = inv_gamma_cdf(y, gamma_shape, gamma_scale, .true., .false., 0.0_r8, missing_r8)
+   y = gamma_cdf(x, gamma_shape, gamma_scale)
+   inv = inv_gamma_cdf(y, gamma_shape, gamma_scale)
    max_diff = max(abs(x-inv), max_diff)
 end do
 
 write(*, *) '----------------------------'
 write(*, *) 'max difference in inversion is ', max_diff
 write(*, *) 'max difference should be less than 1e-11'
-
 if(max_diff < 1e-11_r8) then
    write(*, *) 'Inversion Tests: PASS'
 else
    write(*, *) 'Inversion Tests: FAIL'
 endif
 
-
 end subroutine test_gamma
 
 !-----------------------------------------------------------------------
@@ -110,31 +110,33 @@ function inv_gamma_cdf_params(quantile, p) result(x)
 real(r8),                       intent(in) :: quantile
 type(distribution_params_type), intent(in) :: p
 
-! Could do error checks for gamma_shape and gamma_scale values here
+! Only standard gamma currently supported 
+if(p%lower_bound /= 0.0_r8 .and. p%upper_bound /= missing_r8) then 
+   errstring = 'Only standard gamma distribution with lower bound of 0 is supported' 
+   call error_handler(E_ERR, 'inv_gamma_cdf_params', errstring, source)
+end if
+
 x = inv_cdf(quantile, gamma_cdf_params, inv_gamma_first_guess_params, p)
 
 end function inv_gamma_cdf_params
 !-----------------------------------------------------------------------
 
-function inv_gamma_cdf(quantile, gamma_shape, gamma_scale, &
-               bounded_below, bounded_above, lower_bound, upper_bound) result(x)
+function inv_gamma_cdf(quantile, gamma_shape, gamma_scale) result(x)
 
 real(r8)             :: x
 real(r8), intent(in) :: quantile
 real(r8), intent(in) :: gamma_shape
 real(r8), intent(in) :: gamma_scale
-logical,  intent(in) :: bounded_below, bounded_above
-real(r8), intent(in) :: lower_bound,   upper_bound
 
 ! Given a quantile q, finds the value of x for which the gamma cdf
 ! with shape and scale has approximately this quantile
 
 type(distribution_params_type) :: p
 
 ! Load the p type for the generic cdf calls
-p%params(1) = gamma_shape; p%params(2) = gamma_scale
-p%bounded_below = bounded_below;      p%bounded_above = bounded_above
-p%lower_bound   = lower_bound;        p%upper_bound   = upper_bound
+p%params(1)     = gamma_shape; p%params(2)     = gamma_scale
+p%bounded_below = .true.;      p%bounded_above = .false.
+p%lower_bound   = 0.0_r8;      p%upper_bound   = missing_r8
 
 x = inv_gamma_cdf_params(quantile, p)
 
@@ -174,23 +176,28 @@ function gamma_cdf_params(x, p)
 
 real(r8) :: gamma_shape, gamma_scale
 
+! Only standard gamma currently supported 
+if(p%lower_bound /= 0.0_r8 .and. p%upper_bound /= missing_r8) then 
+   errstring = 'Only standard gamma distribution with lower bound of 0 is supported' 
+   call error_handler(E_ERR, 'gamma_cdf_params', errstring, source)
+end if
+
 gamma_shape = p%params(1);     gamma_scale = p%params(2)
-gamma_cdf_params = gamma_cdf(x, gamma_shape, gamma_scale, &
-                     p%bounded_below, p%bounded_above, p%lower_bound, p%upper_bound)
+gamma_cdf_params = gamma_cdf(x, gamma_shape, gamma_scale)
 
 end function gamma_cdf_params
 
 !---------------------------------------------------------------------------
 
-function gamma_cdf(x, gamma_shape, gamma_scale, bounded_below, bounded_above, lower_bound, upper_bound)
+function gamma_cdf(x, gamma_shape, gamma_scale)
 
 ! Returns the cumulative distribution of a gamma function with shape and scale
 ! at the value x
 
+! Returns a failed_value if called with illegal values
+
 real(r8) :: gamma_cdf
 real(r8), intent(in) :: x, gamma_shape, gamma_scale
-logical,  intent(in) :: bounded_below, bounded_above
-real(r8), intent(in) :: lower_bound,   upper_bound
 
 ! All inputs must be nonnegative
 if(x < 0.0_r8 .or. gamma_shape < 0.0_r8 .or. gamma_scale < 0.0_r8) then
@@ -360,7 +367,7 @@ function random_gamma(r, rshape, rscale)
 ! Draw from U(0, 1) to get a quantile
 quantile = random_uniform(r)
 ! Invert cdf to get a draw from gamma
-random_gamma = inv_gamma_cdf(quantile, rshape, rscale, .true., .false., 0.0_r8, missing_r8)
+random_gamma = inv_gamma_cdf(quantile, rshape, rscale)
 
 end function random_gamma
 
@@ -425,7 +432,7 @@ function inv_gamma_first_guess_params(quantile, p)
 ! A translation routine that is required to use the generic first_guess for
 ! the cdf  optimization routine.
 ! Extracts the appropriate information from the distribution_params_type that is needed
-! for a call to the function approx_inv_normal_cdf below (which is nothing).
+! for a call to the function inv_gamma_first_guess below.
 
 real(r8) :: gamma_shape, gamma_scale
 
@@ -446,26 +453,25 @@ function inv_gamma_first_guess(quantile, gamma_shape, gamma_scale)
 ! For starters, take the mean for this shape and scale
 inv_gamma_first_guess = gamma_shape * gamma_scale
 ! Could use info about sd to further refine mean and reduce iterations
-!!!sd = sqrt(gamma_shape * gamma_scale**2)
 
 end function inv_gamma_first_guess
 
 !---------------------------------------------------------------------------
 
-subroutine set_gamma_params_from_ens(ens, num, bounded_below, bounded_above, &
-                                     lower_bound, upper_bound, p)
+subroutine set_gamma_params_from_ens(ens, num, p)
 
 integer,                        intent(in)    :: num
 real(r8),                       intent(in)    :: ens(num)
-logical,                        intent(in)    :: bounded_below, bounded_above
-real(r8),                       intent(in)    :: lower_bound,   upper_bound
-type(distribution_params_type), intent(inout) :: p
+type(distribution_params_type), intent(out)   :: p
 
 real(r8) :: gamma_shape, gamma_scale
 
+! Set up the description of the gamma distribution defined by the ensemble
+p%distribution_type = GAMMA_DISTRIBUTION
+
 ! Set the bounds info
-p%bounded_below = bounded_below;  p%bounded_above = bounded_above
-p%lower_bound   = lower_bound;    p%upper_bound   = upper_bound
+p%bounded_below = .true.;    p%bounded_above = .false.
+p%lower_bound   = 0.0_r8;    p%upper_bound   = missing_r8
 
 ! Get shape and scale
 call gamma_shape_scale(ens, num, gamma_shape, gamma_scale)

assimilation_code/modules/assimilation/probit_transform_mod.f90

@@ -278,8 +278,7 @@ subroutine to_probit_gamma(ens_size, state_ens, p, probit_ens, use_input_p, &
       call error_handler(E_ERR, 'to_probit_gamma', errstring, source)
    endif
 
-   call set_gamma_params_from_ens(state_ens, ens_size, bounded_below, bounded_above, &
-      lower_bound, upper_bound, p)
+   call set_gamma_params_from_ens(state_ens, ens_size, p)
 endif
 
 do i = 1, ens_size