77 | | == `PDAF_put_state_lestkf` == |
78 | | |
79 | | When the 'flexible' implementation variant is chosen for the assimilation system, the routine `PDAF_put_state_lestkf` has to be used instead of `PDAF_assimilate_lestkf`. The general aspects of the filter specific routines `PDAF_put_state_*` have been described on the page [ModifyModelforEnsembleIntegration Modification of the model code for the ensemble integration]. The interface of the routine is identical with that of `PDAF_assimilate_lestkf` with the exception the specification of the user-supplied routines `U_distribute_state` and `U_next_observation` are missing. |
80 | | |
81 | | The interface when using the LESTKF algorithm is the following: |
82 | | {{{ |
83 | | SUBROUTINE PDAF_put_state_lestkf(U_collect_state, U_init_dim_obs_f, U_obs_op_f, U_init_obs_f, & |
84 | | U_init_obs_l, U_prepoststep, U_prodRinvA_l, U_init_n_domains, & |
85 | | U_init_dim_l, U_init_dim_obs_l, & |
86 | | U_g2l_state, U_l2g_state, U_g2l_obs, & |
87 | | U_init_obsvar, U_init_obsvar_l, status) |
| 72 | == `PDAF_put_state_local` == |
| 73 | |
| 74 | When the 'flexible' implementation variant is chosen for the assimilation system, the routine `PDAF_put_state_local` has to be used instead of `PDAF_assimilate_local`. The general aspects of the filter specific routines `PDAF_put_state_*` have been described on the page [ModifyModelforEnsembleIntegration Modification of the model code for the ensemble integration]. The interface of the routine is identical with that of `PDAF_assimilate_local` with the exception the specification of the user-supplied routines `U_distribute_state` and `U_next_observation` are missing. |
| 75 | |
| 76 | The interface when using one of the local filters is the following: |
| 77 | {{{ |
| 78 | SUBROUTINE PDAF_put_state_lestkf(U_collect_state, & |
| 79 | U_init_dim_obs_f, U_obs_op_f, & |
| 80 | U_prepoststep, U_init_n_domains, U_init_dim_l, & |
| 81 | U_init_dim_obs_l, U_g2l_state, U_l2g_state, & |
| 82 | status) |
146 | | |
147 | | === `U_init_obs_f` (init_obs_f_pdaf.F90) === |
148 | | |
149 | | This routine is used by all filter algorithms with domain-localization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm. |
150 | | The routine is only called if the globally adaptive forgetting factor is used (`type_forget=1` in the example implementation). For the local filters there is also the alternative to use locally adaptive forgetting factors (`type_forget=2` in the example implementation) |
151 | | |
152 | | The interface for this routine is: |
153 | | {{{ |
154 | | SUBROUTINE init_obs_f(step, dim_obs_f, observation_f) |
155 | | |
156 | | INTEGER, INTENT(in) :: step ! Current time step |
157 | | INTEGER, INTENT(in) :: dim_obs_f ! Dimension of full observation vector |
158 | | REAL, INTENT(out) :: observation_f(dim_obs_f) ! Full observation vector |
159 | | }}} |
160 | | |
161 | | The routine is called during the analysis step before the loop over the local analysis domains is entered. It has to provide the full vector of observations in `observation_f` for the current time step. The caller is the routine that computes an adaptive forgetting factor (PDAF_set_forget). |
162 | | |
163 | | Hints: |
164 | | * As for the other 'full' routines: While the global counterpart of this routine (`init_obs`) has to initialize the observation vector only for the local model sub-domain, the 'full' routine has to include observations that spatially belong to neighboring model sub-domains. As an easy choice one can simply initialize a vector of all globally available observations. |
165 | | * If the adaptive forgetting factor is not used, this routine only has to exist. However, no functionality is required. |
166 | | |
167 | | |
168 | | === `U_init_obs_l` (init_obs_l_pdaf.F90) === |
169 | | |
170 | | This routine is used by all filter algorithms with domain-localization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm. |
171 | | |
172 | | The interface for this routine is: |
173 | | {{{ |
174 | | SUBROUTINE init_obs_l(domain_p, step, dim_obs_l, observation_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 ! Local dimension of observation vector |
179 | | REAL, INTENT(out) :: observation_l(dim_obs_l) ! Local observation vector |
180 | | }}} |
181 | | |
182 | | The routine is called during the analysis step during the loop over the local analysis domain. |
183 | | It has to provide the vector of observations for the analysis in the local analysis domain with index `domain_p` in `observation_l` for the current time step. |
184 | | |
185 | | Hints: |
186 | | * For parallel efficiency, the LESTKF algorithm is implemented in a way that first the full vectors are initialized. These are then restricted to the local analysis domain during the loop over all local analysis domains. Thus, if the full vector of observations has been initialized before `U_init_obs_l` is executed (e.g. by `U_init_dim_obs_f`), the operations performed in this routine will be to select the part of the full observation vector that is relevant for the current local analysis domain. |
187 | | * The routine `U_init_dim_obs_l` is executed before this routine. Thus, if that routine already prepares the information which elements of `observation_f` are needed for `observation_l`, this information can be used efficiently here. |
226 | | |
227 | | === `U_prodRinvA_l` (prodrinva_l_pdaf.F90) === |
228 | | |
229 | | This routine is used by the local filters (LSEIK, LETKF, LESTKF). There is a slight difference between LESTKF and LETKF for this routine, which is described below. |
230 | | |
231 | | The interface for this routine is: |
232 | | {{{ |
233 | | SUBROUTINE prodRinvA_l(domain_p, step, dim_obs_l, rank, obs_l, A_l, C_l) |
234 | | |
235 | | INTEGER, INTENT(in) :: domain_p ! Current local analysis domain |
236 | | INTEGER, INTENT(in) :: step ! Current time step |
237 | | INTEGER, INTENT(in) :: dim_obs_l ! Dimension of local observation vector |
238 | | INTEGER, INTENT(in) :: rank ! Rank of initial covariance matrix |
239 | | REAL, INTENT(in) :: obs_l(dim_obs_l) ! Local vector of observations |
240 | | REAL, INTENT(inout) :: A_l(dim_obs_l, rank) ! Input matrix from analysis routine |
241 | | REAL, INTENT(out) :: C_l(dim_obs_l, rank) ! Output matrix |
242 | | }}} |
243 | | |
244 | | 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`. |
245 | | |
246 | | 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. |
247 | | |
248 | | Hints: |
249 | | * 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. |
250 | | * 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. |
251 | | * 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. |
252 | | * 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.) |
253 | | |
254 | | |
361 | | === `U_g2l_obs` (g2l_obs_pdaf.F90) === |
362 | | |
363 | | This routine is used by all filter algorithms with domain-localization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm. |
364 | | |
365 | | The interface for this routine is: |
366 | | {{{ |
367 | | SUBROUTINE g2l_obs(domain_p, step, dim_obs_f, dim_obs_l, mstate_f, mstate_l) |
368 | | |
369 | | INTEGER, INTENT(in) :: domain_p ! Current local analysis domain |
370 | | INTEGER, INTENT(in) :: step ! Current time step |
371 | | INTEGER, INTENT(in) :: dim_obs_f ! Dimension of full observation vector for model sub-domain |
372 | | INTEGER, INTENT(in) :: dim_obs_l ! Dimension of observation vector for local analysis domain |
373 | | REAL, INTENT(in) :: mstate_f(dim_obs_f) ! Full observation vector for model sub-domain |
374 | | REAL, INTENT(out) :: mstate_l(dim_obs_l) ! Observation vector for local analysis domain |
375 | | }}} |
376 | | |
377 | | The routine is called during the loop over the local analysis domains in the analysis step. It has to provide a local observation vector `mstate_l` for the observation domain that corresponds to the local analysis domain with index `domain_p`. Provided to the routine is the full observation vector `mstate_f` from which the local part has to be extracted. |
378 | | |
379 | | Hints: |
380 | | * The vector `mstate_f` that is provided to the routine is one of the observed state vectors that are produced by `U_obs_op_f`. |
381 | | * Some operations performed here are analogous to those required to initialize a local vector of observations in `U_init_obs_l`. If that routine reads first a full vector of observations (e.g. in `U_init_dim_obs_f`), this vector has to be restricted to the relevant observations for the current local analysis domain. For this operation, one can for example initialize an index array when `U_init_dim_obs_l` is executed. (Which happens before `U_g2l_obs`) |
382 | | |
383 | | |
384 | | === `U_init_obsvar` (init_obsvar_pdaf.F90) === |
385 | | |
386 | | This routine is used by the global filter algorithms SEIK, ETKF, and ESTKF as well as the local filters LSEIK, LETKF, ad LESTKF. The routine is only called if the adaptive forgetting factor is used (`type_forget=1` in the example implementation). The difference in this routine between global and local filters is that the global filters use 'global' while the local filters use 'full' quantities. |
387 | | |
388 | | The interface for this routine is: |
389 | | {{{ |
390 | | SUBROUTINE init_obsvar(step, dim_obs_f, obs_f, meanvar_f) |
391 | | |
392 | | INTEGER, INTENT(in) :: step ! Current time step |
393 | | INTEGER, INTENT(in) :: dim_obs_f ! Full dimension of observation vector |
394 | | REAL, INTENT(in) :: obs_f(dim_obs_f) ! Full observation vector |
395 | | REAL, INTENT(out) :: meanvar_f ! Mean observation error variance |
396 | | }}} |
397 | | |
398 | | The routine is called in the local filters before the loop over all local analysis domains is entered. The call is by the routine that computes an adaptive forgetting factor (`PDAF_set_forget`). |
399 | | The routine has to initialize an average full observation error variance, which should be consistent with the observation vector initialized in `U_init_ob_full`. |
400 | | |
401 | | |
402 | | Hints: |
403 | | * For a model with domain-decomposition one might use the mean variance for the model sub-domain of the calling process. Alternatively one can compute a mean variance for the full model domain using MPI communication (e.g. the function `MPI_allreduce`). |
404 | | * The observation vector `obs_p` is provided to the routine for the case that the observation error variance is relative to the value of the observations. |
405 | | * If the adaptive forgetting factor is not used, this routine has only to exist for the compilation, but it does not need functionality. |
406 | | |
407 | | |
408 | | === `U_init_obsvar_l` (init_obsvar_l_pdaf.F90) === |
409 | | |
410 | | This routine is used by all filter algorithms with domain-localization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm. The routine is only called if the local adaptive forgetting factor is used (`type_forget=2` in the example implementation). |
411 | | |
412 | | The interface for this routine is: |
413 | | {{{ |
414 | | SUBROUTINE init_obsvar_l(domain_p, step, dim_obs_l, obs_l, meanvar_l) |
415 | | |
416 | | INTEGER, INTENT(in) :: domain_p ! Current local analysis domain |
417 | | INTEGER, INTENT(in) :: step ! Current time step |
418 | | INTEGER, INTENT(in) :: dim_obs_l ! Local dimension of observation vector |
419 | | REAL, INTENT(in) :: obs_l(dim_obs_l) ! Local observation vector |
420 | | REAL, INTENT(out) :: meanvar_l ! Mean local observation error variance |
421 | | }}} |
422 | | |
423 | | The routine is called in the local filters during the loop over all local analysis domains by the routine that computes a local adaptive forgetting factor (`PDAF_set_forget_l`). The routine has to initialize a local mean observation error variance for all observations used for the analysis in the current local analysis domain. |
424 | | |
425 | | Hints: |
426 | | * If the local adaptive forgetting factor is not used, this routine has only to exist for the compilation, but it does not need functionality. |
427 | | |
428 | | |