Changes between Version 2 and Version 3 of ImplementAnalysisLocal
 Timestamp:
 Nov 16, 2020, 3:22:39 PM (2 years ago)
Legend:
 Unmodified
 Added
 Removed
 Modified

ImplementAnalysisLocal
v2 v3 29 29 With Version 1.16 of PDAF we introduced PDAFOMI (observation module infrastructure). With OMI we provide generic routines for the analysis step, which only distinguish global and local filters. This page describes the implementation of the analysis step for domainlocal filters (LESTKF, LETKF, LNETF, LSEIK). 30 30 31 For the analysis step of the local filters, several operations related to the observations are needed. These operations are requested by PDAF by call s to routines supplied by the user and provided in the OMI structure. The names of the required routines are specified in the call to the routine `PDAF_assimilate_local` in the fullyparallel implementation (or `PDAF_put_state_local` for the 'flexible' implementation) described below. With regard to the parallelization, all these routines (except `U_collect_state`, `U_distribute_state`, and `U_next_observation`) are executed by the filter processes (`filterpe=.true.`) only.32 33 For completeness we discuss here all usersupplied routines that are specified in the interface to `PDAF_assimilate_local `. Many of the routines are identical to those used for the global filters. Hence, when the usersupplied routines for the global filters have been already implemented, one can base on these routines to speed up the implementation. Due to this, it can also be reasonable to first fully implement a global filter version and subsequently implement the corresponding localized filter by modifying and extending the global routines.34 35 36 == `PDAF_assimilate_local ` ==31 For the analysis step of the local filters, several operations related to the observations are needed. These operations are requested by PDAF by callback routines supplied by the user and provided in the OMI structure. The names of the required routines are specified in the call to the routine `PDAF_assimilate_local_omi` in the fullyparallel implementation (or `PDAF_put_state_local_omi` for the 'flexible' implementation) described below. With regard to the parallelization, all these routines (except `U_collect_state`, `U_distribute_state`, and `U_next_observation`) are executed by the filter processes (`filterpe=.true.`) only. 32 33 For completeness we discuss here all usersupplied routines that are specified in the interface to `PDAF_assimilate_local_omi`. Many of the routines are identical to those used for the global filters. Hence, when the usersupplied routines for the global filters have been already implemented, one can base on these routines to speed up the implementation. Due to this, it can also be reasonable to first fully implement a global filter version and subsequently implement the corresponding localized filter by modifying and extending the global routines. 34 35 36 == `PDAF_assimilate_local_omi` == 37 37 38 38 The general aspects of the filterspecific routines `PDAF_assimilate_*` have been described on the page [ModifyModelforEnsembleIntegration Modification of the model code for the ensemble integration] and its subpage on [InsertAnalysisStep inserting the analysis step]. The routine is used in the fullyparallel implementation variant of the data assimilation system. When the 'flexible' implementation variant is used, the routines `PDAF_put_state_*' is used as described further below. 39 The interface for the routine `PDAF_assimilate_local` contains several routine names for routines that operate on the local analysis domains (marked by `_l` at the end of the routine name). In addition, there are names for routines that consider all available observations required to perform local analyses with LESTKF within some subdomain of a domaindecomposed model (marked by `_f` at the end of the routine name). In case of a serial execution of the assimilation program, these will be all globally available observations. However, if the program is executed with parallelization, this might be a smaller set of observations. 40 41 To explain the difference, it is assumed, for simplicity, that a local analysis domain consists of a single vertical column of the model grid. In addition, we assume that the domain decomposition splits the global model domain by vertical boundaries as is typical for ocean models and that the observations are spatially distributed observations of model fields that are part of the state vector. Under these assumptions, the situation is the following: When a model uses domain decomposition, the LESTKF is executed such that for each model subdomain a loop over all local analysis domains (e.g. vertical columns) that belong to the model subdomain is performed. As each model subdomain is treated by a different process, all loops are executed parallel to each other. 42 43 For the update of each single vertical column, observations from some larger domain surrounding the vertical column are required. If the influence radius for the observations is sufficiently small there will be vertical columns for which the relevant observations reside completely inside the model subdomain of the process. However, if a vertical column is considered that is located close to the boundary of the model subdomain, there will be some observations that don't belong spatially to the local model subdomain, but to a neighboring model subdomain. Nonetheless, these observations are required on the local model subdomain. A simple way to handle this situation is to initialize for each process all globally available observations, together with their coordinates. While this method is simple, it can also become inefficient if the assimilation program is executed using a large number of processes. As for an initial implementation it is usually not the concern to obtain the highest parallel efficiency, the description below assumes that 'full' refers to the global model domain. 39 40 The interface for the routine `PDAF_assimilate_local_omi` contains two routine names for routines that operate on the local analysis domains (marked by `_l` at the end of the routine name). Further there are routines that convert between a local and a global model state vector (`U_g2l_state` and `U_l2g_state`). 41 Here, we list the full interface of the routine. Subsequently, the usersupplied routines specified in the call is explained. 44 42 45 43 The interface when using one of the local filters is the following: 46 44 {{{ 47 SUBROUTINE PDAF_assimilate_l estkf(U_collect_state, U_distribute_state, &48 U_init_dim_obs _f, U_obs_op_f, &45 SUBROUTINE PDAF_assimilate_local_omi(U_collect_state, U_distribute_state, & 46 U_init_dim_obs, U_obs_op, & 49 47 U_prepoststep, U_init_n_domains, U_init_dim_l, & 50 48 U_init_dim_obs_l, U_g2l_state, U_l2g_state, & … … 54 52 * [#U_collect_statecollect_state_pdaf.F90 U_collect_state]: The name of the usersupplied routine that initializes a state vector from the array holding the ensemble of model states from the model fields. This is basically the inverse operation to `U_distribute_state` used in [ModifyModelforEnsembleIntegration#PDAF_get_state PDAF_get_state] and also here. 55 53 * [#U_distribute_statedistribute_state_pdaf.F90 U_distribute_state]: The name of a user supplied routine that initializes the model fields from the array holding the ensemble of model state vectors. 56 * [#U_init_dim_obs _finit_dim_obs_f_pdaf.F90 U_init_dim_obs_f]: The name of the usersupplied routine that provides the size of the fullobservation vector57 * [#U_obs_op _fobs_op_f_pdaf.F90 U_obs_op_f]: The name of the usersupplied routine that acts as the fullobservation operator on some state vector54 * [#U_init_dim_obscallback_obs_pdafomi.F90 U_init_dim_obs]: The name of the usersupplied routine that initializes the observation information and provides the size of observation vector 55 * [#U_obs_opcallback_obs_pdafomi.F90 U_obs_op]: The name of the usersupplied routine that acts as the observation operator on some state vector 58 56 * [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep]: The name of the pre/poststep routine as in `PDAF_get_state` 59 57 * [#U_init_n_domainsinit_n_domains_pdaf.F90 U_init_n_domains]: The name of the routine that provides the number of local analysis domains … … 61 59 * [#U_init_dim_obs_linit_dim_obs_l_pdaf.F90 U_init_dim_obs_l]: The name of the routine that initializes the size of the observation vector for a local analysis domain 62 60 * [#U_g2l_stateg2l_state_pdaf.F90 U_g2l_state]: The name of the routine that initializes a local state vector from the global state vector 63 * [#U_l2g_statel2g_state_pdaf.F90 U_l2g_state]: The name of the routine that initializes the corresponding part of the global state vector from the theprovided local state vector61 * [#U_l2g_statel2g_state_pdaf.F90 U_l2g_state]: The name of the routine that initializes the corresponding part of the global state vector from the provided local state vector 64 62 * [#U_next_observationnext_observation.F90 U_next_observation]: The name of a user supplied routine that initializes the variables `nsteps`, `timenow`, and `doexit`. The same routine is also used in `PDAF_get_state`. 65 * `status`: The integer status flag. It is zero, if `PDAF_assimilate_l estkf` is exited without errors.63 * `status`: The integer status flag. It is zero, if `PDAF_assimilate_local_omi` is exited without errors. 66 64 67 65 Note: … … 70 68 71 69 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 usersupplied routines `U_distribute_state` and `U_next_observation` are missing.70 == `PDAF_put_state_local_omi` == 71 72 When the 'flexible' implementation variant is chosen for the assimilation system, the routine `PDAF_put_state_local_omi` has to be used instead of `PDAF_assimilate_local_omi`. 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_omi` with the exception the specification of the usersupplied routines `U_distribute_state` and `U_next_observation` are missing. 75 73 76 74 The interface when using one of the local filters is the following: 77 75 {{{ 78 SUBROUTINE PDAF_put_state_l estkf(U_collect_state, &79 U_init_dim_obs _f, U_obs_op_f, &76 SUBROUTINE PDAF_put_state_local_omi(U_collect_state, & 77 U_init_dim_obs, U_obs_op, & 80 78 U_prepoststep, U_init_n_domains, U_init_dim_l, & 81 79 U_init_dim_obs_l, U_g2l_state, U_l2g_state, & … … 85 83 == Usersupplied routines == 86 84 87 Here, all usersupplied routines are described that are required in the call to `PDAF_assimilate_local ` or `PDAF_put_state_local`. For some of the generic routines, we link to the page on [ModifyModelforEnsembleIntegration modifying the model code for the ensemble integration].85 Here, all usersupplied routines are described that are required in the call to `PDAF_assimilate_local_omi` or `PDAF_put_state_local_omi`. For some of the generic routines, we link to the page on [ModifyModelforEnsembleIntegration modifying the model code for the ensemble integration]. 88 86 89 87 To indicate usersupplied routines we use the prefix `U_`. In the template directory `templates/` as well as in the example implementation in `testsuite/src/dummymodel_1D` these routines exist without the prefix, but with the extension `_pdaf.F90`. In the section titles below we provide the name of the template file in parentheses. … … 93 91 === `U_collect_state` (collect_state_pdaf.F90) === 94 92 95 This routine is independent fromthe filter algorithm used.96 See the mape on [InsertAnalysisStep#U_collect_statecollect_state_pdaf.F90 inserting the analysis step] for the description of this routine.93 This routine is independent of the filter algorithm used. 94 See the page on [InsertAnalysisStep#U_collect_statecollect_state_pdaf.F90 inserting the analysis step] for the description of this routine. 97 95 98 96 === `U_distribute_state` (distribute_state_pdaf.F90) === … … 102 100 103 101 104 === `U_init_dim_obs_f` (init_dim_obs_f_pdaf.F90) === 105 106 This routine is used by all filter algorithms with domainlocalization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm. 107 108 The interface for this routine is: 109 {{{ 110 SUBROUTINE init_dim_obs_f(step, dim_obs_f) 102 === `U_init_dim_obs` (callback_obs_pdafomi.F90) === 103 104 The interface for this routine is: 105 {{{ 106 SUBROUTINE init_dim_obs(step, dim_obs_f) 111 107 112 108 INTEGER, INTENT(in) :: step ! Current time step … … 114 110 }}} 115 111 116 The routine is called at the beginning of each analysis step, before the loop over all local analysis domains is entered. It has to initialize the size `dim_obs_f` of the full observation vector according to the current time step. For simplicity, `dim_obs_f` can be the size for the global model domain. 117 118 Some hints: 119 * 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`. 120 * 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_f`. 121 122 === `U_obs_op_f` (obs_op_f_pdaf.F90) === 123 124 This routine is used by all filter algorithms with domainlocalization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm. 125 126 The interface for this routine is: 127 {{{ 128 SUBROUTINE obs_op_f(step, dim_p, dim_obs_f, state_p, m_state_f) 112 The routine is called at the beginning of each analysis step. For PDAF, it has to initialize the size `dim_obs_f` of the observation vector according to the current time step. `dim_obs_f` is usually be the size for the full model domain. When a domaindecomposed model is used, `dim_obs_f` can be reduced to those observations relevant for the local analysis loop in a process domain. 113 114 With PDAFOMI, the routine just calls a routine from the observation module for each observation type. 115 116 117 === `U_obs_op` (callback_obs_pdafomi.F90) === 118 119 The interface for this routine is: 120 {{{ 121 SUBROUTINE obs_op(step, dim_p, dim_obs_p, state_p, m_state_f) 129 122 130 123 INTEGER, INTENT(in) :: step ! Currrent time step 131 124 INTEGER, INTENT(in) :: dim_p ! PElocal dimension of state 132 INTEGER, INTENT(in) :: dim_obs_f ! Dimension of the fullobserved state125 INTEGER, INTENT(in) :: dim_obs_f ! Dimension of observed state 133 126 REAL, INTENT(in) :: state_p(dim_p) ! PElocal model state 134 127 REAL, INTENT(out) :: m_state_f(dim_obs_f) ! Full observed state 135 128 }}} 136 129 137 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, which 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. 138 139 Hint: 140 * 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. 130 The routine is called during the analysis step. 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`. 131 132 With PDAFOMI, the routine just calls a routine from the observation module for each observation type. 141 133 142 134 … … 179 171 === `U_init_n_domains` (init_n_domains_pdaf.F90) === 180 172 181 This routine is used by all filter algorithms with domainlocalization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm.182 183 173 The interface for this routine is: 184 174 {{{ … … 197 187 198 188 === `U_init_dim_l` (init_dim_l_pdaf.F90) === 199 200 This routine is used by all filter algorithms with domainlocalization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm.201 189 202 190 The interface for this routine is: … … 217 205 218 206 === `U_init_dim_obs_l` (init_dim_obs_l_pdaf.F90) === 219 220 This routine is used by all filter algorithms with domainlocalization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm.221 207 222 208 The interface for this routine is: … … 233 219 It has to initialize in `dim_obs_l` the size of the observation vector used for the local analysis domain with index `domain_p`. 234 220 235 Some hints: 236 * Usually, the observations to be considered for a local analysis are those which reside within some distance from the local analysis domain. Thus, if the local analysis domain is a single vertical column of the model grid and if the model grid is a regular ijgrid, then one could use some range of i/j indices to select the observations and determine the local number of them. More generally, one can compute the physical distance of an observation from the local analysis domain and decide on this basis, if the observation has to be considered. 237 * In the loop over the local analysis domains, the routine is always called before `U_init_obs_l` is executed. Thus, as `U_init_dim_obs_l` has to check which observations should be used for the local analysis domain, one can already initialize an integer array that stores the index of observations to be considered. This index should be the position of the observation in the array `observation_f`. With this, the initialization of the local observation vector in `U_init_obs_l` can be sped up. 238 * For PDAF, we could not join the routines `U_init_dim_obs_l` and `U_init_obs_l`, because the array for the local observations is allocated internally to PDAF after its size has been determined in `U_init_dim_obs_l`. 221 With PDAFOMI, the routine just calls a routine from the observation module for each observation type. PDAFOMI will perform the necessary intializations based on the coordinates of the observations. 239 222 240 223 241 224 === `U_g2l_state` (g2l_state_pdaf.F90) === 242 243 This routine is used by all filter algorithms with domainlocalization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm.244 225 245 226 The interface for this routine is: … … 259 240 Hints: 260 241 * In the simple case that a local analysis domain is a single vertical column of the model grid, the operation in this routine would be to take out of `state_p` the data for the vertical column indexed by `domain_p`. 242 * Usually, one can initialize the indices of the local state vector elements in the global state vector in `U_init_dim_l` and just use these here. 261 243 262 244 263 245 === `U_l2g_state` (l2g_state_pdaf.F90) === 264 265 This routine is used by all filter algorithms with domainlocalization (LSEIK, LETKF, LESTKF) and is independent of the particular algorithm.266 246 267 247 The interface for this routine is: … … 281 261 Hints: 282 262 * In the simple case that a local analysis domain is a single vertical column of the model grid, the operation in this routine would be to write into `state_p` the data for the vertical column indexed by `domain_p`. 263 * Usually, one can initialize the indices of the local state vector elements in the global state vector in `U_init_dim_l` and just use these here. 283 264 284 265 … … 291 272 == Execution order of usersupplied routines == 292 273 293 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 the efficient implementation strategy are given with the descriptions of the routine interfaces above.274 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_l` can prepare an index array that provides the information how to localize a global state vector. Some hints one the efficient implementation strategy are given with the descriptions of the routine interfaces above. 294 275 295 276 Before the analysis step is called the following is executed: … … 299 280 1. [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep] (Call to act on the forecast ensemble, called with negative value of the time step) 300 281 1. [#U_init_n_domainsinit_n_domains_pdaf.F90 U_init_n_domains] 301 1. [#U_init_dim_obs _finit_dim_obs_f_pdaf.F90 U_init_dim_obs_f]302 1. [#U_obs_op _fobs_op_f_pdaf.F90 U_obs_op_f] (Called `dim_ens` times; once for each ensemble member)282 1. [#U_init_dim_obsback_obs_pdafomi.F90 U_init_dim_obs] 283 1. [#U_obs_opcallback_obs_pdafomi.F90 U_obs_op] (Called `dim_ens` times; once for each ensemble member) 303 284 304 285 In the loop over all local analysis domains, it is executed for each local analysis domain: … … 306 287 1. [#U_init_dim_obs_linit_dim_obs_l_pdaf.F90 U_init_dim_obs_l] 307 288 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) 308 1. [#U_prodRinvA_lprodrinva_l_pdaf.F90 U_prodRinvA_l]309 289 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) 310 290 … … 312 292 1. [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep] (Call to act on the analysis ensemble, called with (positive) value of the time step) 313 293 314 In case of the routine `PDAF_assimilate_local `, the following routines are executed after the analysis step:294 In case of the routine `PDAF_assimilate_local_omi`, the following routines are executed after the analysis step: 315 295 1. [#U_distribute_statedistribute_state_pdaf.F90 U_distribute_state] 316 296 1. [#U_next_observationnext_observation_pdaf.F90 U_next_observation]