| 474 | |
| 475 | === `U_likelihood_l` (likelihood_l_pdaf.F90) === |
| 476 | |
| 477 | This routine is used by the LNETF and LKNETF filters. |
| 478 | |
| 479 | The interface for this routine is: |
| 480 | {{{ |
| 481 | SUBROUTINE 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 | |
| 491 | 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))'''. |
| 492 | |
| 493 | 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. |
| 494 | |
| 495 | Hints: |
| 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 | |
| 506 | This routine is used by the local hybrid filter LKNETF. It is a variant of `U_likelihood_l` accounting for hybridization. |
| 507 | |
| 508 | |
| 509 | The interface for this routine is: |
| 510 | {{{ |
| 511 | SUBROUTINE 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 | |
| 522 | This routine is a variant for `U_likelihood_l`. See the description of this routine for its functionality. The |
| 523 | 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))'''. |
| 524 | |
| 525 | 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. |
| 526 | |
| 527 | The 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 | |
| 529 | Hints: |
| 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 | |
504 | | After the loop over all local analysis domains, it is executed: |
| 564 | For `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 | |
| 586 | In 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) |
| 588 | is 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] |
| 591 | are called. |
| 592 | |
| 593 | After the loop(s) over all local analysis domains, it is executed: |