Changes between Initial Version and Version 1 of OMI_observation_operators_PDAF3


Ignore:
Timestamp:
May 27, 2025, 3:39:35 PM (6 days ago)
Author:
lnerger
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • OMI_observation_operators_PDAF3

    v1 v1  
     1= PDAF-OMI Observation Operators =
     2
     3{{{
     4#!html
     5<div class="wiki-toc">
     6<h4>PDAF-OMI Guide</h4>
     7<ol><li><a href="PDAF_OMI_Overview_PDAF3">Overview</a></li>
     8<li><a href="OMI_Callback_obs_pdafomi_PDAF3">callback_obs_pdafomi.F90</a></li>
     9<li><a href="OMI_observation_modules_PDAF3">Observation Modules</a></li>
     10<li>Observation operators</li>
     11<li><a href="ImplementationofAnalysisStep_PDAF3">Implementing the analysis step</a></li>
     12<li><a href="OMI_error_checking_PDAF3">Checking error status</a></li>
     13<li><a href="OMI_debugging_PDAF3">Debugging functionality</a></li>
     14<li><a href="nondiagonal_observation_error_covariance_matrices_PDAF3">Using nondiagonal R-matrices</a></<li>
     15<li><a href="OMI_additional_functionality_PDAF3">Additional OMI Functionality</a></li>
     16<li><a href="Porting_to_OMI_PDAF3">Porting an existing implemention to OMI</a></li>
     17</ol>
     18</div>
     19}}}
     20
     21[[PageOutline(2-3,Contents of this page)]]
     22
     23An observation operator routine is called for each observation type in the routine `obs_op_pdafomi` in the file `callback_obs_pdafomi.F90`.
     24
     25== Observation operators ==
     26
     27OMI currently provides 3 observation operators:
     28 - '''[wiki:PDAFomi_obs_op_gridpoint]'''[[br]]
     29  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`.
     30 - '''[wiki:PDAFomi_obs_op_gridavg]'''[[br]]
     31  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`.
     32 - '''[wiki:PDAFomi_obs_op_interp_lin]'''[[br]]
     33  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.
     34
     35The arguments of the observation operators are
     36{{{
     37  SUBROUTINE PDAFomi_obs_op_X (thisobs,[nrows,] state_p, ostate)
     38
     39     TYPE(obs_f), INTENT(inout) :: thisobs  ! Data type with full observation
     40   [ INTEGER, INTENT(in) :: nrows           ! Number of values to be averaged ]
     41     REAL, INTENT(in)    :: state_p(:)      ! Process-local model state (provided by PDAF)
     42     REAL, INTENT(inout) :: obs_f_all(:)    ! Full observed state for all observation types (array provided by PDAF)
     43}}}
     44Where `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.
     45
     46
     47== Initializing interpolation coefficients ==
     48
     49The 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:
     50 - '''PDAFomi_get_interp_coeff_lin1D'''[[br]]
     51 Simplified initialization for 1-dimensional models
     52 - '''PDAFomi_get_interp_coeff_lin'''[[br]]
     53 Determine interpolation coefficients based on the coordinates of grid points and the observation for a rectangular grid in 1, 2, or 3 dimensions.
     54 - '''PDAFomi_get_interp_coeff_tri'''[[br]]
     55 Determine barycentric interpolation coefficients for triangular grids based on the coordinates of grid points and the observation
     56
     57An 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`.
     58
     59
     60=== PDAFomi_get_interp_coeff_lin ===
     61
     62This routine computes interpolation coefficients for a rectangular grid in 1, 2, or 3 dimensions.
     63
     64The call to this routine is
     65{{{
     66  CALL PDAFomi_get_interp_coeff_lin(num_gp, n_dim, gcoords, ocoord, icoeff)
     67
     68    INTEGER, INTENT(in) :: num_gp         ! Length of thisobs%icoeff_p(:,i)
     69    INTEGER, INTENT(in) :: n_dim          ! Number of dimensions in interpolation
     70    REAL, INTENT(in)    :: gcoords(:,:)   ! Coordinates of grid points
     71    REAL, INTENT(in)    :: ocoord(:)      ! Coordinates of observation (one column ocoord_p(:,i))
     72    REAL, INTENT(inout) :: icoeff(:)      ! Interpolation coefficients (one column thisobs%icoeff_p(:,i))
     73
     74}}}
     75Here 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.
     76
     77In 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.
     78
     79`ocoord(:,i)` holds the list of the coordinates for the observation with index i (different from the use in `gcoords`)
     80
     81The order of the coordinates and coefficients is the following:
     82{{{
     83                       (7)------(8)
     84                       /|       /|    with
     85                     (5)+-----(6)|       - column 1
     86                      | |      | |       / column 2
     87                      |(3)-----+(4)      | column 3
     88                      |/       |/
     89                     (1) ---- (2)
     90
     91   thus gcoords(1,1)/=gcoords(2,1),        but  gcoords(1,1)=gcoords(3,1)=gcoords(5,1),
     92        gcoords(1,2)/=gcoords(3,2),             gcoords(1,2)=gcoords(2,2)=gcoords(5,2),
     93        gcoords(1,3)/=gcoords(5,3)              gcoords(1,3)=gcoords(2,3)=gcoords(3,3)
     94}}}
     95
     96'''Notes:'''
     97 - For 1D linear interpolation (n_dim=1) only the coordinates for grid points 1 and 2 are used to compute the coefficients
     98 - For bi-linear interpolation (n_dim=2) only the coordinates for grid points 1, 2, and 3 are used to compute the coefficients
     99 - For tri-linear interpolation (n_dim=3) only the coordinates for grid points 1, 2, 3, and 5 are used to compute the coefficients
     100
     101
     102=== PDAFomi_get_interp_coeff_lin1D ===
     103
     104This routine is a simplified variant of `PDAFomi_get_interp_coeff_lin`. It computes interpolation coefficients in 1 dimension only
     105
     106The call to this routine is
     107{{{
     108  CALL PDAFomi_get_interp_coeff_lin1D(gcoords, ocoord, icoeff)
     109
     110    REAL, INTENT(in)    :: gcoords(:)  ! Coordinates of grid points; dim=2
     111    REAL, INTENT(in)    :: ocoord      ! Coordinates of observation
     112    REAL, INTENT(inout) :: icoeff(:)   ! Interpolation coefficients; dim=2
     113}}}
     114
     115
     116
     117=== PDAFomi_get_interp_coeff_tri ===
     118
     119For triangular model grid interpolation coefficients are determined as barycentric coordinates. This is performed in this routine.
     120
     121The call to this routine is
     122{{{
     123  CALL PDAFomi_get_interp_coeff_tri(gcoords, ocoord, icoeff)
     124
     125    REAL, INTENT(in)    :: gcoords(:)      ! Coordinates of grid points; dim=(3,2)
     126    REAL, INTENT(in)    :: ocoord(:)       ! Coordinates of observation; dim=2
     127    REAL, INTENT(inout) :: icoeff(:)       ! Interpolation coefficients; dim=3
     128}}}
     129
     130'''Notes:'''
     131 - 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.
     132 - The order of the grid points in `gcoords` has to be consistent with the order of the indices specified in `thisobs%id_obs_p`
     133
     134
     135== Adjoint observation operators ==
     136
     137For the application of 3D-Var, adjoint observation operators are required. These perform the operation of the transposed linear observation operator.
     138
     139OMI provides the adjoint observation operators corresponding to the forward observation operators:
     140 - '''[wiki:PDAFomi_obs_op_adj_gridpoint]'''[[br]]
     141  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`.
     142 - '''[wiki:PDAFomi_obs_op_adj_gridavg]'''[[br]]
     143  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`.
     144 - '''[wiki:PDAFomi_obs_op_adj_interp_lin]'''[[br]]
     145  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.
     146
     147The arguments of the observation operators are
     148{{{
     149  CALL PDAFomi_obs_op_adj_X (thisobs,[nrows,] state_p, ostate)
     150
     151     TYPE(obs_f), INTENT(inout) :: thisobs  ! Data type with full observation
     152     INTEGER, INTENT(in) :: nrows           ! Number of values to be averaged
     153     REAL, INTENT(in)    :: obs_f_all(:)    ! Full observed state for all observation types (array provided by PDAF)
     154     REAL, INTENT(inout) :: state_p(:)      ! Process-local model state (provided by PDAF)
     155}}}
     156Where `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.
     157
     158
     159== Implementing your own observation operator ==
     160
     161The 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.
     162
     163|| 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. ||
     164
     165Each observation operator include the following functionality:
     166 1. Check whether `thisobs%doassim==1`, which indicates that the observation is assimilated
     167 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)
     168 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:
     169{{{
     170  SUBROUTINE PDAFomi_gather_obsstate(thisobs, ostate_p, obs_f_all)
     171
     172    TYPE(obs_f), INTENT(inout) :: thisobs  ! Data type with full observation
     173    REAL, INTENT(in) :: ostate_p(:)        ! Vector of process-local observed state
     174    REAL, INTENT(inout) :: obs_f_all(:)    ! Full observed vector for all types
     175}}}
     176
     177Generally 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.
     178
     179Dependent 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.
     180
     181In the template we also included outputs for the [wiki:OMI_debugging_PDAF3 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.
     182
     183||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.||