User-provided routine for initialization of local observations

PDAF-OMI provides the routine PDAFomi_init_dim_obs_l to initialize local observations in a routine 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
• 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 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