Implementation of the Analysis step for the ENSRF/EAKF (Ensemble Square Root and Ensemble Adjustment Filters)

This page describes the implementation of the analysis step without using PDAF-OMI. Please see the page on the analysis with OMI for the more modern and efficient implementation variant using PDAF-OMI.

The ENSRF/EAKF filters were added with verson 3.0 of PDAF.

Overview

The ENSRF and EAKF are ensemble Kalman filter variants using serial observation processing. The implementation in PDAF follows Houtekamer and Hamill (2002) for the ENSRF and Anderson (2003) for the EAKF variant using local least squares regression. The parallelization follows Anderson and Collins (2007), where the 'high-latency' variant is implemented, since their 'low-latency' variant appears to be overly optimistic while still leading to a partial serialization of the execution.

The variant of the serial-observation processing filter is selected by the subtype in the call to PDAF_init as follows:

subtype	Filter variant
0	ENSRF (Whitaker & Hamill, 2002)
1	EAKF linear regression (Anderson, 2003)

For the analysis step of the ENSRF and EAKF, different operations related to the observations are needed. These operations are requested by PDAF by calling user-supplied routines. Intentionally, the operations are split into separate routines in order to keep the operations rather elementary. This procedure should simplify the implementation. The names of the required routines are specified in the call to the routine PDAF_assimilate_ensrf for the fully-parallel and flexible parallelization implementations (alternatively PDAF_put_state_lenkf for the 'flexible' implementation and for the offline-coupled assimilation). With regard to the parallelization, all these routines are executed by the filter processes (filterpe=.true.) only.

For completeness we discuss here all user-supplied routines that are specified in the interface to PDAF_assimilate_ensrf. Thus, some of the user-supplied routines that are explained on the page explaining the modification of the model code for the ensemble integration are repeated here.

In our study Nerger et al. (2015), we discussed that applying localization can lead to stability issues of the ENSRF. The filter performs a loop over all single observations and with localization the assimilation result depends on the order in which the observations are assimilated. This actually leads to the effect that the assimilation result at some grid point does not only depend on the observations with the localization radius r, but also on observations further away. This happens if these observation at larger distance are assimilated before the observations within the radius r and if they influence the state close to distance r. In this case the innovation for the state close to distance r is already change when assimilating the observation within the radius and hence lead to a different result. This effect has implications on the parallelization since keeping the observation order constant over the full model domain leads to a partial serialization of the algorithm. In the implementation in PDAF, we use the parallelization approach that does not guarantee the same order of the observations. Usually, the differences when changing the observation order are small, but the benefit is a better scaling since the serialization is avoided. Nonetheless, we generally recommend using LESTKF or LETKF, or their global variants ESTKF or ETKF, since they no not depend explicitly on the observation order, and they allow for non-diagonal observation error covariance matrices. However, the ENSRF/EAKF might have a good compute performance.

PDAF_assimilate_ensrf

The general aspects of the filter specific routines PDAF_assimilate_* have been described on the page Modification of the model code for the ensemble integration and its sub-page on inserting the analysis step. The routine is used in the fully-parallel and the flexible implementation variants of the data assimilation system. When the offline model is used, the routines PDAF_put_state_* are used. These have also been used in previous PDAF releases for the 'flexible' implementation variant. Here, we list once more the full interface of the routine. Subsequently, the full set of user-supplied routines specified in the call to PDAF_assimilate_ensrf is explained.

The interface when using the LEnKF is the following:

  SUBROUTINE PDAF_assimilate_ensrf(U_collect_state, U_distribute_state, &
       U_init_dim_obs_f, U_obs_op_f, U_init_obs_f, U_init_obsvars_f, &
       U_localize_covar_serial, &
       U_prepoststep, U_next_observation, outflag)

with the following arguments:

• U_collect_state:
  The name of the user-supplied routine that initializes a state vector from the array holding the ensemble of model states from the model fields. This is the inverse operation to U_distribute_state used in PDAF_get_state as well as here.
• U_distribute_state:
  The name of a user supplied routine that initializes the model fields from the array holding the ensemble of model state vectors. This is the inverse operation to U_collect_state
• U_init_dim_obs_f:
  The name of the user-supplied routine that provides the size of the full observation vector
• U_obs_op_f:
  The name of the user-supplied routine that acts as the full observation operator on some state vector
• U_init_obs_f:
  The name of the user-supplied routine that initializes the full vector of observations
• U_init_obsvars_f:
  The name of the user-supplied routine that initializes the vector of observation error variances.
• U_localize_covar_serial:
  The name of the routine that applies the covariance localization for a single observation
• U_prepoststep:
  The name of the user-supplied pre/poststep routine as in PDAF_get_state
• U_next_observation:
  The name of a user supplied routine that initializes the variables nsteps, timenow, and doexit. The same routine is also used in PDAF_get_state.
• status_pdaf:
  The integer status flag. It is zero, if the routine is exited without errors.

PDAF_put_state_ensrf

For the offline mode of PDAF, the routine PDAF_put_state_ensrf has to be used instead of PDAF_assimilate_ensrf. This routine can also be used when the 'flexible' implementation variant is chosen for the assimilation system, The general aspects of the filter specific routines PDAF_put_state_* have been described on the page Modification of the model code for the ensemble integration. The interface of the routine is identical with that of PDAF_assimilate_ensrf with the exception that the arguments of the user-supplied routines U_distribute_state and U_next_observation are missing.

The interface is the following:

SUBROUTINE PDAF_put_state_ensrf(U_collect_state, &
     U_init_dim_obs_f, U_obs_op_f,  U_init_obs_f, U_init_obsvars_f, &
     U_localize_covar_serial, &
     U_prepoststep, outflag)

User-supplied routines

Here all user-supplied routines are described that are required in the call to PDAF_assimilate_ensrf. For some of the generic routines, we link to the page on modifying the model code for the ensemble integration.

To indicate user-supplied routines we use the prefix U_. In the tutorials in tutorial/ and in the template directory templates/ these routines exist without the prefix, but with the extension _pdaf. The files are named correspondingly. In the section titles below we provide the name of the template file in parentheses.

In the subroutine interfaces some variables appear with the suffix _p. This suffix indicates that the variable is particular to a model sub-domain, if a domain decomposed model is used. Thus, the value(s) in the variable will be different for different model sub-domains.

U_collect_state (collect_state_pdaf.F90)

This routine is independent of the filter algorithm used. See the page on inserting the analysis step for the description of this routine.

U_distribute_state (distribute_state_pdaf.F90)

This routine is independent of the filter algorithm used. See the page on inserting the analysis step for the description of this routine.

U_init_dim_obs_f (init_dim_obs_f_pdaf.F90)

This routine is used by the ENSRF and by all filter algorithms with domain-localization (LSEIK, LETKF, LNETF, LKNETF) and is independent of the particular algorithm.

The interface for this routine is:

SUBROUTINE init_dim_obs_f(step, dim_obs_f)

  INTEGER, INTENT(in)  :: step       ! Current time step
  INTEGER, INTENT(out) :: dim_obs_f  ! Dimension of full observation vector

For the ENSRF the routine is called at the beginning of each analysis step, before the loop over all single observations is entered. It has to initialize the size dim_obs_f of the full observation vector according to the current time step. For simplicity, dim_obs_f can be the size for the global model domain.

Some hints:

• We recommend to not only determine the size of the observation vector at this point. The routine is a good place to also already gather information about the corresponding indices of the state vector needed later to implement the observation operator. In addition, one can already prepare an array that holds the full observation vector, an array storing the coordinates of the observations and possible an array storing observation error variances (if the observation error covariance matrix is diagonal). The required arrays can be defined in a module like mod_assimilation. The information can be used later in U_localize_covar_serial.

Note: PDAF-OMI provides a structured approach for implementing the observation functionality. Because of this we generally recommend to ue the PDAF3 interface that uses PDAF-OMI.

• The routine is similar to init_dim_obs used in the global filters. However, if the global filter is used with a domain-decomposed model, it only initializes the size of the observation vector for the local model sub-domain. This is different for the local filters, as the local analysis also requires observational data from neighboring model sub-domains. Nonetheless, one can base on an implemented routine init_dim_obs to implement init_dim_obs_f.

U_obs_op_f (obs_op_f_pdaf.F90)

This routine is used by the ENSRF and by all filter algorithms with domain-localization (LSEIK, LETKF, LNETF, LKNETF) and is independent of the particular algorithm.

The interface for this routine is:

SUBROUTINE obs_op_f(step, dim_p, dim_obs_f, state_p, m_state_f)

  INTEGER, INTENT(in) :: step               ! Current time step
  INTEGER, INTENT(in) :: dim_p              ! PE-local dimension of state
  INTEGER, INTENT(in) :: dim_obs_f          ! Dimension of the full observed state
  REAL, INTENT(in)    :: state_p(dim_p)     ! PE-local model state
  REAL, INTENT(out) :: m_state_f(dim_obs_f) ! Full observed state

The routine is called during the analysis step, before the loop over the single observations is entered. It has to perform the operation of the observation operator acting on a state vector, which is provided as state_p. The observed state has to be returned in m_state_f. It is the observed state corresponding to the 'full' observation vector.

Hint:

• The routine is similar to init_dim_obs used for the global filters. However, with a domain-decomposed model m_state_f will need to contain parts of the state vector from neighboring model sub-domains. Thus, one needs to collect this information which resides in the memory of other processes. PDAF provides the routine PDAF_gather_obs_f for this task. The example implementation in tutorial/classical/online_2D_parallelmodel shows the use of PDAF_gather_obs_f.

U_init_obs_f (init_obs_f_pdaf.F90)

This routine is used by the ENSRF and by all filter algorithms with domain-localization (LSEIK, LETKF, LNETF, LKNETF) and is independent of the particular algorithm.

The interface for this routine is:

SUBROUTINE init_obs_f(step, dim_obs_f, observation_f)

  INTEGER, INTENT(in) :: step                     ! Current time step
  INTEGER, INTENT(in) :: dim_obs_f                ! Dimension of full observation vector
  REAL, INTENT(out)   :: observation_f(dim_obs_f) ! Full observation vector

The routine is called during the analysis step before the loop over the single observations is entered. It has to provide the full vector of observations in observation_f for the current time step.

Hints:

• As for the other 'full' routines: While the global counterpart of this routine (init_obs) has to initialize the observation vector only for the local model sub-domain, the 'full' routine has to include observations that spatially belong to neighboring model sub-domains. As an easy choice one can simply initialize a vector of all globally available observations.

U_init_obsvars_f (init_obsvars_f_pdaf.F90)

This routine is only used by the ENSRF/EAKF.

The interface for this routine is:

SUBROUTINE init_obsvars_f(step, dim_obs_f, var_f)

  INTEGER, INTENT(in) :: step           ! Current time step
  INTEGER, INTENT(in) :: dim_obs_f      ! Dimension of full observation vector
  REAL, INTENT(out) :: var_f(dim_obs_f) ! vector of observation error variances

The routine is called during the analysis step before the loop over the single observations is entered. It has to provide in var_f a vector of observation error variances corresponding to the full vector of observations.

U_localize_covar_serial (localize_covar_serial_pdaf.F90)

This routine is only used by the ENSRF/EAKF.

The interface for this routine is:

  SUBROUTINE U_localize_covar_serial(iobs, dim_p, dim_obs, HP_p, HXY_p)

  INTEGER, INTENT(in) :: iobs                  !< Index of the assimilated single observation
  INTEGER, INTENT(in) :: dim_p                 !< Process-local state dimension
  INTEGER, INTENT(in) :: dim_obs_f             !< Number of full observations
  REAL, INTENT(inout) :: HP_p(dim_p)           !< Process-local part of matrix HP for observation iobs
  REAL, INTENT(inout) :: HXY_p(dim_obs_F)      !< Process-local part of matrix HX(HX_full) for full observations}}}

The routine is called during the loop over all single observations. The purpose of the routine is to apply covariance localization to the vectors Hi P and Hi PH^T for the assimilation of a single observation (determined by the index iobs related to the observation operator Hi). Here Hi PH^T is the vector relating to the observed covariance matrix for the full observation vector. This vector is required for parallelization.

Hints:

• To compute the localization one can use the routine PDAF_local_weight after computing the distance between two elements in the vector Hi P or Hi PH^T.

U_prepoststep (prepoststep_ens_pdaf.F90)

The general aspects of this routines have already been described on the page on modifying the model code for the ensemble integration for the SEIK filter. For completeness, the description is repeated specifically for the EnKF:

The interface of the routine is identical for all filters, but sizes can vary. Also, the particular operations that are performed in the routine can be specific for each filter algorithm.

The interface for this routine is for the LEnKF

SUBROUTINE prepoststep(step, dim_p, dim_ens, dim_ens_p, dim_obs_p, &
                       state_p, Uinv, ens_p, flag)

  INTEGER, INTENT(in) :: step        ! Current time step
                         ! (When the routine is called before the analysis -step is provided.)
  INTEGER, INTENT(in) :: dim_p       ! PE-local state dimension
  INTEGER, INTENT(in) :: dim_ens     ! Size of state ensemble
  INTEGER, INTENT(in) :: dim_ens_p   ! PE-local size of ensemble
  INTEGER, INTENT(in) :: dim_obs_p   ! PE-local dimension of observation vector
  REAL, INTENT(inout) :: state_p(dim_p) ! PE-local forecast/analysis state
                                     ! The array 'state_p' is not generally not initialized in the case of SEIK/EnKF/ETKF.
                                     ! It can be used freely in this routine.
  REAL, INTENT(inout) :: Uinv(1, 1)  ! Not used not LEnKF
  REAL, INTENT(inout) :: ens_p(dim_p, dim_ens)      ! PE-local state ensemble
  INTEGER, INTENT(in) :: flag        ! PDAF status flag

The routine U_prepoststep is called once at the beginning of the assimilation process. In addition, it is called during the assimilation cycles before the analysis step and after the ensemble transformation. The routine is called by all filter processes (that is filterpe=1).

The routine provides for the user the full access to the ensemble of model states. Thus, user-controlled pre- and post-step operations can be performed. For example the forecast and the analysis states and ensemble covariance matrix can be analyzed, e.g. by computing the estimated variances. In addition, the estimates can be written to disk.

Hint:

• If a user considers to perform adjustments to the estimates (e.g. for balances), this routine is the right place for it.
• Only for the SEEK filter the state vector (state_p) is initialized. For all other filters, the array is allocated, but it can be used freely during the execution of U_prepoststep.
• The array Uinv is not used in the EnKF. Internally to PDAF, it is allocated to be of size (1,1).
• Apart from the size of the array Uinv, the interface is identical for all ensemble filters (SEIK/ETKF/EnKF/LSEIK/LETKF/LEnKF). In general it should be possible to use an identical pre/poststep routine for all these filters.
• The interface through which U_prepoststep is called does not include the array of smoothed ensembles. In order to access the smoother ensemble array one has to set a pointer to it using a call to the routine PDAF_get_smootherens (see page on auxiliary routines)

U_next_observation (next_observation_pdaf.F90)

This routine is independent of the filter algorithm used. See the page on inserting the analysis step for the description of this routine.

Execution order of user-supplied routines

For the ENSRF/EAKF methods, the user-supplied routines are essentially executed in the order they are listed in the interface to PDAF_assimilate_ensrf. The order can be important as some routines can perform preparatory work for later routines. For example, U_init_dim_obs can prepare an index array that provides the information for executing the observation operator in U_obs_op.

Before the analysis step is called the following routine is executed:

1. U_collect_state

The analysis step is executed when the ensemble integration of the forecast is completed. During the analysis step the following routines are executed in the given order:

1. U_prepoststep (Call to act on the forecast ensemble, called with negative value of the time step)
2. U_init_dim_obs
3. U_obs_op (dim_ens calls: one call for each ensemble member; one more call if also applied to ensemble mean)
4. U_init_obs
5. U_init_obsvars

In the loop over all single observations, only one call-back routine is executed:

1. U_localize_covar_serial (once for each single observation)

After the loop over all local analysis domains, it is executed:

1. U_prepoststep (Call to act on the analysis ensemble, called with (positive) value of the time step)

In case of the routine PDAF_assimilate_ensrf, the following routines are executed after the analysis step:

1. U_distribute_state
2. U_next_observation