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. The required extensions are described below.

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: A user supplied routine that initializes the variables nsteps, timenow, and doexit
• U_distributed_state: A user supplied routine that initializes the model fields from the array holding the ensemble of moel state vectors
• U_prepoststep: 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: An integer status flag