Inserting the analysis step into the time stepping for 'fully parallel' data assimilation
Implementation Guide
Contents of this page
Overview
For the fully-parallel variant of the assimilation system, only a single routine has to added to the time stepper of the numerical model. For each filter algorithm there is a specific routine PDAF_assimilate_X
, with X
being the name of the filter. This routine counts the time steps during the model integration and computes the analysis step of the filter when the observation time is reached. It can be convenient to not directly insert PDAF_assimilate_X
into the code but to use an interface routine assimilate_pdaf
that does not need any arguments.
In addition, the initialization routine init_pdaf
discussed before on the page on initializing PDAF, is extended by one additional subroutine call to the PDAF routine PDAF_get_state
.
An example for this implementation variant can be found in the tutorial implementations in the directory tutorial/
. Here online_2D_serialmodel
shows the implementation with a serial (i.e. non-parallel) model, while online_2D_parallelmodel
shows the implementation with a parallel model using domain-decomposition.
PDAF_get_state
The routine PDAF_get_state
has to be called once at the end of the initialization of PDAF. Usually, the call will be added to the routine init_pdaf
that was discussed on the page on initializing PDAF.
The routine PDAF_get_state
has the purpose to initialize the model fields to be propagates from the array holding the ensemble states. In addition, the routine can initialized the information, whether further model integrations have to be computed and how many time steps have to be performed. For the fully-parallel implementation variant, the number of time steps is used inside PDAF, while the flag on further model integrations is not used.
The interface of PDAF_get_state
is the following:
SUBROUTINE PDAF_get_state(nsteps, timenow, doexit, U_next_observation, U_distribute_state, & U_prepoststep, status)
with the following arguments:
nsteps
: An integer specifying upon exit the number of time steps to be performedtimenow
: A real specifying upon exit the current model time. (This value is usually not used in the fully-parallel implemenation variant)doexit
: An integer variable defining whether the assimilation process is completed. For compatibility 1 should be used for exit, 0 for continuing. (This value is not used in the fully-parallel implemenation variant)- U_next_observation: The name of a user supplied routine that initializes the variables
nsteps
,timenow
, anddoexit
- 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
- 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.
status
: The integer status flag. It is zero, ifPDAF_get_state
is exited without errors.
PDAF also has a 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 documentation of PDAF's simplified interface.
Inserting assimilate_pdaf
The right place to insert the interface routine assimilate_pdaf
into the model code, is at the end of the time stepping loop. Thus, usually this is directly before the END DO of the time stepping loop.
The purpose of assimilate_pdaf
is to call the filter-specific PDAF-core routine PDAF_assimilate_X
, with X
being the name of the filter method. It is also possible to insert the call to PDAF_assimilate_X
directly into the model code. However, using the additional interface routine yield usually cleaner source code. This is because of subroutine names that are specified in the call the PDAF_assimilate_X
or when more than one filter are implemented.
One has to ensure that PDAF_assimilate_X
is called at each time step, so that PDAF can count the time steps until the next analysis time.
PDAFomi_assimilate_X
There are a separate routines PDAFomi_assimilate_global
and PDAFomi_assimilate_local
for global and local filter algorithms, respectively. The purpose of the PDAFomi_assimilate_X
routines is to count the time steps. When the forecast phase is complete the routine writes back the forecast model fields into the array holding the ensemble of model state vectors and executes the analysis step of the chosen filter algorithm. The interface to each 'assimilate' 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.
For example, the interface when using the global filter is the following:
SUBROUTINE PDAFomi_assimilate_global(U_distribute_state, U_collect_state, U_init_dim_obs, U_obs_op, & U_prepoststep, U_next_observation, status)
At this state of the implementation only these arguments are relevant:
- 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. (This routine is also used in
PDAF_get_state
) - 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
. - 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.
- U_next_observation: The name of a user supplied routine that initializes the variables
nsteps
,timenow
, anddoexit
. (This routine is also used inPDAF_get_state
) status
: The integer status flag. It is zero, if PDAF_get_state is exited without errors.
The other arguments are names of user-supplied subroutines that are only executed if the analysis step is executed (See the section 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 (Implementation of the Analysis step) separately for each available filter algorithm.
User-supplied routines
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 Compilation and testing), all routines need to exist, but only those described here in detail need to be implemented with functionality.
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.
The user-supplied routines are in general identical for the 'fully parallel' and 'flexible' implementation variants. The only difference is in U_next_observation
as is described below.
U_next_observation
(next_observation_pdaf.F90)
The interface for this routine is
SUBROUTINE next_observation(stepnow, nsteps, doexit, timenow) INTEGER, INTENT(in) :: stepnow ! Number of the current time step INTEGER, INTENT(out) :: nsteps ! Number of time steps until next obs INTEGER, INTENT(out) :: doexit ! Whether to exit forecasting (1 for exit) REAL, INTENT(out) :: timenow ! Current model (physical) time
The routine is called once at the beginning of each forecast phase. It is executed by all processes that participate in the model integrations.
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. For the 'fully parallel' data assimilation variant the flag doexit
is not used and does not need to be set. timenow
is the current model time. However, for the 'fully parallel' data assimilation variant, this value is not relevant.
Some hints:
- 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 - At the first call to
U_next_obs
the variabletimenow
can be initialized with the current model time. At the next call a forecast phase has been completed. Thus, the new value oftimenow
follows from the timer interval for the previous forecast phase. doexit
is not relevant for the fully-parallel implementation. It is recommended to setdoexit=0
in all cases.- If
nsteps=0
ordoexit=1
is set, the ensemble state will not be distributed by PDAF (thusdistribute_state
is not called). If one intends to proceed with ensemble forecasting, one has to set nsteps to a value >0 anddoexit=0
. If nsteps is set to a value larger than the last time step of the model no further analysis step will be performed.
U_distribute_state
(distribute_state_pdaf.F90)
The interface for this routine is
SUBROUTINE distribute_state(dim_p, state_p) INTEGER, INTENT(in) :: dim_p ! State dimension for PE-local model sub-domain REAL, INTENT(inout) :: state_p(dim_p) ! State vector for PE-local model sub-domain
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.
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.
Some hints:
- 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
.
U_prepoststep
(prepoststep_ens_pdaf.F90)
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 ESTKF and LESKTF filters.
The interface for this routine is
SUBROUTINE prepoststep(step, dim_p, dim_ens, dim_ens_p, dim_obs_p, & state_p, Uinv, ens_p, flag) INTEGER, INTENT(in) :: step ! Current time step ! (When the routine is called before the analysis -step is provided.) INTEGER, INTENT(in) :: dim_p ! PE-local state dimension INTEGER, INTENT(in) :: dim_ens ! Size of state ensemble INTEGER, INTENT(in) :: dim_ens_p ! PE-local size of ensemble INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of observation vector REAL, INTENT(inout) :: state_p(dim_p) ! PE-local forecast/analysis state ! The array 'state_p' is not generally not initialized. ! It can be used freely in this routine. REAL, INTENT(inout) :: Uinv(dim_ens-1, dim_ens-1) ! Inverse of matrix U REAL, INTENT(inout) :: ens_p(dim_p, dim_ens) ! PE-local state ensemble INTEGER, INTENT(in) :: flag ! PDAF status flag
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
).
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.
Hint:
- If a user considers to perform adjustments to the estimates (e.g. for balances), this routine is the right place for it.
- 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 ofU_prepoststep
. - 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 routinePDAF_get_smootherens
(see page on auxiliary routines)
U_collect_state
(collect_state_pdaf.F90)
The interface for this routine is
SUBROUTINE collect_state(dim_p, state_p) INTEGER, INTENT(in) :: dim_p ! State dimension for PE-local model sub-domain REAL, INTENT(inout) :: state_p(dim_p) ! State vector for PE-local model sub-domain
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.
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.
Some hints:
- 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
.
Simulating model errors
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, or into assimilate_pdaf
. As this procedure is specific to each model, the is no routine provided by PDAF for this.
Compilation and testing
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.
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_assimilate_X
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.
This test allows to check the following:
- Is
U_prepoststep
working correctly? - Does
U_next_observation
work correctly and is the information from this routine used correctly for the model integration - Do
U_distribute_state
andU_collect_state
work correctly?
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.