= User-provided routine for initialization of local observations =

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

PDAF-OMI provides the routine [wiki:PDAFomi_init_dim_obs_l] to initialize local observations in a routine [wiki:OMI_observation_modules#init_dim_obs_l_OBSTYPE init_dim_obs_p_OBSTYPE] of an PDAF-OMI observation module. `PDAFomi_init_dim_obs_l` provides a convenient functionality with observation search for all isotropic and non-isotropic localization variants in different dimensions. However, because the routine covers all cases with various IF-checks it can be relatively slow. 

While we improved the search performance in PDAF V2.3, we also found that a user-provided routine that is designed for the particular case of the user's application can be significantly faster. Thus, if you have the impression that `PDAFomi_init_dim_obs_l` is 'too slow' in your application, you can implement 'your own' case -specific variant of `PDAFomi_init_dim_obs_l`. In PDAF V2.3 we added routines that help in the implementation, by initialized parameters and allocating arrays of PDAF-OMI, so that one does not need to perform such oeprations in the user code.

== Structure of a user-provided routine ==

An **example** of a user-provided initialization routine for local observations is the file `init_dim_obs_l_pdafomi_user.F90` in `/tutorial/offline_2D_parallel`.

The **calling interface** of the user-provided routine should usually be the same as that of `PDAFomi_init_dim_obs_l`. This is for isotropic localization:
{{{
  SUBROUTINE init_dim_obs_l_pdafomi_user(thisobs_l, thisobs, coords_l, locweight, cradius, &
       sradius, cnt_obs_l_all)

    USE PDAFomi, &            ! Include observation type definitions and OMI-provided routines
         ONLY: obs_f, obs_l, &
         PDAFomi_set_dim_obs_l, PDAFomi_set_localization

    IMPLICIT NONE

    TYPE(obs_f), INTENT(inout) :: thisobs    !< Data type with full observation
    TYPE(obs_l), TARGET, INTENT(inout) :: thisobs_l  !< Data type with local observation
    REAL, INTENT(in) :: coords_l(thisobs%ncoord)     !< Coordinates of current analysis domain
    INTEGER, INTENT(in) :: locweight         !< Type of localization function
    REAL, INTENT(in) :: cradius              !< Vector of localization cut-off radii
    REAL, INTENT(in) :: sradius              !< Vector of support radii of localization function
    INTEGER, INTENT(inout) :: cnt_obs_l_all  !< Local dimension of observation vector over all observation types
}}}

The **structure of the body** of the routine needs to be as follows
{{{

    IF (thisobs%doassim == 1) THEN

       CALL PDAFomi_set_localization(thisobs_l, cradius, sradius, locweight)

       ... Loop over full observations and count number 'cnt_obs' of observations within localization radius ...

       CALL PDAFomi_set_dim_obs_l(thisobs_l, thisobs, cnt_obs_l_all, cnt_obs)

       ... Initialize local observation arrays ...

           thisobs_l%id_obs_l(cnt_l) = i                        ! index of local obs. in full obs. vector
           thisobs_l%distance_l(cnt_l) = SQRT(distance2)        ! distance
           thisobs_l%cradius_l(cnt_l) = thisobs_l%cradius(1)    ! cut-off radius
           thisobs_l%sradius_l(cnt_l) = thisobs_l%sradius(1)    ! support radius

    ENDIF
}}}

**Notes:**
 * The check `thisobs%doassim==1` is mandatory because init_dim_obs_l_OBSTYPE is called for all observations.
 * `cnt_l` is the counter of the valid local observations inside the loop, while `i` is the index of a local observation on the full observation vector
 * The initialization of the local observation arrays of `thisobs_l` can be done by the routine `PDAFomi_store_obs_l_index` provided by PDAF-OMI (see below). However, in our tests this always led to a slower program. Thus, for performance reasons we recommend to use the direct initialization if possible.
 * In the tutorial example routine we use a module to enclose the user routine. This allows for syntax checking, but is not mandatory.
 * `cnt_obs_l_all` should only be passed through from PDAF to the PDAFomi routines. If it is changed by the user, the consistency of the observation-handling inside PDAF-OMI will be broken.

== Recommendations for optimizing performance ==

Based on our experience with different application cases we recommend the following:
 * Implement as specific as possible for your case. IF-checks always deteriorate performance and should be limited as much as possible. However, also costly operations like sine or square root should be minimized. 
 * Compute the distance checks using the square distance `distance2` and only compute the actual distance as `SQRT(distance2)` only for those observations that are within the localization cut-off radius. Since compute a square-root is a costly operation, one should only compute it if needed.
 * For localization in 2 and 3 dimensions with geographic coordinates it can be useful to first check for the distance in one direction. If a vertical localization is used it can be good to check the vertical distance first and only if this is within the cut-off radius proceed to the horizontal distance. One might also check in the order height, latitude, longitude to first determine a box (or rectangle in 2 dimensional cases, omitting height) before calculating the full distance. This can be particularly useful when non-isotropic localization is used. (In PDAF-OMI we strictly first check whether an observation is within a box or rectangle and only for these observations we compute the actual squared distance.)
 * Non-isotropic localization might be implemented by scaling the direction and then computing the radii in this scaled isotropic coordinates. However, this approach can be problematic in case of periodicity. (E.g. computing the haversine function for geographic distances doesn't seem to be possible then). To this end, for non-isotropic localization we implemented in PDAF-OMI a localization ellipse (in 2D) or ellipsoid (in 3D). This avoids scaling in any direction. 
 * As a reference for possible operations you can check the routines `PDAFomi_check_dist2_loop` (isotropic) and `PDAFomi_check_dist2_noniso_loop` (non-isotropic) in `/src/PDAFomi_dim_obs_l`. Generally you could copy one of these routines, paste it into the user code and then adapt it by removing those IF-cases that don't hold for your case.


== PDAF-OMI provided helper routines ==

=== PDAFomi_set_localization ===

This routine stores the localization parameters (cradius, sradius, locweight). It further allocates localization parameter arrays of `thisobs_l`. The inputs to the routine are values provided as arguments in the call to `init_dim_obs_l_pdafomi_user`.

The full interface of the routine is
{{{
  SUBROUTINE PDAFomi_set_localization(thisobs_l, cradius, sradius, locweight)

    TYPE(obs_l), INTENT(inout) :: thisobs_l  ! Data type with local observation
    REAL, INTENT(in) :: cradius              ! Localization cut-off radius
    REAL, INTENT(in) :: sradius              ! Support radius of localization function
    INTEGER, INTENT(in) :: locweight         ! Type of localization function
}}}

For non-isotropic localization there is the alternative routine [wiki:PDAFomi_set_localization_noniso].

=== PDAFomi_set_dim_obs_l ===

This routine stores the local number of observations `cnt_obs` and uses the value of `cnt_obs_l_all` for OMI-internal initializations. The routine also does further initializations for the OMI-internal handling of observations.

The full interface is
{{{
  SUBROUTINE PDAFomi_set_dim_obs_l(thisobs_l, thisobs, cnt_obs_l_all, cnt_obs_l)

    TYPE(obs_f), INTENT(inout) :: thisobs    ! Data type with full observation
    TYPE(obs_l), INTENT(inout) :: thisobs_l  ! Data type with local observation
    INTEGER, INTENT(inout) :: cnt_obs_l_all  ! Local dimension of observation vector over all obs. types
    INTEGER, INTENT(inout) :: cnt_obs_l      ! Local dimension of single observation type vector
}}}
here `thisobs`, `thisobs_l`, and `cnt_obs_l_all` are just passed to the routine.

== Additional helper routines ==

We recommend to perform the initialization of the vectors in `thisobs_l` at the last step of `init_dim_obs_l_pdafomi_user` in this user-provided routine. For convenience we, however, provide a routine for this step. 

=== PDAFomi_store_obs_l_index ===

This routine takes the index 'cnt_l' in the example above, the distance as well as the local cut-off radius (`cradius_l`) and support radius (`sradius_l`) and stores the values in the vectors of `thisobs_l`.

{{{
  SUBROUTINE PDAFomi_store_obs_l_index(thisobs_l, cnt_l, id_obs_l, distance, &
       cradius_l, sradius_l)

    TYPE(obs_l), INTENT(inout) :: thisobs_l  ! Data type with local observation
    INTEGER, INTENT(in) :: cnt_l             ! Element of local observation array to be filled
    INTEGER, INTENT(in) :: id_obs_l          ! Index of local observation in full observation array
    REAL, INTENT(in) :: distance             ! Distance between local analysis domain and observation
    REAL, INTENT(in) :: cradius_l            ! cut-off radius for this local observation
                                             !  (directional radius in case of non-isotropic localization)
    REAL, INTENT(in) :: sradius_l            ! support radius for this local observation
                                             !  (directional radius in case of non-isotropic localization)
}}}

=== PDAFomi_store_obs_l_index_vdist ===

For 2+1D factorized localization the distance in the vertical direction has to be stored separately. For this the alternativve routine `PDAFomi_store_obs_l_index_vdist` is used.

{{{
  SUBROUTINE PDAFomi_store_obs_l_index_vdist(thisobs_l, cnt_l, id_obs_l, distance, &
       cradius_l, sradius_l, vdist)

    TYPE(obs_l), INTENT(inout) :: thisobs_l  ! Data type with local observation
    INTEGER, INTENT(in) :: cnt_l             ! Element of local observation array to be filled
    INTEGER, INTENT(in) :: id_obs_l          ! Index of local observation in full observation array
    REAL, INTENT(in) :: distance             ! Distance between local analysis domain and observation
    REAL, INTENT(in) :: cradius_l            ! cut-off radius for this local observation
                                             !  (directional radius in case of non-isotropic localization)
    REAL, INTENT(in) :: sradius_l            ! support radius for this local observation
                                             !  (directional radius in case of non-isotropic localization)
    REAL, INTENT(in) :: vdist               ! support radius in vertical direction for 2+1D factorized localization
}}}