Merge pull request #706 from iangrooms/kqcef

KQCEF

.gitignore

@@ -199,6 +199,7 @@ test_quad_irreg_interp
 test_quad_reg_interp
 test_table_read
 test_ran_unif
+test_kde_dist
 
 # 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,14 @@ individual files.
 
 The changes are now listed with the most recent at the top.
 
+**August 26 2024 :: KQCEF. Tag 11.7.0**
+
+- Adds a Quantile-Conserving Ensemble Filter Based on Kernel-Density Estimation to DART.
+- New distribution module kde_distribution_mod.
+- New module rootfinding_mod that provides a different implementation of inv_cdf.
+
+*Contributed by Ian Grooms*
+
 **August 15 2024 :: WRF fwd operator bug fixes. Tag v11.6.1**
 
 WRF-DART bug-fixes:

assimilation_code/modules/assimilation/algorithm_info_mod.f90

@@ -22,7 +22,7 @@ module algorithm_info_mod
 
 use distribution_params_mod, only : NORMAL_DISTRIBUTION, BOUNDED_NORMAL_RH_DISTRIBUTION, &
    GAMMA_DISTRIBUTION, BETA_DISTRIBUTION, LOG_NORMAL_DISTRIBUTION, UNIFORM_DISTRIBUTION,  &
-   PARTICLE_FILTER_DISTRIBUTION
+   PARTICLE_FILTER_DISTRIBUTION, KDE_DISTRIBUTION
 
 implicit none
 private
@@ -40,12 +40,13 @@ module algorithm_info_mod
 integer, parameter :: OBS_PARTICLE       = 4
 integer, parameter :: UNBOUNDED_RHF      = 8
 integer, parameter :: GAMMA_FILTER       = 11
-integer, parameter :: BOUNDED_NORMAL_RHF = 101 
+integer, parameter :: BOUNDED_NORMAL_RHF = 101
+integer, parameter :: KDE_FILTER        = 102
 
 public :: obs_error_info, probit_dist_info, obs_inc_info, &
           init_algorithm_info_mod, end_algorithm_info_mod, &
           EAKF, ENKF, BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, &
-          GAMMA_FILTER, KERNEL, OBS_PARTICLE
+          GAMMA_FILTER, KERNEL, OBS_PARTICLE, KDE_FILTER
 
 ! type definitions for the QCF table
 type obs_error_info_type
@@ -223,6 +224,8 @@ subroutine read_qceff_table(qceff_table_filename)
          qceff_table_data(row)%probit_inflation%dist_type = LOG_NORMAL_DISTRIBUTION
       case ('UNIFORM_DISTRIBUTION')
          qceff_table_data(row)%probit_inflation%dist_type = UNIFORM_DISTRIBUTION
+      case ('KDE_DISTRIBUTION')
+         qceff_table_data(row)%probit_inflation%dist_type = KDE_DISTRIBUTION
       !!!case ('PARTICLE_FILTER_DISTRIBUTION')
          !!!qceff_table_data(row)%probit_inflation%dist_type = PARTICLE_FILTER_DISTRIBUTION
       case default
@@ -246,6 +249,8 @@ subroutine read_qceff_table(qceff_table_filename)
          qceff_table_data(row)%probit_state%dist_type = LOG_NORMAL_DISTRIBUTION
       case ('UNIFORM_DISTRIBUTION')
          qceff_table_data(row)%probit_state%dist_type = UNIFORM_DISTRIBUTION
+      case ('KDE_DISTRIBUTION')
+         qceff_table_data(row)%probit_state%dist_type = KDE_DISTRIBUTION
       !!!case ('PARTICLE_FILTER_DISTRIBUTION')
          !!!qceff_table_data(row)%probit_state%dist_type = PARTICLE_FILTER_DISTRIBUTION
       case default
@@ -268,6 +273,8 @@ subroutine read_qceff_table(qceff_table_filename)
          qceff_table_data(row)%probit_extended_state%dist_type = LOG_NORMAL_DISTRIBUTION
       case ('UNIFORM_DISTRIBUTION')
          qceff_table_data(row)%probit_extended_state%dist_type = UNIFORM_DISTRIBUTION
+      case ('KDE_DISTRIBUTION')
+         qceff_table_data(row)%probit_extended_state%dist_type = KDE_DISTRIBUTION
       !!!case ('PARTICLE_FILTER_DISTRIBUTION')
          !!!qceff_table_data(row)%probit_extended_state%dist_type = PARTICLE_FILTER_DISTRIBUTION
       case default
@@ -290,6 +297,8 @@ subroutine read_qceff_table(qceff_table_filename)
          qceff_table_data(row)%obs_inc_info%filter_kind = GAMMA_FILTER
       case ('BOUNDED_NORMAL_RHF')
          qceff_table_data(row)%obs_inc_info%filter_kind = BOUNDED_NORMAL_RHF
+      case ('KDE_FILTER')
+         qceff_table_data(row)%obs_inc_info%filter_kind = KDE_FILTER
       case default
          write(errstring, *) 'Invalid filter kind: ', trim(filter_kind_string(row))
          call error_handler(E_ERR, 'read_qceff_table:', errstring, source)

assimilation_code/modules/assimilation/assim_tools_mod.f90

@@ -78,13 +78,17 @@ module assim_tools_mod
 
 use algorithm_info_mod, only : probit_dist_info, obs_inc_info, EAKF, ENKF, &
                                BOUNDED_NORMAL_RHF, UNBOUNDED_RHF, GAMMA_FILTER, &
-                               KERNEL, OBS_PARTICLE
+                               KERNEL, OBS_PARTICLE, KDE_FILTER
 
 use gamma_distribution_mod, only : gamma_cdf, inv_gamma_cdf, gamma_mn_var_to_shape_scale, &
                                    gamma_gamma_prod
 
 use bnrh_distribution_mod, only   :  inv_bnrh_cdf, bnrh_cdf, inv_bnrh_cdf_like
 
+use kde_distribution_mod, only : kde_cdf_params, inv_kde_cdf_params, obs_dist_types,      &
+                                 pack_kde_params, likelihood_function, separate_ensemble, &
+                                 obs_increment_kde
+
 use distribution_params_mod, only : distribution_params_type, deallocate_distribution_params
 
 
@@ -1025,6 +1029,10 @@ subroutine obs_increment(ens_in, ens_size, obs, obs_var, obs_kind, obs_inc, &
       !!!endif
 
    !--------------------------------------------------------------------------
+
+   else if(filter_kind == KDE_FILTER) then
+      call obs_increment_kde(ens, ens_size, obs, obs_var, bounded_below, &
+         bounded_above, lower_bound, upper_bound, obs_inc)
    else
       call error_handler(E_ERR,'obs_increment', &
               'Illegal value of filter_kind', source)
@@ -1166,9 +1174,6 @@ subroutine obs_increment_bounded_norm_rhf(ens, ens_like, ens_size, prior_var, &
 
 end subroutine obs_increment_bounded_norm_rhf
 
-
-
-
 ! Computes a normal or truncated normal (above and/or below) likelihood.
 function get_truncated_normal_like(x, obs, obs_var, &
    bounded_below, bounded_above, lower_bound, upper_bound)
@@ -1420,6 +1425,8 @@ subroutine obs_increment_kernel(ens, ens_size, obs, obs_var, obs_inc)
 endif
 
 ! Generate a uniform random number and a Gaussian for each new member
+!! IG: This could be done with a single random number and in a way that
+!! preserves the order of the ensemble members
 do i = 1, ens_size
    unif = random_uniform(inc_ran_seq)
    ! Figure out which kernel it's in

assimilation_code/modules/assimilation/distribution_params_mod.f90

@@ -25,10 +25,12 @@ module distribution_params_mod
 integer, parameter :: LOG_NORMAL_DISTRIBUTION        = 5
 integer, parameter :: UNIFORM_DISTRIBUTION           = 6
 integer, parameter :: PARTICLE_FILTER_DISTRIBUTION   = 7
+integer, parameter :: KDE_DISTRIBUTION               = 8
 
 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
+      LOG_NORMAL_DISTRIBUTION, UNIFORM_DISTRIBUTION, PARTICLE_FILTER_DISTRIBUTION, &
+      KDE_DISTRIBUTION
 
 contains