Version 24 (modified by 6 weeks ago) (diff)  ,

PDAFOMI Observation Operators
PDAFOMI Guide
 Overview
 callback_obs_pdafomi.F90
 Observation Modules
 Observation Operators
 Checking error status
 Debugging functionality
 Implementing the analysis step with OMI
 Porting an existing implemention to OMI
 Using domainlimited observations
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:
 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 arraythisobs%id_obs_p
initialized ininit_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 arraythisobs%id_obs_p
initialized ininit_dim_obs_OBSTYPE
.  PDAFomi_obs_op_interp_lin
This observation operator computes the observation by linear interpolation. It uses the index arraythisobs%id_obs_p
and the arraythisobs%icoeff_p
holding interpolation coefficients initialized ininit_dim_obs_OBSTYPE
. To use this observation operator, one has to allocate and initializethisobs%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(:) ! Processlocal 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
. PDAFOMI provides three routines for this task:
 PDAFomi_get_interp_coeff_lin1D
Simplified initialization for 1dimensional 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 bilinear interpolation (n_dim=2) only the coordinates for grid points 1, 2, and 3 are used to compute the coefficients
 For trilinear 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, thusgcoords(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 inthisobs%id_obs_p
Adjoint observation operators
For the application of 3DVar, 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 arraythisobs%id_obs_p
initialized ininit_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 arraythisobs%id_obs_p
initialized ininit_dim_obs_OBSTYPE
.  PDAFomi_obs_op_adj_interp_lin
This observation operator is used for the case of linear interpolation. It uses the index arraythisobs%id_obs_p
and the arraythisobs%icoeff_p
holding interpolation coefficients initialized ininit_dim_obs_OBSTYPE
. To use this observation operator, one has to allocate and initializethisobs%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(:) ! Processlocal 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 PDAFOMI 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:
 Check whether
thisobs%doassim==1
, which indicates that the observation is assimilated  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 processlocal domain (for a model without parallelization this is the global domain)  Call
PDAFomi_gather_obsstate
to gather the full observation vectorobs_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 processlocal 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 multistep 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 PDAFOMI by providing new observation operators for inclusion into PDAFOMI under the LGPL license. 