Modification of the model code for the ensemble integration

Overview

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.

Some operations that are specific to the model and the observations that are assimilated are performed by routines that are supplied by the user and that 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.

External ensemble loop

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.

The extended model code can look like this for the SEIK filter:

 pdaf_modelloop: DO  

    CALL PDAF_get_state(nsteps, ..., doexit, ...)

    ! Check whether forecast has to be performed
    ifcontrol: IF (doexit /= 1) THEN
     
       IF (nsteps > 0) THEN

         ... Time stepping code of the model ...         

       END IF

       CALL PDAF_put_state_seik(...)

    ELSE ifcontrol

       ! No more assimilation work; exit loop
       EXIT pdaf_modelloop

    END IF ifcontrol

 END DO pdaf_modelloop

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).

PDAF_get_state

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.

The interface of PDAF_get_state is the following:

  SUBROUTINE PDAF_get_state(nsteps, timenow, doexit, U_next_obs, U_distribute_state, &
                            U_prepoststep, status)

with the following arguments:

• nsteps: An integer specifying upon exit the number of time steps to be performed
• timenow: A real specifying upon exit the current model time.
• 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.
• U_next_obs: The name of a user supplied routine that initializes the variables nsteps, timenow, and doexit
• U_distributed_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, if PDAF_get_state is existed without errors.

PDAF_put_state_*

There is a separate routine PDAF_put_state_* for each of the filter algorithms. The name of the routine includes the name of the filter at its end. The purpose of the PDAF_put_state_* 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. For example, the interface when using the SEIK filter is the following:

  SUBROUTINE PDAF_put_state_seik(U_collect_state, U_init_dim_obs, U_obs_op, &
                                 U_init_obs, U_prepoststep, U_prodRinvA, U_init_obsvar, status)

with the following arguments:

• U_collect_state: The name of the user-supplied routine that initializes a state vector from the array holding the ensembel of model states from the model fields. This is basically the inverse operation to U_distribute_state used in PDAF_get_state
• U_init_dim_obs: The name of the user-supplied routine that provides the size of observation vector
• U_obs_op: The name of the user-supplied routine that acts as the observation operator on some state vector
• U_init_obs: The name of the user-supplied routine that initializes the vector of observations
• U_prepoststep: The name of the pre/poststep routine as in PDAF_get_state
• U_prodRinvA: 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. This operation occurs during the analysis step of the SEIK filter.
• U_init_obsvar: The name of the user-supplied routine that provides a mean observation error variance to PDAF (This routine will only be executed, if an adaptive forgetting factor is used)
• status: The integer status flag. It is zero, if PDAF_get_state is exited without errors.

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. 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.

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.

This test allows to check the following:

• Is U_prepoststep working correctly?
• Does U_next_obs work correctly and is the information from this routine used correctly for the model integration
• Are U_distribute_state and U_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.

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.