= User-provided routine for initialization of local observations =

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_obs_l) = i                        ! index of local obs. in full obs. vector
           thisobs_l%distance_l(cnt_obs_l) = SQRT(distance2)        ! distance
           thisobs_l%cradius_l(cnt_obs_l) = thisobs_l%cradius(1)    ! cut-off radius
           thisobs_l%sradius_l(cnt_obs_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_obs_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
 * We recommend to compute use 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.
 * In the tutorial example routine we use a module to enclose the user routine. This allows for syntax checking, but is not mandatory.
 * 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, for performance reasons we recommend to use the direct initialization if possible.
 * `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.


== 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_obs_l' in the example above (`idx` below), 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, idx, id_obs_l, distance, &
       cradius_l, sradius_l)

    TYPE(obs_l), INTENT(inout) :: thisobs_l  ! Data type with local observation
    INTEGER, INTENT(in) :: idx               ! 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, idx, id_obs_l, distance, &
       cradius_l, sradius_l, vdist)

    TYPE(obs_l), INTENT(inout) :: thisobs_l  ! Data type with local observation
    INTEGER, INTENT(in) :: idx               ! 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
}}}