Changes between Version 7 and Version 8 of ExternalModelLoop


Ignore:
Timestamp:
May 24, 2025, 11:29:55 AM (8 days ago)
Author:
lnerger
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • ExternalModelLoop

    v7 v8  
    2020== Overview ==
    2121
    22 For the flexible implementation variant of the assimilation system, one has to modify the model code so that it is possible to compute several integrations of model states successively. This is possible by adding an additional loop outside of the regular time-stepping loop of the model. This strategy has the potential to the required chances in the model code to the minimum, while allowing for the flexibility. 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.
     22For the ''flexible'' implementation variant of the assimilation system, one has to modify the model code so that each model task can compute several integrations of model states successively. The required extensions are described below.
     23
     24Because the flexible implementation variant is more complex to implement it can be easier if one first implements the ''fully parallel'' variant and hen modifies this to obtain the code structure for the ''flexible'' variant.
    2325
    2426
    2527== External ensemble loop ==
    2628
    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 `PDAFomi_put_state_global` for global filters. Both routines are described in sections below.
     29Figure 1 compares the code structure of the ''fully parallel'' and ''flexible'' implementation variants. While for the ''fully parallel' variant only the routine `assimilate_pdaf` has to be added in the model time stepping loop, the additional ensemble loop, a call to `PDAF_get_fcst_info`, and a check of the exit flag have to be added.
    2830
    29 The extended model code can look like this for the SEIK filter:
    30  {{{
     31[[Image(//pics/DAextension_flexible_PDAF3.png)]]
     32[[BR]]**Figure 1:** Comparison of the code structure for the ''fully parallel'' and ''flexible'' implementation variants. The ''flexible'' variant required more changes to the code to enable that each model task can integrate more than one ensemble state. Note that the ''fully parallel'' variant defines the number of time steps as `isteps`, which is the total number of steps in the program, while the ''flexible'' variant uses `nsteps`, which is the number of steps in one forecast phase.
     33
     34In the source code, the modification can be realized in the following way:
     35{{{
    3136  pdaf_modelloop: DO 
    3237
    33      CALL PDAF_get_state(nsteps, ..., doexit, ...)
     38     ! PDAF: Get forecast information
     39     CALL PDAF_get_fcst_info(nsteps, timenow, doexit)
    3440
    35      ! Check whether forecast has to be performed
    36      checkforeast: IF (doexit /= 1 .AND. status_pdaf == 0) THEN
    37      
    38         IF (nsteps > 0) THEN
     41     ! Check exit flag
     42     IF (doexit==1) EXIT pdaf_modelloop
    3943
    40           ... Time stepping code of the model ...         
     44     ! Model time stepping loop, for example
     45     DO i = 1, nsteps
    4146
    42         END IF
     47        ... Time stepping code of the model ...
    4348
    44         CALL PDAFomi_put_state_global(...)
    45 
    46      ELSE checkforecast
    47 
    48         ! No more assimilation work; exit loop
    49         EXIT pdaf_modelloop
    50 
    51      END IF checkforecast
     49        ! Let PDAF check forecast progress and perform analysis
     50        CALL assimilate_pdaf()
     51     END DO 
    5252
    5353  END DO pdaf_modelloop
    5454}}}
    55 In this example, which is taken from the example implementation in `templates/online_omi_flexible/`, 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).
     55This example is taken from the example implementation in `templates/online_flexible/`. The unconditional DO loop (while loop) allows for an arbirary number of repetitions. The exit of the loop is controlled by the exit flag `doexit`, which is obtained from `PDAF_get_fcst_info`. The variable `nsteps`, also obtained from PDAF_get_fcst_info, defines the number of time steps to be computed during the current forecast phase. This value has to be used in the model time stepping loop.
    5656
    57 == `PDAF_get_state` ==
     57The code changes for the ''flexible'' parallelization variant can usually be obtained by adapting the code for the ''fully parallel'' implementation variant.
    5858
    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.
     59Apart from the changes described above we need to initialize the ensemble forecasting as in the ''fully parallel'' variant. In particular , we need to write the ensemble state vector values into the model fields. We need to set for PDAF how many time steps need to be done until the first observations are assimilated.  For these operations, we call the routine `PDAF_init_forecast` in the routine `init_pdaf`. This call is identical to that used in the ''fully parallel' variant.
    7460
    7561
    76 == `PDAFomi_put_state_X` ==
     62== Calling `PDAF_init_forecast` ==
    7763
    78 There are separate routines `PDAF_put_state_global` and `PDAF_put_state_local`, respeectively for global and local filter algorithms. 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.
     64The call to `PDAF_init_forecast` is identical for the ''fully parallel'' and ''flexible'' implementation variants.
    7965
    80 For example, the interface when using global filters is the following:
     66The routine is called at the end of the routine `init_pdaf` that was discussed on the [InitPdaf page on initializing PDAF]. The main purpose of this routine is to initialize the model fields for all model tasks from the ensemlbe of state vectors. This is done in the call-back routine `distribute_state_pdaf`. The routine also calls the call-back routine `next_observation_pdaf` to set the number of time steps for the initial forecast phase and an exit flag. These values used are internally by PDAF to control the forecast phase. Further, the routine calls the call-back-routine `prepoststep_pdaf`. This pre/postep routine provides the user access to the initial ensemble.
     67
     68The interface of `PDAF_init_forecast` is:
    8169{{{
    82   SUBROUTINE PDAFomi_put_state_global(U_collect_state, U_init_dim_obs, U_obs_op, &
    83                                  U_prepoststep, status)
     70  SUBROUTINE PDAF_init_forecast(next_observation_pdaf, distribute_state_pdaf, &
     71                                prepoststep_pdaf, status)
    8472}}}
    85 At this state of the implementation only these arguments are relevant:
    86  * [#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`
    87  * `status`: The integer status flag. It is zero, if PDAF_get_state is exited without errors.
     73with the arguments (where the bold names show the arguments relevant to the user for the ''fully-parallel'' variant):
     74 * [#next_observation_pdafnext_observation_pdaf.F90 next_observation_pdaf]:[[BR]] The name of a user supplied routine that initializes the variables `nsteps`, `timenow`, and `doexit`
     75 * [#distribute_state_pdafdistribute_state_pdaf.F90 distribute_state_pdaf]:[[BR]] The name of a user supplied routine that initializes the model fields from the array holding the ensemble of model state vectors
     76 * [#prepoststep_pdafprepoststep_ens_pdaf.F90 prepoststep_pdaf]:[[BR]] 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.
     77 * `status`, `integer, intent(inout)`:[[BR]] The status flag. It is zero, if the routine is exited without errors. We recommend to check the value.
    8878
    89 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.
     79The user-supplied routines are described further below.
     80
     81
     82== `PDAF_get_fcst_info` ==
     83
     84The interface is the following:
     85{{{
     86  SUBROUTINE PDAF_get_fcst_info(nsteps, timenow, doexit)
     87
     88      INTEGER, INTENT(out) :: stepnow  ! The current time step
     89      INTEGER, INTENT(out) :: nsteps   ! The number of time steps of the next forecast phase
     90      INTEGER, INTENT(out) :: doexit   ! Exit flag: (1) exit, (0) continue data assimilation
     91}}}
     92
     93These variables are intialized by `next_observation_pdaf` that was called by `PDAF_init_forecast`. At later forecast phases in the assimilation processing, the analysis routine `PDAF3_assimilation` (or more specifalize routines) will case `next_observation_pdaf`.
     94
     95== Inserting `assimilate_pdaf` into the model code ==
     96
     97The right place to insert the routine `assimilate_pdaf` into the model code is the same as for the ''fully parallel'' implementation variant. Usually this is at the end of the model time stepping loop, thus when the model completed the computation of a single time step. In most cases, this is just before the 'END DO' in the model source code. However, there might be special cases, where the model does some additional operations so that the call to `assimilate_pdaf` should be insered somewhat earlier.
     98
     99== Using `assimilate_pdaf` ==
     100
     101The purpose of `assimilate_pdaf` is to call the universal PDAF-core routine `PDAF3_assimilate` (or a more specific varant of this routine).  The arguments of `PDAF3_assimilate` are mainly the names of user-supplied call-back routines, except from an argument for the status flag. These names are specified in `assimilate_pdaf` as 'external'.
     102
     103The routine `assimilate_pdaf` is called at each time step of the model forecast. This allows to, e.g., apply [wiki:IncrementalAnalysisUpdates incremental analysis updates].
     104
     105Details on the implementations of the user-routines for the analysis step ar provided in the [wiki:ImplementationofAnalysisStep Page on implementating the analysis step].
    90106
    91107
    92108== User-supplied routines ==
    93109
    94 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.
     110Here, we discuss the user-supplied routines that are arguments of `PDAF_init_forecast`.
    95111
    96 To indicate user-supplied routines we use the prefix `U_`. In the tutorials in `tutorial/` and in the template directory `templates/` these routines exist without the prefix, but with the extension `_pdaf`. The files are named correspondingly. In the section titles below we provide the name of the template file in parentheses.
     112In the section titles below we provide the name of the template file in parentheses.
    97113
    98 === `U_next_observation` (next_observation_pdaf.F90) ===
    99114
    100 The interface for this routine is
     115
     116=== `next_observation_pdaf` (next_observation_pdaf.F90) ===
     117
     118The interface is:
    101119{{{
    102 SUBROUTINE next_observation(stepnow, nsteps, doexit, timenow)
     120SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, timenow)
    103121
    104   INTEGER, INTENT(in)  :: stepnow  ! Number of the current time step
    105   INTEGER, INTENT(out) :: nsteps   ! Number of time steps until next obs
    106   INTEGER, INTENT(out) :: doexit   ! Whether to exit forecasting (1 for exit)
    107   REAL, INTENT(out)    :: timenow  ! Current model (physical) time
     122  INTEGER, INTENT(in)  :: stepnow  ! The current time step provided by PDAF
     123  INTEGER, INTENT(out) :: nsteps   ! The number of time steps of the next forecast phase
     124  INTEGER, INTENT(out) :: doexit   ! Exit flag: (1) exit, (0) continue data assimilation
     125  REAL, INTENT(out)    :: timenow  ! Current (physical) model (physical) time. Not used by PDAF itself.
    108126}}}
    109127
    110 The routine is called once at the beginning of each forecast phase. It is executed by all processes that participate in the model integrations.
     128The routine is called by `PDAF_init_forecast` and later at the beginning of each forecast phase by `PDAF3_assimilate`, or similar routines. It is executed by all processes that participate in the model integrations.
    111129
    112 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.
     130Based on the information of the current time step, the routine has to define the number of time steps `nsteps` for the next forecast phase and the exit flag `doexit`. In the user code, one accesses these variables with `PDAF_get_fcst_info` as explained above.
    113131
    114132Some hints:
    115  * 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
    116  * `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.
    117  * 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.
    118  * If `nsteps=0` or `doexit=1` is set, the ensemble state will not be distributed by PDAF (thus `distribute_state` is not called). If one intends to proceed with ensemble forecasting, one has to set nsteps to a value >0 and `doexit=0`.
     133 * We assume that the time interval between successive observations is known. Then, `nsteps` can be simply initialized by dividing the time interval by the size of the time step.
     134 * It is up to the user to define `timenow` because it is only used in the user-code. The user can either ignore it (setting to to 0.0), or could use it freely to indicate the model time. At the first call to `next_observation_pdaf` the variable `timenow` can 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 and can be incremented accordingly. In the user code, one can access `timenow` with a call to [wiki:PDAF_get_fcst_info].
     135 * Inside PDAF, `doexit` is used as follows: If `nsteps=0` or `doexit=1` is set, the ensemble state will not be distributed by PDAF (thus `distribute_state`, see below, is not called). If one intends to proceed with ensemble forecasting, one has to set `nsteps` to a value >0 and `doexit=0`. 
     136 * The total number of time steps set for the assimilation experiment (variable `isteps` in Fig. 1) defines when the model integrations end. If `nsteps` is set so that it specifies a time step larger than `isteps` not all of these steps will be computed.
    119137
    120 === `U_distribute_state` (distribute_state_pdaf.F90) ===
     138=== `distribute_state_pdaf` (distribute_state_pdaf.F90) ===
    121139
    122140The interface for this routine is
    … …  
    124142SUBROUTINE distribute_state(dim_p, state_p)
    125143
    126   INTEGER, INTENT(in) :: dim_p           ! State dimension for PE-local model sub-domain
    127   REAL, INTENT(inout) :: state_p(dim_p)  ! State vector for PE-local model sub-domain
     144  INTEGER, INTENT(in) :: dim_p           ! Size of state vector for process-local sub-domain
     145  REAL, INTENT(inout) :: state_p(dim_p)  ! State vector for process-local sub-domain
    128146}}}
    129147
    130 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.
     148The routine is called by `PDAF_init_forecast` and later at the beginning of each forecast phase by `PDAF3_assimilate`, or similar routines. It is executed by all processes that participate in the model integrations.
    131149
    132 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.
     150PDAF calls this routine providing the state vector `state_p` and its size `dim_p`. The routine has to write the information from the state vector into the model field arrays.
    133151
    134152Some hints:
    135  * 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`.
     153 * The user has defined the setup of the state vector, so that it is know which model fields are stored at which position in the state vector and in which order, for example, a 3-dimensional model field is stored. This was already used in `init_ens_pdaf`.
     154 * The code usually consist of (nested) loops through each of the model fields, one field at a time. In the loop each element of the state vector is written into an element of a model field array.
     155 * 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.
    136156
    137157
    138 === `U_prepoststep` (prepoststep_ens_pdaf.F90) ===
    139 
    140 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.
     158=== `prepoststep_pdaf` (prepoststep_ens_pdaf.F90) ===
    141159
    142160The interface for this routine is
    143161{{{
    144162SUBROUTINE prepoststep(step, dim_p, dim_ens, dim_ens_p, dim_obs_p, &
    145                        state_p, Uinv, ens_p, flag)
     163                       state_p, Ainv, ens_p, flag)
    146164
    147165  INTEGER, INTENT(in) :: step        ! Current time step
    148166                         ! (When the routine is called before the analysis -step is provided.)
    149   INTEGER, INTENT(in) :: dim_p       ! PE-local state dimension
     167  INTEGER, INTENT(in) :: dim_p       ! Process-local state dimension
    150168  INTEGER, INTENT(in) :: dim_ens     ! Size of state ensemble
    151   INTEGER, INTENT(in) :: dim_ens_p   ! PE-local size of ensemble
    152   INTEGER, INTENT(in) :: dim_obs_p   ! PE-local dimension of observation vector
    153   REAL, INTENT(inout) :: state_p(dim_p) ! PE-local forecast/analysis state
    154                                      ! The array 'state_p' is not generally not initialized in the case of SEIK/EnKF/ETKF.
     169  INTEGER, INTENT(in) :: dim_ens_p   ! Process-local ensemble size (not relevant for fully parallel)
     170  INTEGER, INTENT(in) :: dim_obs_p   ! Process-local dimension of observation vector
     171  REAL, INTENT(inout) :: state_p(dim_p) ! Process-local state vector
     172                                     ! The array 'state_p' is not generally not initialized.
    155173                                     ! It can be used freely in this routine.
    156   REAL, INTENT(inout) :: Uinv(dim_ens-1, dim_ens-1) ! Inverse of matrix U
    157   REAL, INTENT(inout) :: ens_p(dim_p, dim_ens)      ! PE-local state ensemble
     174  REAL, INTENT(inout) :: Ainv(dim_ens-1, dim_ens-1) ! Inverse of matrix A )only for special cases)
     175  REAL, INTENT(inout) :: ens_p(dim_p, dim_ens)      ! State ensemble for process-local sub-domain
    158176  INTEGER, INTENT(in) :: flag        ! PDAF status flag
    159177}}}
    160178
    161 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`).
     179The routine `prepoststep_pdaf` 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=.true.`).
    162180
    163 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.
     181With this routine, PDAF provides 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.
    164182
    165 Hint:
    166  * If a user considers to perform adjustments to the estimates (e.g. for balances), this routine is the right place for it.
    167  * 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`.
    168  * 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])
     183The routine also allows to perform adjustments to the state vectors, e.g. for balances, in the ensemble.
    169184
     185Hints:
     186 * The state vector, `state_p`, is usually not initialized because the ensemble DA methods focus on computing the ensemble and not its mean state. Thus the user would need to compute the ensemble mean, if this is required. This array can be used freely in `prepoststep_pdaf`.
     187 * The interface of `prepoststep_pdaf` 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]).
    170188
    171 === `U_collect_state` (collect_state_pdaf.F90) ===
    172 
    173 The interface for this routine is
    174 {{{
    175 SUBROUTINE collect_state(dim_p, state_p)
    176 
    177   INTEGER, INTENT(in) :: dim_p           ! State dimension for PE-local model sub-domain
    178   REAL, INTENT(inout) :: state_p(dim_p)  ! State vector for PE-local model sub-domain
    179 }}}
    180 
    181 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.
    182 
    183 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.
    184 
    185 Some hints:
    186  * 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`.
    187 
    188 == Simulating model errors ==
    189 
    190 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.
    191189
    192190== Compilation and testing ==
    193191
    194 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.
     192One can test the program without hacing implemented the actual assimilation step. In the template code `template/online` all required user-supplied routines are included, but they don't contain functionality. However, one can use them to be able to compile and run the assimilation program for testing. In particular one can check if the ensemble forecasting works correctly.
    195193
    196 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.
    197 
    198 This test allows to check the following:
    199  * Is `U_prepoststep` working correctly?
    200  * Does `U_next_observation` work correctly and is the information from this routine used correctly for the model integration
    201  * Are `U_distribute_state` and `U_collect_state` work correctly?
    202 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.
    203 
    204 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.
     194However, one can only check this in the model code at this point because the forecasted model fields are not yet written back into the state vectors.