Changes between Version 2 and Version 3 of ImplementAnalysislknetf

Timestamp:

Feb 19, 2023, 10:28:49 AM (21 months ago)

Author:

lnerger

Comment:

--

ImplementAnalysislknetf

@@ -272 +272 @@
 === `U_prodRinvA_hyb_l` (prodrinva_hyb_l_pdaf.F90) ===
 
-This routine is used by the local hybrid filter LKNETF. 
+This routine is used by the local hybrid filter LKNETF. It is a variant of `U_proRinvA_l` accounting for hybridization.
 
 The interface for this routine is:
@@ -292 +292 @@
 This routine is also the place to perform observation localization. To initialize a vector of weights, the routine `PDAF_local_weight` can be called. The procedure is used in the example implementation and also demonstrated in the template routine.
 
-The routine also has to apply the hybrid weight `gamma`. This is a simply multiplication with the input value in the loop where `C_l` is initialized.
+The routine also has to apply the hybrid weight `gamma`. This is a simple multiplication with the input value in the loop where `C_l` is initialized.
 
 Hints:
  * This routine is a simple extension of `prodRinvA_l. One can implement the hybrid variant by copying this routine and adapting it. `gamma` is computed inside PDAF and provided to the routine. 
+
 
 === `U_init_n_domains` (init_n_domains_pdaf.F90) ===
@@ -471 +472 @@
 
 
+
+=== `U_likelihood_l` (likelihood_l_pdaf.F90) ===
+
+This routine is used by the LNETF and LKNETF filters. 
+
+The interface for this routine is:
+{{{
+SUBROUTINE likelihood_l(domain_p, step, dim_obs_l, obs_l, resid_l, likely_l)
+
+  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
+  INTEGER, INTENT(in) :: step                 ! Current time step
+  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
+  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
+  REAL, INTENT(inout) :: resid_l(dim_obs_l)   ! Input vector holding the local residual y-Hx
+  REAL, INTENT(out)   :: likely_l(dim_obs_l)  ! Output value of the likelihood
+}}}
+
+The routine is called during the loop over the local analysis domains. The likelihood of the local observations has to be computed for each ensemble member. The likelihood is computed from the observation-state residual according to the assumed observation error distribution. Commonly, the observation errors are assumed to be Gaussian distributed. In this case, the likelihood is '''exp(-0.5*(y-Hx)^T^*R^-1^*(y-Hx))'''.
+
+This routine is also the place to perform observation localization. To initialize a vector of weights, the routine `PDAF_local_weight` can be called. The procedure is used in the example implementation and also demonstrated in the template routine.
+
+Hints:
+ * The routine is a local variant of the routine `U_likelihood`. Thus if that routine has been implemented before, it can be adapted here for the local filter.
+ * The routine is very similar to the routine [wiki:U_prodRinvA_l]. The main addition is the computation of the likelihood after computing '''R^-1^*(y-Hx)''', which corresponds to '''R^-1^*A_p''' in [wiki:U_prodRinvA_l].
+ * The information about the inverse observation error covariance matrix has to be provided by the user. Possibilities are to read this information from a file, or to use a Fortran module that holds this information, which one could already prepare in init_pdaf.
+ * The routine does not require that the product is implemented as a real matrix-vector product. Rather, the product can be implemented in its most efficient form. For example, if the observation error covariance matrix is diagonal, only the multiplication of the inverse diagonal with the vector `resid` has to be implemented.
+ * The observation vector `obs_l` is provided through the interface for cases where the observation error variance is relative to the actual value of the observations.
+
+
+
+=== `U_likelihood_hyb_l` (likelihood_hyb_l_pdaf.F90) ===
+
+This routine is used by the local hybrid filter LKNETF. It is a variant of `U_likelihood_l` accounting for hybridization.
+
+
+The interface for this routine is:
+{{{
+SUBROUTINE likelihood_hyb_l(domain_p, step, dim_obs_l, obs_l, resid_l, gamma, likely_l)
+
+  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
+  INTEGER, INTENT(in) :: step                 ! Current time step
+  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
+  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
+  REAL, INTENT(inout) :: resid_l(dim_obs_l)   ! Input vector holding the local residual y-Hx
+  REAL, INTENT(in)    :: gamma                ! Hybrid weight
+  REAL, INTENT(out)   :: likely_l(dim_obs_l)  ! Output value of the likelihood
+}}}
+
+This routine is a variant for `U_likelihood_l`. See the description of this routine for its functionality.  The 
+The routine is called during the loop over the local analysis domains. The likelihood of the local observations has to be computed for each ensemble member. The likelihood is computed from the observation-state residual according to the assumed observation error distribution. Commonly, the observation errors are assumed to be Gaussian distributed. In this case, the likelihood is '''exp(-0.5*(y-Hx)^T^*R^-1^*(y-Hx))'''.
+
+This routine is also the place to perform observation localization. To initialize a vector of weights, the routine `PDAF_local_weight` can be called. The procedure is used in the example implementation and also demonstrated in the template routine.
+
+The routine also has to apply the hybrid weight `gamma`. This is a simple multiplication with `1-gamma` in the loop where `Rinvresid_l` is initialized.
+
+Hints:
+ * This routine is a simple extension of `U_likelihood_l. One can implement the hybrid variant by copying this routine and adapting it. `gamma` is computed inside PDAF and provided to the routine. 
+
+
 === `U_next_observation` (next_observation_pdaf.F90) ===
 
@@ -478 +538 @@
 == Execution order of user-supplied routines ==
 
-The user-supplied routines are executed in the order listed below. The order can be important as some routines can perform preparatory work for routines executed later on during the analysis. For example, `U_init_dim_obs_l` can prepare an index array that provides the information how to localize a 'full' vector of observations. Some hints one the efficient implementation strategy are given with the descriptions of the routine interfaces above.
+The executation order and how ofter the user routines are called depends on the chosen hybrid filter variant. The two-step variants HNK (subtype=0) and HKN (subtype=1) perform two local analysis loops (one for LNETF and one for LETKF), while the synchronous variance (Hsync, subtype=4) perform only a single loop and computes the LETKF and NETF updates synchronously (which still requires multiple calls to user routines). The order can be important as some routines can perform preparatory work for routines executed later on during the analysis. For example, `U_init_dim_obs_l` is often used to prepare an index array that provides the information how to localize a 'full' vector of observations. Some hints one the efficient implementation strategy are given with the descriptions of the routine interfaces above.
 
 Before the analysis step is called the following is executed:
  1. [#U_collect_statecollect_state_pdaf.F90 U_collect_state] (called once for each ensemble member)
 
-When the ensemble integration of the forecast is completed, the analysis step is executed. Before the loop over all local analysis domains, the following routines are executed:
+When the ensemble integration of the forecast is completed, the analysis step is executed. Before the loop over all local analysis domains, the following routines are executed in the same way for all three hybrid filter variants:
  1. [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep] (Call to act on the forecast ensemble, called with negative value of the time step)
  1. [#U_init_n_domainsinit_n_domains_pdaf.F90 U_init_n_domains]
@@ -491 +551 @@
  1. [#U_init_obsvarinit_obsvar_pdaf.F90 U_init_obsvar] (Only executed, if the global adaptive forgetting factor is used (`type_forget=1` in the example implementation))
 
-In the loop over all local analysis domains, it is executed for each local analysis domain:
+For `Hsync`: In the loop over all local analysis domains, it is executed for each local analysis domain:
  1. [#U_init_dim_linit_dim_l_pdaf.F90 U_init_dim_l]
  1. [#U_init_dim_obs_linit_dim_obs_l_pdaf.F90 U_init_dim_obs_l]
  1. [#U_g2l_stateg2l_state_pdaf.F90 U_g2l_state] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)
- 1. [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs] (A single call to localize the mean observed state)
+ 1. [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)
  1. [#U_init_obs_linit_obs_l_pdaf.F90 U_init_obs_l]
- 1. [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs] (`dim_ens` calls: one call to localize the observed part of each ensemble member) 
  1. [#U_init_obsvar_linit_obsvar_l_pdaf.F90 U_init_obsvar_l] (Only called, if the local adaptive forgetting factor is used (`type_forget=2` in the example implementation))
  1. [#U_prodRinvA_lprodrinva_l_pdaf.F90 U_prodRinvA_l]
+ 1. [#U_likelihood_llikelihood_l_pdaf.F90 U_likelihood_l] (Calls `dim_ens` times, once for each ensemble state)
  1. [#U_l2g_statel2g_state_pdaf.F90 U_l2g_state] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)
 
-After the loop over all local analysis domains, it is executed:
+For `HNK` and `HKN` two local analysis loops are performed with additional initialization of observation information in between:  
+ 1. First local analysis loop:
+  1. [#U_init_dim_linit_dim_l_pdaf.F90 U_init_dim_l]
+  1. [#U_init_dim_obs_linit_dim_obs_l_pdaf.F90 U_init_dim_obs_l]
+  1. [#U_g2l_stateg2l_state_pdaf.F90 U_g2l_state] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)
+  1. [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)
+  1. [#U_init_obs_linit_obs_l_pdaf.F90 U_init_obs_l]
+  1. [#U_likelihood_llikelihood_l_pdaf.F90 U_likelihood_l] (Called `dim_ens` times to determine hybrid weight `gamma`)
+  1. Execute LNETF (for `HNK`) or LETKF (for `HKN`)
+  1. [#U_l2g_statel2g_state_pdaf.F90 U_l2g_state] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)
+ 1. In between both local analysis loops:
+  1. [#U_obs_op_fobs_op_f_pdaf.F90 U_obs_op_f] (Called `dim_ens` times; once for each ensemble member)
+ 1. Second local analysis loop:
+  1. [#U_init_dim_linit_dim_l_pdaf.F90 U_init_dim_l]
+  1. [#U_init_dim_obs_linit_dim_obs_l_pdaf.F90 U_init_dim_obs_l]
+  1. [#U_g2l_stateg2l_state_pdaf.F90 U_g2l_state] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)
+  1. [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)
+  1. [#U_init_obs_linit_obs_l_pdaf.F90 U_init_obs_l]
+  1. [#U_likelihood_llikelihood_l_pdaf.F90 U_likelihood_l] (Called `dim_ens` times to determine hybrid weight `gamma`)
+  1. Execute LETKF (for `HNK`) or LNETF (for `HKN`)
+  1. [#U_l2g_statel2g_state_pdaf.F90 U_l2g_state] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)
+
+In step 7: When the LNETF is executed the routine
+ 1. [#U_likelihood_hyb_llikelihood_hyb_l_pdaf.F90 U_likelihood_hyb_l] (Called `dim_ens` times, once for each ensemble member)
+is called. In contrast when the LETKF is computed the routines
+ 1. [#U_init_obsvar_linit_obsvar_l_pdaf.F90 U_init_obsvar_l] (Only called, if the local adaptive forgetting factor is used (`type_forget=2` in the example implementation))
+ 1. [#U_prodRinvA_hyb_lprodrinva_hyb_l_pdaf.F90 U_prodRinvA_hyb_l]
+are called.
+
+After the loop(s) over all local analysis domains, it is executed:
  1. [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep] (Call to act on the analysis ensemble, called with (positive) value of the time step)