Preserve the ERA-5 surface geopotential and start the model-level integration from it

[s-sasaki-earthsea-wizard]

Description

Before this change, _get_from_cds computed the full-model-level geopotential cube from the downloaded surface geopotential and then discarded that surface field. _load_model_level therefore had to start the hydrostatic integration from the lowest full model level (z[0]), where _calculategeoh's contract is the true surface geopotential, so every reconstructed column sat a near-uniform +9.4 to +10.2 m too high (the height of the lowest full level above ground).

This PR:

• keeps the downloaded surface geopotential in the fetched file as its own z_surface variable, aligned on the level coordinate exactly like lnsp (fetch side);
• returns it from _makeDataCubes (None when absent, or when the stored variable does not match the layout the fetch writes) and prefers it in _load_model_level, falling back to the previous z[0] behavior for files fetched before this change — HRES files, which share this loader, continue through the same fallback path (HRES itself was not independently exercised).

Motivation and Context

Fixes #812; follow-up to #811 as discussed there and in #806.

Relation to #805: that rework of the ML fetch also replaces the downloaded surface z with the recomputed cube. If it lands first, the fetch-side change here is a small rebase onto its fetch path (the acceptance criterion is the same: fetched files retain the true surface geopotential as a separate variable); the loader-side change is independent of it.

How Has This Been Tested?

Measured on real archived CDS responses (the raw lnsp/z and t/q files CDS served for a Kanto scene, 2026-03-01 00:00 UTC), replayed through _get_from_cds with a mocked cdsapi.Client so the exact post-processing of this branch runs on the exact bytes CDS delivered. Full report with commit-pinned RAiDER builds and package versions: https://github.com/s-sasaki-earthsea-wizard/raider-arco-era5/blob/8576fbb0585e42ee05f9a206f16b6dd139d61562/reports/ml_offset_regression_812_fix.md

	z_surface passed to _calculategeoh (bitwise identity)	column offset	columns
dev@c80cf1e, pre-fix file	lowest full level	+9.46 to +10.19 m	n=575
fix, pre-fix file	lowest full level	+9.46 to +10.19 m (identical to the last digit)	n=575
fix, file written by the fixed fetch	the surface z CDS delivered	+0.04 to +0.35 m	n=441

(The last row uses the replayed file's own grid, hence the different column count.) The remaining +0.04 to +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.

Unit tests (all pass locally, pytest test/test_ecmwf_fetch.py test/test_ecmwf_levels.py):

• test_ecmwf_fetch.py: the fetched file carries z_surface bitwise in the same layout as lnsp; the loader reproduces the fetch-time lowest-level height (0.02 m tolerance, which the ~10 m pre-fix behavior exceeds with ample margin); a file with z_surface stripped falls back to the pre-fix behavior, sitting higher by the lowest level's height above ground (to within 0.01 m); a z_surface stored in a foreign layout is ignored rather than misread.
• test_ecmwf_levels.py: the descending-latitude fixture now runs in both layouts (with and without z_surface), guarding the orientation of both the z cube and the stored surface field; a new per-column check pins the reconstruction to the stored cube at 1e-3 m.

Type of change

• Bug fix (non-breaking change which fixes an issue)

Checklist:

• I have added an explanation of what your changes do and why you'd like us to include them.
• I have written new tests for your core changes, as applicable.
• I have successfully ran tests with your changes locally.
• My change requires a change to the documentation.
• I have updated the documentation accordingly.