Changes between Version 2 and Version 3 of ImplementAnalysislknetf


Ignore:
Timestamp:
Feb 19, 2023, 10:28:49 AM (13 months ago)
Author:
lnerger
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • ImplementAnalysislknetf

    v2 v3  
    272272=== `U_prodRinvA_hyb_l` (prodrinva_hyb_l_pdaf.F90) ===
    273273
    274 This routine is used by the local hybrid filter LKNETF.
     274This routine is used by the local hybrid filter LKNETF. It is a variant of `U_proRinvA_l` accounting for hybridization.
    275275
    276276The interface for this routine is:
     
    292292This 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.
    293293
    294 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.
     294The 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.
    295295
    296296Hints:
    297297 * 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.
     298
    298299
    299300=== `U_init_n_domains` (init_n_domains_pdaf.F90) ===
     
    471472
    472473
     474
     475=== `U_likelihood_l` (likelihood_l_pdaf.F90) ===
     476
     477This routine is used by the LNETF and LKNETF filters.
     478
     479The interface for this routine is:
     480{{{
     481SUBROUTINE likelihood_l(domain_p, step, dim_obs_l, obs_l, resid_l, likely_l)
     482
     483  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
     484  INTEGER, INTENT(in) :: step                 ! Current time step
     485  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
     486  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
     487  REAL, INTENT(inout) :: resid_l(dim_obs_l)   ! Input vector holding the local residual y-Hx
     488  REAL, INTENT(out)   :: likely_l(dim_obs_l)  ! Output value of the likelihood
     489}}}
     490
     491The 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))'''.
     492
     493This 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.
     494
     495Hints:
     496 * 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.
     497 * 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].
     498 * 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.
     499 * 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.
     500 * 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.
     501
     502
     503
     504=== `U_likelihood_hyb_l` (likelihood_hyb_l_pdaf.F90) ===
     505
     506This routine is used by the local hybrid filter LKNETF. It is a variant of `U_likelihood_l` accounting for hybridization.
     507
     508
     509The interface for this routine is:
     510{{{
     511SUBROUTINE likelihood_hyb_l(domain_p, step, dim_obs_l, obs_l, resid_l, gamma, likely_l)
     512
     513  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
     514  INTEGER, INTENT(in) :: step                 ! Current time step
     515  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
     516  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
     517  REAL, INTENT(inout) :: resid_l(dim_obs_l)   ! Input vector holding the local residual y-Hx
     518  REAL, INTENT(in)    :: gamma                ! Hybrid weight
     519  REAL, INTENT(out)   :: likely_l(dim_obs_l)  ! Output value of the likelihood
     520}}}
     521
     522This routine is a variant for `U_likelihood_l`. See the description of this routine for its functionality.  The
     523The 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))'''.
     524
     525This 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.
     526
     527The 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.
     528
     529Hints:
     530 * 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.
     531
     532
    473533=== `U_next_observation` (next_observation_pdaf.F90) ===
    474534
     
    478538== Execution order of user-supplied routines ==
    479539
    480 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.
     540The 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.
    481541
    482542Before the analysis step is called the following is executed:
    483543 1. [#U_collect_statecollect_state_pdaf.F90 U_collect_state] (called once for each ensemble member)
    484544
    485 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:
     545When 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:
    486546 1. [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep] (Call to act on the forecast ensemble, called with negative value of the time step)
    487547 1. [#U_init_n_domainsinit_n_domains_pdaf.F90 U_init_n_domains]
     
    491551 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))
    492552
    493 In the loop over all local analysis domains, it is executed for each local analysis domain:
     553For `Hsync`: In the loop over all local analysis domains, it is executed for each local analysis domain:
    494554 1. [#U_init_dim_linit_dim_l_pdaf.F90 U_init_dim_l]
    495555 1. [#U_init_dim_obs_linit_dim_obs_l_pdaf.F90 U_init_dim_obs_l]
    496556 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)
    497  1. [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs] (A single call to localize the mean observed state)
     557 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)
    498558 1. [#U_init_obs_linit_obs_l_pdaf.F90 U_init_obs_l]
    499  1. [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs] (`dim_ens` calls: one call to localize the observed part of each ensemble member)
    500559 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))
    501560 1. [#U_prodRinvA_lprodrinva_l_pdaf.F90 U_prodRinvA_l]
     561 1. [#U_likelihood_llikelihood_l_pdaf.F90 U_likelihood_l] (Calls `dim_ens` times, once for each ensemble state)
    502562 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)
    503563
    504 After the loop over all local analysis domains, it is executed:
     564For `HNK` and `HKN` two local analysis loops are performed with additional initialization of observation information in between: 
     565 1. First local analysis loop:
     566  1. [#U_init_dim_linit_dim_l_pdaf.F90 U_init_dim_l]
     567  1. [#U_init_dim_obs_linit_dim_obs_l_pdaf.F90 U_init_dim_obs_l]
     568  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)
     569  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)
     570  1. [#U_init_obs_linit_obs_l_pdaf.F90 U_init_obs_l]
     571  1. [#U_likelihood_llikelihood_l_pdaf.F90 U_likelihood_l] (Called `dim_ens` times to determine hybrid weight `gamma`)
     572  1. Execute LNETF (for `HNK`) or LETKF (for `HKN`)
     573  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)
     574 1. In between both local analysis loops:
     575  1. [#U_obs_op_fobs_op_f_pdaf.F90 U_obs_op_f] (Called `dim_ens` times; once for each ensemble member)
     576 1. Second local analysis loop:
     577  1. [#U_init_dim_linit_dim_l_pdaf.F90 U_init_dim_l]
     578  1. [#U_init_dim_obs_linit_dim_obs_l_pdaf.F90 U_init_dim_obs_l]
     579  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)
     580  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)
     581  1. [#U_init_obs_linit_obs_l_pdaf.F90 U_init_obs_l]
     582  1. [#U_likelihood_llikelihood_l_pdaf.F90 U_likelihood_l] (Called `dim_ens` times to determine hybrid weight `gamma`)
     583  1. Execute LETKF (for `HNK`) or LNETF (for `HKN`)
     584  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)
     585
     586In step 7: When the LNETF is executed the routine
     587 1. [#U_likelihood_hyb_llikelihood_hyb_l_pdaf.F90 U_likelihood_hyb_l] (Called `dim_ens` times, once for each ensemble member)
     588is called. In contrast when the LETKF is computed the routines
     589 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))
     590 1. [#U_prodRinvA_hyb_lprodrinva_hyb_l_pdaf.F90 U_prodRinvA_hyb_l]
     591are called.
     592
     593After the loop(s) over all local analysis domains, it is executed:
    505594 1. [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep] (Call to act on the analysis ensemble, called with (positive) value of the time step)
    506595