Changes between Version 32 and Version 33 of ImplementAnalysislseik
 Timestamp:
 Sep 6, 2010, 8:52:59 AM (10 years ago)
Legend:
 Unmodified
 Added
 Removed
 Modified

ImplementAnalysislseik
v32 v33 64 64 == Usersupplied routines == 65 65 66 Here all usersupplied routines are described that are required in the call to `PDAF_put_state_lseik`. For some of the generic routines, we link to the page on [ModifyModelforEnsembleIntegration modifying the model code for the ensemble integration].66 Here, all usersupplied routines are described that are required in the call to `PDAF_put_state_lseik`. For some of the generic routines, we link to the page on [ModifyModelforEnsembleIntegration modifying the model code for the ensemble integration]. 67 67 68 68 To indicate usersupplied routines we use the prefix `U_`. In the template directory `templates/` these routines are provided in files with the routines name without this prefix. In the example implementation in `testsuite/src/dummymodel_1D` the routines exist without the prefix, but with the extension `_dummy_D.F90`. In the section titles below we provide the name of the template file in parentheses. … … 89 89 90 90 Some hints: 91 * It can be useful to not only determine the size of the observation vector at this point. One can also already gather information about the location of the observations, which will be used later, e.g. to implement the observation operator. In addition, one can already prepare an array that holds the full observation vector. This can be used later by `U_init_obs_l` to initialize a local vector of observations by selecting the relevant parts of the full observation vector. Arrays for the locations can be defined in a module like `mod_assimilation`.92 * The routine is similar to `init_dim_obs` used in the global filters. However, if the global filter is used with a domaindecomposed model, it only initializes the size of the observation vector for the local model subdomain. This is different for the local filters, as the local analysis also requires observational data from neighboring model subdomains. Anyway, one can base on an implemented routine `init_dim_obs` to implement `init_dim_obs_full`.91 * It can be useful to not only determine the size of the observation vector at this point. One can also already gather information about the location of the observations, which can be used later, e.g. to implement the observation operator. In addition, one can already prepare an array that holds the full observation vector. This can be used later by `U_init_obs_l` to initialize a local vector of observations by selecting the relevant parts of the full observation vector. The required arrays can be defined in a module like `mod_assimilation`. 92 * The routine is similar to `init_dim_obs` used in the global filters. However, if the global filter is used with a domaindecomposed model, it only initializes the size of the observation vector for the local model subdomain. This is different for the local filters, as the local analysis also requires observational data from neighboring model subdomains. Nonetheless, one can base on an implemented routine `init_dim_obs` to implement `init_dim_obs_full`. 93 93 94 94 === `U_obs_op_full` (obs_op_full.F90) === … … 99 99 {{{ 100 100 SUBROUTINE obs_op_full(step, dim_p, dim_obs_f, state_p, m_state_f) 101 101 102 INTEGER, INTENT(in) :: step ! Currrent time step 102 103 INTEGER, INTENT(in) :: dim_p ! PElocal dimension of state … … 108 109 The routine is called during the analysis step, before the loop over the local analysis domain is entered. It has to perform the operation of the observation operator acting on a state vector that is provided as `state_p`. The observed state has to be returned in `m_state_f`. It is the observed state corresponding to the 'full' observation vector. 109 110 110 Hint: 111 * If the observation operator involves a global operation, e.g. some global integration, while using domaindecompostion one has to gather the information from the other model domains using MPI communication. 112 * Analogously to the situation with `init_dim_obs_full`, the routine is similar to `init_dim_obs` used for the global filters. However, with a domaindecompoared model also here `m_state_f` will contain parts of the state vector from neighboring model subdomains. To make these parts accessible, some parallel communication will be necessary (The state information for a neighboring model subdomain, will be in the memory of the process that handles that subdomain). The example implementation in `testsuite/dummymodel_1d` uses the function `MPI_AllGatherV` for this communication. 111 Hint: 112 * Analogously to the situation with `init_dim_obs_full`, the routine is similar to `init_dim_obs` used for the global filters. However, with a domaindecomposed model `m_state_f` will contain parts of the state vector from neighboring model subdomains. To make these parts accessible, some parallel communication will be necessary (The state information for a neighboring model subdomain, will be in the memory of the process that handles that subdomain). The example implementation in `testsuite/dummymodel_1d` uses the function `MPI_AllGatherV` for this communication. 113 113 114 114 === `U_init_obs_full` (init_obs_full.F90) === … … 121 121 SUBROUTINE init_obs_full(step, dim_obs_f, observation_f) 122 122 123 INTEGER, INTENT(in) :: step ! Current time step124 INTEGER, INTENT(in) :: dim_obs_f ! Dimension of full observation vector123 INTEGER, INTENT(in) :: step ! Current time step 124 INTEGER, INTENT(in) :: dim_obs_f ! Dimension of full observation vector 125 125 REAL, INTENT(out) :: observation_f(dim_obs_f) ! Full observation vector 126 126 }}} … … 129 129 130 130 Hints: 131 * As for the other 'full' routines: While the global counterpart of this routine (`init_obs`) has to initialize the observation vector for the local model subdomain, the 'full' routine has to include observations that spatially belong to neighboring model subdomains. As an easy choice one can simply initialize a vector of all globally available observations.131 * 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 subdomain, the 'full' routine has to include observations that spatially belong to neighboring model subdomains. As an easy choice one can simply initialize a vector of all globally available observations. 132 132 * If the adaptive forgetting factor is not used, this routine only has to exist. However, no functionality is required. 133 133 … … 148 148 149 149 The routine is called during the analysis step during the loop over the local analysis domain. 150 It has to provide the vector of observations for analysis of the local analysis domain ofindex `domain_p` in `observation_l` for the current time step.151 152 Hints: 153 * For parallel efficiency the LSEIK algorithm is implemented in a way that first the full vectors are initialized. Thus, if `observation_f` has been initialized before `U_init_obs_local` is executed (e.g. by `U_init_dim_obs_full`), 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.150 It has to provide the vector of observations for the analysis in the local analysis domain owith index `domain_p` in `observation_l` for the current time step. 151 152 Hints: 153 * For parallel efficiency, the LSEIK algorithm is implemented in a way that first the full vectors are initialized. These are then restricted to the local analysis domain when the loop over all local analysis domains is executed. Thus, if `observation_f` has been initialized before `U_init_obs_local` is executed (e.g. by `U_init_dim_obs_full`), 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. 154 154 * The routine `U_init_dim_obs_local` is executed before this routine. Thus, if that routine already prepares the information which elements of `observation_f` are need for `observation_l` this information can be used efficiently here. 155 155 … … 178 178 }}} 179 179 180 The routine is called during the loop over the local analysis domains in the analysis step. 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`.180 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`. 181 181 182 182 This routine is also the place to perform observation localization. To initialize a vector of weights, the routine `PDAF_local_weights` can be called. The procedure is used in the example implementation and also demonstrated in the template routine. … … 201 201 202 202 The routine is called during the analysis step before the loop over the local analysis domains is entered. 203 It has to provide the number of local analysis domains. In case of a domaindecomposed model the number of local analysis domain for the model subd main of the calling process has to be initialized.203 It has to provide the number of local analysis domains. In case of a domaindecomposed model the number of local analysis domain for the model subdomain of the calling process has to be initialized. 204 204 205 205 Hints: … … 288 288 }}} 289 289 290 The routine is called during the loop over the local analysis domains in the analysis step. It has to provide the local state vector `state_l` that corresponds to the local analysis domain with index `domain_p`. Provided to the routine is the state vector `state_p`. With a domain decomposed model, this is the state for the local model subdomain.291 292 290 The routine is called during the loop over the local analysis domains in the analysis step. It has to initialize the part of the global state vector `state_p` that corresponds to the local analysis domain with index `domain_p`. Provided to the routine is the state vector `state_l` for the local analysis domain. 293 291 … … 312 310 }}} 313 311 314 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` ofwhich the local part has to be extracted.312 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. 315 313 316 314 Hints: … … 321 319 === `U_init_obsvar` (init_obsvar.F90) === 322 320 323 This routine is used by the global filter algorithms SEIK and ETKF as well as the local filters LSEIK and LETKF. The routine is only called if the adaptive forgetting factor is used (`type_forget=1` in the example impementation). 324 325 The interface for this routine is: 326 {{{ 327 SUBROUTINE init_obsvar(step, dim_obs_p, obs_p, meanvar) 328 329 INTEGER, INTENT(in) :: step ! Current time step 330 INTEGER, INTENT(in) :: dim_obs_p ! PElocal dimension of observation vector 331 REAL, INTENT(in) :: obs_p(dim_obs_p) ! PElocal observation vector 332 REAL, INTENT(out) :: meanvar ! Mean observation error variance 333 }}} 334 335 The routine is called in the global filters during the analysis or 336 by the routine that computes an adaptive forgetting factor (PDAF_set_forget). 337 The routine has to initialize the mean observation error variance. 338 For the global filters this should be the global mean. 321 This routine is used by the global filter algorithms SEIK and ETKF as well as the local filters LSEIK and LETKF. 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. 322 323 The interface for this routine is: 324 {{{ 325 SUBROUTINE init_obsvar(step, dim_obs_f, obs_f, meanvar_f) 326 327 INTEGER, INTENT(in) :: step ! Current time step 328 INTEGER, INTENT(in) :: dim_obs_f ! Full dimension of observation vector 329 REAL, INTENT(in) :: obs_f(dim_obs_f) ! Full observation vector 330 REAL, INTENT(out) :: meanvar_f ! Mean observation error variance 331 }}} 332 333 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). 334 The routine has to initialize a mean full observation error variance, which should be consistent with the observation vector initialized in `U_init_ob_full`. 335 339 336 340 337 Hints: 341 338 * For a model with domaindecomposition one might use the mean variance for the model subdomain 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`). 342 * The observation vector `obs_p` is provided to the ro tine for the case that the observation error variance is relative to the value of the observations.339 * 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. 343 340 * If the adaptive forgetting factor is not used, this routine has only to exist for the compilation, but it does not need functionality. 344 341 … … 346 343 === `U_init_obsvar_local` (init_obsvar_local.F90) === 347 344 348 This routine is used by the local filters LSEIK and LETKF. The routine is only called if the local adaptive forgetting factor is used (`type_forget=2` in the example imp ementation).345 This routine is used by the local filters LSEIK and LETKF. The routine is only called if the local adaptive forgetting factor is used (`type_forget=2` in the example implementation). 349 346 350 347 The interface for this routine is: … … 359 356 }}} 360 357 361 The routine is called in the local filters during the loop over all local analysis domains by the by theroutine that computes a local adaptive forgetting factor (PDAF_set_forget_local). The routine has to initialize a local mean observation error variance for all observations used for the analysis in the current local analysis domain.358 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_local). The routine has to initialize a local mean observation error variance for all observations used for the analysis in the current local analysis domain. 362 359 363 360 Hints: … … 366 363 == Execution order of usersupplied routines == 367 364 368 The usersupplied routines are executed in the order listed below. The order can be important as some routines can perform preparatory work for later routines. For example, `U_init_dim_obs_local` can prepare an index array that provides the information how to localize a 'full' vector of observations. Some hints one this are given with the descriptions of the routine interfaces above.365 The usersupplied routines are executed in the order listed below. The order can be important as some routines can perform preparatory work for routines executed later on during the analysis. For example, `U_init_dim_obs_local` can prepare an index array that provides the information how to localize a 'full' vector of observations. Some hints one this are given with the descriptions of the routine interfaces above. 369 366 370 367 Before the analysis step is called the following is executed: 371 368 1. [#U_collect_statecollect_state.F90 U_collect_state] (called once for each ensemble member) 372 369 373 When the ensemble integration of the forecast is completed the analysis step is executed. Before the loop over all local analysis domains, the following routines are executed:370 When the ensemble integration of the forecast is completed, the analysis step is executed. Before the loop over all local analysis domains, the following routines are executed: 374 371 1. [#U_prepoststepprepoststep_seik.F90 U_prepoststep] (call to handle the forecast, called with negative value of the time step) 375 372 1. [#U_init_n_domainsinit_n_domains.F90 U_init_n_domains] … … 377 374 1. [#U_obs_op_fullobs_op_full.F90 U_obs_op_full] (Called `dim_ens` times; once for each ensemble member) 378 375 1. [#U_init_obs_fullinit_obs_full.F90 U_init_obs_full] (Only executed, if the global adaptive forgetting factor is used (`type_forget=1` in the example implemention)) 379 1. [#U_init_obsvarinit_obsvar.F90 U_init_obsvar] (Only execute , if the global adaptive forgetting factor is used (`type_forget=1` in the example implemention))376 1. [#U_init_obsvarinit_obsvar.F90 U_init_obsvar] (Only executed, if the global adaptive forgetting factor is used (`type_forget=1` in the example implemention)) 380 377 381 378 The the loop over all local analysis domains, it is executed: … … 386 383 1. [#U_init_obs_localinit_obs_local.F90 U_init_obs_local] 387 384 1. [#U_global2local_obsglobal2local_obs.F90 U_global2local_obs] (`dim_ens` calls; one call to localize the observed part of each ensemble member) 388 1. [#U_init_obsvar_localinit_obsvar_local.F90 U_init_obsvar_local] (Only , if the local adaptive forgetting factor is used (`type_forget=2` in the example implementation))385 1. [#U_init_obsvar_localinit_obsvar_local.F90 U_init_obsvar_local] (Only called, if the local adaptive forgetting factor is used (`type_forget=2` in the example implementation)) 389 386 1. [#U_prodRinvA_localprodrinva_local.F90 U_prodRinvA_local] 390 387 1. [#U_local2global_statelocal2global_state.F90 U_local2global_state] (Called `dim_ens+1` times: Once for each ensemble member and once for the mean state estimate)