diff --git a/data/pipelines/era5_zonal_monthly.py b/data/pipelines/era5_zonal_monthly.py new file mode 100644 index 0000000..97d55d7 --- /dev/null +++ b/data/pipelines/era5_zonal_monthly.py @@ -0,0 +1,153 @@ +# Copyright (C) 2026 University Corporation for Atmospheric Research +# SPDX-License-Identifier: Apache-2.0 +"""ERA5 -> the committed zonal-mean monthly climatology (task 1.1, SCIENCE-1, BLOCKING-5). + +Produces ``studio/science/data/era5_zonal_monthly_v1.npz`` and its manifest. The product is +(month x level x latitude) for temperature, specific humidity and geopotential height -- a few MB, +committed to the repository with a checksum, so no run ever fetches anything (ADR-006: the raw ERA5 +fields are NOT redistributed; the Copernicus licence permits derived products with attribution). + +Conventions, each a decision recorded in ASSUMPTIONS.md (ASSUMPTION-9): + +* **Years 1991-2020** -- the WMO standard climate normal, so "the June climatology" cites a + standard rather than a habit. +* **Zonal mean, monthly mean** -- SCIENCE-1's answer (issue #53), decided 2026-08-18. +* **2.5-degree grid** at retrieval. The zonal mean over 144 longitudes is insensitive to this, and + it keeps the download ~10^2 MB instead of gigabytes. Latitude resolution of the product is 2.5 + degrees; the lookup interpolates. +* **15 pressure levels, 300-5 hPa** -- the lower-to-middle stratosphere the box lives in. The + lookup refuses to extrapolate outside this range. + +Usage (this script does NOT run in Studio's locked venv -- it needs ``cdsapi``; make a scratch env): + + python -m venv /tmp/cds && /tmp/cds/bin/pip install cdsapi netCDF4 numpy + /tmp/cds/bin/python data/pipelines/era5_zonal_monthly.py fetch /path/to/raw/ + /tmp/cds/bin/python data/pipelines/era5_zonal_monthly.py reduce /path/to/raw/ + +``fetch`` needs ``~/.cdsapirc``. ``reduce`` writes the product and prints the sha256 that goes into +the manifest; commit both files together. +""" + +from __future__ import annotations + +import argparse +import datetime +import hashlib +import json +from pathlib import Path + +#: Pressure levels [hPa], descending altitude. The box's native coordinate is pressure (ADR-003). +LEVELS_HPA = [5, 7, 10, 20, 30, 50, 70, 100, 125, 150, 175, 200, 225, 250, 300] + +#: WMO standard climate normal. +YEARS = list(range(1991, 2021)) + +VARIABLES = {"temperature": "t", "specific_humidity": "q", "geopotential": "z"} + +#: Identifier recorded in every run's provenance when the ERA5 dataset is selected. +DATASET_ID = "era5_zonal_monthly_v1" + +#: Standard gravity [m s^-2], to turn geopotential into geopotential height (WMO value). +STANDARD_GRAVITY = 9.80665 + +PRODUCT = Path(__file__).resolve().parents[2] / "studio" / "science" / "data" + + +def fetch(raw_dir: Path) -> None: + """One CDS request per variable: each queues independently, and one failure loses one file.""" + import cdsapi + + client = cdsapi.Client(quiet=True) + for variable in VARIABLES: + target = raw_dir / f"era5_{variable}.nc" + print(f"requesting {variable} -> {target}") + client.retrieve( + "reanalysis-era5-pressure-levels-monthly-means", + { + "product_type": ["monthly_averaged_reanalysis"], + "variable": [variable], + "pressure_level": [str(level) for level in LEVELS_HPA], + "year": [str(year) for year in YEARS], + "month": [f"{month:02d}" for month in range(1, 13)], + "time": ["00:00"], + "data_format": "netcdf", + "grid": [2.5, 2.5], + }, + str(target), + ) + + +def reduce(raw_dir: Path) -> None: + """Raw monthly fields -> (month x level x lat) climatology, zonal mean then multi-year mean.""" + import netCDF4 + import numpy as np + + arrays: dict[str, "np.ndarray"] = {} + latitude = level = None + for variable, short in VARIABLES.items(): + with netCDF4.Dataset(raw_dir / f"era5_{variable}.nc") as ds: + data = np.asarray(ds[short][:], dtype=np.float64) # (time, level, lat, lon) + times = netCDF4.num2date(ds["valid_time"][:], ds["valid_time"].units) + months = np.asarray([t.month for t in times]) + years = np.asarray([t.year for t in times]) + latitude = np.asarray(ds["latitude"][:], dtype=np.float64) + level = np.asarray(ds["pressure_level"][:], dtype=np.float64) + expected = len(YEARS) * 12 + if data.shape[0] != expected: + raise SystemExit(f"{variable}: {data.shape[0]} months, expected {expected} -- refetch") + if sorted(set(years)) != YEARS: + raise SystemExit(f"{variable}: years {min(years)}-{max(years)} != {YEARS[0]}-{YEARS[-1]}") + zonal = data.mean(axis=3) # over longitude + monthly = np.stack([zonal[months == m].mean(axis=0) for m in range(1, 13)]) + arrays[short] = monthly.astype(np.float32) # (12, level, lat) + + assert latitude is not None and level is not None + # BOTH axes ascending, so the lookup's searchsorted works without a flip at read time. CDS + # delivers latitude 90..-90 and pressure levels 300..5; trusting either order as-is is how the + # first cut of this pipeline put the 300 hPa row under every lookup -- caught by the + # stratosphere-shaped physics tests (T(55 hPa) came back 240 K, a tropospheric number). + lat_order = np.argsort(latitude) + latitude = latitude[lat_order] + level_order = np.argsort(level) + level = level[level_order] + for short in arrays: + arrays[short] = arrays[short][:, level_order][:, :, lat_order] + + PRODUCT.mkdir(parents=True, exist_ok=True) + out = PRODUCT / f"{DATASET_ID}.npz" + np.savez_compressed( + out, + latitude_deg=latitude.astype(np.float32), + level_hpa=level.astype(np.float32), + month=np.arange(1, 13, dtype=np.int16), + temperature_k=arrays["t"], + specific_humidity_kg_kg=arrays["q"], + geopotential_height_m=(arrays["z"] / STANDARD_GRAVITY).astype(np.float32), + ) + digest = hashlib.sha256(out.read_bytes()).hexdigest() + manifest = { + "dataset_id": DATASET_ID, + "source": "ERA5 monthly averaged reanalysis on pressure levels (Copernicus CDS)", + "doi": "10.24381/cds.6860a573", + "licence": ( + "Contains modified Copernicus Climate Change Service information (1991-2020). " + "Derived product; raw ERA5 fields are not redistributed." + ), + "convention": "zonal mean, then 1991-2020 mean, per calendar month (SCIENCE-1, issue #53)", + "years": [YEARS[0], YEARS[-1]], + "grid_deg": 2.5, + "levels_hpa": LEVELS_HPA, + "created": datetime.date.today().isoformat(), + "sha256": digest, + } + (PRODUCT / f"{DATASET_ID}.json").write_text(json.dumps(manifest, indent=2) + "\n") + size_mb = out.stat().st_size / 1e6 + print(f"wrote {out} ({size_mb:.1f} MB), sha256={digest}") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("command", choices=["fetch", "reduce"]) + parser.add_argument("raw_dir", type=Path) + args = parser.parse_args() + (fetch if args.command == "fetch" else reduce)(args.raw_dir) diff --git a/docs/studio/ASSUMPTIONS.md b/docs/studio/ASSUMPTIONS.md index 0824bc8..b249520 100644 --- a/docs/studio/ASSUMPTIONS.md +++ b/docs/studio/ASSUMPTIONS.md @@ -191,3 +191,25 @@ observations in ways this repository does not quantify. Per-dataset caveat that each background's diameters are **dry or ambient** is implicit in the dataset (`AMBIENT_BACKGROUNDS` in `coupled/backgrounds.py`) rather than a declared field. +## ASSUMPTION-9 — Climatology product conventions: 1991–2020, 2.5°, 15 levels, no extrapolation + +**Made:** 2026-08-19 · **Affects:** `era5_zonal_monthly_v1` (the committed product), and every run +with `site.dataset = ERA5` · **Tracks:** SCIENCE-1 (answered), #94 (other platforms) + +Four choices inside SCIENCE-1's answer that the answer itself did not fix: + +- **Years 1991–2020** — the WMO standard climate normal, so "the June climatology" cites a standard + rather than a habit. +- **2.5° retrieval grid.** The zonal mean over 144 longitudes is insensitive to this; it keeps the + raw download ~10² MB instead of gigabytes. The product's latitude resolution is therefore 2.5°, + interpolated linearly at lookup. +- **15 pressure levels, 300–5 hPa** — the lower-to-middle stratosphere the box lives in. Linear + interpolation in **log-pressure** between levels. +- **No extrapolation, ever.** Outside 300–5 hPa the derivation raises (ADR-005). A plausible + temperature for 400 hPa from a stratospheric product is exactly the fabricated number this + project forbids. + +**How to revisit.** Regenerate with `data/pipelines/era5_zonal_monthly.py` after editing its +constants; the product version in the filename and the manifest's sha256 change together, and +provenance pins which product each run used, so old runs stay attributable. + diff --git a/docs/studio/PROGRESS.md b/docs/studio/PROGRESS.md index c0a309d..028a8ad 100644 --- a/docs/studio/PROGRESS.md +++ b/docs/studio/PROGRESS.md @@ -29,6 +29,44 @@ derivations to resolve rather than fixtures. --- +### 2026-08-19 — Task 1.1: ERA5 is in the product, the schema, and the wizard + +Raised in use: *"I do not see ERA5 data being utilized"* — correct; the decisions existed and the +pipeline did not. Now: + +**The committed product.** `data/pipelines/era5_zonal_monthly.py` fetches ERA5 monthly means +(1991–2020, the WMO normal; 15 levels, 300–5 hPa; 2.5° grid — ASSUMPTION-9) and reduces 343 MB of +raw fields to a **0.1 MB** `(month × level × lat)` npz, committed at +`studio/science/data/era5_zonal_monthly_v1.npz` with a sha256 manifest. Raw ERA5 is not +redistributed; the derived product is, with attribution (BLOCKING-5). + +**The schema (0.4.0).** `site.dataset` (USER | ERA5) with the emission-system structure: entered +`given_temperature_k` / `given_h2o_ppmv` kept-but-inert, `temperature_k` / `h2o_ppmv` derived. +Pressure stays entered — it places the box, so it is the lookup's coordinate. Default USER, so every +existing config resolves byte-identically. Provenance records the dataset id + checksum +automatically for ERA5 runs (`ProvenanceRecord.datasets`, empty since Phase 0, now used). + +**The wizard.** Stage 1 gains the ERA5 temperature-profile panel with the box marked on it; under +USER it is context, under ERA5 the marker is the run's value. + +**The physics check that paid for itself twice.** The acceptance tests assert the stratosphere looks +like the stratosphere (tropical cold point 185–205 K, summer pole warmer than winter, H₂O single-digit +ppmv). The first cut of the pipeline shipped CDS's *descending* pressure levels into an interpolation +that assumed ascending — every lookup silently read the 300 hPa row, and T(55 hPa) came back 240 K +with 414 ppmv of water. Plausible numbers at a glance; impossible as a stratosphere. Fixed in the +pipeline (both axes sorted ascending) and refused in the loader (a non-monotonic axis raises rather +than being reordered, since the checksum pins what was read). + +**The headline.** ERA5 says **209.42 K** at the ensemble's site (30°N, June, 55 hPa) where the +ensemble typed **210 K** — the paper's number confirmed to half a kelvin by 30 years of reanalysis. +H₂O comes out 4.15 ppmv against the typed 6.91 (RH = 3% assumption), consistent with the known dry +bias (MLS is issue #94). Tropical cold point: 191.4 K. + +357 Python Tier-A (+11, all running in CI — the product is committed, so the climatology is the +first substantive science whose tests CI executes), 55 vitest. + +--- + ### 2026-08-18 — SCIENCE-2 and SCIENCE-3 answered; ERA5 confirmed Three decisions in one message, closing every science question that gated Phases 1–4 (only diff --git a/docs/studio/plan/PHASE_1.md b/docs/studio/plan/PHASE_1.md index 69c99aa..f8e6c61 100644 --- a/docs/studio/plan/PHASE_1.md +++ b/docs/studio/plan/PHASE_1.md @@ -29,7 +29,7 @@ All 42 schema fields are placed and rendered; `test_layout.py` fails if a new on - **1.0 Wizard shell** — *done (2026-08-17)*. Layout manifest, generated form, override/stale actions, review with diff, submit. React + Vite + TS per ADR-007. -- **1.1 Climatology product** — *unblocked (2026-08-18).* SCIENCE-1 is answered: **zonal-mean, +- **1.1 Climatology product** — *done (2026-08-19).* SCIENCE-1 is answered: **zonal-mean, monthly**, so the reduced product is `(lat × month × level)` — of order 86k values per field, a few MB, small enough to **commit with a checksum** rather than fetch at run time. Longitude stays an input for the solar zenith angle but does not select the meteorology, and the convention goes into diff --git a/studio/modelio/preview.py b/studio/modelio/preview.py index 7c99bfb..3552985 100644 --- a/studio/modelio/preview.py +++ b/studio/modelio/preview.py @@ -226,9 +226,12 @@ def concentration_sensitivity(config: RunConfig) -> dict[str, Any]: resolved = resolve(config).config volume = resolved.injection.plume_volume_cm3 - if volume is None or volume <= 0.0: - raise ValueError("plume_volume_cm3 is unresolved; resolve the config before previewing") - air = air_number_density(resolved.site.pressure_mbar, resolved.site.temperature_k) + temperature = resolved.site.temperature_k + if volume is None or volume <= 0.0 or temperature is None: + raise ValueError( + "plume_volume_cm3 / temperature_k are unresolved; resolve the config before previewing" + ) + air = air_number_density(resolved.site.pressure_mbar, temperature) decades = 3.0 volumes = np.logspace( math.log10(volume) - decades, math.log10(volume) + decades, _CURVE_POINTS // 3 @@ -239,7 +242,7 @@ def concentration_sensitivity(config: RunConfig) -> dict[str, Any]: molar_mass_g_per_mol=SO2_MOLAR_MASS_G_PER_MOL, volume_cm3=float(v), pressure_mbar=resolved.site.pressure_mbar, - temperature_k=resolved.site.temperature_k, + temperature_k=temperature, ) for v in volumes ] @@ -253,10 +256,55 @@ def concentration_sensitivity(config: RunConfig) -> dict[str, Any]: } +def climatology_profile(config: RunConfig) -> dict[str, Any]: + """The ERA5 zonal-mean monthly profile at this latitude and month (stage 1). + + Drawn whatever the dataset selection: under USER it is context ("here is what ERA5 thinks this + place looks like") with the typed values marked against it; under ERA5 the marker IS the derived + value. Pure product read -- no model, no JAX. + """ + from studio.science.climatology import DATASET_ID, geopotential_height_m, load, profile + + resolved = resolve(config).config + month = resolved.schedule.month + latitude = resolved.site.latitude_deg + data = profile(month=month, latitude_deg=latitude) + + # Altitude as a second LABELING of the pressure axis, not a second scale: round-kilometre ticks + # placed at the pressures where those altitudes actually sit at this latitude and month, from + # the product's own geopotential. Interpolated in log-p against z, which is near-linear. + z_km = np.asarray(data["geopotential_height_m"], dtype=float) / 1000.0 + log_p = np.log(np.asarray(data["level_hpa"], dtype=float)) + ascending = np.argsort(z_km) # z falls as pressure rises, so sort once for interp + altitude_ticks = [] + for km in range(int(np.ceil(z_km.min())), int(np.floor(z_km.max())) + 1): + if km % 5 == 0: + pressure = float(np.exp(np.interp(km, z_km[ascending], log_p[ascending]))) + altitude_ticks.append({"km": km, "pressure_hpa": pressure}) + + box_altitude_m = geopotential_height_m( + month=month, latitude_deg=latitude, pressure_mbar=resolved.site.pressure_mbar + ) + return { + **data, + "altitude_ticks": altitude_ticks, + "box_altitude_km": box_altitude_m / 1000.0, + "dataset_id": DATASET_ID, + "sha256_12": load().sha256[:12], + "month": month, + "latitude_deg": latitude, + "selected_dataset": resolved.site.dataset.value, + "box_pressure_mbar": resolved.site.pressure_mbar, + "box_temperature_k": resolved.site.temperature_k, + "box_h2o_ppmv": resolved.site.h2o_ppmv, + } + + #: Panel name -> builder. The API exposes exactly these, so a typo in a panel name is a 404 naming #: the ones that exist rather than an empty chart. PANELS = { "sza": sza_diurnal, + "climatology": climatology_profile, "dilution": dilution_curve, "size-distribution": size_distribution, "bins": bin_grid, diff --git a/studio/modelio/provenance.py b/studio/modelio/provenance.py index 2959e95..147f842 100644 --- a/studio/modelio/provenance.py +++ b/studio/modelio/provenance.py @@ -185,6 +185,18 @@ def record_for( """ config.require_consistent() root = repo_root or repository_root() + + # Datasets are recorded automatically from the config, not trusted to the caller: a run whose + # temperature came from the ERA5 climatology but whose provenance does not say so would defeat + # the point of ADR-006. The import is lazy and the load cached; under USER nothing is touched. + recorded_datasets = dict(datasets or {}) + from studio.schema.enums import ClimatologyDataset + + if config.config.site.dataset is ClimatologyDataset.ERA5: + from studio.science.climatology import load + + product = load() + recorded_datasets[product.dataset_id] = product.sha256 submodules = {} for name in MODEL_SUBMODULES: path = root / name @@ -201,7 +213,7 @@ def record_for( studio_version=studio.__version__, sandbox=describe_checkout(root), submodules=submodules, - datasets=dict(datasets or {}), + datasets=recorded_datasets, resolved_config=config.config.model_dump(mode="json"), overrides={path: record.value for path, record in sorted(config.overrides.items())}, ) diff --git a/studio/modelio/scenario.py b/studio/modelio/scenario.py index 5cf3846..a021473 100644 --- a/studio/modelio/scenario.py +++ b/studio/modelio/scenario.py @@ -74,6 +74,16 @@ def to_scenario(config: RunConfig | ResolvedConfig) -> CoupledScenario: "start with no SO2." ) + temperature = run.site.temperature_k + h2o = run.site.h2o_ppmv + if temperature is None or h2o is None: + raise ValueError( + "site.temperature_k / site.h2o_ppmv are unresolved. They are DERIVED fields since the " + "ambient state gained a dataset selector (SCIENCE-1): run the config through " + "studio.resolve.resolve() before converting it, rather than letting the model run at " + "a temperature nobody chose." + ) + day_of_year = run.schedule.day_of_year if day_of_year is None: raise ValueError( @@ -89,9 +99,9 @@ def to_scenario(config: RunConfig | ResolvedConfig) -> CoupledScenario: } return CoupledScenario( - T=run.site.temperature_k, + T=temperature, P=run.site.pressure_mbar, - WTR=run.site.h2o_ppmv, + WTR=h2o, latitude=run.site.latitude_deg, longitude=run.site.longitude_deg, day_of_year=day_of_year, diff --git a/studio/resolve/registry.py b/studio/resolve/registry.py index 6f72d97..80a1ade 100644 --- a/studio/resolve/registry.py +++ b/studio/resolve/registry.py @@ -22,7 +22,7 @@ from dataclasses import dataclass from typing import Any, Final -from studio.schema.enums import EmissionInput +from studio.schema.enums import ClimatologyDataset, EmissionInput from studio.science import ( SO2_MOLAR_MASS_G_PER_MOL, initial_mixing_ratio_pptv, @@ -64,6 +64,36 @@ def compute(self, values: Mapping[str, Any]) -> Any: return self.fn({path: values[path] for path in self.inputs}) +def _site_temperature(values: Mapping[str, Any]) -> float: + """The temperature the run uses: entered, or the ERA5 climatology at (lat, month, p). + + The import is inside the branch on purpose: under USER (the default) the climatology product is + never touched, so a checkout without it still resolves every existing config. + """ + if values["site.dataset"] != ClimatologyDataset.ERA5: + return float(values["site.given_temperature_k"]) + from studio.science.climatology import temperature_k + + return temperature_k( + month=values["schedule.month"], + latitude_deg=values["site.latitude_deg"], + pressure_mbar=values["site.pressure_mbar"], + ) + + +def _site_h2o(values: Mapping[str, Any]) -> float: + """The water vapour the run uses: entered, or ERA5 specific humidity as ppmv.""" + if values["site.dataset"] != ClimatologyDataset.ERA5: + return float(values["site.given_h2o_ppmv"]) + from studio.science.climatology import h2o_ppmv + + return h2o_ppmv( + month=values["schedule.month"], + latitude_deg=values["site.latitude_deg"], + pressure_mbar=values["site.pressure_mbar"], + ) + + def _day_of_year(values: Mapping[str, Any]) -> int: """(month, day) -> day of year on the fixed non-leap calendar (SCIENCE-1).""" return calendar_day_of_year( @@ -164,6 +194,29 @@ def _so2_initial_pptv(values: Mapping[str, Any]) -> float: #: Derived field path -> how to compute it. Completeness against the schema is enforced by test. DERIVATIONS: Final[dict[str, Derivation]] = { + "site.temperature_k": Derivation( + inputs=( + "site.dataset", + "site.given_temperature_k", + "site.latitude_deg", + "site.pressure_mbar", + "schedule.month", + ), + fn=_site_temperature, + summary="the entered temperature, or the ERA5 zonal-mean monthly climatology at " + "(latitude, month, pressure)", + ), + "site.h2o_ppmv": Derivation( + inputs=( + "site.dataset", + "site.given_h2o_ppmv", + "site.latitude_deg", + "site.pressure_mbar", + "schedule.month", + ), + fn=_site_h2o, + summary="the entered water vapour, or ERA5 specific humidity converted to ppmv", + ), "schedule.day_of_year": Derivation( inputs=("schedule.month", "schedule.day_of_month"), fn=_day_of_year, diff --git a/studio/schema/config.py b/studio/schema/config.py index 9e0b955..513aaf1 100644 --- a/studio/schema/config.py +++ b/studio/schema/config.py @@ -38,6 +38,7 @@ from studio.schema.enums import ( BackgroundAerosol, + ClimatologyDataset, DilutionRegime, EmissionInput, PhotolysisMode, @@ -48,7 +49,7 @@ #: The schema's own version. Bumped on any change to field names, semantics or defaults, because #: those change the config hash and therefore run identity (ADR-006). Old configs are never silently #: reinterpreted under new semantics. -SCHEMA_VERSION = "0.3.0" +SCHEMA_VERSION = "0.4.0" #: Stratospheric background gas composition [pptv] used by the 810-run ensemble #: (``run_ensemble.py:56-57``). Module-level so the default is one object with one source, and so a @@ -80,9 +81,12 @@ class SchemaModel(BaseModel): class Site(SchemaModel): """Where the box is, and the thermodynamic state it sits in. - Phase 0 takes T and p as user input. Phase 1 derives them from climatology at a chosen - (lat, lon, altitude-or-tropopause-relative) point -- which is why they are primary fields now - and become ``derived_from`` targets later, not the other way round. + The ambient state has a selectable source (Phase 1): under ``dataset = USER`` the given values + are used as typed; under ``ERA5`` temperature and water vapour are derived from the committed + zonal-mean monthly climatology (SCIENCE-1) at this latitude, month and pressure. Pressure is + always entered -- it is what places the box, so it is the coordinate of the lookup rather than + a result of it. Same structure as the emission system: entered ``given_*`` values that are kept + but inert when unselected, and derived used-values the model reads. """ latitude_deg: float = SciField( @@ -110,15 +114,28 @@ class Site(SchemaModel): provenance=Provenance.PAPER_ENSEMBLE, source="coupled/paper_ensemble/run_ensemble.py:98", ) - temperature_k: float = SciField( + dataset: ClimatologyDataset = SciField( + default=ClimatologyDataset.USER, + unit=Unit.DIMENSIONLESS, + label="Ambient state from", + description=( + "Source of temperature and water vapour. USER takes the given values as typed; ERA5 " + "derives both from the committed zonal-mean monthly climatology (1991-2020) at this " + "latitude, month and pressure. Pressure is always entered: it places the box, so it " + "is the lookup's coordinate, not its result." + ), + provenance=Provenance.CONVENTION, + source="SCIENCE-1 (issue #53), answered 2026-08-18: zonal-mean monthly, ERA5", + ) + given_temperature_k: float = SciField( default=210.0, unit=Unit.KELVIN, gt=0.0, - label="Temperature", + label="Temperature (entered)", description=( - "Box temperature. The box is isobaric and ISOTHERMAL: the temperature feedback is " - "refused while the radiative calculation has no longwave component, so this value " - "holds for the whole run. See switches.heating_to_t and SCIENCE-4 (issue #56)." + "Box temperature as typed. Used when dataset is USER; under ERA5 the temperature is " + "derived and this value is inert. The box is isobaric and ISOTHERMAL either way " + "(SCIENCE-4): whatever the source, the value holds for the whole run." ), provenance=Provenance.PAPER_ENSEMBLE, source="coupled/paper_ensemble/TABLE_microphysics_parameters.md (Site: 210 K, 55 hPa)", @@ -131,27 +148,63 @@ class Site(SchemaModel): label="Pressure", description=( "Box pressure, the model's native pressure unit (mbar == hPa). Also sets the box " - "altitude used to place the aerosol in the TUV-x radiation column." + "altitude used to place the aerosol in the TUV-x radiation column, and is the " + "vertical coordinate of the climatology lookup when dataset is ERA5." ), provenance=Provenance.PAPER_ENSEMBLE, source="coupled/paper_ensemble/TABLE_microphysics_parameters.md (Site: 210 K, 55 hPa)", examples=[55.0, 120.0], ) - h2o_ppmv: float = SciField( + given_h2o_ppmv: float = SciField( default=6.9104, unit=Unit.PPMV, ge=0.0, - label="Water vapour", + label="Water vapour (entered)", description=( - "Water vapour mixing ratio. The ensemble value is RH = 3% precomputed at 210 K / " - "55 hPa; it is a decimal in native units on purpose (ADR-003) -- do not round-trip it " - "through SI." + "Water vapour mixing ratio as typed. Used when dataset is USER; under ERA5 it is " + "derived and this value is inert. The ensemble value is RH = 3% precomputed at " + "210 K / 55 hPa, a decimal in native units on purpose (ADR-003)." ), provenance=Provenance.PAPER_ENSEMBLE, source="coupled/paper_ensemble/run_ensemble.py:62 (LAT_ALT WTR column, RH = 3%)", + ) + temperature_k: float | None = SciField( + default=None, + unit=Unit.KELVIN, + label="Temperature (used)", + description=( + "The temperature the run uses: the entered value, or the ERA5 zonal-mean monthly " + "climatology interpolated to this latitude, month and pressure. Isothermal for the " + "whole run (SCIENCE-4)." + ), + provenance=Provenance.DERIVED, + derived_from=[ + "site.dataset", + "site.given_temperature_k", + "site.latitude_deg", + "site.pressure_mbar", + "schedule.month", + ], + ) + h2o_ppmv: float | None = SciField( + default=None, + unit=Unit.PPMV, + label="Water vapour (used)", + description=( + "The water vapour the run uses: the entered value, or ERA5 specific humidity " + "converted to a volume mixing ratio at this latitude, month and pressure." + ), + provenance=Provenance.DERIVED, + derived_from=[ + "site.dataset", + "site.given_h2o_ppmv", + "site.latitude_deg", + "site.pressure_mbar", + "schedule.month", + ], caveat=( - "Reanalysis stratospheric water vapour is biased dry, so ERA5 is not the recommended " - "source for this field when Phase 1 lands (BLOCKING-5)." + "Reanalysis stratospheric water vapour is biased dry (BLOCKING-5), so the ERA5 value " + "is a floor more than an estimate; MLS as a dedicated H2O source is issue #94." ), ) @@ -853,7 +906,7 @@ class RunConfig(SchemaModel): physics, and two runs whose only difference is a name are the same computation. """ - schema_version: Literal["0.3.0"] = SciField( + schema_version: Literal["0.4.0"] = SciField( default=SCHEMA_VERSION, unit=Unit.DIMENSIONLESS, label="Schema version", diff --git a/studio/schema/enums.py b/studio/schema/enums.py index 0450194..ad10ed7 100644 --- a/studio/schema/enums.py +++ b/studio/schema/enums.py @@ -123,6 +123,24 @@ class EmissionInput(StrEnum): EMISSION_DURATION = "emission_duration" +class ClimatologyDataset(StrEnum): + """Where stage 1's ambient state comes from. + + USER keeps the values typed -- the paper ensemble's parameterisation, and the default, so + existing configs keep their meaning. ERA5 derives temperature and water vapour from the + committed zonal-mean monthly climatology (SCIENCE-1) at this config's latitude, month and + pressure; pressure itself stays entered, because it is what places the box. + + One member per product that EXISTS. MERRA-2 and MLS are issue #94, and adding an enum member + before its product would let a config claim a derivation that cannot run (ADR-005). + """ + + #: Ambient state typed directly (default; the ensemble's own values). + USER = "user" + #: ERA5 zonal-mean monthly climatology, 1991-2020 (era5_zonal_monthly_v1, SCIENCE-1). + ERA5 = "era5" + + class AxisKind(StrEnum): """How a ``RunSet`` axis combines with the others. See ``studio/schema/runset.py``.""" @@ -139,6 +157,7 @@ class AxisKind(StrEnum): __all__ = [ "AxisKind", "BackgroundAerosol", + "ClimatologyDataset", "DilutionRegime", "EmissionInput", "PhotolysisMode", diff --git a/studio/schema/layout.py b/studio/schema/layout.py index ee2bc23..623b421 100644 --- a/studio/schema/layout.py +++ b/studio/schema/layout.py @@ -91,10 +91,22 @@ class Stage: ), Section( title="Ambient state", - fields=("site.temperature_k", "site.pressure_mbar", "site.h2o_ppmv"), - note="Entered directly today. SCIENCE-1 is answered (zonal-mean monthly, ERA5), so " - "the climatology-derived p/T arrives with the reduced product (task 1.1); until " - "then these are typed.", + fields=( + "site.dataset", + "site.pressure_mbar", + "site.given_temperature_k", + "site.given_h2o_ppmv", + ), + note="Pressure is always entered -- it places the box, so it is the climatology " + "lookup's coordinate rather than its result. The entered T and H2O are used under " + "USER and kept but inert under ERA5.", + ), + Section( + title="And the run uses", + fields=("site.temperature_k", "site.h2o_ppmv"), + note="Entered values, or the ERA5 zonal-mean monthly climatology (1991-2020, " + "SCIENCE-1) interpolated to this latitude, month and pressure. ERA5 water vapour " + "is biased dry (BLOCKING-5); MLS as a dedicated H2O source is issue #94.", ), Section( title="Date and time", diff --git a/studio/science/climatology.py b/studio/science/climatology.py new file mode 100644 index 0000000..356cfc2 --- /dev/null +++ b/studio/science/climatology.py @@ -0,0 +1,242 @@ +# Copyright (C) 2026 University Corporation for Atmospheric Research +# SPDX-License-Identifier: Apache-2.0 +"""Reading the committed ERA5 zonal-mean monthly climatology (SCIENCE-1). + +The product is built by ``data/pipelines/era5_zonal_monthly.py`` and committed at +``studio/science/data/era5_zonal_monthly_v1.npz`` with a manifest carrying its sha256. This module +is the ONLY reader: it verifies the checksum once per process, interpolates, and refuses to +extrapolate. No fetching happens here or anywhere at run time (ADR-006). + +Interpolation choices, stated because they are conventions: + +* **Linear in latitude** between the product's 2.5-degree rows. Zonal-mean fields are smooth at + this scale; the residual against the 0.25-degree native grid is far below interannual spread. +* **Linear in log-pressure** between levels, because temperature is closer to linear in log-p than + in p through the stratosphere (pressure falls exponentially with height). +* **No extrapolation.** Outside 300-5 hPa the derivation raises rather than guessing (ADR-005): a + plausible temperature for 400 hPa from a stratospheric product is exactly the fabricated number + this project forbids. +""" + +from __future__ import annotations + +import hashlib +import json +from functools import lru_cache +from pathlib import Path +from typing import Final, NamedTuple + +import numpy as np +import numpy.typing as npt + +#: Molar masses [g mol^-1]: dry air (US Standard Atmosphere 1976), water (CODATA). +MOLAR_MASS_DRY_AIR_G_PER_MOL: Final = 28.9644 +MOLAR_MASS_WATER_G_PER_MOL: Final = 18.0153 + +#: The committed product this module reads. +DATASET_ID: Final = "era5_zonal_monthly_v1" + +_DATA_DIR = Path(__file__).resolve().parent / "data" + + +class Climatology(NamedTuple): + """The product in memory. Arrays are (month, level, latitude); axes ascend except level.""" + + latitude_deg: npt.NDArray[np.floating] + level_hpa: npt.NDArray[np.floating] + temperature_k: npt.NDArray[np.floating] + specific_humidity_kg_kg: npt.NDArray[np.floating] + geopotential_height_m: npt.NDArray[np.floating] + dataset_id: str + sha256: str + + +class ClimatologyUnavailableError(FileNotFoundError): + """The committed product is missing -- regenerate it, never substitute for it.""" + + +@lru_cache(maxsize=1) +def load() -> Climatology: + """Load and checksum-verify the product, once per process. + + Raises: + ClimatologyUnavailableError: If the product or its manifest is absent. The message carries + the regeneration commands, because the reader who hits this is holding a fresh clone. + ValueError: If the file does not match the manifest's sha256 -- an edited or corrupted + product must not silently feed runs (ADR-006). + """ + npz_path = _DATA_DIR / f"{DATASET_ID}.npz" + manifest_path = _DATA_DIR / f"{DATASET_ID}.json" + if not npz_path.is_file() or not manifest_path.is_file(): + raise ClimatologyUnavailableError( + f"the ERA5 climatology product is not present at {npz_path}. It is committed to the " + f"repository, so this means an incomplete checkout or a build that stripped data " + f"files. To regenerate: data/pipelines/era5_zonal_monthly.py (fetch, then reduce; " + f"needs ~/.cdsapirc). Never substitute a fabricated profile." + ) + manifest = json.loads(manifest_path.read_text()) + digest = hashlib.sha256(npz_path.read_bytes()).hexdigest() + if digest != manifest["sha256"]: + raise ValueError( + f"{npz_path.name} does not match its manifest: sha256 {digest} != " + f"{manifest['sha256']}. The product or the manifest was modified; regenerate both " + f"together with data/pipelines/era5_zonal_monthly.py." + ) + with np.load(npz_path) as data: + # Refuse, not fix: the interpolation's searchsorted requires ascending axes, and the first + # cut of the pipeline shipped descending levels. A product violating this is regenerated, + # never silently reordered here -- the checksum pins WHAT was read, so the reader must not + # change what it means. + for axis in ("latitude_deg", "level_hpa"): + values = data[axis] + if not np.all(np.diff(values) > 0): + raise ValueError( + f"{npz_path.name}: {axis} is not strictly ascending; the product predates the " + f"axis-order fix -- regenerate it with data/pipelines/era5_zonal_monthly.py" + ) + return Climatology( + latitude_deg=data["latitude_deg"].astype(np.float64), + level_hpa=data["level_hpa"].astype(np.float64), + temperature_k=data["temperature_k"].astype(np.float64), + specific_humidity_kg_kg=data["specific_humidity_kg_kg"].astype(np.float64), + geopotential_height_m=data["geopotential_height_m"].astype(np.float64), + dataset_id=str(manifest["dataset_id"]), + sha256=digest, + ) + + +def _interpolate( + field: npt.NDArray[np.floating], + *, + month: int, + latitude_deg: float, + pressure_mbar: float, + climatology: Climatology, +) -> float: + """Bilinear in (latitude, log-pressure) at an exact month. Refuses to extrapolate.""" + if not 1 <= month <= 12: + raise ValueError(f"month must be 1-12, got {month}") + lats = climatology.latitude_deg + if not lats[0] <= latitude_deg <= lats[-1]: + raise ValueError(f"latitude {latitude_deg} outside the product's {lats[0]}..{lats[-1]}") + levels = climatology.level_hpa # strictly ascending; enforced at load + p_lo, p_hi = float(levels.min()), float(levels.max()) + if not p_lo <= pressure_mbar <= p_hi: + raise ValueError( + f"pressure {pressure_mbar} mbar is outside the climatology's {p_lo:g}-{p_hi:g} hPa. " + f"The product covers the lower-to-middle stratosphere; extrapolating a stratospheric " + f"profile would fabricate a value (ADR-005)." + ) + plane = field[month - 1] # (level, lat) + + # latitude bracket + j = int(np.searchsorted(lats, latitude_deg)) + j0, j1 = max(j - 1, 0), min(j, len(lats) - 1) + w_lat = 0.0 if j0 == j1 else (latitude_deg - lats[j0]) / (lats[j1] - lats[j0]) + + # log-pressure bracket (ascending, enforced at load) + logp = np.log(levels) + target = float(np.log(pressure_mbar)) + i = int(np.searchsorted(logp, target)) + i0, i1 = max(i - 1, 0), min(i, len(levels) - 1) + w_p = 0.0 if i0 == i1 else (target - logp[i0]) / (logp[i1] - logp[i0]) + + v00, v01 = plane[i0, j0], plane[i0, j1] + v10, v11 = plane[i1, j0], plane[i1, j1] + return float( + (1 - w_p) * ((1 - w_lat) * v00 + w_lat * v01) + w_p * ((1 - w_lat) * v10 + w_lat * v11) + ) + + +def temperature_k(*, month: int, latitude_deg: float, pressure_mbar: float) -> float: + """Climatological temperature at (month, latitude, pressure).""" + c = load() + return _interpolate( + c.temperature_k, + month=month, + latitude_deg=latitude_deg, + pressure_mbar=pressure_mbar, + climatology=c, + ) + + +def specific_humidity_to_ppmv(q_kg_kg: float) -> float: + """Specific humidity (kg water / kg moist air) -> volume mixing ratio [ppmv]. + + Exact mole-fraction form, not the dilute approximation: x = (q/Mw) / (q/Mw + (1-q)/Md). + At stratospheric q ~ 3e-6 the two differ negligibly, but the exact form costs nothing and + removes an approximation nobody needs to remember. + """ + if q_kg_kg < 0.0: + raise ValueError(f"specific humidity cannot be negative, got {q_kg_kg}") + moles_water = q_kg_kg / MOLAR_MASS_WATER_G_PER_MOL + moles_dry = (1.0 - q_kg_kg) / MOLAR_MASS_DRY_AIR_G_PER_MOL + return moles_water / (moles_water + moles_dry) * 1e6 + + +def h2o_ppmv(*, month: int, latitude_deg: float, pressure_mbar: float) -> float: + """Climatological water vapour at (month, latitude, pressure), in the model's ppmv. + + Carries the same caveat as the schema field: reanalysis stratospheric water vapour is biased + dry (BLOCKING-5); MLS as a dedicated H2O source is issue #94. + """ + c = load() + q = _interpolate( + c.specific_humidity_kg_kg, + month=month, + latitude_deg=latitude_deg, + pressure_mbar=pressure_mbar, + climatology=c, + ) + return specific_humidity_to_ppmv(q) + + +def geopotential_height_m(*, month: int, latitude_deg: float, pressure_mbar: float) -> float: + """Climatological geopotential height at (month, latitude, pressure). + + The product's own z field (ERA5 geopotential / g0), so pressure-to-altitude is the atmosphere's + actual relation at this latitude and month rather than a standard-atmosphere approximation -- + the 55 hPa surface is ~500 m higher in the tropics than at the pole, and this carries that. + """ + c = load() + return _interpolate( + c.geopotential_height_m, + month=month, + latitude_deg=latitude_deg, + pressure_mbar=pressure_mbar, + climatology=c, + ) + + +def profile(*, month: int, latitude_deg: float) -> dict[str, list[float]]: + """T and H2O against pressure at (month, latitude) -- the stage-1 preview panel's data.""" + c = load() + out_t, out_q, out_z = [], [], [] + for level in c.level_hpa: + pressure = float(level) + out_t.append(temperature_k(month=month, latitude_deg=latitude_deg, pressure_mbar=pressure)) + out_q.append(h2o_ppmv(month=month, latitude_deg=latitude_deg, pressure_mbar=pressure)) + out_z.append( + geopotential_height_m(month=month, latitude_deg=latitude_deg, pressure_mbar=pressure) + ) + return { + "level_hpa": [float(level) for level in c.level_hpa], + "temperature_k": out_t, + "h2o_ppmv": out_q, + "geopotential_height_m": out_z, + } + + +__all__ = [ + "DATASET_ID", + "MOLAR_MASS_DRY_AIR_G_PER_MOL", + "MOLAR_MASS_WATER_G_PER_MOL", + "Climatology", + "ClimatologyUnavailableError", + "geopotential_height_m", + "h2o_ppmv", + "load", + "profile", + "specific_humidity_to_ppmv", + "temperature_k", +] diff --git a/studio/science/data/era5_zonal_monthly_v1.json b/studio/science/data/era5_zonal_monthly_v1.json new file mode 100644 index 0000000..8db98c4 --- /dev/null +++ b/studio/science/data/era5_zonal_monthly_v1.json @@ -0,0 +1,31 @@ +{ + "dataset_id": "era5_zonal_monthly_v1", + "source": "ERA5 monthly averaged reanalysis on pressure levels (Copernicus CDS)", + "doi": "10.24381/cds.6860a573", + "licence": "Contains modified Copernicus Climate Change Service information (1991-2020). Derived product; raw ERA5 fields are not redistributed.", + "convention": "zonal mean, then 1991-2020 mean, per calendar month (SCIENCE-1, issue #53)", + "years": [ + 1991, + 2020 + ], + "grid_deg": 2.5, + "levels_hpa": [ + 5, + 7, + 10, + 20, + 30, + 50, + 70, + 100, + 125, + 150, + 175, + 200, + 225, + 250, + 300 + ], + "created": "2026-08-19", + "sha256": "617617884a1a7bab5700e0561b90bf35eebdf0fb12703049ebcae66d143fc078" +} diff --git a/studio/science/data/era5_zonal_monthly_v1.npz b/studio/science/data/era5_zonal_monthly_v1.npz new file mode 100644 index 0000000..e9acf01 Binary files /dev/null and b/studio/science/data/era5_zonal_monthly_v1.npz differ diff --git a/studio/tests/unit/test_api.py b/studio/tests/unit/test_api.py index 0d8f309..43de7c5 100644 --- a/studio/tests/unit/test_api.py +++ b/studio/tests/unit/test_api.py @@ -71,7 +71,7 @@ def test_resolve_returns_the_derived_values_and_the_identity(api: Any) -> None: @pytest.mark.parametrize( "config", [ - {"site": {"temperature_k": -5}}, + {"site": {"given_temperature_k": -5}}, {"microphysics": {"n_bins": 100}}, {"switches": {"heating_to_t": True}}, {"site": {"temprature_k": 210}}, diff --git a/studio/tests/unit/test_cli.py b/studio/tests/unit/test_cli.py index 5d71281..0c0d447 100644 --- a/studio/tests/unit/test_cli.py +++ b/studio/tests/unit/test_cli.py @@ -117,7 +117,7 @@ def test_planned_runs_have_distinct_identities(sweep_file: Path) -> None: [ ("not: [valid", "not valid"), ("- a\n- b\n", "must contain a mapping"), - ("site:\n temperature_k: -5\n", "not a valid configuration"), + ("site:\n given_temperature_k: -5\n", "not a valid configuration"), ], ) def test_a_bad_config_file_exits_two_with_a_message( diff --git a/studio/tests/unit/test_climatology.py b/studio/tests/unit/test_climatology.py new file mode 100644 index 0000000..51ead2c --- /dev/null +++ b/studio/tests/unit/test_climatology.py @@ -0,0 +1,203 @@ +# Copyright (C) 2026 University Corporation for Atmospheric Research +# SPDX-License-Identifier: Apache-2.0 +"""The ERA5 zonal-mean monthly climatology: the product, its reader, and the derivations over it. + +Unlike the model-facing tests, these run in CI: the product is COMMITTED (a few MB of npz), not a +gitignored artefact or a submodule. That is deliberate -- it means the climatology path is the first +scientifically substantive code whose tests CI can actually execute. + +The physics assertions use wide, cited bounds. They are not testing ERA5 (ERA5 is the reference); +they are testing that the reduction and interpolation did not scramble axes -- a transposed +(month, level, lat) array still produces plausible-looking numbers at any single point, but it +cannot produce a colder tropical tropopause than midlatitude one AND a summer-warm polar +stratosphere at once. +""" + +from __future__ import annotations + +import pytest + +from studio.resolve import apply_change, resolve +from studio.schema import RunConfig +from studio.science import climatology + + +def _era5(**site: object) -> RunConfig: + return RunConfig.model_validate({"site": {"dataset": "era5", **site}}) + + +@pytest.mark.tier_a +def test_the_product_is_committed_and_matches_its_manifest() -> None: + """Presence and integrity: the sha256 in the manifest is the one the loader verifies.""" + product = climatology.load() + assert product.dataset_id == climatology.DATASET_ID + assert product.temperature_k.shape == (12, len(product.level_hpa), len(product.latitude_deg)) + assert len(product.sha256) == 64 + + +@pytest.mark.tier_a +def test_the_interpolation_is_exact_at_grid_points() -> None: + """At a node of the product, interpolation must return the stored value bit-for-bit.""" + product = climatology.load() + month, level_index, lat_index = 6, 5, 48 # arbitrary interior node + expected = float(product.temperature_k[month - 1, level_index, lat_index]) + got = climatology.temperature_k( + month=month, + latitude_deg=float(product.latitude_deg[lat_index]), + pressure_mbar=float(product.level_hpa[level_index]), + ) + assert got == expected + + +@pytest.mark.tier_a +def test_no_extrapolation_ever() -> None: + """Outside 300-5 hPa the lookup raises rather than guessing (ADR-005).""" + for pressure in (400.0, 1000.0, 2.0): + with pytest.raises(ValueError, match="outside the climatology"): + climatology.temperature_k(month=6, latitude_deg=30.0, pressure_mbar=pressure) + + +@pytest.mark.tier_a +def test_the_stratosphere_looks_like_the_stratosphere() -> None: + """Axis-scrambling detectors, with generous bounds. + + * The tropical tropopause region (~100 hPa) is the coldest place in the product's domain: + ~190-205 K climatologically (Seidel et al. 2001 put the tropical cold point near 190 K). + * 30N at ~50 hPa in June sits near 210-220 K -- the paper ensemble's 210 K at 55 hPa is + exactly this number, which is the strongest single check available. + * The summer polar stratosphere is WARMER than the winter one at 30 hPa by tens of kelvin + (no polar night jet in summer). + """ + tropical_tropopause = climatology.temperature_k(month=2, latitude_deg=0.0, pressure_mbar=100.0) + assert 185.0 < tropical_tropopause < 205.0 + + ensemble_site = climatology.temperature_k(month=6, latitude_deg=30.0, pressure_mbar=55.0) + assert 205.0 < ensemble_site < 222.0 + + july_arctic = climatology.temperature_k(month=7, latitude_deg=80.0, pressure_mbar=30.0) + january_arctic = climatology.temperature_k(month=1, latitude_deg=80.0, pressure_mbar=30.0) + assert july_arctic > january_arctic + 10.0, "summer pole must be warmer than winter pole" + + +@pytest.mark.tier_a +def test_stratospheric_water_vapour_is_a_few_ppmv() -> None: + """The wet stratosphere is single-digit ppmv (3-8 typical); 100x off means a unit slipped. + + The q -> ppmv conversion is where a silent factor would live, so the bound is deliberately + tight enough to catch Mw/Md swapped (x1.6) as well as kg/g confusions (x1000). + """ + value = climatology.h2o_ppmv(month=6, latitude_deg=30.0, pressure_mbar=55.0) + assert 2.0 < value < 10.0, f"{value} ppmv is not a stratospheric water vapour value" + + +@pytest.mark.tier_a +def test_the_conversion_is_the_exact_mole_fraction_form() -> None: + assert climatology.specific_humidity_to_ppmv(0.0) == 0.0 + # Dilute limit: q * Md/Mw * 1e6, to first order. + approx = 3e-6 * (28.9644 / 18.0153) * 1e6 + assert climatology.specific_humidity_to_ppmv(3e-6) == pytest.approx(approx, rel=1e-5) + with pytest.raises(ValueError, match="negative"): + climatology.specific_humidity_to_ppmv(-1e-6) + + +@pytest.mark.tier_a +def test_era5_dataset_derives_the_ambient_state() -> None: + """The schema path end to end: dataset=ERA5 fills T and H2O from the product.""" + resolved = resolve(_era5()).config + assert resolved.site.temperature_k == pytest.approx( + climatology.temperature_k(month=6, latitude_deg=30.0, pressure_mbar=55.0) + ) + assert resolved.site.h2o_ppmv == pytest.approx( + climatology.h2o_ppmv(month=6, latitude_deg=30.0, pressure_mbar=55.0) + ) + # The given values are kept but inert, exactly as in the emission system. + assert resolved.site.given_temperature_k == 210.0 + + +@pytest.mark.tier_a +def test_user_dataset_is_byte_identical_to_the_old_behaviour() -> None: + """The default path must not so much as touch the product.""" + resolved = resolve(RunConfig()).config + assert resolved.site.temperature_k == 210.0 + assert resolved.site.h2o_ppmv == 6.9104 + assert resolved.injection.so2_initial_pptv == pytest.approx(3309115922.996412, rel=1e-15) + + +@pytest.mark.tier_a +def test_changing_the_month_moves_the_era5_temperature() -> None: + """The wizard's claim, on real data: edit the month, the ambient state follows.""" + june = resolve(_era5()) + december = apply_change(june, "schedule.month", 12) + t_june = june.config.site.temperature_k + t_december = december.config.site.temperature_k + assert t_june is not None and t_december is not None + assert t_june != t_december, "30N at 55 hPa has a real seasonal cycle" + assert december.is_consistent + + +@pytest.mark.tier_a +def test_the_concentration_follows_the_era5_temperature() -> None: + """The chain reaches the model's input: T from ERA5 changes the air density, so the pptv.""" + user = resolve(RunConfig()).config + era5 = resolve(_era5()).config + assert era5.site.temperature_k is not None and user.site.temperature_k is not None + if era5.site.temperature_k != user.site.temperature_k: + assert era5.injection.so2_initial_pptv != user.injection.so2_initial_pptv + # Same number density either way; pptv scales with T at fixed p (ideal gas). + ratio = era5.injection.so2_initial_pptv / user.injection.so2_initial_pptv + assert ratio == pytest.approx(era5.site.temperature_k / user.site.temperature_k, rel=1e-12) + + +@pytest.mark.tier_a +def test_era5_runs_record_the_dataset_in_provenance(repo_root: object) -> None: + """ADR-006: a run whose temperature came from ERA5 must say so, with the checksum.""" + from pathlib import Path + + from studio.modelio.provenance import record_for + + root = repo_root if isinstance(repo_root, Path) else None + if root is None or not (root / "stratchem-jax" / ".git").exists(): + pytest.skip("model submodules not checked out") + record = record_for(resolve(_era5()), repo_root=root) + assert record.datasets == {climatology.DATASET_ID: climatology.load().sha256} + assert record_for(resolve(RunConfig()), repo_root=root).datasets == {} + + +@pytest.mark.tier_a +def test_altitude_and_pressure_agree_with_the_site_name() -> None: + """The ensemble's site is literally labelled ``30N_20km`` at 55 hPa; the product must agree. + + Bounds are generous (19-22 km) because the 55 hPa surface moves with season and latitude -- + but a wrong-axis or wrong-units bug lands kilometres away, not hundreds of metres. + """ + z = climatology.geopotential_height_m(month=6, latitude_deg=30.0, pressure_mbar=55.0) + assert 19_000.0 < z < 22_000.0, f"{z} m is not ~20 km" + + +@pytest.mark.tier_a +def test_height_falls_as_pressure_rises() -> None: + """Monotonicity across the whole level range -- the relation the right axis relies on.""" + heights = [ + climatology.geopotential_height_m(month=6, latitude_deg=30.0, pressure_mbar=float(p)) + for p in (5, 30, 100, 300) + ] + assert heights == sorted(heights, reverse=True) + + +@pytest.mark.tier_a +def test_the_profile_panel_carries_the_altitude_labeling(repo_root: object) -> None: + """Round-km ticks placed at their true pressures, and the box's own altitude.""" + from studio.modelio.preview import climatology_profile + + panel = climatology_profile(RunConfig()) + assert panel["box_altitude_km"] == pytest.approx(20.2, abs=0.5) + ticks = panel["altitude_ticks"] + assert [t["km"] for t in ticks] == sorted({t["km"] for t in ticks}), "ascending, no duplicates" + assert all(t["km"] % 5 == 0 for t in ticks), "round kilometres only" + for tick in ticks: + back = climatology.geopotential_height_m( + month=6, latitude_deg=30.0, pressure_mbar=tick["pressure_hpa"] + ) + # Round-trip through the inverse interpolation: the tick must sit where it claims, within + # the interpolation's own error (log-p linear both ways). + assert back / 1000.0 == pytest.approx(tick["km"], abs=0.1) diff --git a/studio/tests/unit/test_config_hash.py b/studio/tests/unit/test_config_hash.py index bd13377..8e8fd56 100644 --- a/studio/tests/unit/test_config_hash.py +++ b/studio/tests/unit/test_config_hash.py @@ -34,7 +34,10 @@ #: * ...46cbe3 -> ...373ab4 (0.2.0) when the temperature feedback was refused. The VALUE of #: heating_to_t did not change -- False either way -- but schema_version is part of the hashed #: payload, which is what makes "old configs are never silently reinterpreted" true. -#: * ...6d3a73 -> the value below, still 0.3.0 and unmerged, when the release date became a month +#: * 0.3.0 -> 0.4.0 when the ambient state gained its dataset selector (SCIENCE-1 implemented): +#: given_temperature_k / given_h2o_ppmv entered, temperature_k / h2o_ppmv derived. The default RUN +#: is unchanged -- USER dataset, the same 210 K -- but the config has three more fields. +#: * ...6d3a73 -> ...7813ae, still 0.3.0 and unmerged, when the release date became a month #: and a day with day_of_year derived (SCIENCE-1: the climatology is a monthly average, so there #: is no year, and the month must not be enterable twice). #: * ...373ab4 -> ...6d3a73 (0.3.0) when the emission system arrived: a selector for which @@ -44,7 +47,7 @@ #: the last bit -- but the config that describes it now has four more fields, so it is a different #: configuration and must hash differently. A schema that grew a field without moving the hash #: would be one where two different configs could share an identity. -GOLDEN_DEFAULT_HASH = "a201086e2dc9039fe95b0d11bd184a7fccad39b5c9cbca83ee83c766957813ae" +GOLDEN_DEFAULT_HASH = "b8ca6becb22b49668e22d1f2dce8451b79bd3159c82fb0966cce182ccdd35ba6" @pytest.mark.tier_a diff --git a/studio/tests/unit/test_modelio_equivalence.py b/studio/tests/unit/test_modelio_equivalence.py index da7f432..34d8d28 100644 --- a/studio/tests/unit/test_modelio_equivalence.py +++ b/studio/tests/unit/test_modelio_equivalence.py @@ -63,7 +63,9 @@ def studio_scenario() -> Any: def test_the_golden_case_is_the_schemas_default(build_scenario: Any) -> None: """``RunConfig()`` with no arguments IS the golden case; nothing has to be set up to get it.""" reference = build_scenario(GOLDEN_AXES) - config = RunConfig() + # temperature_k is DERIVED since 0.4.0 (dataset selector), so the default config resolves it; + # under the default USER dataset it must equal the entered value, which is the golden case's. + config = resolve(RunConfig()).config assert config.site.temperature_k == reference.T assert config.site.pressure_mbar == reference.P assert config.background.aerosol is BackgroundAerosol.SABR_220 diff --git a/studio/tests/unit/test_resolve_graph.py b/studio/tests/unit/test_resolve_graph.py index b3080a5..3f94639 100644 --- a/studio/tests/unit/test_resolve_graph.py +++ b/studio/tests/unit/test_resolve_graph.py @@ -119,6 +119,8 @@ def test_the_schemas_derived_chain() -> None: assert schema_derived_fields() == ( "injection.emission_duration_s", "schedule.day_of_year", + "site.h2o_ppmv", + "site.temperature_k", "injection.emission_rate_kg_s", "injection.plume_length_m", "injection.plume_volume_cm3", @@ -126,6 +128,8 @@ def test_the_schemas_derived_chain() -> None: ) order = schema_derived_fields() for earlier, later in ( + # The mixing ratio reads the DERIVED temperature (0.4.0), so the site must resolve first. + ("site.temperature_k", "injection.so2_initial_pptv"), ("injection.emission_duration_s", "injection.emission_rate_kg_s"), ("injection.emission_duration_s", "injection.plume_length_m"), ("injection.plume_length_m", "injection.plume_volume_cm3"), diff --git a/studio/tests/unit/test_resolve_registry.py b/studio/tests/unit/test_resolve_registry.py index b3c4919..1ff8c23 100644 --- a/studio/tests/unit/test_resolve_registry.py +++ b/studio/tests/unit/test_resolve_registry.py @@ -53,7 +53,7 @@ def test_declared_inputs_match_derived_from_exactly() -> None: @pytest.mark.tier_a def test_an_unregistered_field_raises_rather_than_returning_none() -> None: with pytest.raises(NotImplementedError, match="no derivation registered"): - derivation_for("site.temperature_k") + derivation_for("site.pressure_mbar") @pytest.mark.tier_a diff --git a/studio/tests/unit/test_resolve_resolver.py b/studio/tests/unit/test_resolve_resolver.py index a81d8a9..83de735 100644 --- a/studio/tests/unit/test_resolve_resolver.py +++ b/studio/tests/unit/test_resolve_resolver.py @@ -71,7 +71,8 @@ def test_chained_derivations_resolve_in_order() -> None: [ ("injection.plume_length_m", 30000.0, {VOLUME, SO2_PPTV}), ("injection.plume_width_m", 20.0, {VOLUME, SO2_PPTV}), - ("site.temperature_k", 213.0, {SO2_PPTV}), + # Editing the ENTERED temperature recomputes the derived one, then the mixing ratio. + ("site.given_temperature_k", 213.0, {"site.temperature_k", SO2_PPTV}), ("site.pressure_mbar", 120.0, {SO2_PPTV}), # Mass reaches the REPORTED emission rate as well since 0.3.0: R = M / t, so twice the mass # over the same track is twice the rate. It does NOT reach the length or the volume, which @@ -105,7 +106,7 @@ def test_an_override_is_never_overwritten_by_a_recomputation() -> None: assert overridden.config.injection.so2_initial_pptv == 5.0e9 assert overridden.is_consistent, "an override anchored to the current inputs is not stale" - after_edit = apply_change(overridden, "site.temperature_k", 213.0) + after_edit = apply_change(overridden, "site.given_temperature_k", 213.0) assert after_edit.config.injection.so2_initial_pptv == 5.0e9 @@ -113,7 +114,7 @@ def test_an_override_is_never_overwritten_by_a_recomputation() -> None: def test_a_moved_input_marks_the_override_stale_with_both_values() -> None: """ "Stale" without "here is what it would be" leaves the user to recompute by hand.""" overridden = set_override(resolve(RunConfig()), SO2_PPTV, 5.0e9) - after_edit = apply_change(overridden, "site.temperature_k", 213.0) + after_edit = apply_change(overridden, "site.given_temperature_k", 213.0) assert after_edit.stale_fields == (SO2_PPTV,) (entry,) = after_edit.stale @@ -130,7 +131,7 @@ def test_an_unrelated_edit_does_not_make_an_override_stale() -> None: """Only a change to one of ITS inputs counts. Flagging on every edit would train users to dismiss the flag.""" overridden = set_override(resolve(RunConfig()), VOLUME, 3.0e12) - after = apply_change(overridden, "site.temperature_k", 213.0) + after = apply_change(overridden, "site.given_temperature_k", 213.0) assert after.is_consistent assert after.config.injection.plume_volume_cm3 == 3.0e12 @@ -147,7 +148,7 @@ def test_a_downstream_auto_field_uses_the_overridden_value() -> None: @pytest.mark.tier_a def test_accept_derived_drops_the_override_and_recomputes() -> None: stale = apply_change( - set_override(resolve(RunConfig()), SO2_PPTV, 5.0e9), "site.temperature_k", 213.0 + set_override(resolve(RunConfig()), SO2_PPTV, 5.0e9), "site.given_temperature_k", 213.0 ) accepted = accept_derived(stale, SO2_PPTV) assert accepted.is_consistent @@ -162,7 +163,7 @@ def test_keep_override_re_anchors_and_clears_staleness() -> None: """The user has said, knowingly, that their value still applies -- that is the difference between this and never having flagged it.""" stale = apply_change( - set_override(resolve(RunConfig()), SO2_PPTV, 5.0e9), "site.temperature_k", 213.0 + set_override(resolve(RunConfig()), SO2_PPTV, 5.0e9), "site.given_temperature_k", 213.0 ) kept = keep_override(stale, SO2_PPTV) assert kept.is_consistent @@ -170,7 +171,7 @@ def test_keep_override_re_anchors_and_clears_staleness() -> None: assert kept.overrides[SO2_PPTV].inputs["site.temperature_k"] == 213.0 # ...and it goes stale again on the NEXT change, rather than being permanently silenced - assert apply_change(kept, "site.temperature_k", 220.0).stale_fields == (SO2_PPTV,) + assert apply_change(kept, "site.given_temperature_k", 220.0).stale_fields == (SO2_PPTV,) @pytest.mark.tier_a @@ -178,7 +179,7 @@ def test_a_stale_config_refuses_to_pass_as_consistent() -> None: """A stale config still has a hash, and that is the trap: a stable identity for numbers that do not follow from each other.""" stale = apply_change( - set_override(resolve(RunConfig()), SO2_PPTV, 5.0e9), "site.temperature_k", 213.0 + set_override(resolve(RunConfig()), SO2_PPTV, 5.0e9), "site.given_temperature_k", 213.0 ) with pytest.raises(InconsistentConfigError, match="stale override"): stale.require_consistent() @@ -190,7 +191,7 @@ def test_the_stale_list_travels_with_the_config() -> None: """Serialising must carry the stale list, or a persisted config loses the fact that it is inconsistent -- exactly what the plan forbids.""" stale = apply_change( - set_override(resolve(RunConfig()), SO2_PPTV, 5.0e9), "site.temperature_k", 213.0 + set_override(resolve(RunConfig()), SO2_PPTV, 5.0e9), "site.given_temperature_k", 213.0 ) restored = type(stale).model_validate_json(stale.model_dump_json()) assert restored.stale_fields == (SO2_PPTV,) @@ -209,9 +210,14 @@ def test_editing_a_derived_field_directly_is_an_override() -> None: @pytest.mark.tier_a def test_overriding_a_primary_field_is_refused() -> None: - """Primary fields have no derivation to be stale against; the concept does not apply.""" + """Primary fields have no derivation to be stale against; the concept does not apply. + + Pressure, deliberately: it is one of the few fields that stayed primary through 0.3.0 (the + emission system) and 0.4.0 (the dataset selector), precisely because everything else is + derived FROM it. + """ with pytest.raises(ValueError, match="not a derived field"): - set_override(resolve(RunConfig()), "site.temperature_k", 999.0) + set_override(resolve(RunConfig()), "site.pressure_mbar", 999.0) @pytest.mark.tier_a diff --git a/studio/tests/unit/test_runner.py b/studio/tests/unit/test_runner.py index b46ad9d..9b691d3 100644 --- a/studio/tests/unit/test_runner.py +++ b/studio/tests/unit/test_runner.py @@ -285,7 +285,7 @@ def test_a_bad_input_file_exits_distinctly_from_a_model_failure( import subprocess bad = tmp_path / "bad.json" - bad.write_text('{"config": {"site": {"temperature_k": -5}}}', encoding="utf-8") + bad.write_text('{"config": {"site": {"given_temperature_k": -5}}}', encoding="utf-8") proc = subprocess.run( [sys.executable, "-m", "studio.cli.run", str(bad), str(tmp_path / "out")], capture_output=True, diff --git a/studio/tests/unit/test_schema_export.py b/studio/tests/unit/test_schema_export.py index a79f20d..8d3788c 100644 --- a/studio/tests/unit/test_schema_export.py +++ b/studio/tests/unit/test_schema_export.py @@ -42,7 +42,7 @@ def test_metadata_survives_the_export() -> None: """Unit, provenance and source must reach the client, or the form cannot explain a field.""" schema = run_config_json_schema() site = schema["$defs"]["Site"]["properties"] - temperature = site["temperature_k"][EXTENSION_KEY] + temperature = site["given_temperature_k"][EXTENSION_KEY] assert temperature["unit"] == "K" assert temperature["provenance"] == "paper_ensemble" assert "TABLE_microphysics_parameters.md" in temperature["source"] diff --git a/studio/tests/unit/test_schema_metadata.py b/studio/tests/unit/test_schema_metadata.py index 8a8c482..adbd586 100644 --- a/studio/tests/unit/test_schema_metadata.py +++ b/studio/tests/unit/test_schema_metadata.py @@ -142,9 +142,9 @@ def test_bounded_quantities_declare_their_range() -> None: must_be_bounded = { "site.latitude_deg", "site.longitude_deg", - "site.temperature_k", + "site.given_temperature_k", "site.pressure_mbar", - "site.h2o_ppmv", + "site.given_h2o_ppmv", "schedule.month", "schedule.day_of_month", "schedule.start_utc_hour", @@ -167,7 +167,7 @@ def test_bounded_quantities_declare_their_range() -> None: ("path", "value"), [ ("site.latitude_deg", 91.0), - ("site.temperature_k", 0.0), + ("site.given_temperature_k", 0.0), ("site.pressure_mbar", -1.0), ("schedule.month", 13), ("schedule.day_of_month", 32), diff --git a/studio/tests/unit/test_web_contract.py b/studio/tests/unit/test_web_contract.py index dc97639..67b89d0 100644 --- a/studio/tests/unit/test_web_contract.py +++ b/studio/tests/unit/test_web_contract.py @@ -42,7 +42,7 @@ def _node(schema: dict[str, Any], path: str) -> dict[str, Any]: ("switches.heating_to_t", "const", "single accepted value -> read-only"), ("switches.sulfur", "type", "boolean -> checkbox"), ("dilution.regime", "$ref", "string enum by reference ->