= PDAF-OMI Observation Modules = 

{{{
#!html
<div class="wiki-toc">
<h4>PDAF-OMI Guide</h4>
<ol><li><a href="PDAF_OMI_Overview">Overview</a></li>
<li><a href="OMI_Callback_obs_pdafomi">callback_obs_pdafomi.F90</a></li>
<li>Observation Modules</li>
<li><a href="OMI_observation_operators">Observation operators</a></li>
<li><a href="OMI_debugging">Debugging functionality</a></li>
<li><a href="OMI_ImplementationofAnalysisStep">Implementing the analysis step with OMI</a></li>
<ol>
<li><a href="ImplementAnalysisGlobal">Implementation for Global Filters</a></li>
<li><a href="ImplementAnalysisLocal">Implementation for Local Filters</a></li>
<li><a href="ImplementAnalysislenkfOmi">Implementation for LEnKF</a></li>
</ol>
<li><a href="Porting_to_OMI">Porting an existing implemention to OMI</a></li>
<li><a href="OMI_use_global_obs">Using domain-limited observations</a></li>
</ol>
</div>
}}}

[[PageOutline(2-3,Contents of this page)]]

== Overview ==

The implementation of the observations with OMI is done in observation modules (obs-modules). For each observation type a separate module should be created. 

Each obs-module contains four routines (where 'TYPE' will be replaced by the name of the observation):

 - `init_dim_obs_TYPE` initializes all variables holding the information about one observation type. The information about the observation type is stored in a data structure (Fortran derived type).
 - `obs_op_TYPE` applies the observation operator to a state vector. One can call an observation operator routine provided by PDAF, or one can to implement a new operator.
 - `init_dim_obs_l_TYPE` calls a PDAF-OMI routine to initialize the observation information corresponding to a local analysis domain. One can set localization parameters, like the localization radius, for each observation type.
 - `localize_covar_TYPE` calls a PDAF-OMI routine to apply covariance localization. One can set localization parameters, like the localization radius, for each observation type. 

The template file `obs_TYPE_pdafomi_TEMPLATE.F90` shows the different steps needed when implementing these routines. The main work is to implement `init_dim_obs`, while the other routines mainly call a subroutine provided by PDAF-OMI. 

In the obs-module the subroutines are named according to the observation type. The template file uses generic names which can be replaced by the user. Having distinct names for each observation type is relevant to include the subroutine from the module in the call-back routine with ‘use’. In the header of each obs-module, the user can declare further variables, e.g. assim_TYPE as a flag to control whether the observation type should be assimilated.

'''Note:''' In contrast to the 'classical' implementation of observation routines for PDAF, the global and local filters use the same routines `init_dim_obs` and `obs_op`. PDAF-OMI recognizes whether a global or local filter is used and does the necessary operations by itself.

== Data type obs_f == 

To ensure the functionality within each obs-module, we rely on a derived data type called `obs_f` that contains all information about the observation. One instance of this data type is allocated in each obs-module with the generic variable name `thisobs`. A few of the elements of `obs_f` are initialized by the user when the observation information is initialized on `init_dim_obs_f`. Further variables is set in a call to the routine `PDAFomi_gather_obs`. This information is then used by all other routines in the obs-module. The template file obs_TYPE_pdafomi_TEMPLATE.F90 shows the different steps needed to initialize thisobs. 

The '''mandatory variables''' in `obs_f` that need to be set by the user are:
{{{
  TYPE obs_f
     ! ---- Mandatory variables to be set in INIT_DIM_OBS ----
     INTEGER :: doassim=0                 !< Whether to assimilate this observation type
     INTEGER :: disttype                  !< Type of distance computation to use for localization
                                          !   (0) Cartesian, (1) Cartesian periodic
                                          !   (2) simplified geographic, (3) geographic haversine function
     INTEGER :: ncoord                    !< Number of coordinates use for distance computation
     INTEGER, ALLOCATABLE :: id_obs_p(:,:) !< Indices of process-local observed field in state vector
     ...
  END TYPE obs_f
}}}

In addition there are '''optional variables''' that the be used:
{{{
  TYPE obs_f
     ...
     ! ---- Optional variables - they can be set in INIT_DIM_OBS ----
     REAL, ALLOCATABLE :: icoeff_p(:,:)   !< Interpolation coefficients for obs. operator (optional)
     REAL, ALLOCATABLE :: domainsize(:)   !< Size of domain for periodicity (<=0 for no periodicity) (optional)

     ! ---- Variables with predefined values - they can be changed in INIT_DIM_OBS  ----
     INTEGER :: obs_err_type=0            !< Type of observation error: (0) Gauss, (1) Laplace
     INTEGER :: use_global_obs=1          !< Whether to use (1) global full obs. 
                                          !< or (0) obs. restricted to those relevant for a process domain
     ...
  END TYPE obs_f
}}}

Apart from these variables, there is a number of variables that are set internally when the routine `PDAFomi_gather_obs` is called. 


Next to the derived data type `obs_f`, there is a derived type `obs_l` for localization. This is only used internally. It will be filled in init_dim_obs_l when calling PDAFomi_init_dim_obs_l.

== `init_dim_obs_TYPE` ==

This is the main routine to initialize observation information. 

||= Please see the template file `templates/omi/obs_TYPE_pdafomi_TEMPLATE.F90`[[br]] for a step-by-step description of the implementation steps. =||

Each observation module uses the generic name '''thisobs''' for the variable with observation type `obs_f`. Elements of `thisobs` are accessed like
`thisobs%doassim`. 

The main variables that the filled in this routine are
 1. [wiki:OMI_observation_modules#thisobsdoassim thisobs%doassim]: Specify whether this observation type is assimilated
 1. [wiki:OMI_observation_modules#thisobsdisttype thisobs%disttype]: Specify the type of distance computation
 1. [wiki:OMI_observation_modules#thisobsncoord thisobs%ncoord]: Specify the number of dimensions used to compute distances 
 1. [wiki:OMI_observation_modules#dim_obs_p dim_obs_p]: Count the number of available observations
 1. [wiki:OMI_observation_modules#obs_p obs_p]: Fill the vector of observations
 1. [wiki:OMI_observation_modules#ocoord_p ocoord_p]: store the coordinates of the observations
 1. [wiki:OMI_observation_modules#ivar_obs_p ivar_obs_p]: store the inverse error variance of each observation
 1. [wiki:OMI_observation_modules#thisobsid_obs_p thisobs%id_obs_p]: store the indices of state vector elements that correspond to an observation (A single value for observation at grid points, or multiple values for derived quantities or interpolation)

When the observation operator performs interpolation, one further needs to initialize an array of interpolation coefficients ([wiki:OMI_observation_modules#thisobsicoeff_p thisobs%icoeff_p]). For Cartesian distance computation with periodicity one also needs to set [wiki:OMI_observation_modules#thisobsdomainsize thisobs%domainsize].

When parallel model with domain decomposition is used, the variables with suffix `_p` need to describe the observation information for a particular process domain. The following routine will perform the necessary operations to ensure that the parallelization is taken into account by PDAF.

After these variables are filled, one calls
{{{
    CALL PDAFomi_gather_obs(thisobs, dim_obs_p, obs_p, ivar_obs_p, ocoord_p, &
         thisobs%ncoord, local_range, dim_obs)
}}}
This routine will complete all required initializations for OMI. As such it is mandatory to call the routine. 

The routine `PDAFomi_gather_obs` returns the number of observations `dim_obs` which is the return variable for PDAF.  

Note: The value is `local_range` is only used if [wiki:OMI_observation_modules#thisobsuse_global_obs thisobs%use_global_obs=0].


== `obs_op_TYPE` ==

This routine applies the observation operator to a state vector. It returns the observed state vector to PDAF. The routine is used by all filters.

PDAF-OMI provides several observation operators. For example the observation operator for observations that are grid point values is called as:
{{{
    CALL PDAFomi_obs_op_gridpoint(thisobs, state_p, ostate)
}}}

Here, `state_p` is the state vector and `ostate` is the observed state vector.

For more information on the available observation operator and on how to implement your own observation operator see the [wiki:OMI_Observation_Operators documentation of observation operators for OMI].


== `init_dim_obs_l_TYPE` ==

This routine initializes local observation information. The routine is only used by the domain-localized filters (LESTKF, LETKF, LSEIK, LNETF).

For the initialization the following routine is called:
{{{
    CALL PDAFomi_init_dim_obs_l(thisobs_l, thisobs, coords_l, &
         locweight, local_range, srange, dim_obs_l)
}}}

Here, `thisobs` and `thisobs_l` are the data-type variables `obs_f` and `obs_l`. `dim_obs_l`, the local size of the observation vector, is the direct output of the routine.

''Implementation steps:''
 - Ensure that `coords_l` is filled in `init_dim_l_pdaf`
 - Specify the localization variables (These variables are usually set in `init_pdaf` and included with `use mod_assimilation`)
   - `locweight`: Type of localization (see init_pdaf)
   - `local_range`: The localization weight radius
   - `srange`: The support radius of the localization
 


== `localize_covar_TYPE` ==

This routine initializes local observation information. The routine is only used by the local EnKF (LEnKF).

For the initialization the following routine is called:
{{{
    CALL PDAFomi_localize_covar(thisobs, dim_p, locweight, local_range, srange, &
         coords_p, HP_p, HPH)
}}}

Here, `thisobs` is the data-type variable `obs_f`. `HP_p` and `HPH` are the covariance matrices projected onto the observations. The localization will be applied to these variables.

''Implementation steps:''
 - Ensure that `coords_p` is filled in `localize_covar_pdafomi`
 - Specify the localization variables (These variables are usually set in `init_pdaf` and included with `use mod_assimilation`)
   - `locweight`: Type of localization (see init_pdaf)
   - `local_range`: The localization weight radius
   - `srange`: The support radius of the localization

== Implementing a new observation type ==

To implement a new observation type, the approach is generally as follows:
1.	Create a copy of obs_TYPE_pdafomi_TEMPLATE.F90
1.	Rename the module and its subroutines according to the observation (replacing ‘TYPE’ by name of observation).
1.	Implement init_dim_obs for the observation type following the instructions in the template
1.	Adapt obs_obs for the observation type
1.	Adapt init_dim_obs_l for the observation type (if using a domain_localized filter)
1.	Adapt localize_covar for the observation type (if using a the local EnKF)
1.	Add subroutine calls for the new observation type into the routines in callback_obs_pdafomi.F90


== Implementation hints for init_dim_obs ==

=== `thisobs%doassim` ===

Set this variable to 1 to let the filter assimilate this observation. The setting is usually conditional on the value of `assim_TYPE` which is set in `init_pdaf`: 
{{{
   IF (assim_TYPE) thisobs%doassim = 1
}}}

=== `thisobs%ncoord` ===

This variable specifies the dimension of the distnace computations. Thus thisobs%ncoord=2 will lead to distance computations in 2 dimensions.


=== `thisobs%disttype` ===

This variable specifies the type of distance computation. Possible choices are
 - 0: Cartesian distance in ncoord dimension
 - 1: Cartesian distance in ncoord dimensions with periodicity (Needs specification of thisobs%domainsize(ncoord))
 - 2: Approximate geographic distance in meters with horizontal coordinates in radians (latitude: -pi/2 to +pi/2; longitude -pi/+pi or 0 to 2pi)
 - 3: Geographic distance computation in meters using haversine formula with horizontal coordinates in radians (latitude: -pi/2 to +pi/2; longitude -pi/+pi or 0 to 2pi)

For 0 and 1 any distance unit can be used. The computed distance will be in the same unit. For 2 and 3 the input coordinates are in radians and the distance is computed in meters.

See `/models/lorenz96_omi/` for an example using case 1 with periodicity in one dimension.


=== `dim_obs_p` ===

This is a single integer value giving the number of observations. With a parallel model using domain-decomposition this will be the number of observations  for the process sub-domain. For observation files holding all observations one can read these and then check which observation redice within the process sub-domain. `dim_obs_p` will be used to allocate further arrays and as input argument to `PDAFomi_gather_obs`.

=== `obs_p` ===

This should be a vector of real values. It will be used as an argument to `PDAFomi_gather_obs`. The order of the entries has to be consistent in the arrays `thisobs%id_obs_p`, `obs_p`, `ivar_obs_p`, and `ocoord_p`.


=== `ocoord_p` ===

This should be a rank-2 array of real values with size (thisobs%ncoord, dim_obs_p). It will be used as an argument to `PDAFomi_gather_obs`. The order of the entries has to be consistent in the arrays `thisobs%id_obs_p`, `obs_p`, `ivar_obs_p`, and `ocoord_p`.

The coordinates of the observation with index `k` are given by `ocoord_p(:,k)`. 


=== `ivar_obs_p` ===

This should be a vector of real values. It will be used as an argument to `PDAFomi_gather_obs`. The order of the entries has to be consistent in the arrays `thisobs%id_obs_p`, `obs_p`, `ivar_obs_p`, and `ocoord_p`.


=== `thisobs%id_obs_p` ===

This array is allocated as
{{{
   ALLOCATE(thisobs%id_obs_p(NROWS, dim_obs_p))
}}}
For a fixed value of the second index the NROWS are the indices of the elements of the state vector that are treated in the observation operator.
The value of NROWS depends on the observation operator used for an observation. Examples:
 - Using observations that are grid points values:
   - NROWS=1
   - The entry is the index of a single element of the state vector
 - Using observations that are determined by bi-linear interpolation of 4 grid points: 
   - NROWS=4
   - The entries are the indices of four elements of the state vector


=== `thisobs%domainsize` ===

This array has to be allocated as
{{{
   ALLOCATE(thisobs%domainsize(thisobs%ncoord))
}}}
Here one has to specify the size of the domain in each of its thisobs%ncoord dimensions. The information is used to compute the Cartesian distance with periodicity.

Setting one dimension to 0 or a negative value indicates that there is no periodicity in this direction. 


=== `thisobs%icoeff_p` ===

This array is allocate the in same way as `thisobs%id_obs_p`:

{{{
   ALLOCATE(thisobs%icoeff_p(NROWS, dim_obs_p))
}}}
The value of NROWS has to be the same as for `thisobs%id_obs_p`. For a fixed value of the second index the NROWS of the array hold the interpolation coefficients corresponding to the indices specified in `thisobs%id_obs_p`.

Please see the [wiki:OMI_observation_operators documentation of OMI observation operators] for information on how to initialize the array `thisobs%icoeff_p` using functions provided by PDAF-OMI.


=== `thisobs%obs_err_type` ===

The particle filter methods NETF, LNETF and PF can handle observations with non-Gaussian errors. PDAF-OMI supports the following two choices:
 - 0: Gaussian errors (''default value'')
 - 1: double-exponential (Laplace) errors


=== `thisobs%use_global_obs` ===

In the domain-localized filters (LESTK, LETKF, LSEIK, LNETF) observations are assimilated that are located within the localization around some grid point. When a model uses parallelization with domain-decomposition some of these observations might belong to a different process-domain. In the default mode (`thisobs%use_global_obs`=1) PDAF-OMI gathers all globally available observations so that each process has access to all observations. It can be more efficient to limit the observations on a process-domain to those observations that are located inside the domain or within the localization radius around it. Then, in the local analyses less observations have to be checked for their distance. Setting `thisobs%use_global_obs=0` activates this feature. However, it needs additional preparations to make PDAF-OMI aware of the limiting coordinates of a process sub-domain. 

The use of this feature is described in the [wiki:OMI_use_global_obs documentation on using domain-limited observations].