| | 1 | = User-provided routine for initialization of local observations = |
| | 2 | |
| | 3 | 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. |
| | 4 | |
| | 5 | 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. |
| | 6 | |
| | 7 | == Structure of a user-provided routine == |
| | 8 | |
| | 9 | 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`. |
| | 10 | |
| | 11 | 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: |
| | 12 | {{{ |
| | 13 | SUBROUTINE init_dim_obs_l_pdafomi_user(thisobs_l, thisobs, coords_l, locweight, cradius, & |
| | 14 | sradius, cnt_obs_l_all) |
| | 15 | |
| | 16 | USE PDAFomi, & ! Include observation type definitions and OMI-provided routines |
| | 17 | ONLY: obs_f, obs_l, & |
| | 18 | PDAFomi_set_dim_obs_l, PDAFomi_set_localization |
| | 19 | |
| | 20 | IMPLICIT NONE |
| | 21 | |
| | 22 | TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation |
| | 23 | TYPE(obs_l), TARGET, INTENT(inout) :: thisobs_l !< Data type with local observation |
| | 24 | REAL, INTENT(in) :: coords_l(thisobs%ncoord) !< Coordinates of current analysis domain |
| | 25 | INTEGER, INTENT(in) :: locweight !< Type of localization function |
| | 26 | REAL, INTENT(in) :: cradius !< Vector of localization cut-off radii |
| | 27 | REAL, INTENT(in) :: sradius !< Vector of support radii of localization function |
| | 28 | INTEGER, INTENT(inout) :: cnt_obs_l_all !< Local dimension of observation vector over all observation types |
| | 29 | }}} |
| | 30 | |
| | 31 | The **structure of the body** of the routine needs to be as follows |
| | 32 | {{{ |
| | 33 | |
| | 34 | IF (thisobs%doassim == 1) THEN |
| | 35 | |
| | 36 | CALL PDAFomi_set_localization(thisobs_l, cradius, sradius, locweight) |
| | 37 | |
| | 38 | ... Loop over full observations and count number 'cnt_obs' of observations within localization radius ... |
| | 39 | |
| | 40 | CALL PDAFomi_set_dim_obs_l(thisobs_l, thisobs, cnt_obs_l_all, cnt_obs) |
| | 41 | |
| | 42 | ... Initialize local observation arrays ... |
| | 43 | |
| | 44 | thisobs_l%id_obs_l(cnt_obs_l) = i ! index of local obs. in full obs. vector |
| | 45 | thisobs_l%distance_l(cnt_obs_l) = SQRT(distance2) ! distance |
| | 46 | thisobs_l%cradius_l(cnt_obs_l) = thisobs_l%cradius(1) ! cut-off radius |
| | 47 | thisobs_l%sradius_l(cnt_obs_l) = thisobs_l%sradius(1) ! support radius |
| | 48 | |
| | 49 | ENDIF |
| | 50 | }}} |
| | 51 | |
| | 52 | **Notes:** |
| | 53 | * The check `thisobs%doassim==1` is mandatory because init_dim_obs_l_OBSTYPE is called for all observations. |
| | 54 | * `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 |
| | 55 | * 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. |
| | 56 | * In the tutorial example routine we use a module to enclose the user routine. This allows for syntax checking, but is not mandatory. |
| | 57 | * 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. |
| | 58 | * `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. |
| | 59 | |
| | 60 | |
| | 61 | == PDAF-OMI provided helper routines == |
| | 62 | |
| | 63 | === PDAFomi_set_localization === |
| | 64 | |
| | 65 | 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`. |
| | 66 | |
| | 67 | The full interface of the routine is |
| | 68 | {{{ |
| | 69 | SUBROUTINE PDAFomi_set_localization(thisobs_l, cradius, sradius, locweight) |
| | 70 | |
| | 71 | TYPE(obs_l), INTENT(inout) :: thisobs_l ! Data type with local observation |
| | 72 | REAL, INTENT(in) :: cradius ! Localization cut-off radius |
| | 73 | REAL, INTENT(in) :: sradius ! Support radius of localization function |
| | 74 | INTEGER, INTENT(in) :: locweight ! Type of localization function |
| | 75 | }}} |
| | 76 | |
| | 77 | For non-isotropic localization there is the alternative routine [wiki:PDAFomi_set_localization_noniso]. |
| | 78 | |
| | 79 | === PDAFomi_set_dim_obs_l === |
| | 80 | |
| | 81 | 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. |
| | 82 | |
| | 83 | The full interface is |
| | 84 | {{{ |
| | 85 | SUBROUTINE PDAFomi_set_dim_obs_l(thisobs_l, thisobs, cnt_obs_l_all, cnt_obs_l) |
| | 86 | |
| | 87 | TYPE(obs_f), INTENT(inout) :: thisobs ! Data type with full observation |
| | 88 | TYPE(obs_l), INTENT(inout) :: thisobs_l ! Data type with local observation |
| | 89 | INTEGER, INTENT(inout) :: cnt_obs_l_all ! Local dimension of observation vector over all obs. types |
| | 90 | INTEGER, INTENT(inout) :: cnt_obs_l ! Local dimension of single observation type vector |
| | 91 | }}} |
| | 92 | here `thisobs`, `thisobs_l`, and `cnt_obs_l_all` are just passed to the routine. |
| | 93 | |
| | 94 | == Additional helper routines == |
| | 95 | |
| | 96 | 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. |
| | 97 | |
| | 98 | === PDAFomi_store_obs_l_index === |
| | 99 | |
| | 100 | 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`. |
| | 101 | |
| | 102 | {{{ |
| | 103 | SUBROUTINE PDAFomi_store_obs_l_index(thisobs_l, idx, id_obs_l, distance, & |
| | 104 | cradius_l, sradius_l) |
| | 105 | |
| | 106 | TYPE(obs_l), INTENT(inout) :: thisobs_l ! Data type with local observation |
| | 107 | INTEGER, INTENT(in) :: idx ! Element of local observation array to be filled |
| | 108 | INTEGER, INTENT(in) :: id_obs_l ! Index of local observation in full observation array |
| | 109 | REAL, INTENT(in) :: distance ! Distance between local analysis domain and observation |
| | 110 | REAL, INTENT(in) :: cradius_l ! cut-off radius for this local observation |
| | 111 | ! (directional radius in case of non-isotropic localization) |
| | 112 | REAL, INTENT(in) :: sradius_l ! support radius for this local observation |
| | 113 | ! (directional radius in case of non-isotropic localization) |
| | 114 | }}} |
| | 115 | |
| | 116 | === PDAFomi_store_obs_l_index_vdist === |
| | 117 | |
| | 118 | 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. |
| | 119 | |
| | 120 | {{{ |
| | 121 | SUBROUTINE PDAFomi_store_obs_l_index_vdist(thisobs_l, idx, id_obs_l, distance, & |
| | 122 | cradius_l, sradius_l, vdist) |
| | 123 | |
| | 124 | TYPE(obs_l), INTENT(inout) :: thisobs_l ! Data type with local observation |
| | 125 | INTEGER, INTENT(in) :: idx ! Element of local observation array to be filled |
| | 126 | INTEGER, INTENT(in) :: id_obs_l ! Index of local observation in full observation array |
| | 127 | REAL, INTENT(in) :: distance ! Distance between local analysis domain and observation |
| | 128 | REAL, INTENT(in) :: cradius_l ! cut-off radius for this local observation |
| | 129 | ! (directional radius in case of non-isotropic localization) |
| | 130 | REAL, INTENT(in) :: sradius_l ! support radius for this local observation |
| | 131 | ! (directional radius in case of non-isotropic localization) |
| | 132 | REAL, INTENT(in) :: vdist ! support radius in vertical direction for 2+1D factorized localization |
| | 133 | }}} |
