From a48daa2a52cc1b3ea56caf7d1e4146ef46692db0 Mon Sep 17 00:00:00 2001 From: Yorck Ewerdwalbesloh Date: Thu, 30 Apr 2026 16:13:58 +0200 Subject: [PATCH 1/2] src: PDAF_set_forget(_local) adapted --- src/PDAF_set_forget.F90 | 34 ++++++++++++++++++++++++++++++---- src/PDAF_set_forget_local.F90 | 28 ++++++++++++++++++++++++---- 2 files changed, 54 insertions(+), 8 deletions(-) diff --git a/src/PDAF_set_forget.F90 b/src/PDAF_set_forget.F90 index 5182bb7d7..1765f12a3 100644 --- a/src/PDAF_set_forget.F90 +++ b/src/PDAF_set_forget.F90 @@ -85,17 +85,24 @@ SUBROUTINE PDAF_set_forget(step, filterstr, dim_obs_p, dim_ens, mens_p, & REAL :: var_resid_p, var_resid ! Variance of residual REAL :: var_obs ! Variance of observation errors REAL :: forget_neg, forget_max, forget_min ! Limiting values of forgetting factor + REAL :: var_target + REAL :: spread_fac ! ********************** ! *** INITIALIZATION *** ! ********************** + ! spread_fac defines the minimum ensemble std relative to obs std: + ! sqrt(var_ens_target) = spread_fac * sqrt(var_obs) + ! hence: var_target = spread_fac^2 * var_obs + spread_fac = 1.0 + ! Define limiting values of forgetting factor - ! These are set very arbitrarily for now + ! The input 'forget' is used as obs-type dependent lower bound forget_neg = forget_in - forget_max = 100.0 - forget_min = 0.01 + forget_max = 1.0 + forget_min = forget_in IF (mype == 0) THEN WRITE (*, '(a, 5x, a)') & @@ -180,8 +187,19 @@ SUBROUTINE PDAF_set_forget(step, filterstr, dim_obs_p, dim_ens, mens_p, & CALL PDAF_timeit(51, 'new') + ! Define target ensemble variance in observation space + var_target = spread_fac * spread_fac * var_obs + ! *** Compute optimal forgetting factor *** - forget_out = var_ens / (var_resid - var_obs) + IF (var_target > 0.0) THEN + IF (var_ens < var_target) THEN + forget_out = var_ens / var_target + ELSE + forget_out = 1.0 + ENDIF + ELSE + forget_out = 1.0 + ENDIF ! Apply special condition if observation variance is larger than residual variance IF (forget_out < 0.0) forget_out = forget_neg @@ -198,6 +216,14 @@ SUBROUTINE PDAF_set_forget(step, filterstr, dim_obs_p, dim_ens, mens_p, & ! ******************** IF (mype == 0) THEN +#ifdef PDAF_DEBUG + WRITE (*, '(a, 9x, a, es10.2)') & + 'PDAF', 'Variance of ensemble', var_ens + WRITE (*, '(a, 9x, a, es10.2)') & + 'PDAF', 'Variance of observation errors', var_obs + WRITE (*, '(a, 9x, a, es10.2)') & + 'PDAF', 'Target variance of ensemble', var_target +#endif WRITE (*, '(a, 9x, a, es10.2)') & 'PDAF', '--> Computed forgetting factor', forget_out ENDIF diff --git a/src/PDAF_set_forget_local.F90 b/src/PDAF_set_forget_local.F90 index 06a595a9c..9f1d920f4 100644 --- a/src/PDAF_set_forget_local.F90 +++ b/src/PDAF_set_forget_local.F90 @@ -77,17 +77,26 @@ SUBROUTINE PDAF_set_forget_local(domain, step, dim_obs_l, dim_ens, mens_l, & INTEGER, SAVE :: domain_save = 1 ! Index of domain from last call to routine INTEGER :: n_domains ! Number of local analysis domains REAL :: forget_neg, forget_max, forget_min ! limiting values of forgetting factor + REAL :: var_target + REAL :: spread_fac ! ********************** ! *** INITIALIZATION *** ! ********************** + ! spread_fac defines the minimum ensemble std relative to obs std: + ! sqrt(var_ens_target) = spread_fac * sqrt(var_obs) + ! hence: var_target = spread_fac^2 * var_obs + spread_fac = 1.0 + ! Define limiting values of forgetting factor - ! These are set very arbitrarily for now + ! The input 'forget' is used as obs-type dependent lower bound forget_neg = forget - forget_max = 100.0 - forget_min = 0.01 + forget_max = 1.0 + forget_min = forget + + ! **************************************************** @@ -150,8 +159,19 @@ SUBROUTINE PDAF_set_forget_local(domain, step, dim_obs_l, dim_ens, mens_l, & CALL U_init_obsvar_l(domain, step, dim_obs_l, obs_l, var_obs) CALL PDAF_timeit(52, 'old') + ! Define target ensemble variance in local observation space + var_target = spread_fac * spread_fac * var_obs + ! *** Compute optimal forgetting factor *** - forget = var_ens / (var_resid - var_obs) + IF (var_target > 0.0) THEN + IF (var_ens < var_target .and. dim_obs_l>3) THEN + forget = var_ens / var_target + ELSE + forget = 1.0 + ENDIF + ELSE + forget = 1.0 + ENDIF ! Apply special condition if observation variance is larger than residual variance IF (forget < 0.0) forget = forget_neg From ed59737ead7374a10ee1b5e5eb014b5d2a8b4422 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 30 Apr 2026 16:14:19 +0200 Subject: [PATCH 2/2] docs: document command line input `type_forget` and adapt documentation of command line input `forget` according to the new adaptive forgetting factor routines. --- .../running_tsmp_pdaf/input_cmd.md | 59 ++++++++++++++++--- 1 file changed, 50 insertions(+), 9 deletions(-) diff --git a/docs/users_guide/running_tsmp_pdaf/input_cmd.md b/docs/users_guide/running_tsmp_pdaf/input_cmd.md index 693570fba..70058b49f 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_cmd.md +++ b/docs/users_guide/running_tsmp_pdaf/input_cmd.md @@ -102,23 +102,40 @@ For flexible time stepping, `delt_obs` must be one. - 2: 1 plus timing output - 3: 2 plus debug output +## type_forget ## + +`type_forget` (integer) Type of forgetting factor. Default: `0`. + +For SEIK / LSEIK / ETKF / LETKF / ESTKF / LESTKF: +- `0`: fixed forgetting factor (value given by `forget`) +- `1`: global adaptive forgetting factor (computed from ensemble + spread and observation error variance; `forget` acts as lower bound) +- `2`: local adaptive forgetting factor, for LSEIK / LETKF / LESTKF + only (same adaptive logic applied per local analysis domain; `forget` + acts as lower bound) + +For NETF / LNETF / PF: +- `0`: apply inflation on forecast ensemble +- `2`: apply inflation on analysis ensemble + ## forget ## -`forget` (real) forgetting factor for filter analysis +`forget` (real) forgetting factor for filter analysis. Default: `1.0`. -Example: `-forget 0.98`. +**Fixed mode (`type_forget = 0`):** `forget` is applied directly as the +forgetting factor. Example: `-forget 0.98`. -General advise: Choose forgetting factor close to one. For values -smaller than 0.95, effects like a splitting of the ensemble have been -observed (compare Amezcua et al., Tellus A 2012, 64, 18039, +General advice: choose `forget` close to one. For values smaller than +0.95, effects like a splitting of the ensemble have been observed +(compare Amezcua et al., Tellus A 2012, 64, 18039, ) For EnKF / LEnKF, the forgetting factor leads to a spreading of the -ensemble through manipulating ensemble member by +ensemble through manipulating each ensemble member by -\begin{align*} -x^{f}_{i} &= \bar{x} + (x_{i}-\bar{x}) \cdot \frac{1}{\mathtt{forget}^2}, -\end{align*} +$$ +x^{f}_{i} = \bar{x} + (x_{i}-\bar{x}) \cdot \frac{1}{\mathtt{forget}^2}, +$$ where $x_{i}$ is the state vector ensemble member $i$ and $\bar{x}$ is the ensemble mean of the state vector. @@ -126,6 +143,30 @@ the ensemble mean of the state vector. For ETKF, see e.g. +**Adaptive mode (`type_forget = 1` or `2`):** `forget` is used as a +lower bound on the computed forgetting factor. The adaptive algorithm +inflates the ensemble only when its spread in observation space falls +below the observation error standard deviation. + +Concretely, the forgetting factor is computed as + +$$ +\mathtt{forget_adaptive} = \begin{cases} +\mathtt{forget} & \text{if } \sigma^2_{\mathrm{ens}} < \mathtt{forget} \cdot \sigma^2_{\mathrm{obs}}, \\[5pt] +\sigma^2_{\mathrm{ens}} / \sigma^2_{\mathrm{obs}} & \text{if } \mathtt{forget} \cdot \sigma^2_{\mathrm{obs}} \leq \sigma^2_{\mathrm{ens}} \leq \sigma^2_{\mathrm{obs}}, \\[5pt] +1 & \text{otherwise.} +\end{cases} +$$ + +Note that the result is clipped to $[{\mathtt{forget}},\, 1]$, so +`forget` provides a lower bound (maximum inflation). + +The resulting `forget` factor is then used in the same way as the +global one. + +A typical choice is `-forget 0.95` to allow moderate inflation while +preventing excessive spread. + ## locweight ## Only for localization.