| 242 | | || |
| 243 | | |
| 244 | | The domain-local filters (`template/analysis_step/global/PDAF_assimilate_LOCALTEMPLATE.F90`) have additional arguments to handle the localization of the state vectors and the localization of obsevations. However, they use the same generic routines |
| | 246 | || U_collect_state || Write model fields into a state vector || |
| | 247 | || U_distribute_state || Write a state vector into model fields || |
| | 248 | || U_init_dim_obs || Initialize observations, return observation dimension to PDAF. With PDAF-OMI, this is `init_dim_obs_pdafomi`. || |
| | 249 | || U_obs_op || Observation operator, return observed model state to PDAF. With PDAF-OMI, this is `obs_op_dafomi` || |
| | 250 | || U_init_obs || Return vector of observations to PDAF. When using PDAF-OMI, this routine is not visible to the user. || |
| | 251 | || U_prodRinvA || Compute product of **R**^-1^ with some matrix **A**. Return **R**^-1^**A** to PDAF. When using PDAF-OMI, this routine is not visible to the user. || |
| | 252 | || U_init_obsvar || Return mean observation error variance to PDAF. When using PDAF-OMI, this routine is not visible to the user. || |
| | 253 | || U_next_observation || Return number of time steps in next forecast phase, current model time and exit flag to PDAF || |
| | 254 | || U_prepoststep || Pre/poststep routine || |
| | 255 | |
| | 256 | The domain-local filters (`template/analysis_step/global/PDAF_assimilate_LOCALTEMPLATE.F90`) have additional arguments to handle the localization of the state vectors and the localization of obsevations. However, they use the same generic routines listed above. |
| | 257 | |
| | 258 | Some routines are hidden from the user when using the advanced interfaces like `PDAF3_assimilate` because they are provided by PDAF-OMI or PDAFlocal. However, the analysis routines are implemented without the assumption that PDAF-OMI or PDAFlocal are used. |
| | 259 | |
| | 260 | To get an overview of the available user-supplied call-back routines, you can looks into the file `src/PDAF_cb_procedures.F90`, which declared the interfaces of all call-back routines. |
| | 275 | This routine is the third infrastructure routine. This routine is usually called by `PDAF_assimilate_X`. Also `PDAF3_put_state` and other of the advanced 'assimilate'-routines call this routine. |
| | 276 | |
| | 277 | It provides the full interface in which all user-supplied routines are specified as arguments. For detailed information on using the routines with the full interface, please see the [wiki:ImplementationofAnalysisStep_noOMI Page on Implementing the Analysis Step using PDAF's full interface]. |
| | 278 | |
| | 279 | The routine is an infrastructure routine which is nearly identical for all DA methods. Its functionality is to write the forecasted fields into a state vector from the ensemble array, and check for the completeness of the forecast phase (particularly relevant for the 'flexible parallelization' variant. When the analysis has to be computed, the routine gathers the ensemble on the filter processes and then calls `PDAFX_update` for the analysis step. |
| | 280 | |
| | 281 | Except for the status flag `outflag`, all arguments of `PDAF_put_state_X` are the names of user-provided call-back routines. The user-supplied call-back routines specified as arguments are the same as for `PDAF_assimilate_X`, except that the routines `U_distribute_state` and `U_next_observation` are not present because these routines are already used in `PDAF_assimilate_X` or in `PDAF_get_state` in implementations for PDAF2. |
| | 282 | |
| | 283 | For implementing a new DA method, one mainly needs to adapt the call to PDAFX_update and the related specific call-back routines. |
| | 284 | |
| | 285 | |
| 253 | | The actual DA analysis update is computed in the routine `PDAF_X_update`. The routine is provided with the names of the call-back routines, the arrays as well of the relevant dimensions as describe above. |
| 254 | | |
| 255 | | The structure of the operations in `PDAF_X_update` can be designed freely when implementing a new DA method. The existing methods follow the recommended structure in which `PDAF_X_update` only performes preparations for the actual analysis update and then calls `PDAF_X_analysis` for the computation of the actual update. Here, the operations are different for global and local methods. |
| | 288 | This is the main routine for the data assimilation update, i.e. the analysis step. The routine initializes the observations by calling `PDAFobs_init`. It also calls `U_prepoststep` before and after the actual analysis update. |
| | 289 | |
| | 290 | For the domain-local filters also the local analysis loop is performed in this routine. In this loop the local state vector and the local observations are intialized before the routine for the actual local analysis update is called. |
| | 291 | |
| | 292 | The arguments of this routine are the call-back routines which are needed for the analysis step, and the framework arrays from PDAF (i.e. `state_p`, `ens_p`, `Ainv`, for smoothers also `sens_p`) that are used in the analysis step. |
| | 293 | |
| | 294 | The structure of the operations in `PDAF_X_update` can be designed freely when implementing a new DA method. To keep the initializations clearly separate of the numical calculations we recommend the provided structure in which the actual analysis update is computed in `PDAFX_analysis`. Note, that it is possible to have multiple PDAFX_analysis routines with unique names. Then, one can call different analysis algorithms according to the `subtype`. This is, for example used for the EnSRF and EAKF, which are implemted as different subtypes in `src/PDAF_assimilate_ensrf.F90` |
| | 295 | |
| | 296 | The template files in `PDAF_GLOBALTEMPLATE_update.F90` and `PDAF_LOCALTEMPLATE_update.F90` explain the typical steps which should be usable for other DA methods. |
| 259 | | |
| 260 | | |
| 261 | | === Global Methods === |
| 262 | | |
| 263 | | For global methods, e.g. ESTKF, ETKF, EnKF, NETF, the general structure of `PDAF_X_update` is: |
| 264 | | |
| 265 | | {{{ |
| 266 | | IF (fixed ensemble, i.e. subtype==2 or subtype==3) THEN |
| 267 | | Add central state and ensemble perturbations |
| 268 | | ENDIF |
| 269 | | |
| 270 | | CALL U_prepoststep # Perform pre/poststep before analysis update |
| 271 | | |
| 272 | | CALL PDAF_X_analysis # Compute global analysis update of ensemble |
| 273 | | |
| 274 | | CALL PDAF_smoother # Compute smoother update (Only for filters for which a smoother is implemented) |
| 275 | | |
| 276 | | CALL U_prepoststep # Perform pre/poststep after analysis update |
| 277 | | }}} |
| 278 | | You will find this structure e.g. in `PDAF_ESTKF_update.F90`. This structure is also valid for the LEnKF. The actual implementation, e.g. in `PDAF_ESTKF_update.F90` contains additional functionality e.g. calls for timers and normal screen output or debug outputs. Also there can be several different calls `PDAF_X_analysis` depending on the subtype of a filter. |
| 279 | | |
| 280 | | This general structure should be widely usable. Thus, when implementing a new global ensemble method one could copy an existing file and adapt it for the new method. The main functionality of the new DA method would then be implemented in `PDAF_X_analysis`. |
| 281 | | |
| 282 | | |
| 283 | | === Local Methods === |
| 284 | | |
| 285 | | For domain-localized Methods, e.g. LESTKF, LETKF or LNETF the general structure of `PDAF_X_update` is different from the global methods because the local analysis loop is contained in `PDAF_X_update` and only a local analysis routine is called. The general structure is: |
| 286 | | |
| 287 | | {{{ |
| 288 | | IF (fixed ensemble, i.e. subtype==2 or subtype==3) THEN |
| 289 | | Add central state and ensemble perturbations |
| 290 | | ENDIF |
| 291 | | |
| 292 | | CALL U_prepoststep # Perform pre/poststep before analysis update |
| 293 | | |
| 294 | | CALL U_init_n_domains # Determine number of local analysis domains (n_domains_p) |
| 295 | | |
| 296 | | # --- Perform global operations on observations and ensemble |
| 297 | | CALL U_init_dim_obs # Determine number of observations (PDAF-OMI routine init_dim_obs_pdafomi) |
| 298 | | DO i = 1, dim_ens |
| 299 | | CALL U_obs_op # Apply observation operator to all ensemble state (PDAF-OMI routine obs_op_pdafomi) |
| 300 | | ENDDO |
| 301 | | Compute mean forecast state (state_p) |
| 302 | | Compute mean observed forecast state (HXbar_f) |
| 303 | | # --- End of global operations |
| 304 | | |
| 305 | | # --- Perform local analysis loop |
| 306 | | DO i = 1, n_domains_p |
| 307 | | CALL U_init_dim_l # Set size of local state vector |
| 308 | | CALL U_init_dim_obs_l # Set number of observations within radius (PDAF-OMI routine init_dim_obs_l_pdafomi) |
| 309 | | DO i = 1, dim_ens |
| 310 | | CALL U_g2l_state # Get local state vector for all ensemble states |
| 311 | | ENDDO |
| 312 | | |
| 313 | | CALL PDAF_X_analysis # Compute local analysis update of ensemble |
| 314 | | |
| 315 | | DO i = 1, dim_ens |
| 316 | | CALL U_l2g_state # Update global state vector for all ensemble states |
| 317 | | ENDDO |
| 318 | | CALL U_l2g_state # Update global state vector for ensemble mean state |
| 319 | | |
| 320 | | CALL PDAF_smoother_local # Only for filters for which a smoother is implemented |
| 321 | | ENDDO |
| 322 | | # --- End of local analysis loop |
| 323 | | |
| 324 | | CALL U_prepoststep # Perform pre/poststep after analysis update |
| 325 | | }}} |
| 326 | | You will find this structure e.g. in `PDAF_LESTKF_update.F90`. This structure is not used for the LEnKF because this filter doe snot use a local analysis loop. The actual implementation, e.g. in `PDAF_LESTKF_update.F90` contains additional functionality e.g. calls for timers and normal screen output or debug outputs. Also there can be different calls `PDAF_X_analysis` depending on the subtype of a filter. |
| 327 | | |
| 328 | | This general structure should be widely usable for DA methods that use domain localization. Accordingly, when implementing a new local ensemble method one could copy an existing file and adapt it for the new method. The main functionality of the new DA method would then be implemented in `PDAF_X_analysis`. |
| | 300 | This routine computes the actual analysis update. thus, the actual numerics are in this routine. For ensemble-based Kalman filters and nonlinear filters, this is mainly linear algebra. For these, we recommend to use BLAS and LAPACK library functions, because they usually lead to optimal compute performance. |
| | 301 | |
| | 302 | The provided templates contain some of the of common parts of the analysis step computation, which are stepwise explained. This is based on the computations in the ETKF and LETKF methods. |
| | 303 | |
| | 304 | Note that for the 3D-Var schemes, different solver algorithms are used. These are organized in further subroutines (see e.g. `src/PAF_3dvar_analysis_cvt.F90`). |
| | 305 | |
| | 306 | |
