Changes between Version 30 and Version 31 of ModifyModelforEnsembleIntegration


Ignore:
Timestamp:
Apr 29, 2014, 4:31:06 PM (6 years ago)
Author:
lnerger
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • ModifyModelforEnsembleIntegration

    v30 v31  
    99<li><a href="InitPdaf">Initialization of PDAF</a></li>
    1010<li>Modifications for ensemble integration</li>
     11  <li>&nbsp;&nbsp;&nbsp;<a href="InsertAnalysisStep">Fully parallel implementation variant</a></li>
     12  <li>&nbsp;&nbsp;&nbsp;<a href="ExternalModelLoop">Flexible implementation variant</a></li>
    1113<li><a href="ImplementationofAnalysisStep">Implementation of the analysis step</a></li>
    1214<li><a href="AddingMemoryandTimingInformation">Memory and timing information</a></li>
     
    1921== Overview ==
    2022
    21 Numerical models are typically implemented for normal integration of some initial state. For the data assimilation with filter algorithm, an ensemble of model states has to be integrated for limited time until observations are available and an analysis step of the filter is computed. Subsequently, the updated ensemble has to be integrated again. To allow for these alternating ensemble integrations and analysis steps the model code has to be extended. The recommended implementation strategy for PDAF is to add an additional loop outside of the regular time-stepping loop of the model. This strategy has the potential to reduce the required chances in the model code to the minimum. In addition, a routine that simulates model errors might be required to be inserted into the time stepping loop of the model. The required extensions are described below.
     23Numerical models are typically implemented for normal integration of some initial state. For the data assimilation with a filter algorithm, an ensemble of model states has to be integrated for limited time until observations are available and an analysis step of the filter is computed. Subsequently, the updated ensemble has to be integrated further. To allow for the interruption of the integrations by the analysis step, the model code has to be extended. As described on the page on the [ImplementationConceptOnline implementation concept of the online mode], there are two options for the ensemble integration:
     24 * '''fully parallel''': For this implementation one needs to use a parallel computed with a sufficient number of processes such that the data assimilation program run be run with a concurrent time stepping of all ensemble states. Thus, if one runs each mode task with '''n''' processes and the ensemble has '''m''' members, the program has to run with '''n''' times '''m''' processes. This parallelism allows for a simplified implementation as each model task integrated only one model state and the model is always going forward in time.
     25 * '''flexible''': This variant allows to run the assimilation program in a way so that a model task (set of processors running one model integration) can propagate several ensemble states successively. In the extreme case, this could mean that one only a a single model task that is successively performing the integration of all ensemble states. The implementation for this variant is a bit more complicated, because one has to ensure that the model can jump back in time.
    2226
    23 Operations that are specific to the model and to the assimilated observations are performed by routines that are supplied by the user. These are called through the defined interface of PDAF. Generally, these user-supplied routines have to provide quite elementary operations, like initializing a model state vector for PDAF from model fields or providing the vector of observations. PDAF provides examples for these routines and templates that can be used as the basis for the implementation. As only the interface of these routines is specified, the user can implement the routines like a routine of the model. Thus, the implementation of these routines should not be difficult.
     27The extension to the model code for both cases is depicted in the figure below (See also the page on the [ImplementationConceptOnline implementation concept of the online mode.]) As at this stage of the implementation the calls to `init_parallel_pdaf` and `init_pdaf` are already inserted into the code, the difference is in the addition of routines for the time stepping.
    2428
    25 == External ensemble loop ==
    26 
    27 The external loop for the ensemble integration has to enclose the time stepping loop of the model. Next to the external loop, a control structure for exiting the external loop as well as two calls to subroutines of PDAF have to be added. These are the calls to `PDAF_get_state` and a filter-specific routine like `PDAF_put_state_seik` for the SEIK filter. Both routines are described in sections below.
    28 
    29 The extended model code can look like this for the SEIK filter:
    30  {{{
    31   pdaf_modelloop: DO 
    32 
    33      CALL PDAF_get_state(nsteps, ..., doexit, ...)
    34 
    35      ! Check whether forecast has to be performed
    36      ifcontrol: IF (doexit /= 1) THEN
    37      
    38         IF (nsteps > 0) THEN
    39 
    40           ... Time stepping code of the model ...         
    41 
    42         END IF
    43 
    44         CALL PDAF_put_state_seik(...)
    45 
    46      ELSE ifcontrol
    47 
    48         ! No more assimilation work; exit loop
    49         EXIT pdaf_modelloop
    50 
    51      END IF ifcontrol
    52 
    53   END DO pdaf_modelloop
    54 }}}
    55 In this example, which is taken from the example implementation in `testsuite/src/dummymodel_1D`, we use an unconditional DO loop (while loop). The exit flag `doexit` for this loop is set within `PDAF_get_state`. In addition, the variable `nsteps` is initialized, which defines the number of time steps to be performed during the current forecast phase. Thus, we only execute the time stepping code if `nsteps>0`. (If this has to be implemented using an IF-clause as in the example should be checked for the particular code).
    56 
    57 == `PDAF_get_state` ==
    58 
    59 The routine `PDAF_get_state` has the purpose to initialize the information, whether further model integrations have to be computed and how many time steps have to be performed. In addition, the model fields to be propagated are initialized from the array holding the ensemble states.
    60 
    61 The interface of `PDAF_get_state` is the following:
    62 {{{
    63   SUBROUTINE PDAF_get_state(nsteps, timenow, doexit, U_next_observation, U_distribute_state, &
    64                             U_prepoststep, status)
    65 }}}
    66 with the following arguments:
    67  * `nsteps`: An integer specifying upon exit the number of time steps to be performed
    68  * `timenow`: A real specifying upon exit the current model time. 
    69  * `doexit`: An integer variable defining whether the assimilation process is completed and the program should exit the while loop. For compatibility 1 should be used for exit, 0 for continuing in the loop.
    70  * [#U_next_observationnext_observation.F90 U_next_observation]: The name of a user supplied routine that initializes the variables `nsteps`, `timenow`, and `doexit`
    71  * [#U_distribute_statedistribute_state.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
    72  * [#U_prepoststepprepoststep_seik.F90 U_prepoststep]: The name of a user supplied routine that is called before and after the analysis step. Here the user has the possibility to access the state ensemble and can e.g. compute estimated variances or can write the ensemble states the state estimate into files.
    73  * `status`: The integer status flag. It is zero, if `PDAF_get_state` is exited without errors.
    74 
    75 PDAF also has a [PdafSimplifiedInterface Simplified Interface] providing the routine `PDAF_get_state_si`. In the simplified interface, the names of all user-supplied call back routines are predefined such that they not appear in the call to `PDAF_get_state_si`. More information on the pre-defined names is provided in the [PdafSimplifiedInterface documentation of PDAF's simplified interface].
    76 
    77 == `PDAF_put_state_X` ==
    78 
    79 There is a separate routine `PDAF_put_state_X` for each of the filter algorithms. The name of the routine includes the name of the filter at its end (instead of `X`). The purpose of the `PDAF_put_state_X` routines is to write back the forecast model fields into the array holding the ensemble of model state vectors. In addition, the routine checks if the current forecast phase is completed. If not, the routine is exited and the next cycle of the ensemble loop is performed. If the current forecast phase is completed, the routine executes the analysis step of the chosen filter algorithm. The interface to each put-state routine is specific for each filter algorithm, because the names of several user-supplied routines have to be specified, which are specific for each filter algorithm. However, at the stage of implementing the ensemble integration only the first and last arguments of the routines are relevant.
    80 
    81 For example, the interface when using the SEIK filter is the following:
    82 {{{
    83   SUBROUTINE PDAF_put_state_seik(U_collect_state, U_init_dim_obs, U_obs_op, &
    84                                  U_init_obs, U_prepoststep, U_prodRinvA, U_init_obsvar, status)
    85 }}}
    86 At this state of the implementation only these arguments are relevant:
    87  * [#U_collect_statecollect_state.F90 U_collect_state]: The name of the user-supplied 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_dist_state` used in `PDAF_get_state`
    88  * `status`: The integer status flag. It is zero, if PDAF_get_state is exited without errors.
    89 
    90 The other arguments are names of user-supplied subroutines that are only executed if the analysis step is executed (See the section [#Compilationandtesting 'Compilation and testing'] for how to provide these routines for compilation at this stage). These routines are explained in the next section of the implementation guide ([ImplementationofAnalysisStep Implementation of the Analysis step]) separately for each available filter algorithm.
    91 
    92 PDAF also has a [PdafSimplifiedInterface Simplified Interface] providing the routine `PDAF_out_state_X_si`. In the simplified interface, the names of all user-supplied call back routines are predefined such that they not appear in the call to `PDAF_put_state_X_si`. More information on the pre-defined names is provided in the [PdafSimplifiedInterface documentation of PDAF's simplified interface].
    93 
    94 == User-supplied routines ==
    95 
    96 Here, only the user-supplied routines are discussed that are required at this stage of the implementation (that is, the ensemble integration). For testing (see [#Compilationandtesting 'Compilation and testing']), all routines need to exist, but only those described here in detail need to be implemented with functionality.
    97 
    98 To 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.
    99 
    100 === `U_next_observation` (next_observation_pdaf.F90) ===
    101 
    102 The interface for this routine is
    103 {{{
    104 SUBROUTINE next_observation(stepnow, nsteps, doexit, timenow)
    105 
    106   INTEGER, INTENT(in)  :: stepnow  ! Number of the current time step
    107   INTEGER, INTENT(out) :: nsteps   ! Number of time steps until next obs
    108   INTEGER, INTENT(out) :: doexit   ! Whether to exit forecasting (1 for exit)
    109   REAL, INTENT(out)    :: timenow  ! Current model (physical) time
    110 }}}
    111 
    112 The routine is called once at the beginning of each forecast phase. It is executed by all processes that participate in the model integrations.
    113 
    114 Based on the information of the current time step, the routine has to define the number of time steps `nsteps` for the next forecast phase. In addition, the flag `doexit` has to be initialized to provide the information if the external ensemble loop can be exited. `timenow` is the current model time. This variable should also be initialized. It is particularly important, if an ensemble task integrates more than one model state. In this case `timenow` can be used to correctly jump back in time.
    115 
    116 Some hints:
    117  * If the time interval between successive observations is known, `nsteps` can be simply initialized by dividing the time interval by the size of the time step
    118  * `doexit` should be 0 to continue the assimilation process. In most cases `doexit` is set to 1, when `PDAF_get_state` is called after the last analysis for which observations are available.
    119  * At the first call to `U_next_obs` the variable `timenow` should be initialized with the current model time. At the next call a forecast phase has been completed. Thus, the new value of `timenow` follows from the timer interval for the previous forecast phase.
    120 
    121 === `U_distribute_state` (distribute_state_pdaf.F90) ===
    122 
    123 The interface for this routine is
    124 {{{
    125 SUBROUTINE distribute_state(dim_p, state_p)
    126 
    127   INTEGER, INTENT(in) :: dim_p           ! State dimension for PE-local model sub-domain
    128   REAL, INTENT(inout) :: state_p(dim_p)  ! State vector for PE-local model sub-domain
    129 }}}
    130 
    131 This routine is called during the forecast phase as many times as there are states to be integrated by a model task. Again, the routine is executed by all processes that belong to model tasks.
    132 
    133 When the routine is called a state vector `state_p` and its size `dim_p` are provided. As the user has defined how the model fields are stored in the state vector, one can initialize the model fields from this information. If the model is not parallelized, `state_p` will contain a full state vector. If the model is parallelized using domain decomposition, `state_p` will contain the part of the state vector that corresponds to the model sub-domain for the calling process.
    134 
    135 Some hints:
    136  * If the state vector does not include all model fields, it can be useful to keep a separate array to store those additional fields. This array has to be kept separate from PDAF, but can be defined using a module like `mod_assimilation`.
     29[[Image(//pics/da_extension2x.png)]]
     30[[BR]]'''Figure 1:''' (left) Generic structure of a model code, (center) extension for ''fully-parallel'' data assimilation system with PDAF, (right) extension for ''flexible'' data assimilation system with PDAF.
    13731
    13832
    139 === `U_prepoststep` (prepoststep_ens_pdaf.F90) ===
    14033
    141 The interface of the routine is identical for all filters. However, the particular operations that are performed in the routine can be specific for each filter algorithm. Here, we exemplify the interface on the example of the SEIK filter.
     34Operations that are specific to the model and to the assimilated observations are performed by call-back routines that are supplied by the user. These are called through the defined interface of PDAF. Generally, these user-supplied routines have to provide quite elementary operations, like initializing a model state vector for PDAF from model fields or providing the vector of observations. PDAF provides examples for these routines and templates that can be used as the basis for the implementation. As only the interface of these routines is specified, the user can implement the routines like a routine of the model. Thus, the implementation of these routines should not be difficult. The implementation of the call-back routines is identical for the fully-parallel and flexible implementation variants.
    14235
    143 The interface for this routine is
    144 {{{
    145 SUBROUTINE prepoststep(step, dim_p, dim_ens, dim_ens_p, dim_obs_p, &
    146                        state_p, Uinv, ens_p, flag)
     36The implementations that are required for the fully-parallel and the flexible implementation variants are described on separate pages:
     37 * '''fully parallel''': [InsertAnalysisStep Inserting the routine for the analysis step]
     38 * '''flexible''': [ExternalModelLoop Modification of the model code for the '''flexible''' implementation variant]
    14739
    148   INTEGER, INTENT(in) :: step        ! Current time step
    149                          ! (When the routine is called before the analysis -step is provided.)
    150   INTEGER, INTENT(in) :: dim_p       ! PE-local state dimension
    151   INTEGER, INTENT(in) :: dim_ens     ! Size of state ensemble
    152   INTEGER, INTENT(in) :: dim_ens_p   ! PE-local size of ensemble
    153   INTEGER, INTENT(in) :: dim_obs_p   ! PE-local dimension of observation vector
    154   REAL, INTENT(inout) :: state_p(dim_p) ! PE-local forecast/analysis state
    155                                      ! The array 'state_p' is not generally not initialized in the case of SEIK/EnKF/ETKF.
    156                                      ! It can be used freely in this routine.
    157   REAL, INTENT(inout) :: Uinv(dim_ens-1, dim_ens-1) ! Inverse of matrix U
    158   REAL, INTENT(inout) :: ens_p(dim_p, dim_ens)      ! PE-local state ensemble
    159   INTEGER, INTENT(in) :: flag        ! PDAF status flag
    160 }}}
    161 
    162 The routine `U_prepoststep` is called once at the beginning of the assimilation process. In addition, it is called during the assimilation cycles before the analysis step and after the ensemble transformation. The routine is called by all filter processes (that is `filterpe=1`).
    163 
    164 The routine provides for the user the full access to the ensemble of model states. Thus, user-controlled pre- and post-step operations can be performed.  For example the forecast and the analysis states and ensemble covariance matrix can be analyzed, e.g. by computing the estimated variances. If the smoother is used, also the smoothed ensembles can be analyzed. In addition, the estimates can be written to disk.
    165 
    166 Hint:
    167  * If a user considers to perform adjustments to the estimates (e.g. for balances), this routine is the right place for it.
    168  * Only for the SEEK filter the state vector (`state_p`) is initialized. For all other filters, the array is allocated, but it can be used freely during the execution of `U_prepoststep`.
    169  * The interface through which `U_prepoststep` is called does not include the array of smoothed ensembles. In order to access the smoother ensemble array one has to set a pointer to it using a call to the routine `PDAF_get_smootherens` (see page on [AuxiliaryRoutines auxiliary routines])
    170 
    171 
    172 === `U_collect_state` (collect_state_pdaf.F90) ===
    173 
    174 The interface for this routine is
    175 {{{
    176 SUBROUTINE collect_state(dim_p, state_p)
    177 
    178   INTEGER, INTENT(in) :: dim_p           ! State dimension for PE-local model sub-domain
    179   REAL, INTENT(inout) :: state_p(dim_p)  ! State vector for PE-local model sub-domain
    180 }}}
    181 
    182 This routine is called during the forecast phase as many times as there are states to be integrated by a model task. It is called at the end of the integration of a member state of the ensemble. The routine is executed by all processes that belong to model tasks.
    183 
    184 When the routine is called, a state vector `state_p` and its size `dim_p` are provided. The operation to be performed in this routine is inverse to that of the routine `U_distribute_state`. That is, the state vector `state_p` has to be initialized from the model fields. If the model is not parallelized, `state_p` will contain a full state vector. If the model is parallelized using domain decomposition, `state_p` will contain the part of the state vector that corresponds to the model sub-domain for the calling process.
    185 
    186 Some hints:
    187  * If the state vector does not include all model fields, it can be useful to keep a separate array to store those additional fields. This array has to be kept separate from PDAF, but can be defined using a module like `mod_assimilation`.
    188 
    189 == Simulating model errors ==
    190 
    191 The implementation of the filter algorithms does not support the specification of a model error covariance matrix. This was left out, because in the SEEK and SEIK filter, the handling can be extremely costly, as the model error covariance matrix has to be projected onto the ensemble space. Instead PDAF support the simulation of model errors by disturbing fields during the model integration. For this, some routine will be required that is inserted into the time stepping loop of the model. As this procedure is specific to each model, the is no routine provided by PDAF for this.
    192 
    193 == Compilation and testing ==
    194 
    195 To compile the extended model code with PDAF, one has to extend the Makefile for the model by adding the additional user-supplied routines. While all of the user-supplied routines need to exist not all of them need to be fully implemented at this time if the following procedure is used. The routines that will not be called are `U_init_dim_obs`, `U_obs_op`, `U_init_obs`, `U_prodRinvA`, `U_init_obsvar`. A simple way to provide them for the compilation could be to copy the corresponding files (i.e. named without `U_`) from the template directory `templates/` and to include these files in the compilation and linking. These templates are simple stubs without any functionality.
    196 
    197 At this implementation stage one can use the preprocessor definition `PDAF_NO_UPDATE` (available from Version 1.6.1). With this, the actual analysis step of the chosen filter algorithm is not executed. Accordingly, only the user-supplied routines used in `PDAF_get_state` as well as the routine `U_collect_state` need to be implemented with functionality. The other routines will not be executed, because they are only called during the analysis step. Generally with `PDAF_NO_UPDATE` the program performs just an ensemble integration. That is, PDAF is initialized by `PDAF_init`. Then a forecast is computed by using `PDAF_get_state` and the chosen `PDAF_put_state_*` routine. At the initial time `U_prepoststep` is executed by `PDAF_get_state`. `U_next_obs` will provide the number of time steps to be computed by the model and `U_distributed_state` will initialize the model fields. Subsequently the ensemble integration is performed and the forecast fields are written back to the ensemble array by `U_collect_state`. Upon completion of the forecast phase, the routine `U_prepoststep` is executed twice. The first time is the regular call before the analysis is executed. Thus, it allows to access the forecast ensemble. If the analysis would not be deactivated, the second call to `U_prepoststep` would be after the analysis allowing access to the ensemble directly after the analysis. As the analysis is deactivated here, the ensemble will be the same as in the first call.
    198 
    199 This test allows to check the following:
    200  * Is `U_prepoststep` working correctly?
    201  * Does `U_next_observation` work correctly and is the information from this routine used correctly for the model integration
    202  * Are `U_distribute_state` and `U_collect_state` work correctly?
    203 One could also comment out the actual time stepping part of the model. This would allow to only test the interfacing between PDAF and the model.
    204 
    205 It is important to ensure that the ensemble integration performs correctly. The simplest case should be a parallel configuration in which the number of model tasks equals the ensemble size as here the model tasks always compute forward in time. If the number of model tasks is smaller than the ensemble size, some model tasks will have to integrate multiple states of the ensemble. If a model task has to integrate two states, the model will have to jump back in time for the integration of the second state. It might be that some arrays of the model need to be re-initialized to ensure that the second integration is consistent. Also, one might need to check if the initialization of forcing fields (e.g. wind stress over the ocean) performs correctly for the second integration. (Sometimes model are implemented with the constraint that the model time always increases, which is the normal case for pure model simulations without assimilation.) A useful test is to initialize an ensemble in which all states are equal. If this ensemble is integrated the forecast states of the ensemble should, of course, still be equal.
    206 
    207