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
• 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 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