Version 3 (modified by 4 weeks ago) (diff) | ,
---|

# User-provided routine for initialization of local observations

#### Contents of this page

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