= Implementation of the Analysis step for the ENSRF/EAKF (Ensemble Square Root and Ensemble Adjustment Filters) = {{{ #!html }}} [[PageOutline(2-3,Contents of this page)]] || This page describes the implementation of the analysis step without using PDAF-OMI. Please see the [wiki:ImplementationofAnalysisStep 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 [ModifyModelforEnsembleIntegration Modification of the model code for the ensemble integration] and its sub-page on [InsertAnalysisStep 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_statecollect_state_pdaf.F90 U_collect_state]: [[BR]]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_statedistribute_state_pdaf.F90 U_distribute_state]: [[BR]]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_finit_dim_obs_f_pdaf.F90 U_init_dim_obs_f]: [[BR]]The name of the user-supplied routine that provides the size of the full observation vector * [#U_obs_op_fobs_op_f_pdaf.F90 U_obs_op_f]: [[BR]]The name of the user-supplied routine that acts as the full observation operator on some state vector * [#U_init_obs_finit_obs_f_pdaf.F90 U_init_obs_f]: [[BR]]The name of the user-supplied routine that initializes the full vector of observations * [#U_init_obsvars_finit_obsvars_pdaf.F90 U_init_obsvars_f]: [[BR]]The name of the user-supplied routine that initializes the vector of observation error variances. * [#U_localize_covar_seriallocalize_covar_serial_pdaf.F90 U_localize_covar_serial]: [[BR]]The name of the routine that applies the covariance localization for a single observation * [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep]: [[BR]]The name of the user-supplied pre/poststep routine as in `PDAF_get_state` * [#U_next_observationnext_observation_pdaf.F90 U_next_observation]: [[BR]]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`: [[BR]]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 [ModifyModelforEnsembleIntegration 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 [ModifyModelforEnsembleIntegration 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 [InsertAnalysisStep#U_collect_statecollect_state_pdaf.F90 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 [InsertAnalysisStep#U_distribute_statedistribute_state_pdaf.F90 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 [wiki:PDAF_gather_obs_f 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 [ModifyModelforEnsembleIntegration#U_prepoststepprepoststep_ens_pdaf.F90 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 [AuxiliaryRoutines auxiliary routines]) === `U_next_observation` (next_observation_pdaf.F90) === This routine is independent of the filter algorithm used. See the page on [InsertAnalysisStep#U_next_observationnext_observation_pdaf.F90 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_statecollect_state_pdaf.F90 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_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep] (Call to act on the forecast ensemble, called with negative value of the time step) 1. [#U_init_dim_obsinit_dim_obs_pdaf.F90 U_init_dim_obs] 1. [#U_obs_opobs_op_pdaf.F90 U_obs_op] (`dim_ens` calls: one call for each ensemble member; one more call if also applied to ensemble mean) 1. [#U_init_obsinit_obs_pdaf.F90 U_init_obs] 1. [#U_init_varsinit_obsvars_pdaf.F90 U_init_obsvars] In the loop over all single observations, only one call-back routine is executed: 1. [#U_localize_covar_seriallocalize_covar_serial_pdaf.F90 U_localize_covar_serial] (once for each single observation) After the loop over all local analysis domains, it is executed: 1. [#U_prepoststepprepoststep_ens_pdaf.F90 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_statedistribute_state_pdaf.F90 U_distribute_state] 1. [#U_next_observationnext_observation_pdaf.F90 U_next_observation]

Implementation Guide