| 117 | | * 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. |
| | 117 | * 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. |
| | 118 | |
| | 119 | === init_obscovar_pdafomi === |
| | 120 | |
| | 121 | The interface for this routine is: |
| | 122 | {{{ |
| | 123 | SUBROUTINE init_obscovar_pdafomi(step, dim_obs, dim_obs_p, covar, m_state_p, & |
| | 124 | isdiag) |
| | 125 | |
| | 126 | INTEGER, INTENT(in) :: step ! Current time step |
| | 127 | INTEGER, INTENT(in) :: dim_obs ! Dimension of observation vector |
| | 128 | INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of observation vector |
| | 129 | REAL, INTENT(out) :: covar(dim_obs, dim_obs) ! Observation error covariance matrix |
| | 130 | REAL, INTENT(in) :: m_state_p(dim_obs_p) ! PE-local observation vector |
| | 131 | LOGICAL, INTENT(out) :: isdiag ! Whether the observation error covar. matrix is diagonal |
| | 132 | }}} |
| | 133 | |
| | 134 | The 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. |
| | 135 | |
| | 136 | The operation is for the global observation space. Thus, it is independent of whether the filter is executed with or without parallelization. |
| | 137 | |
| | 138 | Hints: |
| | 139 | * 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. |
| | 140 | * 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. |
| | 141 | * 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. |
| | 142 | |
| | 143 | === prodRinvA_l_pdafomi === |
| | 144 | |
| | 145 | The interface for this routine is: |
| | 146 | {{{ |
| | 147 | SUBROUTINE prodRinvA_l_pdafomi(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) |
| | 148 | |
| | 149 | INTEGER, INTENT(in) :: domain_p ! Current local analysis domain |
| | 150 | INTEGER, INTENT(in) :: step ! Current time step |
| | 151 | INTEGER, INTENT(in) :: dim_obs_l ! Dimension of local observation vector |
| | 152 | INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix |
| | 153 | REAL, INTENT(in) :: obs_l(dim_obs_l) ! Local vector of observations |
| | 154 | REAL, INTENT(inout) :: A_l(dim_obs_l, rank) ! Input matrix from analysis routine |
| | 155 | REAL, INTENT(out) :: C_l(dim_obs_l, rank) ! Output matrix |
| | 156 | }}} |
| | 157 | |
| | 158 | The 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`. |
| | 159 | |
| | 160 | 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. |
| | 161 | |
| | 162 | Hints: |
| | 163 | * 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. |
| | 164 | * 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`. |
| | 165 | * 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. |
| | 166 | * 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. |
| | 167 | * 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.) |
| | 168 | * 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. |
| | 169 | |
| | 170 | === likelihood_l_pdafomi === |
| | 171 | |
| | 172 | The interface for this routine is: |
| | 173 | {{{ |
| | 174 | SUBROUTINE likelihood_l_pdafomi(domain_p, step, dim_obs_l, obs_l, resid_l, likely_l) |
| | 175 | |
| | 176 | INTEGER, INTENT(in) :: domain_p ! Current local analysis domain |
| | 177 | INTEGER, INTENT(in) :: step ! Current time step |
| | 178 | INTEGER, INTENT(in) :: dim_obs_l ! Dimension of local observation vector |
| | 179 | REAL, INTENT(in) :: obs_l(dim_obs_l) ! Local vector of observations |
| | 180 | REAL, INTENT(inout) :: resid_l(dim_obs_l) ! Input vector holding the local residual y-Hx |
| | 181 | REAL, INTENT(out) :: likely_l(dim_obs_l) ! Output value of the likelihood |
| | 182 | }}} |
| | 183 | |
| | 184 | The 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))'''. |
| | 185 | |
| | 186 | 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. |
| | 187 | |
| | 188 | Hints: |
| | 189 | * 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. |
| | 190 | * 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`. |
| | 191 | * 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]. |
| | 192 | * 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. |
| | 193 | * 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. |
| | 194 | * 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. |
| | 195 | * 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. |
