| 1 | = likelihood_pdaf = |
| 2 | |
| 3 | The page document the user-supplied call-back routine `likelihood_pdaf`. |
| 4 | |
| 5 | The routine `likelihood_pdaf` (called `U_likelihood` inside the PDAF core routines) is a call-back routine that has to be provided by the user. |
| 6 | The routine is used with the nonlinear filter NETF and called during the analysis step. The purpose of the routine is to compute the likelihood of the observation for a given ensemble member. The likelihood depends on the assumed observation error distribution. For a Gaussian observation error, the likelihood is '''exp(-0.5*(y-Hx)^T^*R^-1^*(y-Hx))'''. The vector '''y-HX =''' `resid` is proved as an input argument. The likelihood has to be returned in the variable `likely`. |
| 7 | |
| 8 | The interface is the following: |
| 9 | {{{ |
| 10 | SUBROUTINE likelihood_pdaf(step, dim_obs_p, observation_p, resid, likely) |
| 11 | }}} |
| 12 | with |
| 13 | * `step` : `integer, intent(in)`[[BR]] Current time step |
| 14 | * `dim_obs_p` : `integer, intent(in)`[[BR]] Number of observations at current time step (i.e. the size of the observation vector) |
| 15 | * `obs_p` : `real, intent(in), dimension(dim_obs_p)`[[BR]] Vector of observations |
| 16 | * `resid` : `real, intent(in), dimension(dim_obs_p)`[[BR]] Input vector holding the residual |
| 17 | * `likely` : `real, intent(out)`[[BR]] Output value of the likelihood |
| 18 | |
| 19 | |
| 20 | Notes: |
| 21 | * In case of a parallelization 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 vector. |
| 22 | |
| 23 | Hints: |
| 24 | * The routine is very similar to the routine [wiki:prodRinvA_pdaf]. The main addition is the computation of the likelihood after computing '''R^-1^*(y-Hx)''', which corresponds to '''R^-1^*A_p''' in [wiki:prodRinvA_pdaf]. |
| 25 | * 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. |
| 26 | * 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. |
| 27 | * 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. |
| 28 | |
| 29 | |