Additional Functionality of PDAF-OMI

PDAF-OMI provide some additional functionality which is described here.

The additional functionalities are:

• Non-isotropic localization
• Omitting observations that are potential outliers
• Using domain-limited observations
• 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 init_dim_obs_l_OBSTYPE and 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 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 PDAFomi_set_inno_omit and 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 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 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

• PDAFomi_set_disttype
• PDAFomi_set_doassim
• PDAFomi_set_domainsize
• PDAFomi_set_icoeff_p
• PDAFomi_set_id_obs_p
• PDAFomi_set_inno_omit
• PDAFomi_set_inno_omit_ivar
• PDAFomi_set_ncoord
• PDAFomi_set_obs_err_type
• PDAFomi_set_use_global_obs