= Using non-diagonal R matrices with OMI = [[PageOutline(2-3,Contents of this page)]] This feature was introduced with PDAF V2.3. The default mode of PDAF-OMI is to use a diagonal observation error covariance matrix **R** and specifying the observation error variances, i.e. the diagonalof **R** only. This is in line with the common choice in data assimilation to assume that observation errors are uncorrelated. However, there are also observation types with significant observation error correlations, which should be represented by a non-diagonal observation error covariance matrix. With PDAF V2.3 support for such nondiagonal **R** matrices was added to OMI. PDAF-OMI's support for nondiagonal **R**-matrices consists in given the user access to the routines that perform operations involve **R**. This differs with the filter type, e.g. in LESTKF and LETKF a produce of some matrix with the inverse of **R** ha sto be computed, while in the traditional, perturbed observations, EnKF the matrix **R** has to be added to some other matrix. For the particle filter and the NETF the computation of the likelihood involves **R**. == Routines to perform the analysis step == For 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. ||= '''Global Filter'' =||= **diagonal R** =||= **non-diagonal R** =||= **additional routine(s)** =|| ||= 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 || ||= 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 || ||= 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 || ||= 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 || ||= 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 || ||= 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 || [[BR]] ||= '''Local Filter'' =||= **diagonal R** =||= **non-diagonal R** =||= **additional routine(s)** =|| ||= 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 || ||= 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 || ||= 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 || ||= 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 || ||= 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 || ||= 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 || [[BR]] ||= '''3D-Var'' =||= **diagonal R** =||= **non-diagonal R** =||= **additional routine(s)** =|| ||= 3DVar =|| [wiki:PDAFomi_assimilate_3dvar][[BR]][wiki:PDAFomi_assimilate_3dvar] || [wiki:PDAFomi_put_state_3dvar_nondiagR][[BR]][wiki:PDAFomi_put_state_3dvar_nondiagR ] || prodRinvA_pdafomi || ||= En3DVar ESTKF =|| [wiki:PDAFomi_assimilate_en3dvar_estkf][[BR]][wiki:PDAFomi_assimilate_en3dvar_estkf] || [wiki:PDAFomi_put_state_en3dvar_estkf_nondiagR][[BR]][wiki:PDAFomi_put_state_en3dvar_estkf_nondiagR] || prodRinvA_pdafomi || ||= En3DVar LESTKF =|| [wiki:PDAFomi_assimilate_en3dvar_lestkf][[BR]][wiki:PDAFomi_assimilate_en3dvar_lestkf] || [wiki:PDAFomi_put_state_en3dvar_lestkf_nondiagR][[BR]][wiki:PDAFomi_put_state_en3dvar_lestkf_nondiagR] || prodRinvA_l_pdafomi || ||= hyb3DVar ESTKF =|| [wiki:PDAFomi_assimilate_hyb3dvar_estkf][[BR]][wiki:PDAFomi_assimilate_hyb3dvar_estkf] || [wiki:PDAFomi_put_state_hyb3dvar_estkf_nondiagR][[BR]][wiki:PDAFomi_put_state_hyb3dvar_estkf_nondiagR] || prodRinvA_pdafomi || ||= hyb3DVar LESTKF =|| [wiki:PDAFomi_assimilate_hyb3dvar_lestkf][[BR]][wiki:PDAFomi_assimilate_hyb3dvar_lestkf] || [wiki:PDAFomi_put_state_hyb3dvar_lestkf_nondiagR][[BR]][wiki:PDAFomi_put_state_hyb3dvar_lestkf_nondiagR] || prodRinvA_l_pdafomi || == Call-back routines handling observations == The 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. Note 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. === prodRinvA_pdafomi === The interface for this routine is: {{{ SUBROUTINE prodRinvA_pdafomi(step, dim_obs_p, rank, obs_p, A_p, C_p) INTEGER, INTENT(in) :: step ! Current time step INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of obs. vector INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix REAL, INTENT(in) :: obs_p(dim_obs_p) ! PE-local vector of observations REAL, INTENT(in) :: A_p(dim_obs_p,rank) ! Input matrix from analysis routine REAL, INTENT(out) :: C_p(dim_obs_p,rank) ! Output matrix }}} The 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`. For 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. Hints: * 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`. * 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. * 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. * 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.) * 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. == likelihood_pdafomi == The interface for this routine is: {{{ SUBROUTINE likelihood_pdafomi(step, dim_obs_p, obs_p, residual, likely) INTEGER, INTENT(in) :: step ! Current time step INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of obs. vector REAL, INTENT(in) :: obs_p(dim_obs_p) ! PE-local vector of observations REAL, INTENT(in) :: residual(dim_obs_p) ! Input vector holding the residual y-Hx REAL, INTENT(out) :: likely ! Output value of the likelihood }}} The 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))'''. For 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. Hints: * 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`. * 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]. * 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_p` is provided through the interface for cases where the observation error variance is relative to the actual value of the observations. * 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. === add_obs_err_pdafomi === The interface for this routine is: {{{ SUBROUTINE add_obs_err_pdafomi(step, dim_obs, C) INTEGER, INTENT(in) :: step ! Current time step INTEGER, INTENT(in) :: dim_obs ! Dimension of obs. vector REAL, INTENT(inout) :: C(dim_obs, dim_obs) ! Matrix to that the observation ! error covariance matrix is added }}} During 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`. The operation is for the global observation space. Thus, it is independent of whether the filter is executed with or without parallelization. Hints: * 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. * 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. * 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.