Developer Guide: Adding a New Observation Type to the OMI Interface

Developer Guide: Adding a New Observation Type to the OMI Interface#

This page describes how to add a new observation type — generically called C below — to the TSMP-PDAF OMI (Observation Model Interface) framework. It covers the architecture, the role of each source file, and the changes required, using the existing GRACE-DA and SM-DA implementations as reference.

Architecture Overview#

Observation handling sits between the PDAF library and the eCLM model state in three layers:

PDAF library (src/)
    │  calls callback routines at each analysis step
    ▼
callback_obs_pdafomi.F90        ← single routing hub
    │  dispatches to one routine per observation type
    ▼
obs_GRACE_pdafomi.F90           ← one module per observation type
obs_SM_pdafomi.F90
obs_C_pdafomi.F90  (new)
    │  read observations, build H, apply H to state vector
    ▼
eCLM state vector (via enkf_clm_mod, mod_assimilation, …)

Each observation-type module is self-contained. The only coupling between modules is through the central routing file and the shared mod_assimilation module.

Files Involved#

To add a new observation type C, touch exactly eight things:

File

Action

interface/framework/obs_C_pdafomi.F90

Create (new file, modelled on the template)

interface/framework/Makefile

Add obs_C_pdafomi.o to MOD_USER_PDAFOMI

interface/framework/callback_obs_pdafomi.F90

Add calls for C to every callback routine

interface/framework/mod_assimilation.F90

Declare cradius_C, sradius_C, and any other C-specific parameters

interface/framework/mod_read_obs.F90

Add 'C' case to update_obs_type

interface/model/eclm/enkf_clm_mod_5.F90

Add clmupdate_C flag (parallel to clmupdate_tws, clmupdate_swc)

interface/framework/init_pdaf.F90

Import assim_C and set it from clmupdate_C

interface/framework/init_pdaf_parse.F90

Import rms_obs_C and add a parse call for it

The template at templates/omi/obs_OBSTYPE_pdafomi_TEMPLATE.F90 is a good starting skeleton. The existing SM and GRACE modules show how the template is applied in this codebase.

Step 1 — Create obs_C_pdafomi.F90#

Module-level variables#

At minimum, declare assim_C (a logical switch controlled by callback_obs_pdafomi), rms_obs_C (a constant fallback error standard deviation), and the PDAFomi types thisobs (of type obs_f, thread-safe) and thisobs_l (of type obs_l, declared THREADPRIVATE).

assim_C and rms_obs_C are wired up externally:

  • init_pdaf.F90 imports assim_C and sets it from clmupdate_C (i.e. assim_C = (clmupdate_C /= 0)), so the observation type is enabled automatically when the corresponding model-state update flag is active. Commented-out placeholders mark the insertion point.

  • init_pdaf_parse.F90 imports rms_obs_C and populates it via a parse call that reads rms_obs_C from enkfpf.par. Commented-out placeholders mark the insertion point here as well.

Additional module-level allocatable arrays — e.g. a cached climatological mean needed for an anomaly operator — follow the same pattern as in obs_GRACE_pdafomi.F90.

The obs_f type#

obs_f carries all observation metadata and is populated inside init_dim_obs_C. Mandatory fields:

Field

Type

Meaning

thisobs%doassim

INTEGER

1 = assimilate, 0 = skip

thisobs%disttype

INTEGER

distance metric for localization (see below)

thisobs%ncoord

INTEGER

number of spatial coordinates (usually 2)

thisobs%id_obs_p(:,:)

INTEGER, ALLOCATABLE

state-vector index for each process-local observation

Frequently useful optional fields:

  • icoeff_p: bilinear interpolation weights

  • obs_err_type: 0 = Gaussian, 1 = Laplace

  • inno_omit: innovation outlier rejection threshold

disttype choices:

  • 0 — Cartesian distance; coordinates in any consistent unit (GRACE uses integer grid-cell indices, so cradius_GRACE is in number of grid cells).

  • 3 — geographic distance; coordinates in degrees, radii in km (used by SM).

Use disttype=3 for point observations distributed geographically (C, SM) and disttype=0 for products already mapped to grid-cell indices (e.g. GRACE).

Routine 1 — init_dim_obs_C#

Called once per analysis step before the filter loop. Responsibilities:

  1. Set the mandatory thisobs fields.

  2. Read the global observation vector from the NetCDF file using read_obs_nc_type from mod_read_obs.

  3. Handle the dim_obs == 0 case (no observations scheduled for this step): allocate dummy arrays of size 1 and call PDAFomi_gather_obs before returning.

  4. For each global observation that falls within the PE-local domain, record the state-vector index in thisobs%id_obs_p and populate the local arrays obs_p, ivar_obs_p, and ocoord_p.

  5. Call PDAFomi_gather_obs to finalize observation bookkeeping across MPI ranks.

Reading observations: Both SM and GRACE use read_obs_nc_type, which checks the type_clm NetCDF variable and returns dim_obs = 0 if the file does not contain observations of the requested type.

Who reads the file: SM reads only on mype_filter==0 and broadcasts; GRACE reads on all filter processes because the spatial averaging operator needs each PE to know its contribution. For point observations like C, the SM pattern (read + broadcast) is simpler.

Snapping observations to the state vector: For point observations (SM, C), loop over all global observations and all PE-local gridcells and record the state-vector index for matching pairs. For GRACE the logic is inverted: each observation aggregates multiple gridcells, so the loop iterates over gridcells and maps each one to the observation it contributes to.

Routine 2 — obs_op_C#

Applies the observation operator H, mapping the PE-local model state to the observation space. For a simple point observation (SM pattern): allocate a local output array, extract state_p(obs_index_p(i)) for each local observation, then call PDAFomi_gather_obsstate to assemble the global result in ostate.

For aggregating operators (GRACE pattern), sum all gridcell contributions mapping to the same observation, apply any transformation (e.g. subtract temporal mean for TWS anomaly), and then call PDAFomi_gather_obsstate. Follow the SM pattern for C unless your observations represent spatial averages over a large footprint.

Routines 3 & 4 — Localization (required for LESTKF/LETKF)#

init_dim_obs_l_C counts observations within the localization radius of the current local analysis domain. It is a thin wrapper around PDAFomi_init_dim_obs_l, passing thisobs_l, thisobs, the local domain coordinates coords_l (in the same units and disttype as ocoord_p), and the localization parameters cradius_C and sradius_C.

localize_covar_C is required only for the LEnKF. It is an equally thin wrapper around PDAFomi_localize_covar with the same localization parameters.

Step 2 — Update callback_obs_pdafomi.F90#

This file dispatches PDAF callbacks to each observation module. Add C in the same pattern as GRACE and SM. The commented-out !USE obs_C_pdafomi placeholders show exactly where to insert new entries.

For each of the eight callback subroutines, add a USE obs_C_pdafomi statement and one call to the corresponding C routine:

  • init_dim_obs_pdafomi: set assim_C = .true., call init_dim_obs_C, and add dim_obs_C to the running total. The assim_* flags exist to disable an observation type for debugging; keep assim_C = .true. here because the module itself returns dim_obs = 0 when no observations are present. (This logic is planned to be complemented in the future by setting assim_C from the EnKF input file)

  • obs_op_pdafomi: call obs_op_C. Call order does not affect correctness because offsets in ostate are determined by the order of init_dim_obs_* calls above.

  • init_dim_obs_l_pdafomi: call init_dim_obs_l_C. The variable dim_obs_l is incremented inside each call; PDAF initializes it to 0 before invoking the callback.

  • localize_covar_pdafomi: call localize_covar_C.

  • Remaining callbacks (add_obs_err, init_obscovar, prodRinvA, prodRinvA_l, deallocate_obs): each needs a corresponding C call, all of which are thin wrappers around the matching PDAFomi generic routine, identical in structure to the existing GRACE and SM entries.

Step 3 — Add Parameters to mod_assimilation.F90#

Declare cradius_C and sradius_C next to the existing GRACE and SM radius parameters. These are read from the enkfpf.par input file via the existing namelist mechanism; add corresponding entries to the [DA] section parser in the same module.

Step 4 — Register the New Type in mod_read_obs.F90#

The routine update_obs_type reads the type_clm string from the observation file and sets the clmupdate_* flags for the next analysis cycle. Add a 'C' case that sets the appropriate flags. If no existing clmupdate_* flag covers your new variable, declare a new clmupdate_C in enkf_clm_mod.F90 and wire it to the state-vector fill and scatter logic in the eCLM interface routines (eclm/).

Key Design Differences Between GRACE-DA and SM-DA#

Localization distance units#

Type

disttype

Coordinate arrays

Radius unit

GRACE

0 (Cartesian)

integer grid-cell indices

grid cells

SM

3 (haversine)

degrees longitude/latitude

km

Observation operator complexity#

  • SM: direct extraction from state_p using a precomputed index.

  • GRACE: aggregating — all gridcell values within the GRACE footprint are summed and averaged, then the temporal mean is subtracted to produce a TWS anomaly. thisobs%id_obs_p(1, g) maps each gridcell g to its observation (inverse of the usual convention).

Follow the SM pattern for C point observations; follow the GRACE pattern if observations represent spatial averages over a large footprint.

Coverage checks and error specification#

GRACE discards observations where fewer than a minimum number of gridcells lie within the support radius (important near coastlines). SM and C point observations do not need this, but consider rejecting observations over water bodies or outside the CLM domain.

Error specification is controlled by multierr: constant (multierr=0), per-observation from obserr_clm (multierr=1), or full covariance matrix from obscov_clm (multierr=2).

Assimilating Multiple Observation Types at the Same Timestep#

The framework generates a state vector for each type individually before the assimilation, some things would need to be adapted when mutliple observation types are assimilated at the same timestep. Currently, one observation file only consists of one observation type. As SM observations are usually assimilated at noon and GRACE observations are assimilated at the end of the month at midnight, this should not provide any problems.