nodata marker is -999., which real data can take — but it cannot simply be changed (it feeds the thermodynamics)

[mvreeuwijk]

Summary

modmpi's nodata (the value the slab/column averaging routines write where a
slab has no fluid cells) is -999., and it should not be. It is a value real data
can take — a pressure perturbation or a temperature could be -999 — so missing
data cannot be told apart from a number. It is also the NetCDF _FillValue
(modstat_nc), so it is what readers mask on.

The obvious fix — pick a value nothing physical can reach — does not work as a
drop-in, and this issue exists to record why, because the reason is not obvious
and the naive attempt fails in a confusing way.

Why it cannot simply be swapped

The marker is not confined to output. It reaches the thermodynamics.
-999. survives today only because it is mild enough to produce slightly wrong
numbers instead of a trap.

Measured on tests/cases/091 (raised STL, slabs kb and kb+1 fully solid),
Intel 2021a Debug: setting nodata = -1e30 makes the case die with

Caught signal 8 (Floating point exception: floating-point divide by zero)
  modthermodynamics_mp_fromztop_()   src/modthermodynamics.f90:432
  modthermodynamics_mp_diagfld_()    src/modthermodynamics.f90:339

because at the top ghost level th0av comes out exactly equal and opposite to the
level below:

th0av(65) = -25232932864.0
th0av(64) = +25232932864.0
thetah(65) = (th0av(65)*dzf + th0av(64)*dzf)/(2*dzh) = 0   exactly
thvh = 0  ->  divide by zero

Restoring the kb special case in avexy_ibm does not help — this is the marker
value, not that hack.

So the marker leaks into th0av and propagates into the hydrostatic march. That
leak exists today with -999.; it is simply ~1e-20 relative and therefore
invisible. This is the same contamination BCcleanup_backlog.md describes when it
argues the -999-poisoned march mis-weights buried layers.

A tiny value is worse, not safer. -1e-99 (or anything below ~1e-38) is not an
option: the dumps are single precision, so it underflows to -0. on write and the
marker silently becomes an ordinary zero — undetectable.

What a fix has to do

One of:

1. Mark only at dump time. Keep the averaging routines returning something
   physics-safe and substitute the marker in modstatsdump when writing. Clean
   separation; the physics never sees it.
2. Pick an unmistakable value (e.g. -1e30) and guard every consumer. Truer to
   intent, and it forces the leak into the open instead of leaving it masked — but
   it is a bug hunt: the fromztop divide-by-zero is one consumer, and there may be
   others that merely go quiet.

Option 1 is probably the smaller change. Option 2 fixes an actual latent bug.

Related

• The kb special case in avexy_ibm fills a fully-solid kb slab with
  sum(var)/IId(ke) — an average over no fluid cells, scaled by an unrelated
  slab's count. Its stated reason ("value at kb is used in modthermo") no longer
  holds: modthermodynamics substitutes the base state itself wherever
  IIcs(k)==0 (Remove thls dependency in Exner function calculation in modthermodynamics #302), and it did so even with the hack present, because the hack
  only touched a local copy of IIs. Its one surviving effect is putting a
  fabricated number into the statistics at kb while every other buried slab is
  marked — which is why kb and kb+1 disagree in xytdump for a raised-terrain
  case (290-ish vs -999). Worth removing with this.
• tests/regression/bc_cleanup/run_test.py now asserts the current contract
  (_assert_sentinels_only_in_solid_slabs): tdump/fielddump carry no marker,
  xytdump may carry it but only at fully-solid slabs. That test should be updated
  with whatever this issue decides.

Context

Found while reviewing #332. The marker is now the named parameter modmpi::nodata
rather than a literal repeated in four places, so the value is changeable in one
spot — but deliberately left at -999. there, with the reasoning above recorded in
the comment, precisely so nobody swaps it without doing the work.