= PDAFlocal - the localization module of PDAF =


[[PageOutline(2-3,Contents of this page)]]

== Overview ==

PDAFlocal is optional functionality and was introduced with PDAF V2.3. 

PDAFlocal provides an improved handling of the localization of state vectors as used in the loca filters LESTKF, LETKF, LNETF, LSEIK, and LKNETF. A leading motivation to implement PDAFlocal was that it leads to much better performance when the Python interface pyPDAF is used. However, it is also useful for native implementations in e.g. Fortran because it also provides additional functionality for e.g. vertical localization or weakly coupled assimilation.

Compared to not using PDAFlocal removes the need to implement the routines `g2l_state_pdaf` and 'l2g_state_pdaf' (see e.g. [wiki:ImplementAnalysisLocal]). For this, PDAFlocal provides three initialization routines and a set of alternative routines that are called for the analysis step. 

Please note that if you have an existing implementation there is no direct need to switch to using PDAFlocal. If might be useful if you plan to use the additional functionality of `PDAFlocal_set_increment_weights` (see below).

PDAFlocal provides the initialization routines
 * `PDAFlocal_set indices` : This routine stores the index array that describes the indices of the element of a local state vector in the global state vector. It is used for mapping between both the local and global vectors.
 * `PDAFlocal_set_increment_weights` : This routine allows to provide PDAFlocal with a vector of local increment weights. When the local state vector is written back to the global state vector after a local analysis update, this allows to prescribe increment weights for each element of the local state vector. With this, one can e.g. implement a vertical localization if the local state vector is a column of the model grid.
 * `PDAFlocal_clear_increment_weights` : This routine clears and deallocates increment weights which have been set before with `PDAFlocal_set_increment_weights`. After executing this routine, the weights are reset to one for all elements of the local state vector (i.e. unweighted increment)
These routines are called in `init_dim_l_pdaf`.

== Initialization routines ==

Here we provide an overview of the interfaces of the 3 initialization routines.

=== PDAFlocal_set_indices ===

This routine is called in `init_dim_l_pdaf` when the dimension of the local state vector has to be determined. In addition to this dimension, the implementation guide recommends to also initialize the coordinates of the local analysis domain and the indices of the elements ofthe local state vector in the non-localized (global or domain-decomposed) state vector `state_p`. Thise indices are then used in the mapping from the global to the local state vector before the analysis and back after the analysis update in the local analysis loop. With PDAFlocal, one uses `PDAFlocal_set_indices` to pass the vector of indices to PDAFlocal, which later performs the mapping for all ensemble states.

The interface is:
{{{
SUBROUTINE PDAFlocal_set_indices(dim_l, map)

  INTEGER, INTENT(in) :: dim_l          ! Dimension of local state vector
  INTEGER, INTENT(in) :: map(dim_l)     ! Index array for mapping
}}}

Hints:
 * For complex cases it might feel more intuitive to perform the mapping between the global and local state vectors in an explicit loop. However, this might lead to repeated computation of the indices and can hence be less efficient than computing the indices before. 
 * The initialization of hthe index vector `map` is analogous to a loop that directly performs the initialization of a local state vector. However, here only the indices are stored.


=== PDAFlocal_set_increment_weights ===

This routine provides the optional functionality to prescribe the assimilation increment of each element of the local state vector a different weight. This can be used, e.g. to implement a vertical localization in the case that the local state vector is a full vertical column of the model grid. In this case, one can make the increment weight depending on the height (or depth) of a grid point. Another application is to implement weakly-coupled assimilation in which the local state vector contains all variables, but only a subset of them is updated. This is achieved by givening those element that should not be updated the weight 0.

The routine can be called in `init_dim_l_pdaf` when the dimension of the local state vector has to be determined. In addition to this dimension, the implementation guide recommends to also initialize the coordinates of the local analysis domain and the indices of the elements ofthe local state vector in the non-localized (global or domain-decomposed) state vector `state_p`. In addition one can initialize a vector of increment weights and provide it to PDAFlocal by calling this routine.

The interface is:
{{{
SUBROUTINE PDAFlocal_set_increment_weights(dim_l, weights)

  INTEGER, INTENT(in) :: dim_l          ! Dimension of local state vector
  REAL, INTENT(in) :: weights(dim_l)    ! Weights array
}}}

Hints:
 * The loop to initialize `weights` can be analogous to the loop that initializes the index array for mapping between the global and local state vector.
 * It is possible to clear the weights array by calling `PDAFlocal_clear_increment_weights` (see below). Afterwards, unit weights are used.
 * One can call this routine again to change the weights. One does not need to call `PDAFlocal_clear_increment_weights` before. In any case, one has to be careful to ensure that the size of `weights` remains consistent when the size of the local state vector changes.


=== PDAFlocal_clear_increment_weights ===

If `PDAFlocal_set_increment_weights` was used to prescribe the assimilation increment of each element of the local state vector a different weight this routine can be used to clear and deallocate the weights array. Afterwards a weight of one is used for the increment of all elements of the local state vector when initializing the global state vector after the analysis update in the local analysis loop.

The interface is  without arguments:
{{{
SUBROUTINE PDAFlocal_clear_increment_weights()
}}}

Hint:
 * After calling this routine, the local increment weights are one until one calls `PDAFlocal_set_increment_weights` with prescribed weights again.

== Routines to perform the analysis step with OMI ==

Generally the routines to be called for the analysis step are analogous to those when not used PDAFlocal (e.g. in implementations before the release of PDAF V2.3). The difference is in the name and in the fact that the routines `g2l_state_pdaf` and `l2g_state_pdaf` are not present in the interface of the PDAFlocal routines.

=== Ensemble filters ===

**for diagonal R matrix**

Depending on the fully parallel are flexible implementation choose the routine 
||= **fully parallel** =||= **flexible** =||
|| [wiki:PDAFlocalomi_assimilate] || [wiki:PDAFlocalomi_put_state] || 

**For non-diagonal R matrix**

See [wiki:OMI_nondiagonal_observation_error_covariance_matrices] for information on using non-diagonal **R**-matrices with OMI. The routines are only partly generic depending on the needed observation-specific routine:
||= **filter** =||= **fully parallel** =||= **flexible** =||
||= LESTKF[[BR]]LETKF[[BR]] LSEIK[[BR]] =|| [wiki:PDAFlocalomi_assimilate_nondiagR] ||[wiki:PDAFlocalomi_put_state_nondiagR] ||
||= LNETF =|| [wiki:PDAFlocalomi_assimilate_lnetf_nondiagR] ||[wiki:PDAFlocalomi_put_state_lnetf_nondiagR] ||
||= LLNETF =|| [wiki:PDAFlocalomi_assimilate_lknetf_nondiagR] ||[wiki:PDAFlocalomi_put_state_lknetf_nondiagR] ||

=== 3D-Var ===

Only En3DVar and hybrid 3D-Var use a local filter. The routines are the following.

**for diagonal R matrix**

||= **Method** =||= **fully parallel** =||= **flexible** =||
||= En3DVar =|| [wiki:PDAFlocalomi_assimilate_en3dvar_lestkf] || [wiki:PDAFlocalomi_put_state_en3dvar_lestkf] || 
||= Hyb3DVar =|| [wiki:PDAFlocalomi_assimilate_hyb3dvar_lestkf] || [wiki:PDAFlocalomi_put_state_hyb3dvar_lestkf] || 


**for non-diagonal R matrix**

See [wiki:OMI_nondiagonal_observation_error_covariance_matrices] for information on using non-diagonal **R**-matrices with OMI.

||= **Method** =||= **fully parallel** =||= **flexible** =||
||= En3DVar =|| [wiki:PDAFlocalomi_assimilate_en3dvar_lestkf_nondiagR] || [wiki:PDAFlocalomi_put_state_en3dvar_lestkf_nondiagR] || 
||= Hyb3DVar =|| [wiki:PDAFlocalomi_assimilate_hyb3dvar_lestkf_nondiagR] || [wiki:PDAFlocalomi_put_state_hyb3dvar_lestkf_nondiagR] || 


== Routines to perform the analysis step not using OMI ==

Thesse are the routines for the older observation handling without OMI. There is no distinction with regard to diagonal or non-diagonal **R** because without PDAF-OMI one always needs to implement all observation routines on the user side.

Generally the routines to be called for the analysis step are analogous to those when not used PDAFlocal (e.g. in implementations before the release of PDAF V2.3). The difference is in the name and in the fact that the routines `g2l_state_pdaf` and `l2g_state_pdaf` are not present in the interface of the PDAFlocal routines.

=== Ensemble filters ===

Depending on the fully parallel are flexible implementation choose the routine 
||= **Filter** =||= **fully parallel** =||= **flexible** =||
||= LESTKF =|| [wiki:PDAFlocal_assimilate_lestkf] || [wiki:PDAFlocal_put_state_lestkf] || 
||= LETKF =|| [wiki:PDAFlocal_assimilate_letkf] || [wiki:PDAFlocal_put_state_letkf] || 
||= LSEIK =|| [wiki:PDAFlocal_assimilate_lseik] || [wiki:PDAFlocal_put_state_lseik] || 
||= LNETF =|| [wiki:PDAFlocal_assimilate_lnetf] || [wiki:PDAFlocal_put_state_lnetf] || 
||= LKNETF =|| [wiki:PDAFlocal_assimilate_lknetf] || [wiki:PDAFlocal_put_state_lknetf] || 

=== 3D-Var ===

Only En3DVar and hybrid 3D-Var use a local filter. The routines are the following.

**for diagonal R matrix**

||= **Method** =||= **fully parallel** =||= **flexible** =||
||= En3DVar =|| [wiki:PDAFlocal_assimilate_en3dvar_lestkf] || [wiki:PDAFlocal_put_state_en3dvar_lestkf] || 
||= Hyb3DVar =|| [wiki:PDAFlocal_assimilate_hyb3dvar_lestkf] || [wiki:PDAFlocal_put_state_hyb3dvar_lestkf] || 

== Porting to PDAFlocal ==

If you like to port you existing code to using PDAFlocal, the steps are the following
1. Change the call to the analysis routine from `PDAFomi_` to `PDAFlocalomi_` (or `PDAF_` to `PDAFlocal_` for non-OMI implementations)
1. Remove `g2l_state_pdaf` and `l2g_state_pdaf` from the interface of the `PDAFlocalomi_`
 routine
1. In `init_dim_l_pdaf` ensure that the index array for the mapping between the global and local state vectors is implemented (this should exist if the implementation followed the tutorial and templates where it is the array `id_lstate_in_pstate`).
1. Insert the call to `PDAFlocal_set_indices`.

After these steps it should be possible to recompile and run. If this was successful, one can safely remove the files holding `g2l_state_pdaf` and `l2g_state_pdaf`.