Version 4 (modified by 5 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 still 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, whilei
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 routinePDAFomi_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 asSQRT(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) andPDAFomi_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