Changes between Initial Version and Version 1 of OMI_additional_functionality_PDAF3

Timestamp:

May 27, 2025, 3:53:07 PM (6 days ago)

Author:

lnerger

Comment:

--

OMI_additional_functionality_PDAF3

@@ +1 @@
+= Additional Functionality of PDAF-OMI =
+
+{{{
+#!html
+<div class="wiki-toc">
+<h4>PDAF-OMI Guide</h4>
+<ol><li><a href="PDAF_OMI_Overview_PDAF3">Overview</a></li>
+<li><a href="OMI_Callback_obs_pdafomi_PDAF3">callback_obs_pdafomi.F90</a></li>
+<li><a href="OMI_observation_modules_PDAF3">Observation Modules</a></li>
+<li><a href="OMI_observation_operators_PDAF3">Observation operators</a></li>
+<li><a href="OMI_error_checking_PDAF3">Checking error status</a></li>
+<li><a href="OMI_debugging_PDAF3">Debugging functionality</a></li>
+<li>Additional OMI Functionality</li>
+<li><a href="Porting_to_OMI_PDAF3">Porting an existing implemention to OMI</a></li>
+<br>
+<li>Related pages in Implementation Guide<li>
+<ol>
+<li><a href="ImplementationofAnalysisStep_PDAF3">Implementing the analysis step</a></li>
+<li><a href="nondiagonal_observation_error_covariance_matrices_PDAF3">Using nondiagonal R-matrices</a></<li>
+</ol>
+</ol>
+</div>
+}}}
+
+[[PageOutline(2-3,Contents of this page)]]
+PDAF-OMI provide some additional functionality which is described here.
+
+The additional functionalities are:
+- [wiki:PDAFomi_additional_functionality_PDAF3#Non-isotropiclocalization Non-isotropic localization]
+- [wiki:PDAFomi_additional_functionality_PDAF3#Omittingobservationsthatarepotentialoutliers Omitting observations that are potential outliers]
+- [wiki:PDAFomi_additional_functionality_PDAF3#Usingdomain-limitedobservations Using domain-limited observations]
+- [wiki:PDAFomi_additional_functionality_PDAF3#Initializationroutines Initialization routines] (added in PDAF V2.3.1)
+
+== Non-isotropic localization ==
+
+The localization cut-off radius and the support radius of the localization function are specified in each observation module in the calls to [wiki:OMI_observation_modules_PDAF3#init_dim_obs_l_OBSTYPE init_dim_obs_l_OBSTYPE] and [wiki:OMI_observation_modules_PDAF3#init_dim_obs_l_OBSTYPE localize_covar_OBSTYPE]. 
+
+If the cut-off and support radius are defined as scalar variables, the localization is **isotropic** (same radius in all directions). This was the only option before PDAF V2.2. 
+
+With PDAF V2.2 a **non-isotropic localization** is supported. For this one has to define the variables for the cut-off and support radii as vectors with the same length as the number of coordinates used in the distance calculation of the localization. Then each element of the vectors specifies a directional radius (e.g. in longitude and latitude in geographic coordinates; if 3-dimensional localization is used, one can in addition specify a separate radius for the vertical direction) and a separate radius can be specified for each direction.
+
+The non-isotropic localization is computed as an ellipse in 2 dimensions and as an ellipsoid in 3 dimensions. Thus, the distance between the model grid point and an observation (or between two model grid points in case of the LEnKF) is compared to a directional cut-off radius of the ellipse or ellipsoid that is computed from the elements of the vector specifying the cut-off radii (likewise fo rthe support radii).
+
+
+
+== Omitting observations that are potential outliers ==
+
+If the difference between observations and the ensemble mean state if very large, the observations can have the effect that the data assimilation increment becomes very large and this might lead to unrealistic increments. A large difference might indicate that the observation value is not reliable (thus an 'outlier'), but it can also happen if a forecast state estimate is e.g. biased or if an initial state estimate is far away from the current situation measured by the observations.
+
+PDAF V2.2 introduced an option to omit observations if their difference from the ensemble mean is too large. Here 'too large' is measured as a multiple of the observation error variance which is compared to the square of the difference. 
+
+=== Activating observation omission === 
+
+Whether observations are omitted is controlled by the variable
+{{{
+   REAL :: thisobs%inno_omit    ! Omit observation if absolute innovation is larger this factor times (default=0.0)
+                                !     observation standard deviation (only active for >0)
+}}}
+For the default value `thisobs%inno_omit=0.0` no observation are omitted. 
+To use this omission functionality set in [wiki:OMI_observation_modules#init_dim_obs_OBSTYPE init_dim_obs_OBSTYPE] in the observation modules `thisobs%inno_omit` to a value > 0.0. This activates the omission and specifies the omission limit. If thisobs%inno_omit=2.0, an observation would be omitted if the absolute value of the difference between the observed ensemble mean state and the observation value is larger than 2 times the observation error standard deviation.
+
+To omit an observation, its inverse observation error variance is set by default to the very small value
+{{{
+   REAL :: thisobs%inno_omit_ivar=1.0e-12  ! Value of inverse variance to omit observation
+                                           !     (should be much larger than actual observation error variance)
+}}}
+Thus, the observation is prescribed to have a large error variance and hence has no effect in the assimilation. (The observation is not removed from the assimilation, but just made irrelevant). Depending on the observation, this default value might not be small enough to make an observation irrelevant. If this is the case, one can also specify a different value of `thisobs%inno_omit_ivar` in the observation module. 
+
+|| **Note:** From PDAF V2.3.1 one can also call the routines [wiki:PDAFomi_set_inno_omit] and [wiki:PDAFomi_set_inno_omit_ivar] instead of setting the values directly in the code ||
+
+When `thisobs%inno_omit` is set to >0, you should see a screen output like
+{{{
+PDAFomi     --- Exclude obs. type ID   1  if innovation^2 >     3.00 times obs. error variance
+}}}
+(Note that is line is a bit misleading. It's actually 'if innovation > 3.00 times obs. error standard deviation')
+
+=== Statistics about omitted observations ===
+
+The number of omitted observation can be seen from the output of the observation statistics. For domain-localized filters, like the LESTKF the output looks like the following:
+{{{
+PDAFomi     --- Global statistics for locally omitted observations:
+PDAFomi        Local domains with omitted observations:               607
+PDAFomi        Local domains without omitted observations:             41
+PDAFomi        Maximum number of locally omitted observations:          3
+PDAFomi        Maximum number of locally used observations:            14
+PDAFomi        Total avg. locally omitted observations:      1.5
+PDAFomi        Total avg. locally used observations:         8.1
+PDAFomi        Avg. omitted for domains with omitted observations:          1.6
+PDAFomi        Avg. used for domains with omitted observations:             8.2
+}}}
+(This example uses the tutorial case `tutorial/offline_2D_serial`: We add `thisobs%inno_omit=3.0` in `obs_A_pdafomi.F90` and run as `./PDAF_offline -dim_ens 9 -filtertype 7 -cradius 10`. Thus, we omit observations if the squared difference between model and observations is large than 3 observation error variancse). In this case observations were omitted in 607 local analysis domains, while no observations were omitted in 41 domains. The maximum number of locally omitted observations if 3, while on average 1.5 observations were omitted. 
+
+In case of the LEnKF, the output looks like:
+{{{
+PDAFomi         Global statistics for omitted observations:
+PDAFomi           Global number of omitted observations:          4
+PDAFomi           Global number of used observations:            24
+}}}
+
+
+== Using domain-limited observations == 
+
+=== Overview (thisobs%use_global_obs=0) ===
+
+In the domain-localized filters (LESTK, LETKF, LSEIK, LNETF) observations are assimilated that are located within the localization around some grid point. When a model uses parallelization with domain-decomposition some of these observations might belong to a different process-domain. In the default mode (`thisobs%use_global_obs`=1) PDAF-OMI gathers all globally available observations so that each process has access to all observations. 
+
+It can be more efficient to limit the observations on a process-domain to those observations that are located inside the domain or within the localization radius around it. Then, in the local analyses less observations have to be checked for their distance. Setting `thisobs%use_global_obs=0` (or `CALL PDAFomi_set_use_global_obs(thisobs,0)` from PDAF V2.3.1) activates this feature. However, it needs additional preparations to make PDAF-OMI aware of the limiting coordinates of a process sub-domain.
+
+In order to make use of the domain-limited observations, one has to provide PDAF-OMI with the limiting coordinates of a process-subdomain. There are two routines that can to this task:
+- `PDAFomi_set_domain_limits`
+- `PDAFomi_get_domain_limits_unstrc`
+One of these routines has to be called in `init_pdaf` to set the domain information. Their use is described below.
+
+After this call, OMI is set up to use domain-limited observations with `thisobs%use_global_obs=0`. The necessary operation will be performed by PDAF-OMI inside `PDAFomi_gather_obs` and inside the observation operators. This feature has to be activated separately for each observation type. However, mixing the settings 0 and 1 is possible.
+
+'''Note:''' The domain-limitation is only implemented to work in two dimensions. 
+
+
+=== `PDAFomi_set_domain_limits` ===
+
+In this routine one provide the limiting coordinates of the process domain so that PDAF-OMI can store the information.
+
+The interface is
+{{{
+  SUBROUTINE PDAFomi_set_domain_limits(lim_coords)
+
+    REAL, INTENT(in) :: lim_coords(2,2)     !< coordinate array (1: longitude, 2: latitude)
+                                            !< geographic range: longitude (-pi, pi), latitude (-pi/2, pi/2)
+                                            !< Cartesian range: (x) coordinate grows from left to right; (y) from bottom to top
+}}}
+here `lim_coords` are
+{{{
+           (2,1)
+   (:,1)+---------+
+        |         |           - (1,1) longitude of the western edge of the domain (or x-coordinate of the left edge)
+        |         |           - (1,2) longitude of the eastern edge of the domain (or x-coordinate of the right edge)
+  (1,1) |         | (1,2)     - (2,1) latitude of the northern edge of the domain (or y-coordinate of the upper edge)
+        |         |           - (2,2) latitude of the southern edge of the domain (or y-coordinate of the lower edge)
+        |         |
+        +---------+(:,2)
+           (2,2)
+}}}
+Thus, (:,1) specifies the north-western corner of the sub-domain and (:,2) the souther-estern corner. The first (x) coordinate grows from left to right, while the second coordinate (y) increases from bottom to the top.
+
+If the model grid is not decomposed in cardinal directions, but e.g. rotated, the coordinates should specife the extrema. Thus, lim_coords(1,1) would be the coordinate of the northernmost grid point of a domain.
+
+Note:
+- The domain limits are coded for two dimensions, which are usually the horizontal directions. If the observations have e.g. three dimensions, the first two will be used for the domain-limitation. In case of geographic coordinates these are longituare and latitude.
+
+
+=== `PDAFomi_get_domain_limits_unstrc` ===
+
+This routine is find the extreme coordinates for a model domain. The routine is provided with the coordinates of all grid pints of a domain and then find the limiting coordinates. It is designed for unstructured grid and we have only tested it with the ocean model FESOM. (The tricky part is when a process-domain crosses the date line).
+
+The interface is
+{{{
+  SUBROUTINE PDAFomi_get_domain_limits_unstr(verbose, npoints_p, coords_p)
+
+    INTEGER, INTENT(in) :: verbose          !< verbosity flag (1: write output)
+    INTEGER, INTENT(in) :: npoints_p        !< number of process-local grid points
+    REAL, INTENT(in) :: coords_p(:,:)       !< geographic coordinate array, dimension (2, npoints_p)
+                                            !<   (row 1: longitude, 2: latitude)
+                                            !<   ranges: longitude (-pi, pi), latitude (-pi/2, pi/2)
+}}}
+
+This function only supports geographic coordinates given in radians.
+
+== Determining observation localization weights ==
+
+When using non-diagonal observation error covariance matrices **R**, one has to code the localization weighting in the user-supplied routines. To obtain the localization weight according to the distance between the local analysis domain and the observation, one can utilize the routine `PDAFomi_observation_localization_weights, which ensures that the OMI-internal information is correctly handled.
+
+=== PDAFomi_observation_localization_weights ===
+
+This routine was introduced with PDAF V2.3.
+
+This routine returns a vector of observation localization weights computed according to the localization specifications in `thisobs_l` and observation coordinates in `thisobs`. The routine can be called in `prodRinvA_l_pdafomi` or `likelihood_l_pdafomi` in order to obtain the localization weights for observations. For diagonal observation error covariance matrices, the weighting is computed by OMI internally, but for nondiagonal R-matrices, the use had to implement the application of the localization. 
+
+The interface is the following:
+{{{
+  SUBROUTINE PDAFomi_observation_localization_weights(thisobs_l, thisobs, ncols, &
+             A_l, weight, verbose)
+
+    TYPE(obs_l), INTENT(inout) :: thisobs_l  ! Data type with local observation
+    TYPE(obs_f), INTENT(inout) :: thisobs    ! Data type with full observation
+    INTEGER, INTENT(in) :: ncols             ! Rank of initial covariance matrix
+    REAL, INTENT(in) :: A_l(:, :)            ! Input matrix, size (thisobs_l%dim_obs_l, ncols)
+    INTEGER, INTENT(in) :: verbose           ! Verbosity flag
+    REAL, INTENT(out) :: weight(thisobs_l%dim_obs_l) ! Localization weights
+}}}
+
+
+Notes:
+ * Array `A_l` is usually the input argument to the routine [wiki:prodRinvA_l_pdaf prodRinvA_l_pdaf]. It is only used when the regulated localization is used. 
+ * The regulated localization computes an alternative localization function based on the ration of state error variances to observation error variances. the regulated localization method is described in the paper Nerger et al. (Quarterly Journal of the Royal Meteorological Society, 138 (2012) 802-812; see [wiki:PublicationsandPresentations publications].)
+
+
+== Initialization routines ==
+
+With PDAF V2.3.1 we have introduced a number of routines than can be used to initialize different quantities in `thisobs`. With these routines can can avoid to set these variables in `thisobs` directly in the code. One can also avoid to allocate arrays in `thisobs` in the user code, which can be useful e.g. for programming languages like C. The routines are
+ * [wiki:PDAFomi_set_disttype]
+ * [wiki:PDAFomi_set_doassim]
+ * [wiki:PDAFomi_set_domainsize]
+ * [wiki:PDAFomi_set_icoeff_p]
+ * [wiki:PDAFomi_set_id_obs_p]
+ * [wiki:PDAFomi_set_inno_omit]
+ * [wiki:PDAFomi_set_inno_omit_ivar]
+ * [wiki:PDAFomi_set_ncoord]
+ * [wiki:PDAFomi_set_obs_err_type]
+ * [wiki:PDAFomi_set_use_global_obs]