Initialization of PDAF and the ensemble by PDAF_init

Overview

After the adaption of the parallelization, the initialization of PDAF has to be implemented. PDAF_internally, the initialization is performed by the routine PDAF_init. Typically, the initializations of all variables required for the call to PDAF_init can be collected into a single subroutine. This implementation strategy is useful, as only a single addition subroutine call has to be added to the model source code. In the example in testsuite/src/dummymodel_1D the routine in the file init_pdaf_dummy_D.F90 shows this strategy.

PDAF_init itself calls a user-supplied routine to initialize the ensemble of model states through its interface. In the example there are separate routines for the different filters. They are in the files init_seik_dummy_D.F90 (for SEIK, LSEIK, ETKF, LETKF), init_enkf_dummy_D.F90 (for EnKF), and init_seek_dummy_D.F90 (for SEEK).

Inserting init_pdaf

The right place to insert a routine like init_pdaf to the model code is in between the initialization part of the model and its time stepping loop. In the routine a number of variables have to be defined that are used in the call to PDAF_init as described in 'Required arguments for `PDAF_init`'. (Please note: All names of subroutines that start with PDAF_ are core routines of PDAF, while subroutines whose name end with _pdaf are generally user-supplied interface routines) There are also a few variables that are initialized in init_pdaf but not used in the call to PDAF_init. These are variables that are specific for the data assimilation system, but only shared in between the user-supplied routines. For the example case, these variables are described in 'Other variables for the assimilation'

The example implementation allows to parse all variables through a command line parser. This method provides a convenient way to define an experiment and could also be used for other models.

Required arguments for PDAF_init

The call to PDAF_init has the following structure:

CALL PDAF_init(filtertype, subtype, step_null, &
               filter_param_i, length_filter_param_i, &
               filter_param_r, length_filter_param_r, &
               COMM_model, COMM_filter, COMM_couple, &
               task_id, n_modeltasks, filterpe, &
               U_ensemble_init, screen, status_pdaf)

The required variables are the following:

• filtertype: An integer defining the type of filter algorithm. Available are
  • 0: SEEK
  • 1: SEIK
  • 2: EnKF
  • 3: LSEIK
  • 4: ETKF
  • 5: LETKF
• subtype: An integer defining the sub-type of the filter algorithm (see example code for choices)
• step_null: An integer defining the initial time step. For some cases it can use useful to set step_null larger to 0.
• filter_param_i: Integer array collecting several variables for PDAF. The first two variables are mandatory and equal for all filters. Further variables are optional (see example code). The mandatory variables are in the following order:
  • The size of the local state vector for the current process.
  • The ensemble size for all ensemble-based filters (or the rank of the state covariance matrix for mode-based filters like SEEK)
• length_filter_param_i: An Integer defining the length of the array filter_param_i. The entries in the array are parsed up to this index.
• filter_param_r: Array of reals collecting flowting point variables for PDAF. The first variable is mandatory and equal for all filters. Further variables are optional (see example code). The mandatory variable is:
  • The value of the forgetting factor (required to be larger than zero)
• length_filter_param_r: An Integer defining the length of the array filter_param_r. The entries in the array are parsed up to this index.
• COMM_model: The communicator variable COMM_model as initialized by init_parallel_pdaf. If the model-communicator is named differently in the actual program, the name has to be adapted
• COMM_filter: The communicator variable COMM_filter as initialized by init_parallel_pdaf.
• COMM_couple: The communicator variable COMM_couple as initialized by init_parallel_pdaf.
• task_id: The index of the model tasks as initialized by init_parallel_pdaf.
• n_modeltasks: The number of model tasks as defined before the call to init_parallel_pdaf.
• filterpe: The flag showing if a process belongs to COMM_filter as initialized by init_parallel_pdaf.
• U_ensemble_init: The name of the user-supplied routine that is called by PDAF_init to initialize the ensemble of model states. (See 'User-supplied routine `U_ensemble_init`'
• screen: An integer defining whether information output is written to the screen (i.e. standard output). The following choices are available:
  • 0: quite mode - no information is displayed.
  • 1: Display standard information (recommended)
  • 2: Display detailed information for debugging
• status_pdaf: An integer used as status flag of PDAF. If status_pdaf is zero upon exit from PDAF_init the initialization was successful. An error occurred for non-zero values. (The error codes are documented in the routine PDAF_init.)

Other variables for the assimilation

The routine init_pdaf in the example also initializes several variables that are not used to call PDAF_init. These variables control some functionality of the user-supplied routines for the data assimilation system and are shared with these routines through mod_assimilation. These variables are for example:

• delt_obs: An integer specifying the number of time steps between two analysis steps
• rms_obs: Assumed observation error
• local_range: Localization radius in grid points for the observation domain in LSEIK and LETKF
• locweight: Type of localizing weight

It is useful to define variables like this at this central position. Of course, this definition has to be adapted to the particluar model used.

User-supplied routine U_ensemble_init

The user-supplied routine the we named U_ensemble_init here, is called by PDAF through the defined interface described below. The routine is called by all MPI processes that compute the filter analysis step (i.e. those for which 'filterpe' is set to true. In the standard configuration of init_parallel_pdaf these are all processes of the first model task, i.e. task_id=1.) U_ensemble_init is only called by PDAF_init if no error occurred before; thus the status flag is zero.

The interface is the following:

SUBROUTINE U_ensemble_init(filtertype, dim_p, dim_ens, &
                           state_p, Uinv, ens_p, flag)

with

• filtertype: The integer defining the type of filter algorithm as given in the call to PDAF_init
• dim_p: An integer holding the size of the state dimension for the calling process as specified in the call to PDAF_init
• dim_ens: An integer holding the size of the ensemble (or the rank of the state covariance matrix for SEEK)
• state_p: A real array of size (dim_p) for the local model state of the calling process (Only relevant for mode-based filters)
• Uinv: A real array of size (dim_ens-1, dim_ens-1) for the inverse of matrix U from the decomposition of the state error covariance matrix P = VUV^T (Only relevant for mode-based filters like SEEK.)
• ens_p: A real array of size (dim_p, dim_ens) the has to hold upon exit the ensemble of model states.
• flag: Status flag for PDAF. It is 0 upon entry and can be set by in the user-supplied routine, depending on the success of the ensemble initialization. Preferably, values above 102 should be used for failures to avoid conflicts with the error codes defined within PDAF_init.

In the initialization routine U_ensemble_init one has to distinguish between ensemble-based and mode-based filters. The only mode based filter supplied with PDAF is SEEK, while all other methods are ensemble-based.

Initialization for ensemble-based filters

Generally, the filters SEIK, LSEIK, EnKF, ETKF, and LETKF are ensemble-based filters.For these filters only the array ens_p needs to be initialized by the ensemble of model states. If a model with domain decomposition is used, the full ensemble for the local sub-domain of the MPI process has to be initialized.

The arrays state_p and Uinv are allocated to their correct sizes because they are used during the assimilation cycles. They are not yet initialized and it is allowed to use these arrays in the initialization. An exception from this is EnKF for which Uinv is allocated only with size (1,1), because Uinv is not using for EnKF.

Initialization for mode-based filters

The only mode-based filter supplied with PDAF is currenly the SEEK filter. For this filter the initialization bases on the decomposition of the state error covariance matrix in the form P = VUV^T. According to this decomposition, the array ens_p has to be initialized to hold the modes from matrix V and Uinv holds the inverse of matrix U. In addition state_p has to be initialized with the initial state estimate. If a model with domain decomposition is used, the part of all modes for the local sub-domain of the MPI process and the corresponding part of the state vector has to be initialized.