Changes between Initial Version and Version 1 of PDAFomi_additional_functionality

Timestamp:

Mar 8, 2024, 11:59:02 AM (2 months ago)

Author:

lnerger

Comment:

--

PDAFomi_additional_functionality

@@ +1 @@
+= Additional Functionality of PDAF-OMI =
+
+{{{
+#!html
+<div class="wiki-toc">
+<h4>PDAF-OMI Guide</h4>
+<ol><li><a href="PDAF_OMI_Overview">Overview</a></li>
+<li><a href="OMI_Callback_obs_pdafomi">callback_obs_pdafomi.F90</a></li>
+<li><a href="OMI_observation_modules">Observation Modules</a></li>
+<li><a href="OMI_observation_operators">Observation operators</a></li>
+<li><a href="OMI_error_checking">Checking error status</a></li>
+<li><a href="OMI_debugging">Debugging functionality</a></li>
+<li><a href="OMI_ImplementationofAnalysisStep">Implementing the analysis step with OMI</a></li>
+<ol>
+<li> <a href="ImplementFilterAnalysisOverview"> General overview for ensemble filters</a></li>
+<ol>
+<li><a href="ImplementAnalysisGlobal">Implementation for Global Filters</a></li>
+<li><a href="ImplementAnalysisLocal">Implementation for Local Filters</a></li>
+<li><a href="ImplementAnalysislenkfOmi">Implementation for LEnKF</a></li>
+</ol>
+<li> <a href="Implement3DVarAnalysisOverview"> General overview for 3D-Var methods</a></li>
+<ol>
+<li><a href="ImplementAnalysis_3DVar">Implementation for 3D-Var</a></li>
+<li><a href="ImplementAnalysis_3DEnVar">Implementation for 3D Ensemble Var</a></li>
+<li><a href="ImplementAnalysis_Hyb3DVar">Implementation for Hybrid 3D-Var</a></li>
+</ol>
+</ol>
+<li><a href="Porting_to_OMI">Porting an existing implemention to OMI</a></li>
+<li>Using domain-limited observations</li>
+</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#Non-isotropiclocalization Non-isotropic localization]
+- [wiki:PDAFomi_additional_functionality#Omittingobservationsthatarepotentialoutliers Omitting observations that are potential outliers]
+- [wiki:PDAFomi_additional_functionality#Usingdomain-limitedobservations Using domain-limited observations]
+
+== 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#init_dim_obs_l_OBSTYPE init_dim_obs_l_OBSTYPE] and [wiki:OMI_observation_modules#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 === 
+
+To use this omission functionality one has to augment [wiki:OMI_observation_modules#init_dim_obs_OBSTYPE] in the observation modules as follows: [[BR]]
+Set the variable
+{{{
+   REAL :: thisobs%inno_omit=0.0    ! Omit observation if squared innovation larger this factor times
+                                    !     observation variance (only active for >0)
+}}}
+to a value > 0.0. This activates the omission and specifies the omission limit. If inno_omit=2.0, an observation would be omitted if the squared difference between the observed ensemble mean state and the observation value is larger than 2 times the observation error variance.
+
+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 sufficient 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. 
+
+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
+}}}
+
+=== 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` 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.