PDAF-OMI Observation Operators

The observation operator is called for each observation type in the routine obs_op_pdafomi in the file callback_obs_pdafomi.F90.

OMI currently provides 3 observation operators:

• PDAFomi_obs_op_gridpoint
  This observation operator is used for the case that observations are model variables located at grid points. Thus, the operation is to select single element from the state vector according to the index array thisobs%id_obs_p initialized in init_dim_obs_f_TYPE.
• PDAFomi_obs_op_gridavg
  This observation operator is used for the case that observations are the average of model variables at grid points. The averages are computed according to the number of rows in the index array thisobs%id_obs_p initialized in init_dim_obs_f_TYPE.
• PDAFomi_obs_op_interp_lin
  This observation operator computes the observation by linear interpolation. It uses the index array thisobs%id_obs_p and the array thisobs%icoeff_p holding interpolation coefficients initialized in init_dim_obs_f_TYPE. To use this observation operator, one has to allocate and initialize thisobs%icoeff_p as described below.

The arguments of the observation operators are

       CALL PDAFomi_obs_op_X (thisobs,[nrows,] state_p, ostate)

Where thisobs is the observation type variable, state_p is the input state vector provided by PDAF and ostate is the observed state that will be returned to PDAF. nrows only exists for the observation operators X=gridavg and X=interp_lin and specifies the number of grid points involved in the observation operation. For X=gridpoint, this argument does not exist.

Initializing interpolation coefficients

The observation operator PDAFomi_obs_op_interp_lin requires that the interpolation coefficients have been initialized in the array thisobs%icoeff_p. This initialization is performed in init_dim_obs_TYPE. PDAF-OMI provides three routines for this task:

• PDAFomi_get_interp_coeff_lin1D
  Simplified initialization for 1-dimensional models
• PDAFomi_get_interp_coeff_lin
  Determine interpolation coefficients based on the coordinates of grid points and the observation for a rectangular grid in 1, 2, or 3 dimensions.
• PDAFomi_get_interp_coeff_tri
  Determine barycentric interpolation coefficients for triangular grids based on the coordinates of grid points and the observation

An example of initializing interpolation coefficients with PDAFomi_get_interp_coeff_lin and of using PDAFomi_obs_op_interp_lin is provided in tutorial/online_2D_serialmodel_omi/obs_C_pdafomi.F90.

PDAFomi_get_interp_coeff_lin

The call to this routine is

  CALL PDAFomi_get_interp_coeff_lin(num_gp, n_dim, gcoords, ocoord, icoeff)

    INTEGER, INTENT(in) :: num_gp         !< Length of thisobs%icoeff_p(:,i)
    INTEGER, INTENT(in) :: n_dim          !< Number of dimensions in interpolation
    REAL, INTENT(in)    :: gcoords(:,:)   !< Coordinates of grid points
    REAL, INTENT(in)    :: ocoord(:)      !< Coordinates of observation (one column ocoord_p(:,i))
    REAL, INTENT(inout) :: icoeff(:)      !< Interpolation coefficients (one column thisobs%icoeff_p(:,i))

Here it is required that num_gp=2 for n_dim=1; num_gp=4 for n_dim=2; num_gp=8 for n_dim=3.

In the array gcoords in init_dim_obs_TYPE, the first index specifies the grid point while the second specifies the coordinate, thus gcoords(j,:) is the list of coordinates for grid point j. The coordinates need to be consistent with the indices specified in thisobs%id_obs_p since these specify the elements of the state vector that are interpolated. Only the first n_dim entries of ocoord will be used for the interpolation.

ocoord_p(:,i) holds the list of the coordinates for the observation with index i.

The order of the coordinates and coefficients is the following:

                       (7)------(8) 
                       /|       /|    with
                     (5)+-----(6)|       - column 1
                      | |      | |       / column 2
                      |(3)-----+(4)      | column 3
                      |/       |/
                     (1) ---- (2)

   thus gcoords(1,1)/=gcoords(2,1),        but  gcoords(1,1)=gcoords(3,1)=gcoords(5,1),
        gcoords(1,2)/=gcoords(3,2),             gcoords(1,2)=gcoords(2,2)=gcoords(5,2), 
        gcoords(1,3)/=gcoords(5,3)              gcoords(1,3)=gcoords(2,3)=gcoords(3,3)

Notes:

• For 1D linear interpolation (n_dim=1) only the coordinates for grid points 1 and 2 are used to compute the coefficients
• For bi-linear interpolation (n_dim=2) only the coordinates for grid points 1, 2, and 3 are used to compute the coefficients
• For tri-linear interpolation (n_dim=3) only the coordinates for grid points 1, 2, 3, and 5 are used to compute the coefficients