Merge branch 'main' into qceff_dartlab

Author: hkershaw-brown

.gitignore

@@ -207,6 +207,7 @@ test_normal_dist
 test_beta_dist
 test_kde_dist
 test_window
+test_force_bounds
 
 # Directories to NOT IGNORE ... same as executable names
 # as far as I know, these must be listed after the executables

CHANGELOG.rst

@@ -22,6 +22,27 @@ individual files.
 
 The changes are now listed with the most recent at the top.
 
+**January 23 2025 :: Pangu-DART. Tag v11.9.0**
+
+- Pangu-Weather ML model DART interface
+
+*contributed by Nuo Chen, University of Oklahoma*
+
+**January 22 2025 :: Bug-fix: Gamma and Beta Distributions. Tag v11.8.9**
+
+Bug fixes:
+
+    - Beta distribution only supporting standard Beta, bounded 0-1.
+    - Gamma distribution only supporting standard, lower bound 0.
+    - Beta and Gamma bounds are forced in the QCEFF table.
+
+Updates:
+
+    - Explicitly setting distribution type, now have UNSET.
+    - Message about failing to converge changed to E_ALLMSG to be visible
+      on all mpi ranks.
+    - remove unused test_obs directory
+
 **January 14 2025 :: Bug-fix MOM6 potential temperature. Tag v11.8.8**
 
 - MOM6 model_interpolate for potential temperature

assimilation_code/modules/assimilation/algorithm_info_mod.f90

@@ -217,9 +217,19 @@ subroutine read_qceff_table(qceff_table_filename)
       case ('BOUNDED_NORMAL_RH_DISTRIBUTION')
          qceff_table_data(row)%probit_inflation%dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION
       case ('GAMMA_DISTRIBUTION')
+         ! Force standard Gamma distribution
          qceff_table_data(row)%probit_inflation%dist_type = GAMMA_DISTRIBUTION
+         qceff_table_data(row)%probit_inflation%bounded_above = .false.
+         qceff_table_data(row)%probit_inflation%bounded_below = .true.
+         qceff_table_data(row)%probit_inflation%upper_bound = MISSING_R8
+         qceff_table_data(row)%probit_inflation%lower_bound = 0.0_r8
       case ('BETA_DISTRIBUTION')
+         ! Force standard Beta distribution
          qceff_table_data(row)%probit_inflation%dist_type = BETA_DISTRIBUTION
+         qceff_table_data(row)%probit_inflation%bounded_above = .true.
+         qceff_table_data(row)%probit_inflation%bounded_below = .true.
+         qceff_table_data(row)%probit_inflation%upper_bound = 1.0_r8
+         qceff_table_data(row)%probit_inflation%lower_bound = 0.0_r8
       case ('LOG_NORMAL_DISTRIBUTION')
          qceff_table_data(row)%probit_inflation%dist_type = LOG_NORMAL_DISTRIBUTION
       case ('UNIFORM_DISTRIBUTION')
@@ -242,9 +252,19 @@ subroutine read_qceff_table(qceff_table_filename)
       case ('BOUNDED_NORMAL_RH_DISTRIBUTION')
          qceff_table_data(row)%probit_state%dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION
       case ('GAMMA_DISTRIBUTION')
+         ! Force standard Gamma distribution
          qceff_table_data(row)%probit_state%dist_type = GAMMA_DISTRIBUTION
+         qceff_table_data(row)%probit_state%bounded_above = .false.
+         qceff_table_data(row)%probit_state%bounded_below = .true.
+         qceff_table_data(row)%probit_state%upper_bound = MISSING_R8
+         qceff_table_data(row)%probit_state%lower_bound = 0.0_r8
       case ('BETA_DISTRIBUTION')
+         ! Force standard Beta distribution
          qceff_table_data(row)%probit_state%dist_type = BETA_DISTRIBUTION
+         qceff_table_data(row)%probit_state%bounded_above = .true.
+         qceff_table_data(row)%probit_state%bounded_below = .true.
+         qceff_table_data(row)%probit_state%upper_bound = 1.0_r8
+         qceff_table_data(row)%probit_state%lower_bound = 0.0_r8
       case ('LOG_NORMAL_DISTRIBUTION')
          qceff_table_data(row)%probit_state%dist_type = LOG_NORMAL_DISTRIBUTION
       case ('UNIFORM_DISTRIBUTION')
@@ -266,9 +286,19 @@ subroutine read_qceff_table(qceff_table_filename)
       case ('BOUNDED_NORMAL_RH_DISTRIBUTION')
          qceff_table_data(row)%probit_extended_state%dist_type = BOUNDED_NORMAL_RH_DISTRIBUTION
       case ('GAMMA_DISTRIBUTION')
+         ! Force standard Gamma distribution
          qceff_table_data(row)%probit_extended_state%dist_type = GAMMA_DISTRIBUTION
+         qceff_table_data(row)%probit_extended_state%bounded_above = .false.
+         qceff_table_data(row)%probit_extended_state%bounded_below = .true.
+         qceff_table_data(row)%probit_extended_state%upper_bound = MISSING_R8
+         qceff_table_data(row)%probit_extended_state%lower_bound = 0.0_r8
       case ('BETA_DISTRIBUTION')
+         ! Force standard Beta distribution
          qceff_table_data(row)%probit_extended_state%dist_type = BETA_DISTRIBUTION
+         qceff_table_data(row)%probit_extended_state%bounded_above = .true.
+         qceff_table_data(row)%probit_extended_state%bounded_below = .true.
+         qceff_table_data(row)%probit_extended_state%upper_bound = 1.0_r8
+         qceff_table_data(row)%probit_extended_state%lower_bound = 0.0_r8
       case ('LOG_NORMAL_DISTRIBUTION')
          qceff_table_data(row)%probit_extended_state%dist_type = LOG_NORMAL_DISTRIBUTION
       case ('UNIFORM_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/beta_distribution_mod.f90

@@ -4,6 +4,9 @@
 
 ! Thanks to Chris Riedel who developed the methods in this module.
 
+! This module only supports standard beta distributions for which the 
+! lower bound is 0 and the upper bound is 1. 
+
 module beta_distribution_mod
 
 use types_mod,      only : r8, PI, missing_r8
@@ -12,7 +15,7 @@ module beta_distribution_mod
 
 use random_seq_mod, only : random_seq_type, random_uniform
 
-use distribution_params_mod, only : distribution_params_type
+use distribution_params_mod, only : distribution_params_type, BETA_DISTRIBUTION
 
 use normal_distribution_mod, only : inv_cdf
 
@@ -36,11 +39,11 @@ subroutine test_beta
 
 ! This routine provides limited tests of the numerics in this module. It begins
 ! by comparing a handful of cases of the pdf and cdf to results from Matlab. It
-! then tests the quality of the inverse cdf for a single shape/scale pair. Failing
+! then tests the quality of the inverse cdf for a single alpha/beta pair. Failing
 ! these tests suggests a serious problem. Passing them does not indicate that 
 ! there are acceptable results for all possible inputs. 
 
-real(r8) :: x, y, p, inv
+real(r8) :: x, y, inv
 real(r8) :: alpha, beta, max_diff
 integer :: i
 
@@ -61,41 +64,37 @@ subroutine test_beta
 ! Compare to matlab
 write(*, *) 'Absolute value of differences should be less than 1e-15'
 do i = 1, 7
-   pdf_diff(i) = beta_pdf(mx(i), malpha(i), mbeta(i)) - mpdf(i)
-   cdf_diff(i) = beta_cdf(mx(i), malpha(i), mbeta(i), 0.0_r8, 1.0_r8) - mcdf(i)
+   pdf_diff(i) = abs(beta_pdf(mx(i), malpha(i), mbeta(i)) - mpdf(i))
+   cdf_diff(i) = abs(beta_cdf(mx(i), malpha(i), mbeta(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 Comparison Tests: FAIL'
 endif
 
-
 ! Test many x values for cdf and inverse cdf for a single set of alpha and beta
 alpha = 5.0_r8
 beta = 2.0_r8
 
 max_diff = -1.0_r8
 do i = 0, 1000
    x = i / 1000.0_r8
-   p = beta_pdf(x, alpha, beta)
-   y = beta_cdf(x, alpha, beta, 0.0_r8, 1.0_r8)
-   inv = inv_beta_cdf(y, alpha, beta, 0.0_r8, 1.0_r8)
+   y = beta_cdf(x, alpha, beta)
+   inv = inv_beta_cdf(y, alpha, beta)
    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-14'
-
 if(max_diff < 1e-14_r8) then
    write(*, *) 'Inversion Tests: PASS'
 else
    write(*, *) 'Inversion Tests: FAIL'
 endif
 
-
 end subroutine test_beta
 
 !-----------------------------------------------------------------------
@@ -106,20 +105,25 @@ function inv_beta_cdf_params(quantile, p) result(x)
 real(r8), intent(in)                       :: quantile
 type(distribution_params_type), intent(in) :: p
 
+! Only standard beta is currently supported. Fail if bounds are not 0 and 1
+if(p%lower_bound /= 0.0_r8 .or. p%upper_bound /= 1.0_r8) then
+  errstring = 'Only standard beta distribution with bounds of 0 and 1 is supported'
+  call error_handler(E_ERR, 'inv_beta_cdf_params', errstring, source)
+endif
+
 x = inv_cdf(quantile, beta_cdf_params, inv_beta_first_guess_params, p)
 
 end function inv_beta_cdf_params
 
 !-----------------------------------------------------------------------
 
-function inv_beta_cdf(quantile, alpha, beta, lower_bound, upper_bound) result(x)
+function inv_beta_cdf(quantile, alpha, beta) result(x)
 
 real(r8)             :: x
 real(r8), intent(in) :: quantile
 real(r8), intent(in) :: alpha, beta
-real(r8), intent(in) :: lower_bound,   upper_bound
 
-! Given a quantile, finds the value of x for which the scaled beta cdf
+! Given a quantile, finds the value of x for which the beta cdf
 ! with alpha and beta has approximately this quantile
 
 type(distribution_params_type) :: p
@@ -129,15 +133,13 @@ function inv_beta_cdf(quantile, alpha, beta, lower_bound, upper_bound) result(x)
   call error_handler(E_ERR, 'inv_beta_cdf', errstring, source)
 endif
 
+! Load the p type for the generic cdf calls
 p%params(1) = alpha;  p%params(2) = beta
-! Beta must be bounded on both sides
-p%lower_bound = lower_bound;       p%upper_bound = upper_bound
+p%bounded_below = .true.;          p%bounded_above = .true.
+p%lower_bound   = 0.0_r8;          p%upper_bound   = 1.0_r8
 
 x = inv_beta_cdf_params(quantile, p)
 
-! Undo the scaling 
-x = x * (upper_bound - lower_bound) + lower_bound
-
 end function inv_beta_cdf
 
 !---------------------------------------------------------------------------
@@ -184,25 +186,34 @@ function beta_cdf_params(x, p)
 real(r8), intent(in)                       :: x
 type(distribution_params_type), intent(in) :: p
 
+! A translation routine that is required to use the generic cdf optimization routine
+! Extracts the appropriate information from the distribution_params_type that is needed
+! for a call to the function beta_cdf below. 
+
 real(r8) :: alpha, beta
 
+! Only standard beta is currently supported. Fail if bounds are not 0 and 1
+if(p%lower_bound /= 0.0_r8 .or. p%upper_bound /= 1.0_r8) then
+  errstring = 'Only standard beta distribution with bounds of 0 and 1 is supported'
+  call error_handler(E_ERR, 'beta_cdf_params', errstring, source)
+endif
+
 alpha = p%params(1);    beta = p%params(2)
-beta_cdf_params = beta_cdf(x, alpha, beta, p%lower_bound, p%upper_bound)
+beta_cdf_params = beta_cdf(x, alpha, beta)
 
 end function beta_cdf_params
 
 !---------------------------------------------------------------------------
 
-function beta_cdf(x, alpha, beta, lower_bound, upper_bound)
+function beta_cdf(x, alpha, beta)
 
 ! Returns the cumulative distribution of a beta function with alpha and beta
 ! at the value x
 
-! Returns a large negative value if called with illegal values
+! Returns a failed_value if called with illegal values
 
 real(r8) :: beta_cdf
 real(r8), intent(in) :: x, alpha, beta
-real(r8), intent(in) :: lower_bound,   upper_bound
 
 ! Parameters must be positive
 if(alpha <= 0.0_r8 .or. beta <= 0.0_r8) then
@@ -251,7 +262,7 @@ function random_beta(r, alpha, beta)
 ! Draw from U(0, 1) to get a quantile
 quantile = random_uniform(r)
 ! Invert cdf to get a draw from beta
-random_beta = inv_beta_cdf(quantile, alpha, beta, 0.0_r8, 1.0_r8)
+random_beta = inv_beta_cdf(quantile, alpha, beta)
 
 end function random_beta
 
@@ -339,24 +350,25 @@ function inv_beta_first_guess_params(quantile, p)
 real(r8), intent(in)                       :: quantile
 type(distribution_params_type), intent(in) :: 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 inv_beta_first_guess below.
+
 real(r8) :: alpha, beta
 
 alpha = p%params(1);     beta = p%params(2)
-inv_beta_first_guess_params = inv_beta_first_guess(quantile, alpha, beta, &
-   p%bounded_below, p%bounded_above, p%lower_bound, p%upper_bound)
+inv_beta_first_guess_params = inv_beta_first_guess(quantile, alpha, beta)
 
 end function inv_beta_first_guess_params
 
 !---------------------------------------------------------------------------
 
-function inv_beta_first_guess(x, alpha, beta, &
-   bounded_below, bounded_above, lower_bound, upper_bound)
+function inv_beta_first_guess(x, alpha, beta)
 
 real(r8) :: inv_beta_first_guess
 real(r8), intent(in) :: x
 real(r8), intent(in) :: alpha, beta
-logical,  intent(in) :: bounded_below, bounded_above
-real(r8), intent(in) :: lower_bound,   upper_bound
 
 ! Need some sort of first guess, should be smarter here
 ! For starters, take the mean for this alpha and beta
@@ -389,19 +401,22 @@ end subroutine beta_alpha_beta
 
 !---------------------------------------------------------------------------
 
-subroutine set_beta_params_from_ens(ens, num, lower_bound, upper_bound, p)
+subroutine set_beta_params_from_ens(ens, num, p)
 
-integer,                        intent(in)    :: num
-real(r8),                       intent(in)    :: ens(num)
-real(r8),                       intent(in)    :: lower_bound,   upper_bound
-type(distribution_params_type), intent(inout) :: p
+integer,  intent(in)                        :: num
+real(r8), intent(in)                        :: ens(num)
+type(distribution_params_type), intent(out) :: p
 
 real(r8) :: alpha, beta
 
+! Set up the description of the beta distribution defined by the ensemble
+p%distribution_type = BETA_DISTRIBUTION
+
 ! Set the bounds info
-p%lower_bound   = lower_bound;    p%upper_bound   = upper_bound
+p%bounded_below = .true.;         p%bounded_above = .true.
+p%lower_bound   = 0.0_r8;         p%upper_bound   = 1.0_r8
 
-! Get alpha and beta for the scaled ensemble
+! Get alpha and beta for the ensemble
 call beta_alpha_beta(ens, num, alpha, beta)
 p%params(1) = alpha
 p%params(2) = beta

assimilation_code/modules/assimilation/distribution_params_mod.f90

@@ -7,17 +7,8 @@ module distribution_params_mod
 implicit none
 private
 
-type distribution_params_type
-   integer               :: distribution_type
-   logical               :: bounded_below, bounded_above
-   real(r8)              :: lower_bound,   upper_bound
-   real(r8)              :: params(2)
-   integer               :: ens_size
-   real(r8), allocatable :: ens(:)
-   real(r8), allocatable :: more_params(:)
-end type
-
 ! Defining parameter strings for different prior distributions that can be used for probit transform
+integer, parameter :: UNSET = -1
 integer, parameter :: NORMAL_DISTRIBUTION            = 1
 integer, parameter :: BOUNDED_NORMAL_RH_DISTRIBUTION = 2
 integer, parameter :: GAMMA_DISTRIBUTION             = 3
@@ -27,6 +18,16 @@ module distribution_params_mod
 integer, parameter :: PARTICLE_FILTER_DISTRIBUTION   = 7
 integer, parameter :: KDE_DISTRIBUTION               = 8
 
+type distribution_params_type
+   integer               :: distribution_type = UNSET
+   logical               :: bounded_below, bounded_above
+   real(r8)              :: lower_bound, upper_bound
+   real(r8)              :: params(2)
+   integer               :: ens_size
+   real(r8), allocatable :: ens(:)
+   real(r8), allocatable :: more_params(:)
+end type
+
 public :: distribution_params_type, deallocate_distribution_params, &
    NORMAL_DISTRIBUTION, BOUNDED_NORMAL_RH_DISTRIBUTION, GAMMA_DISTRIBUTION, BETA_DISTRIBUTION, &
       LOG_NORMAL_DISTRIBUTION, UNIFORM_DISTRIBUTION, PARTICLE_FILTER_DISTRIBUTION, &