Changes between Initial Version and Version 1 of nondiagonal_observation_error_covariance_matrices_PDAF3


Ignore:
Timestamp:
May 25, 2025, 7:33:17 PM (8 weeks ago)
Author:
lnerger
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • nondiagonal_observation_error_covariance_matrices_PDAF3

    v1 v1  
     1= The PDAF3 interface for non-diagonal R matrices =
     2
     3{{{
     4#!html
     5<div class="wiki-toc">
     6<h4>PDAF-OMI Guide</h4>
     7<ol><li><a href="PDAF_OMI_Overview">Overview</a></li>
     8<li><a href="OMI_Callback_obs_pdafomi">callback_obs_pdafomi.F90</a></li>
     9<li><a href="OMI_observation_modules">Observation Modules</a></li>
     10<li><a href="OMI_observation_operators">Observation operators</a></li>
     11<li><a href="OMI_error_checking">Checking error status</a></li>
     12<li><a href="OMI_debugging">Debugging functionality</a></li>
     13<li><a href="OMI_ImplementationofAnalysisStep">Implementing the analysis step with OMI</a></li>
     14<ol>
     15<li> <a href="ImplementFilterAnalysisOverview"> General overview for ensemble filters</a></li>
     16<ol>
     17<li>Implementation for Global Filters</li>
     18<li><a href="ImplementAnalysisLocal">Implementation for Local Filters</a></li>
     19<li><a href="ImplementAnalysislenkfOmi">Implementation for LEnKF</a></li>
     20</ol>
     21<li> <a href="Implement3DVarAnalysisOverview"> General overview for 3D-Var methods</a></li>
     22<ol>
     23<li><a href="ImplementAnalysis_3DVar">Implementation for 3D-Var</a></li>
     24<li><a href="ImplementAnalysis_3DEnVar">Implementation for 3D Ensemble Var</a></li>
     25<li><a href="ImplementAnalysis_Hyb3DVar">Implementation for Hybrid 3D-Var</a></li>
     26</ol>
     27</ol>
     28<li>Using nondiagonal R-matrices</li>
     29<li><a href="Porting_to_OMI">Porting an existing implemention to OMI</a></li>
     30<li><a href="PDAFomi_additional_functionality">Additional OMI Functionality</a></li>
     31</ol>
     32</div>
     33}}}
     34
     35[[PageOutline(2-3,Contents of this page)]]
     36
     37The default mode of the PDAF3 assimilation routines are designed to use a diagonal observation error covariance matrix **R** and specifying the observation error variances, i.e. the diagonal of **R** only. This is in line with the common choice in data assimilation to assume that observation errors are uncorrelated.
     38
     39However, there are also observation types with significant observation error correlations, which should be represented by a non-diagonal observation error covariance matrix. The PDAF3 assimilation interface provides support for such nondiagonal **R** matrices.
     40
     41To support nondiagonal **R**-matrices consists in giving the user access to the routines that perform operations involving **R**. The operations dedepend on the filter type, e.g. in LESTKF and LETKF a product of some matrix with the inverse of **R** has to be computed, while in the traditional, perturbed observations, EnKF the matrix **R** has to be added to some other matrix. For the particle filter, the NETF and the hybrid filter LKNETF, the computation of the likelihood involves **R**.
     42
     43
     44
     45== Routines to perform the analysis step ==
     46
     47For using a nondiagonal **R** there is a modified variant of the PDAF-OMI routines `PDAFomi_assimilate_*` or `PDAFomi_put_state_*`. Compared to the case of a diagonal **R**, the routines are slightly less generic. Below we provide an overview of the routines including the links to the interface description of each. The last column gives the name of the additional routine(s) provided in the interface to implement the non-diagonal **R**-matrix.
     48
     49||= '''Global Filter'' =||= **diagonal R** =||= **non-diagonal R** =||= **additional routine(s)** =||
     50||= EnKF =|| [wiki:PDAFomi_assimilate_global][[BR]][wiki:PDAFomi_put_state_global] ||  [wiki:PDAFomi_assimilate_enkf_nondiagR][[BR]][wiki:PDAFomi_put_state_enkf_nondiagR]  || add_obs_err_pdafomi [[BR]]init_obscovar_pdafomi ||
     51||= ESTKF =|| [wiki:PDAFomi_assimilate_global][[BR]][wiki:PDAFomi_put_state_global] ||  [wiki:PDAFomi_assimilate_global_nondiagR][[BR]][wiki:PDAFomi_put_state_global_nondiagR]  || prodRinvA_pdafomi ||
     52||= ETKF =|| [wiki:PDAFomi_assimilate_global][[BR]][wiki:PDAFomi_put_state_global] ||  [wiki:PDAFomi_assimilate_global_nondiagR][[BR]][wiki:PDAFomi_put_state_global_nondiagR]  || prodRinvA_pdafomi ||
     53||= PF =|| [wiki:PDAFomi_assimilate_global][[BR]][wiki:PDAFomi_put_state_global] ||  [wiki:PDAFomi_assimilate_nonlin_nondiagR][[BR]][wiki:PDAFomi_put_state_nonlin_nondiagR]  || likelihood_pdafomi ||
     54||= NETF =|| [wiki:PDAFomi_assimilate_global][[BR]][wiki:PDAFomi_put_state_global] ||  [wiki:PDAFomi_assimilate_nonlin_nondiagR][[BR]][wiki:PDAFomi_put_state_nonlin_nondiagR]  || likelihood_pdafomi ||
     55||= SEIK =|| [wiki:PDAFomi_assimilate_global][[BR]][wiki:PDAFomi_put_state_global] ||  [wiki:PDAFomi_assimilate_global_nondiagR][[BR]][wiki:PDAFomi_put_state_global_nondiagR]  || prodRinvA_pdafomi ||
     56[[BR]]
     57||= '''Local Filter'' =||= **diagonal R** =||= **non-diagonal R** =||= **additional routine(s)** =||
     58||= LEnKF =|| [wiki:PDAFomi_assimilate_lenkf][[BR]][wiki:PDAFomi_put_state_lenkf] ||  [wiki:PDAFomi_assimilate_lenkf_nondiagR][[BR]][wiki:PDAFomi_put_state_lenkf_nondiagR]  || add_obs_err_pdafomi [[BR]]init_obscovar_pdafomi ||
     59||= LESTKF =|| [wiki:PDAFomi_assimilate_local][[BR]][wiki:PDAFomi_put_state_local] ||  [wiki:PDAFomi_assimilate_local_nondiagR][[BR]][wiki:PDAFomi_put_state_local_nondiagR]  || prodRinvA_l_pdafomi ||
     60||= LETKF =|| [wiki:PDAFomi_assimilate_local][[BR]][wiki:PDAFomi_put_state_local] ||  [wiki:PDAFomi_assimilate_local_nondiagR][[BR]][wiki:PDAFomi_put_state_local_nondiagR]  || prodRinvA_l_pdafomi ||
     61||= LSEIK =|| [wiki:PDAFomi_assimilate_local][[BR]][wiki:PDAFomi_put_state_local] ||  [wiki:PDAFomi_assimilate_local_nondiagR][[BR]][wiki:PDAFomi_put_state_local_nondiagR]  || prodRinvA_l_pdafomi ||
     62||= LNETF =|| [wiki:PDAFomi_assimilate_local][[BR]][wiki:PDAFomi_put_state_local] ||  [wiki:PDAFomi_assimilate_lnetf_nondiagR][[BR]][wiki:PDAFomi_put_state_lnetf_nondiagR]  || likelihood_l_pdafomi ||
     63||= LKNETF =|| [wiki:PDAFomi_assimilate_local][[BR]][wiki:PDAFomi_put_state_local] ||  [wiki:PDAFomi_assimilate_lknetf_nondiagR][[BR]][wiki:PDAFomi_put_state_lknetf_nondiagR]  ||likelihood_l_pdafomi[[BR]]likelihood_hyb_l_pdafomi[[BR]] prodRinvA_l_pdafomi[[BR]] prodRinvA_hyb_l_pdafomi  ||
     64[[BR]]
     65||= '''3D-Var'' =||= **diagonal R** =||= **non-diagonal R** =||= **additional routine(s)** =||
     66||= 3DVar =|| [wiki:PDAFomi_assimilate_3dvar][[BR]][wiki:PDAFomi_put_state_3dvar] || [wiki:PDAFomi_assimilate_3dvar_nondiagR][[BR]][wiki:PDAFomi_put_state_3dvar_nondiagR ] || prodRinvA_pdafomi ||
     67||= En3DVar ESTKF =|| [wiki:PDAFomi_assimilate_en3dvar_estkf][[BR]][wiki:PDAFomi_put_state_en3dvar_estkf] || [wiki:PDAFomi_assimilate_en3dvar_estkf_nondiagR][[BR]][wiki:PDAFomi_put_state_en3dvar_estkf_nondiagR] || prodRinvA_pdafomi ||
     68||= En3DVar LESTKF =|| [wiki:PDAFomi_assimilate_en3dvar_lestkf][[BR]][wiki:PDAFomi_put_state_en3dvar_lestkf] || [wiki:PDAFomi_assimilate_en3dvar_lestkf_nondiagR][[BR]][wiki:PDAFomi_put_state_en3dvar_lestkf_nondiagR] ||  prodRinvA_l_pdafomi ||
     69||= hyb3DVar ESTKF =|| [wiki:PDAFomi_assimilate_hyb3dvar_estkf][[BR]][wiki:PDAFomi_put_state_hyb3dvar_estkf] || [wiki:PDAFomi_assimilate_hyb3dvar_estkf_nondiagR][[BR]][wiki:PDAFomi_put_state_hyb3dvar_estkf_nondiagR] ||  prodRinvA_pdafomi ||
     70||= hyb3DVar LESTKF =|| [wiki:PDAFomi_assimilate_hyb3dvar_lestkf][[BR]][wiki:PDAFomi_put_state_hyb3dvar_lestkf] || [wiki:PDAFomi_assimilate_hyb3dvar_lestkf_nondiagR][[BR]][wiki:PDAFomi_put_state_hyb3dvar_lestkf_nondiagR] || prodRinvA_l_pdafomi ||
     71
     72== Call-back routines handling observations ==
     73
     74The tables above show the additional routine(s) that need to be implemented the nondiagonal **R** matrix for  the different filters and 3D-Var variants. Here we describe the interface of these routines. Generally the required routine(s) for the chosen DA method should be added to `callback_obs_pdafomi.F90` and to each observation module, where `thisobs` or `thisobs_l` can be included.
     75
     76Note that a mixed implementation is not possible. If one observation type has a non-diagonal **R** matrix, the additional routine(s) have to be implemented for all observations, even if an observation type uses a diagonal **R**. However, for observation wit a diagonal **R**, the corresponding PDAF-OMI routine to be called in the observation module can be called. Below we provide the information where this routine can be found.
     77
     78=== prodRinvA_pdafomi ===
     79
     80The interface for this routine is:
     81{{{
     82SUBROUTINE prodRinvA_pdafomi(step, dim_obs_p, rank, obs_p, A_p, C_p)
     83
     84  INTEGER, INTENT(in) :: step                ! Current time step
     85  INTEGER, INTENT(in) :: dim_obs_p           ! PE-local dimension of obs. vector
     86  INTEGER, INTENT(in) :: rank                ! Rank of initial covariance matrix
     87  REAL, INTENT(in)    :: obs_p(dim_obs_p)    ! PE-local vector of observations
     88  REAL, INTENT(in)    :: A_p(dim_obs_p,rank) ! Input matrix from analysis routine
     89  REAL, INTENT(out)   :: C_p(dim_obs_p,rank) ! Output matrix
     90}}}
     91
     92The routine computes the product of the inverse of the observation error covariance matrix with matrix `A_p`. For the ESTKF this matrix holds the observed part of the ensemble perturbations. The product has to be given as `C_p`.
     93
     94For a model with domain decomposition, `A_p` contains the part of the matrix that resides on the model sub-domain of the calling process. The product has to be computed for this sub-domain, too.
     95
     96Hints:
     97 * The matrix `A_p` relates to the observations of all observation types. Thus, one needs to take the offset of an observation in the observation vector of all observation types into account. This offset is given by `thisobs%off_obs_f`.
     98 * The routine does not require that the product is implemented as a real matrix-matrix 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 diagonal with matrix `A_p` has to be implemented.
     99 * The observation vector `obs_p` is provided through the interface for cases where the observation error variance is relative to the actual value of the observations.
     100 * The interface has a difference for ESTKF and ETKF: For ETKF the third argument is the ensemble size (`dim_ens`), while for the ESTKF it is the rank (`rank`) of the covariance matrix (usually ensemble size minus one). In addition, the second dimension of `A_p` and `C_p` has size `dim_ens` for ETKF, while it is `rank` for the ESTKF.  (Practically, one can usually ignore this difference as the fourth argument of the interface can be named arbitrarily in the routine.)
     101 * For diagonal **R** PDAF-OMI uses the routine `PDAFomi_prodRinvA` in `/src/PDAFomi_obs_f.F90` as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.
     102
     103=== likelihood_pdafomi ===
     104
     105The interface for this routine is:
     106{{{
     107SUBROUTINE likelihood_pdafomi(step, dim_obs_p, obs_p, residual, likely)
     108
     109  INTEGER, INTENT(in) :: step                ! Current time step
     110  INTEGER, INTENT(in) :: dim_obs_p           ! PE-local dimension of obs. vector
     111  REAL, INTENT(in)    :: obs_p(dim_obs_p)    ! PE-local vector of observations
     112  REAL, INTENT(in)    :: residual(dim_obs_p) ! Input vector holding the residual y-Hx
     113  REAL, INTENT(out)   :: likely              ! Output value of the likelihood
     114}}}
     115
     116The routine computes the likelihood of the observations 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))'''.
     117
     118For a model with domain decomposition, `resid` contains the part of the matrix that resides on the model sub-domain of the calling process. The likelihood has to be computed for the global state vector. Thus some parallel communication might be required to complete the computation.
     119
     120Hints:
     121 * The matrix `residual` relates to the observations of all observation types. Thus, one needs to take the offset of an observation in the observation vector of all observation types into account. This offset is given by `thisobs%off_obs_f`.
     122 * The routine is very similar to the routine [wiki:U_prodRinvA]. 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].
     123 * 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.
     124 * 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.
     125 * The observation vector `obs_p` is provided through the interface for cases where the observation error variance is relative to the actual value of the observations.
     126 * For diagonal **R** PDAF-OMI uses the routine `PDAFomi_likelihood` in `/src/PDAFomi_obs_f.F90` as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.
     127
     128=== add_obs_err_pdafomi ===
     129
     130The interface for this routine is:
     131{{{
     132SUBROUTINE add_obs_err_pdafomi(step, dim_obs, C)
     133
     134  INTEGER, INTENT(in) :: step                ! Current time step
     135  INTEGER, INTENT(in) :: dim_obs             ! Dimension of obs. vector
     136  REAL, INTENT(inout) :: C(dim_obs, dim_obs) ! Matrix to that the observation
     137                                             !    error covariance matrix is added
     138}}}
     139
     140During the analysis step of the EnKF, the projection of the ensemble covariance onto the observation space is computed. This matrix is provided to the routine as `C_p`. The routine has to add the observation error covariance matrix to `C_p`.
     141
     142The operation is for the global observation space. Thus, it is independent of whether the filter is executed with or without parallelization.
     143
     144Hints:
     145 * Matrix C relates to the observations of all types. Thus, for a single obsevation type one has to take the offset of this observtion type in the observation vector into account. For this, PDAF-OMI initializes an array `obsdims`, which can also be used by a user-implemented routine. See the routine `PDAFomi_add_obs_error` in `/src/PDAFomi_obs_f.F90` for how this can be implemented.
     146 * The routine does not require that the observation error covariance matrix is added as a full matrix. If the matrix is diagonal, only the diagonal elements have to be added.
     147 * For diagonal **R**, PDAF-OMI uses the routine `PDAFomi_add_obs_error` in `/src/PDAFomi_obs_f.F90` as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.
     148
     149=== init_obscovar_pdafomi ===
     150
     151The interface for this routine is:
     152{{{
     153SUBROUTINE init_obscovar_pdafomi(step, dim_obs, dim_obs_p, covar, m_state_p, &
     154     isdiag)
     155
     156  INTEGER, INTENT(in) :: step                ! Current time step
     157  INTEGER, INTENT(in) :: dim_obs             ! Dimension of observation vector
     158  INTEGER, INTENT(in) :: dim_obs_p           ! PE-local dimension of observation vector
     159  REAL, INTENT(out) :: covar(dim_obs, dim_obs) ! Observation error covariance matrix
     160  REAL, INTENT(in)  :: m_state_p(dim_obs_p)  ! PE-local observation vector
     161  LOGICAL, INTENT(out) :: isdiag             ! Whether the observation error covar. matrix is diagonal
     162}}}
     163
     164The routine has to initialize the global observation error covariance matrix `covar`. In addition, the flag `isdiag` has to be initialized to provide the information, whether the observation error covariance matrix is diagonal.
     165
     166The operation is for the global observation space. Thus, it is independent of whether the filter is executed with or without parallelization.
     167
     168Hints:
     169 * Matrix `covar` relates to the observations of all types. Thus, for a single obsevation type one has to take the offset of this observtion type in the observation vector into account. For this, PDAF-OMI initializes an array `obsdims`, which can also be used by a user-implemented routine. See the routine `PDAFomi_init_obscovar` in `/src/PDAFomi_obs_f.F90` for how this can be implemented.
     170 * The local observation vector `m_state_p` is provided to the routine for the case that the observation errors are relative to the value of the observation.
     171 * For diagonal **R**, PDAF-OMI uses the routine `PDAFomi_init_obscovar` in `/src/PDAFomi_obs_f.F90` as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.
     172
     173=== prodRinvA_l_pdafomi ===
     174
     175The interface for this routine is:
     176{{{
     177SUBROUTINE prodRinvA_l_pdafomi(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l)
     178
     179  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
     180  INTEGER, INTENT(in) :: step                 ! Current time step
     181  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
     182  INTEGER, INTENT(in) :: rank                 ! Rank of initial covariance matrix
     183  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
     184  REAL, INTENT(inout) :: A_l(dim_obs_l, rank) ! Input matrix from analysis routine
     185  REAL, INTENT(out)   :: C_l(dim_obs_l, rank) ! Output matrix
     186}}}
     187
     188The routine is called during the loop over the local analysis domains. In the algorithm, the product of the inverse of the observation error covariance matrix with some matrix has to be computed. For the SEIK filter this matrix holds the observed part of the ensemble perturbations for the local analysis domain of index `domain_p`. The matrix is provided as `A_l`. The product has to be given as `C_l`.
     189
     190This 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.
     191
     192Hints:
     193 * The routine is a local variant of the routine `U_prodRinvA`. Thus if that routine has been implemented before, it can be adapted here for the local filter.
     194 * The matrix `A_p` relates to the observations of all observation types. Thus, one needs to take the offset of an observation in the observation vector of all observation types into account. This offset is given by `thisobs_l%off_obs_l`.
     195 * The routine does not require that the product is implemented as a real matrix-matrix 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 diagonal with matrix `A_l` has to be implemented.
     196 * 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.
     197 * The interface has a difference for LESTKF and LETKF: For LETKF the third argument is the ensemble size (`dim_ens`), while for LESTKF it is the rank (`rank`) of the covariance matrix (usually ensemble size minus one). In addition, the second dimension of `A_l` and `C_l` has size `dim_ens` for LETKF, while it is `rank` for LESTKF.  (Practically, one can usually ignore this difference as the fourth argument of the interface can be named arbitrarily in the routine.)
     198 * For diagonal **R** PDAF-OMI uses the routine `PDAFomi_prodRinvA_l` in `/src/PDAFomi_obs_k.F90` as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.
     199
     200=== likelihood_l_pdafomi ===
     201
     202The interface for this routine is:
     203{{{
     204SUBROUTINE likelihood_l_pdafomi(domain_p, step, dim_obs_l, obs_l, resid_l, likely_l)
     205
     206  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
     207  INTEGER, INTENT(in) :: step                 ! Current time step
     208  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
     209  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
     210  REAL, INTENT(inout) :: resid_l(dim_obs_l)   ! Input vector holding the local residual y-Hx
     211  REAL, INTENT(out)   :: likely_l(dim_obs_l)  ! Output value of the likelihood
     212}}}
     213
     214The routine is called during the loop over the local analysis domains. In the NETF, as in other particle filters, 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))'''.
     215
     216This 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.
     217
     218Hints:
     219* 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.
     220 * The matrix `resid_l` relates to the observations of all observation types. Thus, one needs to take the offset of an observation in the observation vector of all observation types into account. This offset is given by `thisobs_l%off_obs_l`.
     221 * 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].
     222 * 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.
     223 * 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.
     224 * 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.
     225 * For diagonal **R** PDAF-OMI uses the routine `PDAFomi_likelihood_l` in `/src/PDAFomi_obs_l.F90` as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.
     226
     227=== prodRinvA_hyb_l_pdafomi ===
     228
     229The interface for this routine is:
     230{{{
     231SUBROUTINE prodRinvA_hyb_l_pdafomi(domain_p, step, dim_obs_l, dim_ens, obs_l, gamma, A_l, C_l)
     232
     233  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
     234  INTEGER, INTENT(in) :: step                 ! Current time step
     235  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
     236  INTEGER, INTENT(in) :: dim_ens              ! Ensemble size
     237  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
     238  REAL, INTENT(in)    :: gamma                ! Hybrid weight
     239  REAL, INTENT(inout) :: A_l(dim_obs_l, dim_ens) ! Input matrix from analysis routine
     240  REAL, INTENT(out)   :: C_l(dim_obs_l, dim_ens) ! Output matrix
     241}}}
     242
     243The routine is called during the loop over the local analysis domains. In the algorithm, the product of the inverse of the observation error covariance matrix with some matrix has to be computed. This matrix holds the observed part of the ensemble perturbations for the local analysis domain of index `domain_p`. The matrix is provided as `A_l`. The product has to be given as `C_l`.
     244
     245This 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.
     246
     247The routine is analogous to `prodRinvA_l_pdafomi, but also has to apply the hybrid weight `gamma` of the hybrid filter LKNETF. This is a simple multiplication with the input value in the loop where `C_l` is initialized.
     248
     249Hints:
     250 * This routine is a simple extension of `prodRinvA_l_pdafomi`. One can implement the hybrid variant by copying this routine and adapting it. `gamma` is computed inside PDAF and provided to the routine.
     251 * Please see the further implementation hints for `prodRinvA_l_pdafomi`.
     252 * For diagonal **R** PDAF-OMI uses the routine `PDAFomi_prodRinvA_hyb_l` in `/src/PDAFomi_obs_l.F90` as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.
     253
     254=== likelihood_hyb_l_pdafomi ===
     255
     256The interface for this routine is:
     257{{{
     258SUBROUTINE likelihood_hyb_l_pdafomi(domain_p, step, dim_obs_l, obs_l, resid_l, gamma, likely_l)
     259
     260  INTEGER, INTENT(in) :: domain_p             ! Current local analysis domain
     261  INTEGER, INTENT(in) :: step                 ! Current time step
     262  INTEGER, INTENT(in) :: dim_obs_l            ! Dimension of local observation vector
     263  REAL, INTENT(in)    :: obs_l(dim_obs_l)     ! Local vector of observations
     264  REAL, INTENT(inout) :: resid_l(dim_obs_l)   ! Input vector holding the local residual y-Hx
     265  REAL, INTENT(in)    :: gamma                ! Hybrid weight
     266  REAL, INTENT(out)   :: likely_l(dim_obs_l)  ! Output value of the likelihood
     267}}}
     268
     269This routine is a variant for `U_likelihood_l_pdafomi`. See the description of this routine for its functionality. 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))'''.
     270
     271This 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.
     272
     273The 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.
     274
     275Hints:
     276 * 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.
     277 * Please see the further implementation hints for `prodRinvA_l_pdafomi`.
     278 * For diagonal **R** PDAF-OMI uses the routine `PDAFomi_likelihood_hyb_l` in `/src/PDAFomi_obs_l.F90` as the routine that is called within the observation module. This routine can serve as a template for an implementation by the user.
     279
     280== Hints ==
     281
     282For now, we do not have an example code that we could provide. However, for the implementation the follonig hints might be useful:
     283 * The matrix **R** could be constructed e.g. from decorrelation length scales or by direct implementation.
     284 * The routine [wiki:PDAF_correlation_function] can be used to obtain values of a correlation function.
     285 * Th eroutine [wiki:PDAFomi_observation_localization_weights] can be used to obtain the values of the localization weights.
     286 * The matrix **R** cannot be stored fully for larger numbers of observations. E.g. for 100 000 observations, the full matric **R** would have 10^10^ entries, which corresponds to 80 GB memory at double precision. Thus, one would need compressed ways of storing the matrix or one can implement the products with the inverse of **R** or the addition of **R** in form of oeprations.
     287 * For the localized operations in `prodRinvA_l_pdafomi` and `likelihood_l_pdafomi` it should be possible to handle the part of **R** that relates only to the local observations used in a local analysis domain. However, the inverse of this local matrix part is not identical to the local part of the global inverse. The different might be mitigated by applying the localization weight to **R**.
     288 * Examples of implementations are also available without PDAF-OMI, e.g. in the tutorial and in /models/lorenz96
     289.