ERA-5 model levels: loader passes the lowest full-level z where `_calculategeoh` expects the true surface geopotential (~+10 m column height bias)

[s-sasaki-earthsea-wizard]

Context

Follow-up to #806, as discussed in #811. #811 fixes the north–south mirror of the geopotential cube; this issue tracks the smaller residual item acknowledged there, so it doesn't get lost when #806 is closed by the merge.

The defect

_load_model_level passes z[0] — the lowest full model level after the axis fix — as z_surface to _calculategeoh, whose contract is the true 2D surface geopotential. The hydrostatic integration therefore starts from the lowest-full-level height instead of the terrain surface.

Measured effect (same harness, raw file, and frozen environment as the #806 / #811 measurements, RAiDER at #811 head 919d6b6): every reconstructed column carries a near-uniform vertical offset of +9.44 to +10.20 m (mean +9.81 m) against a true-surface reference — the height of the lowest full level above ground. Within-column drift is ≤ 0.43 m (std), so the error is essentially a rigid upward shift of each column. Details: reports/ml_offset_regression_pr811.md.

Why the loader can't do better today

Since f81225b (#751), the fetch step recomputes the full-level z cube with calcgeoh and overwrites the downloaded surface z in the stored file. The true surface field is discarded at fetch time, so no loader-side change alone can recover it for freshly fetched files.

Proposed fix

1. Fetch: keep the downloaded surface geopotential as its own variable (e.g. z_surface) alongside the recomputed full-level cube, instead of overwriting it.
2. Loader: use z_surface when present; fall back to the current post-Fixes a north–south flip of the geopotential cube in the ECMWF/ERA-5 model-level reader. #811 behavior (z[0], lowest full level) for existing cached files, so nothing breaks.
3. Test: extend test/test_ecmwf_levels.py with a fixture carrying a distinct z_surface, asserting the surface height matches the terrain to well under the ~10 m lowest-level offset.

Interaction with #805

The reworked ML fetch in #805 (checked at its current head 6586c94) also replaces the downloaded surface z with the recomputed cube, in both its single-request and batch paths. So whichever fetch implementation lands first, the acceptance criterion is the same: fetched model-level files retain the true surface geopotential as a separate variable, and the loader prefers it. If #805 merges first, item 1 above moves into that fetch path.

I'm happy to submit a PR for this once #811 is merged.

[jlmaurer]

@s-sasaki-earthsea-wizard thanks! this is interesting, and definitely important to fix. If you want to go ahead with the PR I'll review! @royagrace pinging you here as well since this is affecting the geometric heights calculation relevant to #805.

[s-sasaki-earthsea-wizard]

Thanks @jlmaurer — yes, I'd like to take this on. I'll open a PR.

@royagrace, since #805 touches the same fetch path, I'd rather not duplicate your work. My plan is to keep the downloaded surface geopotential as a separate variable at fetch time and have _load_model_level prefer it, with a fallback to the current z[0] for existing cached files. If that overlaps with what you already have, I'm happy to scope my PR to the loader side only and let the fetch change land with yours — whichever is easier on your end.

I'll follow up here once the PR is up.

[s-sasaki-earthsea-wizard]

The PR is up: #820.

Verification on real archived CDS responses (replayed through _get_from_cds with a mocked client, so the exact post-processing of the branch ran on the exact bytes CDS served): https://github.com/s-sasaki-earthsea-wizard/raider-arco-era5/blob/8576fbb0585e42ee05f9a206f16b6dd139d61562/reports/ml_offset_regression_812_fix.md

• current dev (c80cf1e), pre-fix file: the loader still receives the lowest full level bitwise; columns +9.46 to +10.19 m high (unchanged after Update xarray usage in RAiDER to use lazy loading #808).
• fix branch, pre-fix file: identical to the last digit — the fallback preserves existing behavior for cached files (and HRES, which continues through the same fallback path).
• fix branch, file written by the fixed fetch: the loader receives, bitwise, the surface geopotential CDS delivered; the offset collapses to +0.04–0.35 m.

The remaining +0.04–0.35 m is reproduced quantitatively by a separate pre-existing defect — the fetch stores virtual temperature in t because calcgeoh mutates its argument in place — filed as #819; the surface field itself arrives bitwise intact.

On #805: the fetch-side change here is one added assign variable, so if #805 lands first I'll rebase it onto its fetch path — the loader side is independent either way.