= PDAF-OMI Observation Operators =
{{{
#!html
}}}
[[PageOutline(2-3,Contents of this page)]]
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:
- '''[wiki:PDAFomi_obs_op_gridpoint]'''[[br]]
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`.
- '''[wiki:PDAFomi_obs_op_gridavg]'''[[br]]
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`.
- '''[wiki:PDAFomi_obs_op_interp_lin]'''[[br]]
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'''[[br]]
Simplified initialization for 1-dimensional models
- '''PDAFomi_get_interp_coeff_lin'''[[br]]
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'''[[br]]
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(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(:,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:
- '''[wiki:PDAFomi_obs_op_adj_gridpoint]'''[[br]]
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`.
- '''[wiki:PDAFomi_obs_op_adj_gridavg]'''[[br]]
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`.
- '''[wiki:PDAFomi_obs_op_adj_interp_lin]'''[[br]]
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_X (thisobs,[nrows,] state_p, ostate)
}}}
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 [wiki:OMI_debugging 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.||