PDAF-OMI Observation Operators

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

Observation operators

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_OBSTYPE.
• 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_OBSTYPE.
• 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_OBSTYPE. To use this observation operator, one has to allocate and initialize thisobs%icoeff_p as described below.

The arguments of the observation operators are

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

     TYPE(obs_f), INTENT(inout) :: thisobs  ! Data type with full observation
   [ INTEGER, INTENT(in) :: nrows           ! Number of values to be averaged ]
     REAL, INTENT(in)    :: state_p(:)      ! Process-local model state (provided by PDAF)
     REAL, INTENT(inout) :: obs_f_all(:)    ! Full observed state for all observation types (array provided by PDAF)

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_OBSTYPE. 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/obs_C_pdafomi.F90.

PDAFomi_get_interp_coeff_lin

This routine computes interpolation coefficients for a rectangular grid in 1, 2, or 3 dimensions.

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, the first index specifies the grid point while the second specifies the coordinate. Thus, gcoords(1,1) is the first coordinate for grid point 1, gcoords(1,2) is the second coordiate for grid point 1, while gcoords(2,1) is the first coordiate for grid point 2. 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(:,i) holds the list of the coordinates for the observation with index i (different from the use in gcoords)

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

PDAFomi_get_interp_coeff_lin1D

This routine is a simplified variant of PDAFomi_get_interp_coeff_lin. It computes interpolation coefficients in 1 dimension only

The call to this routine is

  CALL PDAFomi_get_interp_coeff_lin1D(gcoords, ocoord, icoeff)

    REAL, INTENT(in)    :: gcoords(:)  ! Coordinates of grid points; dim=2
    REAL, INTENT(in)    :: ocoord      ! Coordinates of observation
    REAL, INTENT(inout) :: icoeff(:)   ! Interpolation coefficients; dim=2

PDAFomi_get_interp_coeff_tri

For triangular model grid interpolation coefficients are determined as barycentric coordinates. This is performed in this routine.

The call to this routine is

  CALL PDAFomi_get_interp_coeff_tri(gcoords, ocoord, icoeff)

    REAL, INTENT(in)    :: gcoords(:)      ! Coordinates of grid points; dim=(3,2)
    REAL, INTENT(in)    :: ocoord(:)       ! Coordinates of observation; dim=2
    REAL, INTENT(inout) :: icoeff(:)       ! Interpolation coefficients; dim=3

Notes:

• In the array gcoords, 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 order of the grid points in gcoords has to be consistent with the order of the indices specified in thisobs%id_obs_p

Adjoint observation operators

For the application of 3D-Var, adjoint observation operators are required. These perform the operation of the transposed linear observation operator.

OMI provides the adjoint observation operators corresponding to the forward observation operators:

• PDAFomi_obs_op_adj_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_OBSTYPE.
• PDAFomi_obs_op_adj_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_OBSTYPE.
• PDAFomi_obs_op_adj_interp_lin
  This observation operator is used for the case of 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_OBSTYPE. 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_adj_X (thisobs,[nrows,] state_p, ostate)

     TYPE(obs_f), INTENT(inout) :: thisobs  ! Data type with full observation
     INTEGER, INTENT(in) :: nrows           ! Number of values to be averaged
     REAL, INTENT(in)    :: obs_f_all(:)    ! Full observed state for all observation types (array provided by PDAF)
     REAL, INTENT(inout) :: state_p(:)      ! Process-local model state (provided by PDAF)

Where thisobs is the observation type variable, ostate is the input vector in the observation space provided by PDAF and state_p is the output state vector 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.

Implementing your own observation operator

The current set of observation operators provided by PDAF-OMI is rather fundamental. However, there are also observation types which are not variables of state vector, but functions of several values. Likewise one might want to use a more sophisticated interpolation than the linear one or some interpolation that treats the horizontal and vertical directions separately. For these cases you can implement you own operators.

The template observation operator in /template/omi/ob_op_pdafomi_TEMPLATE.F90 can be used as the basis for the implementation. It describes the steps needs in the implementation.

Each observation operator include the following functionality:

1. Check whether thisobs%doassim==1, which indicates that the observation is assimilated
2. initialize the observation vector ostate_p for the observation type for which the routine is called. For a parallel model this is done for for the process-local domain (for a model without parallelization this is the global domain)
3. Call PDAFomi_gather_obsstate to gather the full observation vector obs_f_all. This call is mandatory even if the model is not parallelized. The interface is:

     SUBROUTINE PDAFomi_gather_obsstate(thisobs, ostate_p, obs_f_all)
   
       TYPE(obs_f), INTENT(inout) :: thisobs  ! Data type with full observation
       REAL, INTENT(in) :: ostate_p(:)        ! Vector of process-local observed state
       REAL, INTENT(inout) :: obs_f_all(:)    ! Full observed vector for all types
   

Generally one can implement any function that computes an observation from the elements of the state vector that is provided to the routine as state_p. The coordinate information is provided in thisobs. One can also modify the interface to the obsration operator routine if, e.g., additional information is required.

Dependent on the complexity of an observation operator it might be useful to separate the functional computations from the interpolation. For this one could consider multi-step observation operators in separate routines.

In the template we also included outputs for the debug functionality. These should be adapted to a particular observation operator so that meaningful outputs are generated that help to check whether the operations in the routine are correct.

In case that you implement a new observation operator: The community of PDAF users would be grateful if you can help to advance PDAF-OMI by providing new observation operators for inclusion into PDAF-OMI under the LGPL license.