Changes between Version 1 and Version 2 of ImplementAnalysisLocal


Ignore:
Timestamp:
Nov 16, 2020, 1:46:38 PM (9 days ago)
Author:
lnerger
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • ImplementAnalysisLocal

    v1 v2  
    4545The interface when using one of the local filters is the following:
    4646{{{
    47   SUBROUTINE PDAF_assimilate_lestkf(U_collect_state, U_distribute_state, U_init_dim_obs_f, U_obs_op_f, &
     47  SUBROUTINE PDAF_assimilate_lestkf(U_collect_state, U_distribute_state, &
     48                                  U_init_dim_obs_f, U_obs_op_f, &
    4849                                  U_prepoststep, U_init_n_domains, U_init_dim_l, &
    4950                                  U_init_dim_obs_l, U_g2l_state, U_l2g_state, &
     
    5556 * [#U_init_dim_obs_finit_dim_obs_f_pdaf.F90 U_init_dim_obs_f]: The name of the user-supplied routine that provides the size of the full observation vector
    5657 * [#U_obs_op_fobs_op_f_pdaf.F90 U_obs_op_f]: The name of the user-supplied routine that acts as the full observation operator on some state vector
    57  * [#U_init_obs_finit_obs_f_pdaf.F90 U_init_obs_f]: The name of the user-supplied routine that initializes the full vector of observations
    58  * [#U_init_obs_linit_obs_l_pdaf.F90 U_init_obs_l]: The name of the user-supplied routine that initializes the vector of observations for a local analysis domain
    5958 * [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep]: The name of the pre/poststep routine as in `PDAF_get_state`
    60  * [#U_prodRinvA_lprodrinva_l_pdaf.F90 U_prodRinvA_l]: The name of the user-supplied routine that computes the product of the inverse of the observation error covariance matrix with some matrix provided to the routine by PDAF.
    6159 * [#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
    6260 * [#U_init_dim_linit_dim_l_pdaf.F90 U_init_dim_l]: The name of the routine that provides the state dimension for a local analysis domain
     
    6462 * [#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
    6563 * [#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 the provided local state vector
    66  * [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs]: The name of the routine that initializes a local observation vector from a full observation vector
    67  * [#U_init_obsvarinit_obsvar_pdaf.F90 U_init_obsvar]: The name of the user-supplied routine that provides a global mean observation error variance (This routine will only be executed, if an adaptive forgetting factor is used)
    68  * [#U_init_obsvar_linit_obsvar_l_pdaf.F90 U_init_obsvar_l]: The name of the user-supplied routine that provides a mean observation error variance for the local analysis domain (This routine will only be executed, if a local adaptive forgetting factor is used)
    6964 * [#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`.
    7065 * `status`: The integer status flag. It is zero, if `PDAF_assimilate_lestkf` is exited without errors.
     
    7570
    7671
    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
     74When 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
     76The 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)
    8883}}}
    8984
    9085== User-supplied routines ==
    9186
    92 Here, all user-supplied routines are described that are required in the call to `PDAF_assimilate_lestkf` or `PDAF_put_state_lestkf`. For some of the generic routines, we link to the page on [ModifyModelforEnsembleIntegration modifying the model code for the ensemble integration].
     87Here, all user-supplied 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].
    9388
    9489To indicate user-supplied 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.
     
    144139Hint:
    145140 * The routine is similar to `init_dim_obs` used for the global filters. However, with a domain-decomposed model `m_state_f` will contain parts of the state vector from neighboring model sub-domains. To make these parts accessible, some parallel communication will be necessary (The state information for a neighboring model sub-domain, will be in the memory of the process that handles that sub-domain). The example implementation in `testsuite/dummymodel_1d` uses the function `MPI_AllGatherV` for this communication.
    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.
    188141
    189142
     
    224177
    225178
    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 
    255179=== `U_init_n_domains` (init_n_domains_pdaf.F90) ===
    256180
     
    359283
    360284
    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 
    429285=== `U_next_observation` (next_observation_pdaf.F90) ===
    430286
     
    445301 1. [#U_init_dim_obs_finit_dim_obs_f_pdaf.F90 U_init_dim_obs_f]
    446302 1. [#U_obs_op_fobs_op_f_pdaf.F90 U_obs_op_f] (Called `dim_ens` times; once for each ensemble member)
    447  1. [#U_init_obs_finit_obs_f_pdaf.F90 U_init_obs_f] (Only executed, if the global adaptive forgetting factor is used (`type_forget=1` in the example implemention))
    448  1. [#U_init_obsvarinit_obsvar_pdaf.F90 U_init_obsvar] (Only executed, if the global adaptive forgetting factor is used (`type_forget=1` in the example implemention))
    449303
    450304In the loop over all local analysis domains, it is executed for each local analysis domain:
     
    452306 1. [#U_init_dim_obs_linit_dim_obs_l_pdaf.F90 U_init_dim_obs_l]
    453307 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)
    454  1. [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs] (A single call to localize the mean observed state)
    455  1. [#U_init_obs_linit_obs_l_pdaf.F90 U_init_obs_l]
    456  1. [#U_g2l_obsg2l_obs_pdaf.F90 U_g2l_obs] (`dim_ens` calls: one call to localize the observed part of each ensemble member)
    457  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))
    458308 1. [#U_prodRinvA_lprodrinva_l_pdaf.F90 U_prodRinvA_l]
    459309 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)
     
    462312 1. [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep] (Call to act on the analysis ensemble, called with (positive) value of the time step)
    463313
    464 In case of the routine `PDAF_assimilate_lestkf`, the following routines are executed after the analysis step:
     314In case of the routine `PDAF_assimilate_local`, the following routines are executed after the analysis step:
    465315 1. [#U_distribute_statedistribute_state_pdaf.F90 U_distribute_state]
    466316 1. [#U_next_observationnext_observation_pdaf.F90 U_next_observation]