Adding a Filter Algorithm to PDAF

This text describes the implementation strategy and internal structure of PDAF valid for version 1.7.0 and later. If you use an earlier version of PDAF, we recommend to update to the most recent version. In this text, we assume that the reader is already familiar with PDAF to the extend of experience with the implementation of a model with PDAF that is described in the implementation guide.

The internal structure of PDAF is organized into a generic part providing the infrastructure to perform ensemble forecasts and filter analysis steps. This generic part is independent of the particular filter algorithm and only distinguishes between ensemble based filters (all filters except SEEK) and mode based filters (currently only SEEK). The filter-specific routines are called through an internal interface.

Each filter algorithm consists of 6 routines that are described below. All routines are called trough the internal interface of PDAF, except for the "put state" routine (PDAF_put_state_X where X is the name of the selected filter).

PDAF's Internal Interface

Before explaining the filter-specific routines and the calling interface of each routine, we provide an overview of the internal interface routines of PDAF. The internal interface of PDAF is structured as follows ('X' is the name of the filter algorithm):

The routine PDAF_init calls

• PDAF_init_filters - interface routine to PDAF_X_init
  
  • PDAF_X_init performs the filter-specific initialization of parameters and calls the user-supplied routine that initializes the initial ensemble of model states.
• PDAF_alloc_filters - interface routine to PDAF_X_alloc
  
  • PDAF_X_alloc allocates the filter-specific arrays.
• PDAF_options_filters - interface routine to PDAF_X_options
  • PDAF_X_options is an optional routine. Its purpose is to display an overview of available options for the filter algorithm.

The routine PDAF_get_state calls

• PDAF_initinfo_filters - interface routine to PDAF_X_initinfo
  • PDAF_X_initinfo displays information on the current configuration of the filter algorithm.

The routine PDAF_print_info only includes the interface to PDAF_X_memtime

• PDAF_X_memtime displays information on the execution duration of the different parts of the assimilation process as well as information on the amount of allocated memory. This functionality is optional.

The routine PDAF_put_state_X is called directly from the model code. There is a separate routine for each filter, mainly because of the fact that different user-supplied routines may be needed for the analysis step of the filter. However, the operations performed directly in PDAF_put_state_X are widely generic and the filter-specific analysis step is typically implemented as another subroutine. The standard implementation calls

• PDAF_X_initinfo
  • This is the same routine that is called by PDAF_get_state. Here, it is only called, if the offline mode of PDAF is used.
• PDAF_X_update
  • This routine contains the actual assimilation or analysis step of the filter algorithm.

When PDAF_init is called, the filter algorithm is chosen by its ID number. Internally to PDAF, each filter is identified by a string that is defined in PDAF_init_filters. The interface routines have a very simple structure. In general, they select the filter-specific routine based on the string identifying the filters. When a filter algorithm is added a line for the corresponding filter-specific routine has to be inserted to each of the interface routines.

Filter-specific routines

When a filter algorithm is added, the following filter routines have to be implemented and inserted to each interface routines described above.

• PDAF_X_init
• PDAF_X_alloc
• PDAF_X_options
• PDAF_X_initinfo
• PDAF_X_memtime

In addition, the routine

• PDAF_put_state_X

has to be implemented that is called directly in the model code.

We recommend to base on the routines of an existing filter, as most of the routines can be easily adapted to a new filter method.

PDAF_X_init

The routine PDAF_X_init performs the initialization of filter-specific parameters. In addition, it prints information about the configuration.

The interface is as follows:

  SUBROUTINE PDAF_X_init(subtype, param_int, dim_pint, param_real, dim_preal, &
                                 ensemblefilter, fixedbasis, verbose, outflag)

with the following arguments:

• subtype: The subtype index of the filter algorithm [integer, input].
• param_int: The array of integer parameters [integer(dim_pint), input].
• dim_pint: The number of parameters in param_int [integer, input].
• param_real: The array of real parameters [real(dim_preal), input].
• dim_preal: The number of parameters in param_real [integer, input]
• ensemblefilter: Flag, whether the filter is an ensemble filter or a mode-based filter [logical, output].
• fixedbasis: Flag, whether the chosen subtype is a filter with fixed ensemble, such that only the ensemble mean is integrated by the model [logical, output].
• verbose: Verbosity flag [integer, input]. Valid are the values provided to PDAF_init.
• outflag: Error flag [integer, output]

The required operations are to initialize the PDAF-internal parameter variables from the provided values of subtype, param_int, and param_real. In the addition, the logical flags ensemblefilter and fixedbasis have to be set. The existing implementations also include some screen output about the configuration.

Please note:

• The routine should check, whether the provided value of subtype is a valid choice. If this is not the case, the error flag should be set to 2.
• Only parameters from param_int and param_real up to the value dim_pint and dim_preal should be considered in the initialization. The array may be bigger, but the user defined which parameters are to be used be setting the values of dim_pint and dim_preal.
• The error flag outflag is initial set to 0.
• The internal parameters are declared and stored in the Fortran module PDAF_mod_filter. If a new filter algorithm requires additional parameters, their declaration should be added to the module.

PDAF_X_alloc

The routine PDAF_X_alloc allocates arrays for the data assimilation, like the ensemble array and a state vector. The success of the allocation is checked.

The interface is as follows:

  SUBROUTINE PDAF_X_alloc(subtype, outflag)

with the following arguments:

• subtype: The subtype index of the filter algorithm [integer, input].
• outflag: Error flag [integer, input/output]

All arrays that need to be allocated are declared in the Fortran module PDAF_mod_filter. Here, also the dimensions of the arrays are declared. For the allocation of arrays, one has to distinguish between processes that compute the analysis step and those that only participate in the ensemble forecast. For the processes that compute the analysis (those with filterpe=.true.) it is mandatory to allocate the following two arrays:

• state: The state vector of size dim_p.
• eofV: This is the ensemble matrix in all ensemble-based filters. For SEEK it is the matrix holding eigenvectors. eofV has size (dim_p, dim_ens).

Depending on the filter algorithm some of the following arrays also need to be allocated:

• eofU: This is the eigenvalue matrix U used in the SEIK and SEEK filters (here, its size is (rank,rank)). For ETKF, it is the matrix A of size (dim_ens,dim_ens). The array only needs to be allocated if the algorithm uses such a matrix. (For EnKF, which does not use this matrix, it is allocated with size (1,1).)
• state_inc: The increment to the state vector computed in the analysis step. It only needs to be allocated in this routine, if incremental analysis updating is implemented. Otherwise, it is sufficient to allocate and deallocate state_inc in the routine performing the analysis step. The size of state_inc is dim_p.
• bias: If the filter algorithm is implemented with bias correction, the vector bias with size dim_bias_p is allocated.

Processes that only participate in the computation of the ensemble forecast, but are not involved in computing the analysis step, operate only on a sub-ensemble. Accordingly, an ensemble array for this sub-ensemble has to be allocated. This is:

• eofV: This is the ensemble matrix in all ensemble-based filters. For SEEK it is the matrix holding eigenvectors. For the processes with filterpe=.false., eofV has size (dim_p, dim_ens_l).