Column-based forward operators in DART

[hkershaw-brown]

Column-based forward operators in DART

Current state of https://github.com/NCAR/DART/blob/main/observations/forward_operators/obs_def_*_mod.f90 for forward
operators (FOs) that need

• a whole vertical column of model state
  or
• b the number of vertical levels in the model, obs_defs get this by calling
  model_interpolate repeatedly, one level and one quantity at a time.

As of DART v11.25.
Projects:
DART The Need for Speed #995.
DART Column Operators #366 #754 #438.

---

1. The two idioms in use today

1a. "Count levels by interpolating until it fails"

There is no model_mod interface that reports the number of vertical levels at a
location, so operators discover it by probing VERTISLEVEL 1, 2, 3, … until
model_interpolate returns a non-zero istatus:

Module	Loop	Probe quantity	Cap
obs_def_rttov13_mod.f90:3540	COUNTLEVELS	QTY_PRESSURE	maxlevels = 10000
obs_def_rttov_mod.f90:3539	COUNTLEVELS	QTY_PRESSURE	maxlevels = 10000
obs_def_COSMOS_mod.f90:510	COUNTLEVELS	QTY_GEOPOTENTIAL_HEIGHT	maxlayers = 1000
obs_def_cice_mod.f90:426	CATEGORIES	QTY_SEAICE_CATEGORY	MAX_CATEGORIES = 100
obs_def_cwp_mod.f90:491	LEVELS	QTY_PRESSURE	10000 (literal)
obs_def_tpw_mod.f90:253	LEVELS	QTY_PRESSURE	10000 (literal)
obs_def_CO_Nadir_mod.f90:296, :458	model_levels	QTY_PRESSURE	unbounded do
obs_def_upper_atm_mod.f90:302, :651, :677	LEVELS, FILLINTERFACES, FILLMIDPOINTS	QTY_DENSITY_ION_E, QTY_GEOMETRIC_HEIGHT, QTY_PRESSURE	MAXLEVELS = 300 (compile-time)

Consequences of the idiom itself:

• A deliberate failure per column. Every column costs one extra
  model_interpolate that is expected to fail. Some model_mods log/print on
  failure, and some fail with E_ERR on the quantity being probed — e.g.
  tiegcm/model_mod.f90:308 hard-errors on
  QTY_GEOPOTENTIAL_HEIGHT, which is exactly what
  get_expected_gnd_gps_vtec/get_expected_O_N2_ratio want to probe. The comment
  there says as much.
• The model side has been written to server this idiom. CLM has an explicit early-return
  for the probe: "Some applications just need to know the number of vertical
  levels. This is done by trying to 'interpolate' height on a large number of
  levels. When the interpolation fails, you've gone one level too far."
  (clm/model_mod.f90:1237).
• Compile-time level caps. obs_def_upper_atm_mod requires preprocess +
  recompile if the model has more than 300 levels
  (:308).
• Fragile probe locations. obs_def_cice_mod counts categories at a
  hard-wired (0.0, 0.0) (:419-423),
  which is open ocean in the Gulf of Guinea. It works only because the
  quantity exists on the grid there for global models (not true for regional models?)
• Not always cached. RTTOV caches the count once globally
  (arrays_prealloced), and obs_def_cice_mod caches Ncat once
  (:167). COSMOS,
  cwp, tpw, CO_Nadir and upper_atm re-count on every observation.

1b. "Build the column one model_interpolate call per (level, quantity)"

Then the column data is fetched the same way. Worst offender, RTTOV
GETLEVELDATA (obs_def_rttov13_mod.f90:3603-3822):
42 possible interpolate calls inside the level loop, plus ~16 surface calls
after it. With a 55-level model and clouds enabled that's roughly
55 × 9 ≈ 500 calls per observation; with the aerosol options on it approaches
55 × 42 ≈ 2300.

---

2. Affected forward operators, ranked by cost

Tier 1 — full multi-quantity column per observation

Operator	Routine	Calls per obs (order of)
RTTOV v12 / v13 radiances — obs_def_rttov_mod, obs_def_rttov13_mod	get_expected_radiance	nlev × (3 … 42) + ~16 surface
Cloud water path — obs_def_cwp_mod	get_expected_cwp	nlev × (1 … 6) + 7 surface
IASI / MOPITT CO — obs_def_CO_Nadir_mod	get_expected_iasi_co, get_expected_mopitt_co	nlev (pressure) + n_retrieval_layers × nlev for QTY_CO
TES nadir radiance (Mars) — obs_def_TES_nadir_mod	get_expected_TES_nadir_obs	ens_size × N_layers temperature calls
O/N₂ ratio — obs_def_upper_atm_mod	get_expected_O_N2_ratio	nlev × 5 (interfaces + 4 midpoint quantities)
Ground GPS VTEC — obs_def_upper_atm_mod	get_expected_gnd_gps_vtec	nlev × 2
Total precipitable water — obs_def_tpw_mod	get_expected_tpw	nlev × 2 + 2 surface
COSMOS neutron intensity — obs_def_COSMOS_mod	get_expected_neutron_intensity	nlayer × 3 (count + height + moisture)
Sea ice aggregates — obs_def_cice_mod	4 get_expected_agreg_*	Ncat × (1 … 2), ×2 for freeboard/thickness

Two of these are particularly bad

• obs_def_CO_Nadir_mod — the retrieval-layer loop do i=start_i,end_i wraps
  the model-level loop do j=1,num_levs, and re-interpolates QTY_CO at model
  level j inside both
  (:325-352).
  With mopitt_dim = 10 and 32 CAM levels that is ~320 QTY_CO interpolations
  where 32 would do — each model level is fetched up to 10 times.
• obs_def_TES_nadir_mod — do imem = 1, ens_size / do k = 1, N_layers calls
  interpolate with the full ens_size each time
  (:874-885),
  so with 80 members and N_layers = 30 that is 2400 full-ensemble interpolations
  to fill 30 levels. The @todo right above it names the problem: "Possibly we
  could overload interpolate and set_location to take an array of locations …
  For now this is a loop around ensemble_size calling interpolate —
  overcommunicating."

Tier 2 — many quantities at one location (repeated horizontal search, no column)

These don't need a column, but they pay the horizontal-location cost N times for
one point, so they benefit from the same fix (locate once, read many):

Operator	Quantities at one location
obs_def_radar_mod — get_expected_radar_ref, get_expected_fall_velocity	up to 5
obs_def_upper_atm_mod — get_expected_upper_atm_density, get_expected_oxygen_ion_density	4–6
obs_def_oxygen_ion_density_mod — get_expected_oxygen_ion_val	7
obs_def_land_mod — albedo / biomass / PAR operators	3–4
obs_def_GRACE_mod — get_expected_TWS	6 (4 hard-coded soil depths + snow + aquifer)
obs_def_rel_humidity_mod, obs_def_abs_humidity_mod, obs_def_dew_point_mod, obs_def_altimeter_mod, obs_def_wind_speed_mod	2–3

obs_def_GRACE_mod is a Tier-1 case in waiting: it hard-codes soil depths
0.05/0.25/0.7/1.5 m and thicknesses 100/300/600/1000 mm with
!>@todo Determine number of layers and their thicknesses automatically
(:171) — i.e.
it wants a column interface and currently fakes it. (It also E_ERRs
unconditionally at :168.)

Tier 3 — ray/path integration (many different locations)

obs_def_gps_mod — get_expected_gpsro_ref, non-local operator
(:481-530). The
INTEGRATE loop steps along the ray and calls ref_local twice per step; each
ref_local does 3 interpolate calls (T, Q, P) at the same new lat/lon/height
(:599-607).

This is not a column caching problem,
But the 3 calls at each step are at an identical location, so each step
triggers 3 identical vertical-column constructions in the model (see 3). A
"give me these N quantities at this one location" interface fixes 2/3 of the work,
and a location-array interface would let the whole ray go in one call.
Q. Is a location array and a column array the same thing (interface wise?)

---

3. Which models are expensive?

The per-call cost splits into a horizontal locate and a vertical locate.
Repeating model_interpolate down a column repeats both, even though only the
level index changed.

3a. Horizontal search repeated per level

MPAS-A — models/mpas_atm/model_mod.f90. Worst case. Note users reporting slowness with this: #995
model_interpolate calls cell_ok_to_interpolate
(:1179) → find_closest_cell_center
→ xyz_find_nearest against a get_close structure over all nCells
(:6380-6417). find_triangle then
calls find_closest_cell_center again (:5920).
So an unstructured nearest-neighbour search runs ≥2× per interpolate call. For an
RTTOV column that is ~1000 identical nearest-cell searches for one lat/lon.
highest_obs_pressure_mb > 0 adds a compute_pressure_at_loc on top, which the
code itself flags: "this is expensive - only do it if users want to reject
observations at the top of the model" (:1266).

CAM-SE — models/cam-se/model_mod.f90.
coord_ind_cs searches the cubed-sphere for the enclosing quad on every call
(:592). The code already knows this is
redundant: "SENote2: In the CLASSIC, there is a possibility that the cell_corner
was already found and this call can be skipped. Understand that and implement as
needed." (:587).

CLM — models/clm/model_mod.f90. Worst case for the land operators.
Both compute_gridcell_value and get_grid_vertval do a linear scan over the
entire state range of the variable to find the columns/PFTs in the target
gridcell — ELEMENTS : do indexi = index1, indexN
(:1422,
:1661), and get_grid_vertval scans
twice (once to count, once to gather,
:1625). Cost is O(number of state elements
of that variable) per call. The code flags it:
"this is the loop that can exploit the knowledge of what columnids or pftids are
needed for any particular gridcell" (:1411).
COSMOS drives this with 3 × nlayers calls per observation, and re-counts layers
every observation.

POP / CICE (dipole & tripole grids) — models/POP/model_mod.f90.
Regular-grid lookup table narrows the candidates, then get_dipole_quad does an
in_quad test over every overlapping target quad
(:1329-1340). Bounded but not free, and
repeated per category by the obs_def_cice_mod operators.

MOM6 — models/MOM6/model_mod.f90. The only user of
GRID_QUAD_FULLY_IRREGULAR, whose documented search is "do an exhaustive search
of the target grid quads which overlap that regular grid box"
(quad_utils_mod.f90:40-45).

CAM-FV — models/cam-fv/model_mod.f90. Horizontal locate is cheap
(regular lon/lat), but interp_handle = get_interp_handle(obs_qty)
(:783) *"should this be a pointer?? is
it replicating the internal arrays on assignment?" What is this comment getting at?

WRF — models/wrf/model_mod.f90. Map projection inverse is closed-form and
cheap, but get_domain_info walks nested domains highest resolution (I think) outward until the
point is inside (:6336-6398) — redone per
call, per level.

3b. Vertical column rebuilt per call — the bigger cost when which_vert /= VERTISLEVEL

This is what makes COSMOS (VERTISHEIGHT), TES nadir (VERTISPRESSURE),
GPS non-local (VERTISHEIGHT), and the pressure-level branches of cwp/tpw so
expensive:

• CAM-FV: get_quad_vals loops do icorner = 1,4 calling
  find_vertical_levels (:914) — a full
  pressure or height column built at 4 corners on every single
  model_interpolate.
• CAM-SE: with dry_mass_vertical_coordinate and precise = .true. (the
  default in model_interpolate, :544),
  build_dry_mass_pressure_columns reads specific humidity, cloud liquid and
  cloud ice at every level to build one pressure column
  (:1668-1762) — ~3 × nlev state
  reads per interpolate call. The comment concedes it: "This is expensive but
  probably necessary to get unbiased forward operators."
• WRF: get_model_height_profile_distrib fetches PH at 4 corners for all
  bt levels (:5387-5479); the pressure
  path is equivalent. Rebuilt per call.

So a GPS non-local ray step under CAM-SE costs 3 identical dry-mass column
builds ≈ 3 × 3 × 32 ≈ 288 state reads to produce 3 numbers. Under CAM-FV it is
3 × 4 corners of column construction.

3c. The multiplier that makes it MPI traffic, not just flops

Every leaf state read goes through get_state
(distributed_state_mod.f90:70).
When the owning task is remote it issues a one-sided get_from_fwd / MPI_Get
for that single element, all copies
(:112).
There is no batching across levels. So a column built one level at a time is
nlev separate small RMA transfers where one contiguous fetch of the column would
do. This is the reason the redundancy hurts at scale rather than just being untidy.

3d. Radiance-specific multiplier: per-channel observations

RTTOV brightness temperatures are one observation per channel. Every channel at
the same footprint re-enters get_expected_radiance and re-runs GETLEVELDATA
from scratch — arrays_prealloced caches only the allocation, never the data.
A 10-channel sounder footprint therefore rebuilds the identical 55-level column
10 times: ~5000 interpolate calls where ~500 would do, and there is no guarantee
same-location obs are even adjacent in the obs sequence.

---

4. Policy rejections collide with the probe-until-failure idiom

The probe idiom (1a) reads any non-zero istatus as "the model has no such
level". But several model_mods also return non-zero istatus to mean
"DART policy says don't assimilate here" — the no_obs_assim_above_level /
highest_obs_pressure_mb top cutoffs. The operator cannot tell the two apart, so
a QC setting silently changes the geometry the operator believes the model has.

4a. CAM: the level cutoff is a no-op — accidentally

obs_too_high (cam_common_code_mod.f90:487)
has a VERTISLEVEL branch:

if (which_vert == VERTISLEVEL) then
   ! level 1 is top; watch less than/greater than in tests
   if (vert_value < no_assim_above_level) my_status = 14
   return
endif

no_assim_above_level is never assigned anywhere in DART. It is declared and
initialized to -1.0_r8 at
:144, read once at
:511, and that is the
only other reference in the tree. init_discard_high_obs
(:1266-1332) sets
no_assim_above_pressure, no_assim_above_height and a local
no_assim_above_scaleh — but not the level. So vert_value < -1.0 is never true
and CAM-FV/CAM-SE never reject a VERTISLEVEL location, whatever
no_obs_assim_above_level is set to.

Two consequences, pulling opposite ways:

• Good, by accident? or put in because of model_level counts?: the level-probe loops (RTTOV, CO_Nadir, upper_atm) on CAM
  terminate on check_good_levels returning my_status = 8 when
  vert_value > ref_nlevels
  (:1152), which is the
  genuine model top. They are not truncated by the namelist.
• Bad: a real observation whose vertical coordinate is VERTISLEVEL is never
  discarded, contrary to what no_obs_assim_above_level advertises and what the
  init_discard_high_obs log message claims ("Discarding observations higher than
  model level N"). The namelist is honoured for pressure, height and scale-height
  obs only. Q. Who cares about obs VERTISLEVEL, OSSEs?

Fixing the 'bug; (no_assim_above_level = no_obs_assim_above_level) would
break every column operator on CAM: level 1 is the CAM model top, so
vert_value < no_assim_above_level is true for exactly levels 1 … N-1 so the probe
loop would fail on its first call and every column operator would return
numlevels = 0. That coupling has to be resolved before the assignment is added.

4b. CAM: VERTISHEIGHT / VERTISPRESSURE rejections do fire for fwd ops #438

discarding_high_obs is set true whenever no_obs_assim_above_level > 0
(cam-fv/model_mod.f90:338-341,
cam-se/model_mod.f90:385-388), and the
height/pressure branches of obs_too_high are live. The operator that gets hurt is
nonlocal GPS RO refractivity:

• ref_local builds its location as VERTISHEIGHT
  (obs_def_gps_mod.f90:592).
• obs_too_high returns 14 for all ensemble members at once (istatus(:) = status1,
  cam-fv/model_mod.f90:800).
• track_status therefore sets return_now, and get_expected_gpsro_ref returns
  immediately (obs_def_gps_mod.f90:499-501).

So the entire GPS observation is discarded as soon as any single point on its ray
path rises above the cutoff,
ray_top by construction: the converters default ray_htop to 15,000 m
(convert_cosmic_gps_cdf.f90:91)
and 30,000 m
(convert_gpsro_bufr.f90:147).
With a cutoff below ray_htop, every GPSRO_REFRACTIVITY obs using the GPSEXC
(nonlocal) form is rejected regardless of its own height Q. How easy it is to set ray_htop to "fail" a bunch of obs? Not sure.

The local (GPSREF) form is unaffected: it makes one ref_local call at the obs
location itself.

Same mechanism, same VERTISHEIGHT/VERTISPRESSURE exposure, applies to COSMOS
(soil moisture at VERTISHEIGHT), TES nadir (VERTISPRESSURE), and the
non-model_levels branches of cwp/tpw, but those run on CLM / (and others?), none of
which implement a top cutoff - maybe a problem in strongly coupled DA?

4c. MPAS-A: the cutoff does apply to VERTISLEVEL, and truncates columns silently

MPAS-A checks highest_obs_pressure_mb inside model_interpolate
(mpas_atm/model_mod.f90:1274-1290).
The guard excludes only SURFACE and UNDEFINED:

if(.not.(is_vertical(location, "SURFACE")) .and. (.not.(is_vertical(location, "UNDEFINED")))) then
   call compute_pressure_at_loc(state_handle, ens_size, location, lpres, istatus)
   where (lpres < highest_obs_pressure_mb * 100.0_r8)
      istatus(:) = 201
   end where

@hkershaw-brown check refactored code.

VERTISLEVEL is not excluded, so the RTTOV COUNTLEVELS probe is subject to it.
The module default is highest_obs_pressure_mb = 100.0 (100 mb,
:254); the shipped
mpas_atm/work/input.nml sets 1.0. So:

• With the module default (any run that omits the entry from &model_nml),
  COUNTLEVELS stops at ~100 mb instead of the model top. RTTOV then computes radiances from a
  profile that is not the whole atmosphere levels. The obs is not rejected; it
  returns a wrong value - assuming that you should be using the whole atmospheric column.
• numlevels is computed once and cached globally via
  arrays_prealloced
  (obs_def_rttov13_mod.f90:3537).
  The truncated count from whatever the first radiance observation happened to be
  is then used for every observation for the rest of the run.
• COUNTLEVELS exits on any(this_istatus /= 0), so one ensemble member being
  above the cutoff truncates the column for all of them. Not sure how varied ensemble members are? Could this be several levels different between members.
• Whether you clip off the top or not depends on the input.nml highest_obs_pressure_mb (or default model_mod) value so yo may not be aware that you are setting it.

Cost note: compute_pressure_at_loc for a VERTISLEVEL location calls
convert_vert_distrib(… VERTISPRESSURE …)
(:4374-4385) — a full barycentric
vertical conversion, with its own cell search on every call in the column
loop, purely to evaluate the cutoff. That roughly doubles the per-level cost
already described in 3a.

4d. Other models

• WRF has no top-pressure cutoff; the level probe terminates on the genuine
  boundsCheck/zloc bad-level path. (It does allocate(v_h, v_p) of
  (0:bt, ens_size) unconditionally on every model_interpolate call,
  :1161 a per-call allocate/deallocate
  inside every column loop.) Also see WRF model_mod - Why not allow last level? #366
• pangu carries its own copy of obs_too_high (copy of cam model_mod)
  (:3006-3041) with
  no_assim_above_level = 13_r8 hard-coded at
  :187 and the opposite inequality
  (vert_value > no_assim_above_level, since pangu level 1 is the surface). It is
  currently inert — the obs_too_high call in model_interpolate is commented out
  (:581-582) but it would truncate at
  level 13 if re-enabled, and discarding_high_obs already defaults to .true.
  \

4e. Related: unguarded probe loops

obs_def_CO_Nadir_mod's model_levels: do is unbounded and writes
p_col(:, lev) before testing istatus
(:296-305,
:458-467),
with p_col dimensioned max_model_levs = 100
(:78) and no
bound check. A model with ≥100 levels (WACCM has 70 or 126) overruns the array.
By contrast RTTOV, COSMOS, cice and upper_atm all check their cap and error out.

---

5. Summary of possible capabilities needed for a column interface:

1. Number of levels at a location (1a) — replaces 8 probe-until-failure loops
   and the maxlevels/MAXLEVELS/maxlayers compile-time caps.
2. A whole column of one quantity at a location, all ensemble members
   (1b) — replaces the nlev × 1 loops in tpw, VTEC, cwp.
3. Several quantities at one location (2 Tier 2, 3b) locate once (not just applicable to columns),
   including the vertical column construction, then read N quantities.
   Also fixes GPS ref_local's three identical column builds.
4. A vector of locations in one call (2 Tier 3) the GPS ray (Q. Are people using this version of the GPS fwd op?) , and the
   TES_nadir per-member loop, note @todo in the code about this already.
5. A distinction between "no such level / outside the domain" and "DART policy
   says don't assimilate here" (4). Today both come back as a single non-zero
   istatus, so a QC namelist setting silently redefines how many levels the
   operator thinks the model has (MPAS, 4c), and a cutoff applied to one point on
   a ray discards the whole observation (CAM + nonlocal GPS, 4b). Capability 1
   ("how many levels?") means you can distinguish these.