Skip to content

Commit

Permalink
Merge pull request #316 from NCAR/filter-assim-refactor
Browse files Browse the repository at this point in the history 
Filter assim refactor
  • Loading branch information
@hkershaw-brown
hkershaw-brown committed Dec 8, 2021
2 parents b907819 + 5ec0345 commit 83b239f
Show file tree
Hide file tree
Showing 4 changed files with 493 additions and 720 deletions.
7 changes: 6 additions & 1 deletion CHANGELOG.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,12 @@ individual files.

The changes are now listed with the most recent at the top.

**November 22 2021 :: Bug fix for groups with posterior spatially-varying adaptive inflation. Tag: v9.12.13**
**December 7 2021 :: Refactored filter_assim. Tag: v9.12.0**

- Filter_assim refactored so each process calcuates increments
- Code readability changes

**November 22 2021 :: Bug fix for groups with posterior spatially-varying adaptive inflation. Tag: v9.11.13**

- Removed the additional outlier threshold test for each group when using posterior
spatially-varying adaptive inflation. The outlier test is done for the entire ensemble
Expand Down
133 changes: 123 additions & 10 deletions assimilation_code/modules/assimilation/adaptive_inflate_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,17 +19,18 @@ module adaptive_inflate_mod
implicit none
private

public :: update_inflation, do_obs_inflate, &
do_varying_ss_inflate, do_single_ss_inflate, inflate_ens, &
adaptive_inflate_init, adaptive_inflate_type, &
deterministic_inflate, solve_quadratic, &
log_inflation_info, get_minmax_task_zero, mean_from_restart, &
sd_from_restart, &
output_inf_restart, get_inflate_mean, get_inflate_sd, &
get_is_prior, get_is_posterior, do_ss_inflate, &
set_inflation_mean_copy, set_inflation_sd_copy, get_inflation_mean_copy, &
public :: update_inflation, do_obs_inflate, &
update_single_state_space_inflation, update_varying_state_space_inflation, &
do_varying_ss_inflate, do_single_ss_inflate, inflate_ens, &
adaptive_inflate_init, adaptive_inflate_type, &
deterministic_inflate, solve_quadratic, &
log_inflation_info, get_minmax_task_zero, mean_from_restart, &
sd_from_restart, &
output_inf_restart, get_inflate_mean, get_inflate_sd, &
get_is_prior, get_is_posterior, do_ss_inflate, &
set_inflation_mean_copy, set_inflation_sd_copy, get_inflation_mean_copy, &
get_inflation_sd_copy, do_rtps_inflate, validate_inflate_options, &
print_inflation_restart_filename, &
print_inflation_restart_filename, &
PRIOR_INF, POSTERIOR_INF, NO_INFLATION, OBS_INFLATION, VARYING_SS_INFLATION, &
SINGLE_SS_INFLATION, RELAXATION_TO_PRIOR_SPREAD, ENHANCED_SS_INFLATION

Expand Down Expand Up @@ -95,6 +96,10 @@ module adaptive_inflate_mod
! Flag indicating whether module has been initialized
logical :: initialized = .false.

! Used for precision tests in inflation update routines
real(r8), parameter :: small = epsilon(1.0_r8) ! threshold for avoiding NaNs/Inf


!===============================================================================

contains
Expand Down Expand Up @@ -645,6 +650,114 @@ subroutine update_inflation(inflate_handle, inflate, inflate_sd, prior_mean, pri

end subroutine update_inflation

!-------------------------------------------------------------------------------
!> Computes updated inflation mean and inflation sd for single state space inflation

subroutine update_single_state_space_inflation(inflate, inflate_mean, inflate_sd, &
ss_inflate_base, orig_obs_prior_mean, orig_obs_prior_var, obs, obs_err_var, &
ens_size, inflate_only)

type(adaptive_inflate_type), intent(in) :: inflate
real(r8), intent(inout) :: inflate_mean
real(r8), intent(inout) :: inflate_sd
real(r8), intent(in) :: ss_inflate_base
real(r8), intent(in) :: orig_obs_prior_mean
real(r8), intent(in) :: orig_obs_prior_var
real(r8), intent(in) :: obs
real(r8), intent(in) :: obs_err_var
integer, intent(in) :: ens_size
logical, intent(in) :: inflate_only

real(r8) :: gamma, ens_var_deflate, r_var, r_mean

! If either inflation or sd is not positive, not really doing inflation
if(inflate_mean <= 0.0_r8 .or. inflate_sd <= 0.0_r8) return

! For case with single spatial inflation, use gamma = 1.0_r8
gamma = 1.0_r8
! Deflate the inflated variance; required for efficient single pass
! This is one of many places that assumes linear state/obs relation
! over range of ensemble; Essentially, we are removing the inflation
! which has already been applied in filter to see what inflation should
! have been needed.
ens_var_deflate = orig_obs_prior_var / &
(1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2

! If this is inflate_only (i.e. posterior) remove impact of this obs.
! This is simulating independent observation by removing its impact.
if(inflate_only .and. &
ens_var_deflate > small .and. &
obs_err_var > small .and. &
obs_err_var - ens_var_deflate > small ) then
r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var)
r_mean = r_var *(orig_obs_prior_mean / ens_var_deflate - obs / obs_err_var)
else
r_var = ens_var_deflate
r_mean = orig_obs_prior_mean
endif

! Update the inflation mean value and standard deviation
call update_inflation(inflate, inflate_mean, inflate_sd, &
r_mean, r_var, ens_size, obs, obs_err_var, gamma)

end subroutine update_single_state_space_inflation

!-------------------------------------------------------------------------------
!> Computes updated inflation mean and inflation sd for varying state space inflation

subroutine update_varying_state_space_inflation(inflate, inflate_mean, inflate_sd, &
ss_inflate_base, orig_obs_prior_mean, orig_obs_prior_var, obs, obs_err_var, &
ens_size, reg_factor, correl, inflate_only)

type(adaptive_inflate_type), intent(in) :: inflate
real(r8), intent(inout) :: inflate_mean
real(r8), intent(inout) :: inflate_sd
real(r8), intent(in) :: ss_inflate_base
real(r8), intent(in) :: orig_obs_prior_mean
real(r8), intent(in) :: orig_obs_prior_var
real(r8), intent(in) :: obs
real(r8), intent(in) :: obs_err_var
integer, intent(in) :: ens_size
real(r8), intent(in) :: reg_factor
real(r8), intent(in) :: correl
logical, intent(in) :: inflate_only

real(r8) :: gamma, ens_var_deflate, r_var, r_mean
real(r8) :: diff_sd, outlier_ratio
logical :: do_adapt_inf_update

if(inflate_mean <= 0.0_r8 .or. inflate_sd <= 0.0_r8) return

! Gamma is less than 1 for varying ss, see adaptive inflate module
gamma = reg_factor * abs(correl)

! Remove the impact of inflation to allow efficient single pass with assim.
if ( abs(gamma) > small ) then
ens_var_deflate = orig_obs_prior_var / &
(1.0_r8 + gamma*(sqrt(ss_inflate_base) - 1.0_r8))**2
else
ens_var_deflate = orig_obs_prior_var
endif

! If this is inflate only (i.e. posterior) remove impact of this obs.
if(inflate_only .and. &
ens_var_deflate > small .and. &
obs_err_var > small .and. &
obs_err_var - ens_var_deflate > small ) then
r_var = 1.0_r8 / (1.0_r8 / ens_var_deflate - 1.0_r8 / obs_err_var)
r_mean = r_var *(orig_obs_prior_mean / ens_var_deflate - obs / obs_err_var)
else
r_var = ens_var_deflate
r_mean = orig_obs_prior_mean
endif

! IS A TABLE LOOKUP POSSIBLE TO ACCELERATE THIS?
! Update the inflation values
call update_inflation(inflate, inflate_mean, inflate_sd, &
r_mean, r_var, ens_size, obs, obs_err_var, gamma)

end subroutine update_varying_state_space_inflation


!-------------------------------------------------------------------------------
!> Uses one of 2 algorithms in references on DART web site to update the
Expand Down

0 comments on commit 83b239f

Please sign in to comment.