Version 28 (modified by 9 years ago) (diff)  ,

Adding a Filter Algorithm to PDAF
Contents of this page
This text describes the implementation strategy and internal structure of PDAF valid for version 1.8.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 filterspecific routines are called through an internal interface.
Each filter algorithm consists of 3 mandatory routines plus 2 optional routines. All routines are described below. They are called through 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 filterspecific 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 toPDAF_X_init
PDAF_X_init
performs the filterspecific initialization of parameters and calls the usersupplied routine that initializes the initial ensemble of model states.
PDAF_alloc_filters
 interface routine toPDAF_X_alloc
PDAF_X_alloc
allocates the filterspecific arrays.
PDAF_options_filters
 interface routine toPDAF_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_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 usersupplied 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 filterspecific analysis step is typically implemented as another subroutine.
The standard implementation calls
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 filterspecific routine based on the string identifying the filters. When a filter algorithm is added a line for the corresponding filterspecific routine has to be inserted to each of the interface routines.
Internal dimensions
PDAF internally stores the dimensions of the assimilation system. The dimensions are declared in the Fortran module PDAF_mod_filter
. Important are the following dimensions:
dim_p  The size of the state vector (with parallelization the size of the local state vector for the current process) 
dim_ens  The overall size of the ensemble 
dim_ens_l  If the ensemble integration is distributed over several ensemble tasks, this variable stores the size of the subensemble handled by the current process. (dim_ens_l equals dim_ens if no parallelization or if only a single model task is used.)

rank  The maximum rank of the ensemble covariance matrix. In almost all cases, it is dim_ens1 .

dim_eof  For mode based filters (currently only SEEK), this is the number of modes used in the state covariance matrix. 
Filterspecific 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
(optional)PDAF_X_memtime
(optional)
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 filterspecific 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 inparam_int
[integer, input].param_real
: The array of real parameters [real(dim_preal), input].dim_preal
: The number of parameters inparam_real
[integer, input]ensemblefilter
: Flag, whether the filter is an ensemble filter or a modebased filter [logical, output].fixedbasis
: Flag, whether the chosensubtype
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 toPDAF_init
.outflag
: Error flag [integer, output]
The required operations are to initialize the PDAFinternal 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
andparam_real
up to the valuedim_pint
anddim_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 ofdim_pint
anddim_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 sizedim_p
.eofV
: This is the ensemble matrix in all ensemblebased 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 deallocatestate_inc
in the routine performing the analysis step. The size ofstate_inc
isdim_p
.bias
: If the filter algorithm is implemented with bias correction, the vectorbias
with sizedim_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 subensemble. Accordingly, an ensemble array for this subensemble has to be allocated. This is:
eofV
: This is the ensemble matrix in all ensemblebased filters. For SEEK it is the matrix holding eigenvectors. For the processes withfilterpe=.false.
,eofV
has size (dim_p
,dim_ens_l
).
PDAF_X_options
The optional routine PDAF_X_options
displays information on the available options for the filter algorithm.
The interface is as follows:
SUBROUTINE PDAF_X_options()
It has no arguments!
The following display is recommended:
 Available subtypes (At least '0' for standard implementation; '5' for offline mode)
 Parameters used from the parameter arrays
param_int
andparam_real
.
PDAF_X_memtime
The optional routine PDAF_X_memtime
displays information about allocated memory and the execution time of different parts of the filter algorithm.
The interface is as follows:
SUBROUTINE PDAF_X_memtime(printtype)
with the following argument:
printtype
: The type of the output to be done [integer, input]. For the filter algorithms that are included in the PDAF source code package the following choices are implemented: 1: Display general timers
 2: Display allocated memory
 3: Display detailed timers
The timing operation are implement using the module PDAF_timer
, which provides the function PDAF_timeit
. Memory allocation is computed using PDAF_memcount
, which is provided by the module PDAF_memcounting
.
PDAF_put_state_X
This routine is directly inserted into the model code, if the online mode of PDAF is used. The text on the implementation of the analysis step in the Implementation Guide explains the interface for the algorithms that are included in the PDAF package. Apart from the usual integer status flag, the interface contains the names of the usersupplied routines that are required for the analysis step. Usually, the minimum set of routines are:
 A routine to write model fields back into the ensemble state array (
U_collect_state
).  A routine that determines the size of the observation vector (
U_init_dim_obs
).  A routine that contains the implementation of the observation operator (
U_obs_op
).  A routine that fills the observation vector (
U_init_obs
).  The pre and poststep routine in which the forecast and analysis ensembles can be analyzed (
U_prepoststep
).
Additional routines are possible depending on the requirements of the filter algorithm. When one implements a new filter, one should check, which of the existing routines can be reused. This will facilitate the switching between different filter methods. For example, the SEIK, ETKF, and ESTKF filters (as well as the LSEIK, LETKF, and LESTKF algorithms) use the identical set of usersupplied routines. Thus, switching between these filters is possible without additional implementation work.
With regard to the operations, PDAF_put_state_X
prepares for the actual analysis step, that is called inside PDAF_put_state_X
as a separate routine. The required operations are:
 Write model fields back into the ensemble array (by calling
U_collect_state
)  Increment the counter for the integrated ensemble members (named
counter
and provided by the modulePDAF_mod_filter
.  Check, if the ensemble integration is completed (in that case, it is
member = local_dim_ens + 1
). If not, exit thePDAF_put_state_X
.  If the ensemble integration is completed, the following operations are required:
 If more than one model task is used: Collect the subensembles from all model tasks onto the processes that perform the analysis step. This operation is done by the subroutine
PDAF_gather_ens
.  Call the routine that computes the analysis step for the chosen filter algorithm
 Reset the control variables for the ensemble forecast (
initevol=1
,member=1
,step=step_obs+1
).
 If more than one model task is used: Collect the subensembles from all model tasks onto the processes that perform the analysis step. This operation is done by the subroutine