Refactor adaptive forgetting factor in PDAF_set_forget / PDAF_set_forget_local

docs/users_guide/running_tsmp_pdaf/input_cmd.md

@@ -102,30 +102,71 @@ 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,
 <http://dx.doi.org/10.3402/tellusa.v64i0.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.
 
 For ETKF, see
 e.g. <https://github.com/PDAF/PDAF/blob/ae9545227bd4804469dff389a9baadcc9e31906e/src/PDAF_etkf_analysis.F90#L441-L444>
 
+**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.

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

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