diff --git a/changelog/104.docs.md b/changelog/104.docs.md new file mode 100644 index 00000000..73ea3bec --- /dev/null +++ b/changelog/104.docs.md @@ -0,0 +1 @@ +Added a "Concentration-driven runs" documentation page comparing the MAGICC7, FaIR2 and CICEROSCMPY2 adapters' `RunMode.CONCENTRATION_DRIVEN` support. diff --git a/changelog/104.feature.md b/changelog/104.feature.md new file mode 100644 index 00000000..46a06b80 --- /dev/null +++ b/changelog/104.feature.md @@ -0,0 +1 @@ +Added `RunMode.CONCENTRATION_DRIVEN` support to the MAGICC7 adapter, driving the model from `Atmospheric Concentrations|*` inputs (overlaid on the RCMIP3 baseline) instead of emissions. Covers CO2, CH4 and N2O via MAGICC's per-gas concentration flags, plus the F-gas and Montreal-halocarbon groups via the bundled-array `FGAS_FILES_CONC` / `MHALO_FILES_CONC` flags. The concentration-driven decision is made per scenario, so a gas can be concentration-driven in one scenario and emissions-driven in another within the same batch. Requires `rcmip3_bundle_path`. diff --git a/changelog/104.internal.md b/changelog/104.internal.md new file mode 100644 index 00000000..f5d1540d --- /dev/null +++ b/changelog/104.internal.md @@ -0,0 +1 @@ +Added an opt-in MAGICC7 validation harness against IPCC AR6 Table 7.SM.4: load the AR6 probabilistic drawnset, run concentration-driven SSP scenarios, rebase GSAT to 1995-2014 and compare per-period percentiles to the published MAGICC7 column. Gated on a local drawnset (not vendored) and the MAGICC binary, so skipped by default. diff --git a/docs/source/concentration-driven.md b/docs/source/concentration-driven.md new file mode 100644 index 00000000..29ef2772 --- /dev/null +++ b/docs/source/concentration-driven.md @@ -0,0 +1,74 @@ +# Concentration-driven runs + +OpenSCM-Runner can drive its three modern adapters — **MAGICC7**, **FaIR2** +and **CICEROSCMPY2** — from prescribed atmospheric concentrations instead of +emissions, via {class}`openscm_runner.RunMode`: + +```python +from openscm_runner import RunMode, run + +run( + climate_models_cfgs={"MAGICC7": ({"core_climatesensitivity": 3, + "rcmip3_bundle_path": "/path/to/rcmip3"},)}, + scenarios=scenarios, # may carry ``Atmospheric Concentrations|*`` rows + output_variables=("Surface Air Temperature Change", + "Atmospheric Concentrations|CO2"), + mode=RunMode.CONCENTRATION_DRIVEN, +) +``` + +All three adapters share the same contract: + +- **Inputs.** Concentrations are supplied as ``Atmospheric Concentrations|`` + rows in the scenarios DataFrame, overlaid year-by-year (user wins) on a + baseline loaded from the canonical RCMIP Phase 3 Zenodo bundle (record + `20430630`). Species the user does not supply fall back to the baseline. +- **`rcmip3_bundle_path`.** Required on the cfg for concentration-driven runs, + so the baseline is reproducible across machines. +- **Outputs.** Each adapter can return ``Atmospheric Concentrations|`` + for the driven species, so you can check that what comes out matches what was + prescribed. + +## Capability matrix + +| | **MAGICC7** | **FaIR2** | **CICEROSCMPY2** | +|---|---|---|---| +| CO2 / CH4 / N2O | ✅ | ✅ | ✅ | +| F-gases (PFCs / HFCs / SF6 / NF3) | ✅ | ✅ | ✅ | +| Montreal halocarbons | ✅ | ✅ | ✅ | +| Requires `rcmip3_bundle_path` | ✅ | ✅ (or conc rows in the scenario) | ✅ | +| Hybrid baseline + user overlay | ✅ | ✅ | ✅ | +| Returns concentrations as output | ✅ | ✅ | ✅ | +| **Mode scoping** | per-scenario, per-gas (WMGHGs) / per-group (F-gas, MHalo) | per-calibration instance (auto-detected) + protocol mixed-mode | adapter-wide global flag | +| Per-gas mixed mode | WMGHGs independent; F-gas & MHalo all-or-nothing within group | protocol-strict only (CO2 emissions + non-CO2 concentration) | ❌ all-or-nothing | +| Per-scenario mixing within a batch | ✅ (no batch-consistency constraint) | batch-intersection rule applies in mixed-mode | ❌ global | + +**Species coverage is at parity:** all three drive the full WMGHG + +F-gas + Montreal-halocarbon set. + +## Where the adapters differ + +The remaining differences are in *how finely modes can be mixed*, and they +follow each model's native interface rather than any missing feature: + +- **MAGICC7** writes a separate config and ``CONC.IN`` files per + ``(scenario, model)``, so the decision is genuinely per-scenario and, for + CO2/CH4/N2O, per-gas. The F-gas and Montreal-halocarbon groups each share a + single ``*_SWITCHFROMCONC2EMIS_YEAR`` flag, so within a group driving is + all-or-nothing: the writer fills the whole group from the baseline and warns + on any unexpected empty slot. ``HALON1202`` has no RCMIP3 trajectory and + falls back to MAGICC's built-in default. +- **FaIR2** uses a per-species ``input_mode`` that is shared across the whole + FaIR *instance* (all scenarios in a calibration). It therefore keeps a + batch-consistency rule: its protocol mixed-mode (CO2 emissions-driven, + non-CO2 concentration-driven) only engages when every scenario in the batch + supplies the concentrations. +- **CICEROSCMPY2** exposes a single adapter-level ``conc_run`` flag that + applies uniformly to every scenario and species in the call. + +```{note} +Because MAGICC7 drives each scenario independently, it has **no** +batch-consistency requirement: a gas can be concentration-driven in one +scenario and emissions-driven in another within the same batch. This is +intentionally more permissive than FaIR2's shared-instance constraint. +``` diff --git a/docs/source/index.md b/docs/source/index.md index 641d90d7..9a03d5bd 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -28,6 +28,7 @@ ```{toctree} :caption: Contents :maxdepth: 2 +concentration-driven idealised-experiments notebooks development diff --git a/src/openscm_runner/adapters/magicc7/_concentrations_translator.py b/src/openscm_runner/adapters/magicc7/_concentrations_translator.py new file mode 100644 index 00000000..61bc8998 --- /dev/null +++ b/src/openscm_runner/adapters/magicc7/_concentrations_translator.py @@ -0,0 +1,493 @@ +""" +Translate RCMIP3 + user-overlay concentrations into MAGICC7 inputs. + +This module is the conc-driven sibling to the SCEN7-writing path in +:mod:`magicc7`. The strategy mirrors the FaIR2 / CICEROSCMPY2 +adapters: load a baseline from the canonical RCMIP3 Zenodo bundle +(record `20430630`_) and overlay any ``Atmospheric Concentrations|*`` +rows the caller supplies on top, then write one MAGICC ``CONC.IN`` +file per (scenario, species) and emit cfg patches pointing MAGICC at +those files. + +Per-scenario mixed-mode: a gas is concentration-driven for a given +scenario iff that scenario supplies an +``Atmospheric Concentrations|`` trajectory (either directly +in the user ScmRun, or — failing that — via the RCMIP3 baseline). +Because MAGICC writes an independent cfg + ``CONC.IN`` files per +``(scenario, model)``, the decision is genuinely per-scenario: a gas +can be conc-driven in one scenario and emissions-driven in another in +the same batch. This differs from the FaIR2 adapter, whose +``input_mode`` is a per-species flag shared across the whole FaIR +*instance* (all scenarios), forcing a batch-consistency rule that +MAGICC does not need. Species not supplied for a scenario fall back +silently to that scenario's SCEN7 emissions-driven path. + +.. _20430630: https://zenodo.org/records/20430630 +""" +from __future__ import annotations + +import logging +from collections.abc import Iterable +from pathlib import Path + +import pandas as pd + +from ._compat import pymagicc + +LOGGER = logging.getLogger(__name__) + + +# RCMIP3-side short species name -> MAGICC7-side short species name. +# Most of CO2/CH4/N2O/SF6/CF4 round-trip unchanged; this map only +# covers cases where the canonical RCMIP3 CSV's leaf species (after +# stripping intermediate IAMC categories via +# :func:`openscm_runner.io.rcmip3.canonicalise_rcmip3_variable`) does +# NOT match MAGICC's species name. +# +# Halons: RCMIP3 uses ``H-1211``; MAGICC uses ``HALON1211``. +# HFC4310mee: this repo's existing magicc7 ``_VARIABLE_MAP`` already +# folds the ``mee`` suffix away on the emissions path, so we mirror +# that here for the concentrations path. +RCMIP_TO_MAGICC_SPECIES: dict[str, str] = { + "H-1202": "HALON1202", + "H-1211": "HALON1211", + "H-1301": "HALON1301", + "H-2402": "HALON2402", + "HFC4310mee": "HFC4310", + "cC4F8": "CC4F8", +} + + +# MAGICC's namelist binds the bundled-array concentration flags +# (``FGAS_FILES_CONC`` / ``MHALO_FILES_CONC``) *positionally* to these +# species-name lists, taken verbatim from ``MAGCFG_DEFAULTALL.CFG``. +# They are hardcoded (rather than read from the run dir) so this +# translator stays binary-free and unit-testable; a +# ``@pytest.mark.magicc`` cross-check test asserts they still match the +# installed binary's ``fgas_names`` / ``mhalo_names``. Order matters: +# the writer fills one array slot per name below. +FGAS_NAMES: tuple[str, ...] = ( + "CF4", "C2F6", "C3F8", "C4F10", "C5F12", "C6F14", "C7F16", "C8F18", + "CC4F8", "HFC23", "HFC32", "HFC4310", "HFC125", "HFC134A", "HFC143A", + "HFC152A", "HFC227EA", "HFC236FA", "HFC245FA", "HFC365MFC", "NF3", + "SF6", "SO2F2", +) +MHALO_NAMES: tuple[str, ...] = ( + "CFC11", "CFC12", "CFC113", "CFC114", "CFC115", "HCFC22", "HCFC141B", + "HCFC142B", "CH3CCL3", "CCL4", "CH3CL", "CH2CL2", "CHCL3", "CH3BR", + "HALON1211", "HALON1301", "HALON2402", "HALON1202", +) +_FGAS_SET: frozenset[str] = frozenset(FGAS_NAMES) +_MHALO_SET: frozenset[str] = frozenset(MHALO_NAMES) + + +def to_magicc_species(species_short: str) -> str: + """ + Map an RCMIP3 leaf species name to its MAGICC-side name. + + Applies the explicit :data:`RCMIP_TO_MAGICC_SPECIES` renames first + (halons ``H-1211`` -> ``HALON1211``, ``HFC4310mee`` -> ``HFC4310``, + ``cC4F8`` -> ``CC4F8``), then falls back to an upper-casing rule for + pure case differences: RCMIP3 ships mixed-case leaves (``HFC134a``, + ``CCl4``, ``CH3Cl``, ``Halon1211``) whereas MAGICC's + ``FGAS_NAMES`` / ``MHALO_NAMES`` are upper-case. Anything else + round-trips unchanged. + """ + if species_short in RCMIP_TO_MAGICC_SPECIES: + return RCMIP_TO_MAGICC_SPECIES[species_short] + upper = species_short.upper() + if upper in _FGAS_SET or upper in _MHALO_SET: + return upper + return species_short + + +# Default fall-back concentration units, used when +# ``pymagicc.definitions.MAGICC7_CONCENTRATIONS_UNITS`` does not list +# a species (e.g. on older pymagicc versions). The triplet ppm/ppb/ppt +# is the long-standing MAGICC convention: CO2 -> ppm, CH4 / N2O -> +# ppb, all halocarbons / PFCs / SF6 -> ppt. +_FALLBACK_CONC_UNITS: dict[str, str] = { + "CO2": "ppm", + "CH4": "ppb", + "N2O": "ppb", + # All F-gases / PFCs / SF6 / Montreal halocarbons report in ppt. + **{sp: "ppt" for sp in FGAS_NAMES}, + **{sp: "ppt" for sp in MHALO_NAMES}, +} + + +# MAGICC7 exposes per-gas concentration-driving flags +# (``FILE__CONC`` + ``_SWITCHFROMCONC2EMIS_YEAR``) only for +# the three main WMGHGs. F-gases and Montreal halocarbons instead share +# bundled-array flags (``FGAS_FILES_CONC`` indexed positionally by +# :data:`FGAS_NAMES`, ``MHALO_FILES_CONC`` by :data:`MHALO_NAMES`, each +# with a single shared ``*_SWITCHFROMCONC2EMIS_YEAR``). :func:`classify_conc_species` +# routes a species to the right path; everything else is unsupported and +# falls back to the SCEN7 emissions path. +PER_GAS_CONC_SPECIES: frozenset[str] = frozenset({"CO2", "CH4", "N2O"}) +SUPPORTED_PER_GAS_CONC_SPECIES: frozenset[str] = ( + PER_GAS_CONC_SPECIES | _FGAS_SET | _MHALO_SET +) + + +def classify_conc_species(magicc_species: str) -> str | None: + """ + Classify a MAGICC-side species into its conc-driving cfg mechanism. + + Returns ``"per_gas"`` for CO2/CH4/N2O (per-gas ``FILE__CONC`` + flags), ``"fgas"`` / ``"mhalo"`` for species in the bundled-array + groups, or ``None`` if MAGICC has no concentration-driving path for + it (caller should fall back to SCEN7 emissions). + """ + if magicc_species in PER_GAS_CONC_SPECIES: + return "per_gas" + if magicc_species in _FGAS_SET: + return "fgas" + if magicc_species in _MHALO_SET: + return "mhalo" + return None + + +def cfg_keys_for_species(magicc_species: str) -> tuple[str, str]: + """ + Return the ``(file__conc, _switchfromconc2emis_year)`` + MAGICC cfg flag pair for a given MAGICC-side species name. + + Only valid for the per-gas WMGHGs (``CO2`` / ``CH4`` / ``N2O``); + raises :class:`ValueError` otherwise. F-gases / Montreal halocarbons + use bundled-array flags handled separately (see + :func:`classify_conc_species`). + """ + if magicc_species not in PER_GAS_CONC_SPECIES: + raise ValueError( + f"cfg_keys_for_species only handles per-gas species " + f"{sorted(PER_GAS_CONC_SPECIES)}; got {magicc_species!r}. " + f"F-gases / Montreal halocarbons use the bundled-array path." + ) + name = magicc_species.lower() + return (f"file_{name}_conc", f"{name}_switchfromconc2emis_year") + + +def _conc_units_lookup() -> dict[str, str]: + """ + Build an openscm-variable -> magicc concentration-unit lookup. + + Prefer pymagicc's shipped ``MAGICC7_CONCENTRATIONS_UNITS`` table + where available (so unit names track the installed pymagicc + version); fall back to the CO2/CH4/N2O defaults if the table is + missing or doesn't list a species. Callers convert per-species + values to these units before writing the ``.IN`` file. + """ + if pymagicc is None: + return dict(_FALLBACK_CONC_UNITS) + table = getattr( + pymagicc.definitions, "MAGICC7_CONCENTRATIONS_UNITS", None, + ) + if table is None: + return dict(_FALLBACK_CONC_UNITS) + lookup: dict[str, str] = {} + for _, row in table.iterrows(): + magicc_var = row.get("magicc_variable") + unit = row.get("concentration_unit") or row.get("unit") + if not magicc_var or not unit: + continue + if magicc_var.endswith("_CONC"): + magicc_var = magicc_var[: -len("_CONC")] + lookup[magicc_var] = unit + return lookup + + +def build_concentrations_overlay( + scenario_run, + rcmip3_bundle_path, + scenario_names: Iterable[str], + nystart: int = 1750, + nyend: int = 2500, +) -> pd.DataFrame: + """ + Build per-(scenario, MAGICC-species) overlay rows for conc-driven mode. + + Reads ``rcmip_phase3_concentrations_v2.0.0.csv`` for the baseline + and overlays any ``Atmospheric Concentrations|*`` rows the caller + supplied in ``scenario_run``. The decision is per-scenario: every + ``(scenario, species)`` for which a trajectory exists (baseline or + overlay) and which MAGICC can concentration-drive + (:func:`classify_conc_species`) survives. There is no + batch-consistency requirement — a gas can be conc-driven in one + scenario and emissions-driven in another in the same batch. + + Returns a DataFrame with one row per (scenario, MAGICC-species): + + scenario | magicc_species | unit | series + + where ``series`` is a :class:`pandas.Series` indexed by integer + year holding the concentration trajectory. Returns an empty + DataFrame if nothing qualifies (e.g. the caller's batch covers a + scenario the RCMIP3 baseline does not). + """ + scenario_names = list(scenario_names) + if not scenario_names: + return pd.DataFrame() + + # Caller overlay first (so it can shadow the RCMIP3 baseline + # year-by-year if both are present for the same species). + user_rows: list[dict] = [] + if scenario_run is not None and not scenario_run.empty: + conc_run = scenario_run.filter( + variable="Atmospheric Concentrations|*", log_if_empty=False, + ) + if not conc_run.empty: + ts = conc_run.timeseries(time_axis="year") + for index_tuple, values in ts.iterrows(): + meta = dict(zip(ts.index.names, index_tuple)) + scenario = meta.get("scenario") + if scenario not in scenario_names: + continue + variable = meta["variable"] + species_short = variable.split("|", 1)[1] + magicc_species = to_magicc_species(species_short) + series = pd.Series({ + int(year): float(val) + for year, val in values.items() + if nystart <= int(year) <= nyend + and pd.notna(val) + }) + if series.empty: + continue + user_rows.append({ + "scenario": scenario, + "magicc_species": magicc_species, + "unit": meta.get("unit"), + "series": series, + }) + + user_df = pd.DataFrame(user_rows) + + # RCMIP3 baseline second. Cells the caller also supplies are + # merged year-by-year (user values win), so a sparse user override + # (e.g. CO2 for 2050+2100 only) does NOT replace the full + # historical trajectory with a 2-year stub. + base_df = _load_rcmip3_baseline( + Path(rcmip3_bundle_path), scenario_names, nystart, nyend, + ) + + if user_df.empty and base_df.empty: + return pd.DataFrame() + + if user_df.empty: + merged = base_df + elif base_df.empty: + merged = user_df + else: + base_by_key = { + (r["scenario"], r["magicc_species"]): r + for _, r in base_df.iterrows() + } + merged_rows: list[dict] = [] + user_keys: set[tuple[str, str]] = set() + for _, urow in user_df.iterrows(): + key = (urow["scenario"], urow["magicc_species"]) + user_keys.add(key) + brow = base_by_key.get(key) + if brow is None: + merged_rows.append(dict(urow)) + continue + combined = brow["series"].copy() + combined.update(urow["series"]) + merged_rows.append({ + "scenario": urow["scenario"], + "magicc_species": urow["magicc_species"], + "unit": urow["unit"] or brow["unit"], + "series": combined.sort_index(), + }) + for key, brow in base_by_key.items(): + if key not in user_keys: + merged_rows.append(dict(brow)) + merged = pd.DataFrame(merged_rows) + + # Keep only species MAGICC can concentration-drive (per-gas WMGHGs + # plus the bundled-array F-gas / Montreal-halocarbon groups). + # Anything else falls through to the SCEN7 emissions-driven path. + unsupported = ( + set(merged["magicc_species"].unique()) + - SUPPORTED_PER_GAS_CONC_SPECIES + ) + if unsupported: + LOGGER.info( + "MAGICC7 conc-driven: dropping %d species MAGICC cannot " + "concentration-drive; they will be driven by SCEN7 emissions " + "instead: %s", + len(unsupported), sorted(unsupported), + ) + merged = merged[ + merged["magicc_species"].isin(SUPPORTED_PER_GAS_CONC_SPECIES) + ].copy() + if merged.empty: + return pd.DataFrame() + + # No batch-consistency filter: every surviving (scenario, species) + # row is conc-driven for that scenario. MAGICC writes an independent + # cfg + CONC.IN per (scenario, model), so unlike FaIR2 a gas need + # not be supplied by every scenario in the batch. + for scenario in scenario_names: + species = sorted( + merged.loc[merged["scenario"] == scenario, "magicc_species"] + ) + if species: + LOGGER.info( + "MAGICC7 conc-driven: scenario %s — %d species driven by " + "concentration: %s", + scenario, len(species), species, + ) + + # Apply MAGICC's canonical concentration units. Conversion from + # the RCMIP3 / user-supplied unit is currently a no-op assumption + # (RCMIP3 ships ppm/ppb/ppt directly); revisit if a bundle ever + # ships something else for a given species. + unit_lookup = _conc_units_lookup() + merged["unit"] = merged["magicc_species"].apply( + lambda sp: unit_lookup.get(sp, _FALLBACK_CONC_UNITS.get(sp, "1")), + ) + + return merged.reset_index(drop=True) + + +def _load_rcmip3_baseline( + rcmip3_bundle_path: Path, + scenario_names: list[str], + nystart: int, + nyend: int, +) -> pd.DataFrame: + """ + Read the RCMIP3 baseline for the batch's scenarios. + + Returns a frame with the same column layout as the caller-overlay + half of :func:`build_concentrations_overlay`. Empty frame if the + bundle has no rows for any of the requested scenarios. + """ + from ...io.rcmip3 import ( + canonicalise_rcmip3_variable, + load_rcmip3_concentrations, + ) + + df = load_rcmip3_concentrations( + rcmip3_bundle_path, scenarios=scenario_names, + ) + if df.empty: + return pd.DataFrame() + + year_cols = [ + c + for c in df.columns + if isinstance(c, str) and c.isdigit() + and nystart <= int(c) <= nyend + ] + + rows: list[dict] = [] + for _, csv_row in df.iterrows(): + variable = canonicalise_rcmip3_variable(csv_row["Variable"]) + if "|" not in variable: + continue + species_short = variable.split("|", 1)[1] + magicc_species = to_magicc_species(species_short) + series = pd.Series({ + int(year): float(csv_row[year]) + for year in year_cols + if pd.notna(csv_row[year]) + }) + if series.empty: + continue + rows.append({ + "scenario": csv_row["Scenario"], + "magicc_species": magicc_species, + "unit": csv_row["Unit"], + "series": series, + }) + return pd.DataFrame(rows) + + +# MAGICC7 reads every ``*_CONC.IN`` file as a contiguous *annual* +# series spanning its full internal integration window: the binary's +# own shipped concentration files (e.g. ``SSP245_CO2_CONC.IN``) all +# carry ``THISFILE_ANNUALSTEPS = 1`` over 1700-2500 (801 rows). Feeding +# it a sparse / shorter series makes the Fortran ``readdata`` routine +# hit end-of-file. We therefore resample whatever (possibly sparse) +# trajectory we have onto this grid before writing. +MAGICC_CONC_FIRSTYEAR = 1700 +MAGICC_CONC_LASTYEAR = 2500 + + +def _to_annual_magicc_grid(series: pd.Series) -> pd.Series: + """ + Resample a (year-indexed) concentration series onto MAGICC's + annual 1700-2500 grid. + + Values are linearly interpolated between supplied years and held + constant (the nearest endpoint) outside the supplied range, so a + bundle that only covers, say, 1750-2100 still yields a file MAGICC + can read end to end. + """ + s = series.copy() + s.index = s.index.astype(int) + s = s.sort_index() + annual = pd.Index(range(MAGICC_CONC_FIRSTYEAR, MAGICC_CONC_LASTYEAR + 1)) + return ( + s.reindex(s.index.union(annual)) + .interpolate(method="index") + .reindex(annual) + .ffill() + .bfill() + ) + + +def write_conc_in_file( + out_path: str, + scenario: str, + species: str, + unit: str, + series: pd.Series, + magicc_version: str, +) -> str: + """ + Write a single MAGICC ``CONC.IN`` file via pymagicc and return its path. + + Wraps :class:`pymagicc.io.MAGICCData` with the right per-gas + metadata. The trajectory is first resampled onto MAGICC's annual + 1700-2500 grid (see :func:`_to_annual_magicc_grid`); pymagicc then + derives ``THISFILE_FIRSTYEAR`` / ``LASTYEAR`` / ``ANNUALSTEPS`` from + that index and sets ``THISFILE_REGIONMODE`` / ``THISFILE_DATTYPE`` + automatically for concentration files; we only need to set the + variable name, unit and region. + """ + if series.empty: + raise ValueError( + f"Cannot write concentration .IN for {scenario}/{species}: " + "trajectory is empty.", + ) + series = _to_annual_magicc_grid(series) + # pymagicc cross-checks the data variable against the one it parses + # back out of the filename, and it uses its own openscm naming + # (``HFC134a``, ``Halon1211``, ``CCl4``) rather than MAGICC's + # upper-case species token. Derive that canonical name so the check + # passes; for CO2/CH4/N2O it is identical to ``species``. + openscm_variable = pymagicc.definitions.convert_magicc7_to_openscm_variables( + f"{species}_CONC", + ) + frame = pd.DataFrame({ + "model": ["unspecified"], + "scenario": [scenario], + "region": ["World"], + "variable": [openscm_variable], + "unit": [unit], + "todo": ["SET"], + **{int(year): [float(value)] for year, value in series.items()}, + }) + writer = pymagicc.io.MAGICCData(frame) + writer.metadata = { + "header": ( + f"CONC.IN file written by openscm_runner for the " + f"{scenario} scenario, species {species}" + ), + } + writer.write(out_path, magicc_version=magicc_version) + return out_path diff --git a/src/openscm_runner/adapters/magicc7/_run_magicc_parallel.py b/src/openscm_runner/adapters/magicc7/_run_magicc_parallel.py index 10efdc6b..39aee249 100644 --- a/src/openscm_runner/adapters/magicc7/_run_magicc_parallel.py +++ b/src/openscm_runner/adapters/magicc7/_run_magicc_parallel.py @@ -107,6 +107,9 @@ def _run_func( scenario = cfg.pop("scenario") model = cfg.pop("model") output_config = cfg.pop("output_config") + # Adapter-internal dispatch keys (not MAGICC namelist entries): + cfg.pop("magicc_conc_driven", None) + cfg.pop("rcmip3_bundle_path", None) res = magicc.run(**cfg) if res.metadata["stderr"]: diff --git a/src/openscm_runner/adapters/magicc7/magicc7.py b/src/openscm_runner/adapters/magicc7/magicc7.py index e79d1622..08dd644f 100644 --- a/src/openscm_runner/adapters/magicc7/magicc7.py +++ b/src/openscm_runner/adapters/magicc7/magicc7.py @@ -7,10 +7,19 @@ from scmdata import ScmRun, run_append +from ..._run_mode import RunMode from ...progress import progress from ...settings import config from ..base import _Adapter from ._compat import pymagicc +from ._concentrations_translator import ( + FGAS_NAMES, + MHALO_NAMES, + build_concentrations_overlay, + cfg_keys_for_species, + classify_conc_species, + write_conc_in_file, +) from ._run_magicc_parallel import run_magicc_parallel LOGGER = logging.getLogger(__name__) @@ -42,6 +51,18 @@ def _convert_to_pymagicc_var(in_var): return out +def _with_mode_applied(cfg, mode): + """ + Inject the adapter-level ``mode`` into a cfg as the internal + ``magicc_conc_driven`` flag. + + Returns a new dict; the caller's cfg is not mutated. + """ + new = dict(cfg) + new["magicc_conc_driven"] = mode == RunMode.CONCENTRATION_DRIVEN + return new + + class MAGICC7(_Adapter): """ Adapter for running MAGICC7 @@ -51,6 +72,9 @@ class MAGICC7(_Adapter): """ model_name = "MAGICC7" + supported_modes = frozenset( + {RunMode.EMISSIONS_DRIVEN, RunMode.CONCENTRATION_DRIVEN} + ) def _init_model(self): # pylint:disable=arguments-differ if pymagicc is None: @@ -97,16 +121,35 @@ def _run(self, scenarios, cfgs, output_variables, output_config): # TODO: add use of historical data properly # pylint:disable=fixme LOGGER.warning("Historical data has not been checked") - magicc_df = scenarios.timeseries().reset_index() - magicc_df["variable"] = magicc_df["variable"].apply( + cfgs = [_with_mode_applied(cfg, self.mode) for cfg in cfgs] + + emissions_df = scenarios.timeseries().reset_index() + # Filter to emissions before the variable rename: the rename's + # ``HFC4310mee`` -> ``HFC4310`` substring substitution would + # otherwise mangle ``Atmospheric Concentrations|HFC4310mee`` + # rows too. Concentration rows are handled separately below + # via the conc-driven path. + emissions_df = emissions_df[ + emissions_df["variable"].str.startswith("Emissions|") + ].copy() + emissions_df["variable"] = emissions_df["variable"].apply( lambda x: x.replace("Sulfur", "SOx") .replace("HFC4310mee", "HFC4310") .replace("VOC", "NMVOC") ) - magicc_scmdf = self._convert_to_magicc_units(magicc_df) + magicc_scmdf = self._convert_to_magicc_units(emissions_df) full_cfgs = self._write_scen_files_and_make_full_cfgs(magicc_scmdf, cfgs) + if self.mode == RunMode.CONCENTRATION_DRIVEN: + conc_patches = self._write_conc_in_files_and_cfg_updates( + scenarios=scenarios, cfgs=cfgs, + ) + full_cfgs = [ + self._merge_conc_patch(cfg, conc_patches) + for cfg in full_cfgs + ] + pymagicc_vars = [_convert_to_pymagicc_var(v) for v in output_variables] res = run_magicc_parallel(full_cfgs, pymagicc_vars, output_config) @@ -146,6 +189,171 @@ def _fix_pint_incompatible_units(inp): return out + @staticmethod + def _merge_conc_patch(cfg, conc_patches): + """ + Layer the auto-generated conc cfg patch under the caller cfg. + + Caller-supplied ``*_switchfromconc2emis_year`` keys win over + the auto-generated default of 9999, mirroring the precedence + rule used elsewhere in the adapter (caller cfgs > built-in + scenario setup). + """ + patch = conc_patches.get((cfg["scenario"], cfg["model"]), {}) + if not patch: + return cfg + overrides = { + k: v for k, v in cfg.items() + if k.endswith("_switchfromconc2emis_year") + } + return {**cfg, **patch, **overrides} + + def _write_conc_in_files_and_cfg_updates( + self, scenarios, cfgs, out_directory=None, + ): + """ + Write per-(scenario, gas) concentration ``.IN`` files and + return per-(scenario, model) cfg patches that point MAGICC at + them. + + The decision is per-scenario: a gas is concentration-driven for + a scenario iff that scenario supplies a trajectory for it + (baseline or overlay); gases it doesn't supply remain + emissions-driven via the existing SCEN7 path. CO2/CH4/N2O use + per-gas ``file__conc`` + ``_switchfromconc2emis_year`` + flags (switch year defaults to 9999). F-gases / Montreal + halocarbons share the bundled-array flags ``fgas_files_conc`` / + ``mhalo_files_conc`` (positional, one slot per + :data:`FGAS_NAMES` / :data:`MHALO_NAMES` entry) with a single + shared ``*_switchfromconc2emis_year`` (defaults to 10000) — so + within a group conc-driving is all-or-nothing. Callers can + override any ``*_switchfromconc2emis_year`` per cfg (see + :meth:`_merge_conc_patch`). + """ + rcmip3_bundle_path = self._resolve_rcmip3_bundle_path(cfgs) + + if out_directory is None: + out_directory = os.path.join(self._run_dir(), "openscm-runner") + os.makedirs(out_directory, exist_ok=True) + + scenario_model_pairs = list( + scenarios.meta[["scenario", "model"]] + .drop_duplicates() + .itertuples(index=False, name=None) + ) + scenario_names = sorted({s for s, _ in scenario_model_pairs}) + overlay = build_concentrations_overlay( + scenarios, rcmip3_bundle_path, scenario_names, + ) + if overlay.empty: + return {} + + magicc_version = self.get_version()[1] + patches: dict = {} + for scenario, model in scenario_model_pairs: + scen_overlay = overlay[overlay["scenario"] == scenario] + if scen_overlay.empty: + continue + cfg_patch: dict = {} + fgas_written: dict[str, str] = {} + mhalo_written: dict[str, str] = {} + for _, row in scen_overlay.iterrows(): + species = row["magicc_species"] + file_name = ( + f"{scenario}_{model}_{species}_CONC.IN".upper() + .replace("/", "-") + .replace("\\", "-") + .replace(" ", "-") + ) + out_path = os.path.join(out_directory, file_name) + write_conc_in_file( + out_path=out_path, + scenario=scenario, + species=species, + unit=row["unit"], + series=row["series"], + magicc_version=magicc_version, + ) + kind = classify_conc_species(species) + if kind == "per_gas": + file_key, switch_key = cfg_keys_for_species(species) + cfg_patch[file_key] = out_path + cfg_patch[switch_key] = 9999 + elif kind == "fgas": + fgas_written[species] = out_path + elif kind == "mhalo": + mhalo_written[species] = out_path + + if fgas_written: + cfg_patch["fgas_files_conc"] = self._build_bundled_conc_array( + FGAS_NAMES, fgas_written, "F-gas", scenario, + ) + cfg_patch["fgas_switchfromconc2emis_year"] = 10000 + if mhalo_written: + cfg_patch["mhalo_files_conc"] = self._build_bundled_conc_array( + MHALO_NAMES, mhalo_written, "Montreal-halocarbon", + scenario, + ) + cfg_patch["mhalo_switchfromconc2emis_year"] = 10000 + + patches[(scenario, model)] = cfg_patch + return patches + + # Bundled-array slots MAGICC ships with no concentration file of its + # own (it falls back to a built-in default); leaving these empty is + # expected, not a sign of a missing trajectory. + _EXPECTED_EMPTY_CONC_SLOTS = frozenset({"HALON1202"}) + + @classmethod + def _build_bundled_conc_array(cls, names, written, group_label, scenario): + """ + Build a positional ``*_files_conc`` array over ``names``. + + ``written`` maps MAGICC species name -> CONC.IN path. Slots with + no written file become ``""`` (MAGICC uses its built-in default + for that species). Because the group shares a single switch year, + an empty slot is concentration-driven from MAGICC's default + rather than from emissions, so warn on any unexpected gap. + """ + array = [written.get(name, "") for name in names] + unexpected = [ + name + for name in names + if name not in written + and name not in cls._EXPECTED_EMPTY_CONC_SLOTS + ] + if unexpected: + LOGGER.warning( + "MAGICC7 conc-driven: scenario %s — %s group is " + "concentration-driven but %d species have no trajectory " + "(baseline or overlay) and will use MAGICC's built-in " + "default concentrations, not emissions: %s", + scenario, group_label, len(unexpected), sorted(unexpected), + ) + return array + + @staticmethod + def _resolve_rcmip3_bundle_path(cfgs): + """ + Pull the (required) ``rcmip3_bundle_path`` cfg key from the + first cfg; raise a clear error if no cfg supplies it. + + Mirrors the FaIRv2 / CICEROSCMPY2 adapter convention so + conc-driven runs are reproducible across machines (the + baseline does not depend on the binary's bundled + historical data). + """ + for cfg in cfgs: + if cfg.get("rcmip3_bundle_path"): + return cfg["rcmip3_bundle_path"] + raise ValueError( + "MAGICC7 concentration-driven mode requires the " + "``rcmip3_bundle_path`` cfg key. Pass the path to a local " + "copy of the RCMIP Phase 3 Zenodo bundle " + "(https://zenodo.org/records/20430630) or the bundle root " + "directory." + ) + def _write_scen_files_and_make_full_cfgs(self, scenarios, cfgs, out_directory=None): full_cfgs = [] run_id_block = 0 diff --git a/tests/integration/_ar6_validation.py b/tests/integration/_ar6_validation.py new file mode 100644 index 00000000..89d8f11f --- /dev/null +++ b/tests/integration/_ar6_validation.py @@ -0,0 +1,108 @@ +"""Helpers for validating MAGICC7 against IPCC AR6 Table 7.SM.4. + +The "Future Warming (GSAT)" SSP block of AR6 WG1 Table 7.SM.4 is a +concentration-driven, probabilistic reproduction target: each emulator +column is the 5th / 50th / 95th percentile of a calibrated ensemble, +driven from the SSP concentrations, expressed as GSAT relative to the +1995-2014 mean. For MAGICC7 the ensemble is the AR6 probabilistic +drawnset (600 parameter sets) run with MAGICC v7.5.x. + +These helpers wire that up so a test (or a notebook) can: + +- load the drawnset JSON into openscm-runner MAGICC7 cfgs, +- load the assessed reference table, and +- reduce a run result to per-period GSAT percentiles on the AR6 grid. + +The drawnset is licensed (see its README) and is **not** vendored into +this repo; point the loader at a local copy. Reproducing the published +numbers also needs the full RCMIP3 concentration bundle; with only the +mini test bundle this exercises the pipeline (a "smoke" run), not the +published values. +""" +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np +import pandas as pd + +GSAT_VARIABLE = "Surface Air Temperature Change" +AR6_REFERENCE_PERIOD = range(1995, 2015) # 1995-2014 inclusive +AR6_PERIODS: dict[str, range] = { + "2021-2040": range(2021, 2041), + "2041-2060": range(2041, 2061), + "2081-2100": range(2081, 2101), +} +AR6_QUANTILES = (0.05, 0.50, 0.95) # -> (lower, central, upper) + + +def load_ar6_drawnset(path, n: int | None = None) -> list[dict]: + """Load the AR6 probabilistic drawnset JSON into MAGICC7 cfgs. + + Each drawnset member is ``{"nml_allcfgs": {...}, "paraset_id": int}``; + we lower-case the namelist keys (the adapter's convention) and use + ``paraset_id`` as ``run_id`` so ensemble members stay distinguishable + through :func:`openscm_runner.utils.calculate_quantiles`. + + Parameters + ---------- + path + Path to the drawnset ``.json``. + n + If given, take only the first ``n`` members (for smoke runs). + """ + members = json.loads(Path(path).read_text())["configurations"] + if n is not None: + members = members[:n] + cfgs = [] + for member in members: + cfg = {k.lower(): v for k, v in member["nml_allcfgs"].items()} + cfg["run_id"] = member["paraset_id"] + cfgs.append(cfg) + return cfgs + + +def load_table_7sm4(path) -> pd.DataFrame: + """Load the cleaned AR6 Table 7.SM.4 reference fixture.""" + return pd.read_csv(path) + + +def assessed_row(table: pd.DataFrame, metric: str, source: str, period_startswith: str = ""): + """Return ``(lower, central, upper)`` for one metric/source/period.""" + sel = table[(table["metric"] == metric) & (table["source"] == source)] + if period_startswith: + sel = sel[sel["period"].fillna("").str.startswith(period_startswith)] + if len(sel) != 1: + raise LookupError( + f"expected exactly one row for metric={metric!r} source={source!r} " + f"period~{period_startswith!r}, got {len(sel)}" + ) + row = sel.iloc[0] + return float(row["lower"]), float(row["central"]), float(row["upper"]) + + +def gsat_percentiles_by_period( + res, + periods: dict[str, range] = AR6_PERIODS, + reference_period: range = AR6_REFERENCE_PERIOD, + quantiles=AR6_QUANTILES, +) -> dict[str, tuple[float, float, float]]: + """Reduce a MAGICC run to per-period GSAT percentiles on the AR6 grid. + + Rebases ``Surface Air Temperature Change`` to the ``reference_period`` + mean (matching AR6's "relative to 1995-2014"), takes each ensemble + member's mean warming over each period, then the cross-member + percentiles. Returns ``{period_label: (lower, central, upper)}``. + """ + rebased = res.filter(variable=GSAT_VARIABLE).relative_to_ref_period_mean( + year=list(reference_period), + ) + ts = rebased.timeseries() + out: dict[str, tuple[float, float, float]] = {} + for label, years in periods.items(): + year_cols = [c for c in ts.columns if int(getattr(c, "year", c)) in years] + per_member = ts[year_cols].mean(axis=1).to_numpy() + lo, ce, up = np.percentile(per_member, [q * 100 for q in quantiles]) + out[label] = (float(lo), float(ce), float(up)) + return out diff --git a/tests/integration/test_magicc7.py b/tests/integration/test_magicc7.py index 3df86cf6..4f3e04c0 100644 --- a/tests/integration/test_magicc7.py +++ b/tests/integration/test_magicc7.py @@ -1,14 +1,26 @@ import os.path +from pathlib import Path +import pandas as pd import pymagicc.io import pytest from scmdata import ScmRun import openscm_runner.run +from openscm_runner import RunMode from openscm_runner.adapters import MAGICC7 from openscm_runner.testing import _AdapterTester from openscm_runner.utils import calculate_quantiles +RCMIP3_MINI_BUNDLE = ( + Path(__file__).parent.parent / "test-data" / "rcmip3-mini" +) +# Variant bundle that also ships a few F-gas / Montreal-halocarbon ssp245 +# concentration trajectories, for exercising the bundled-array path. +RCMIP3_MINI_HALO_BUNDLE = ( + Path(__file__).parent.parent / "test-data" / "rcmip3-mini-halo" +) + @pytest.mark.magicc class TestMagicc7Adapter(_AdapterTester): @@ -209,6 +221,243 @@ def test_variable_naming(self, test_scenarios): raise AssertionError(missing_vars) +@pytest.mark.magicc +def test_conc_driven_writes_conc_in_files_and_patches_cfgs( + tmp_path, monkeypatch, +): + """Focused test of the conc-file-writing path. + + Builds an adapter in CONCENTRATION_DRIVEN mode, runs + ``_write_conc_in_files_and_cfg_updates`` directly with + ``out_directory=tmp_path`` (so we don't need ``_run_dir()`` and a + real MAGICC install), and verifies that: + + - The per-(scenario, model) patch sets ``file_co2_conc``, + ``file_ch4_conc``, ``file_n2o_conc`` to paths of written + ``.IN`` files (RCMIP3-mini's ssp245 baseline supplies all + three). + - The matching ``*_switchfromconc2emis_year`` flags default to + 9999. + - User-supplied ``Atmospheric Concentrations|CO2`` rows shadow + the baseline CO2 trajectory in the written file. + """ + monkeypatch.setattr( + MAGICC7, "get_version", classmethod(lambda cls: "v7.5.3"), + ) + + user_run = ScmRun(pd.DataFrame({ + "model": ["test-model"], + "scenario": ["ssp245"], + "region": ["World"], + "variable": ["Atmospheric Concentrations|CO2"], + "unit": ["ppm"], + 2050: [999.0], + 2100: [999.0], + })) + + adapter = MAGICC7( + cfgs=[ + { + "core_climatesensitivity": 3, + "rcmip3_bundle_path": str(RCMIP3_MINI_BUNDLE), + "scenario": "ssp245", + "model": "test-model", + }, + ], + mode=RunMode.CONCENTRATION_DRIVEN, + output_variables=("Surface Air Temperature Change",), + ) + + patches = adapter._write_conc_in_files_and_cfg_updates( + scenarios=user_run, + cfgs=adapter.cfgs, + out_directory=str(tmp_path), + ) + + assert ("ssp245", "test-model") in patches + patch = patches[("ssp245", "test-model")] + for gas in ("co2", "ch4", "n2o"): + assert patch[f"{gas}_switchfromconc2emis_year"] == 9999 + path = patch[f"file_{gas}_conc"] + assert os.path.exists(path), path + assert path.endswith(f"_{gas.upper()}_CONC.IN") + + # User CO2 trajectory survived into the written file. + co2_data = pymagicc.io.MAGICCData(patch["file_co2_conc"]) + co2_ts = co2_data.timeseries(time_axis="year") + assert co2_ts.iloc[0].loc[2050] == pytest.approx(999.0) + assert co2_ts.iloc[0].loc[2100] == pytest.approx(999.0) + + # MAGICC reads CONC.IN as a contiguous annual series over its full + # internal window; a sparse / truncated file trips a Fortran + # end-of-file error at runtime. Assert the writer resampled the + # (here sparse) bundle onto the annual 1700-2500 grid. + years = co2_ts.columns.astype(int) + assert years.min() == 1700 + assert years.max() == 2500 + assert list(years) == list(range(1700, 2501)) + # Values are held constant beyond the last supplied year (2100). + assert co2_ts.iloc[0].loc[2500] == pytest.approx(999.0) + + +@pytest.mark.magicc +def test_conc_driven_writes_fgas_mhalo_bundled_arrays(tmp_path, monkeypatch): + """F-gases / Montreal halocarbons drive via the bundled-array flags. + + The halo bundle supplies ssp245 trajectories for HFC134a (F-gas), + SF6 (F-gas) and CFC12 (Montreal halocarbon). Verify the per-(scenario, + model) patch builds the positional ``fgas_files_conc`` / + ``mhalo_files_conc`` arrays in ``FGAS_NAMES`` / ``MHALO_NAMES`` order, + fills only the supplied slots, leaves the rest (incl. HALON1202) + empty, and sets the shared group switch year to 10000. + """ + from openscm_runner.adapters.magicc7._concentrations_translator import ( + FGAS_NAMES, + MHALO_NAMES, + ) + + monkeypatch.setattr( + MAGICC7, "get_version", classmethod(lambda cls: "v7.5.3"), + ) + + adapter = MAGICC7( + cfgs=[ + { + "core_climatesensitivity": 3, + "rcmip3_bundle_path": str(RCMIP3_MINI_HALO_BUNDLE), + "scenario": "ssp245", + "model": "test-model", + }, + ], + mode=RunMode.CONCENTRATION_DRIVEN, + output_variables=("Surface Air Temperature Change",), + ) + + # No user overlay: drive purely from the halo bundle baseline. + empty_run = ScmRun(pd.DataFrame({ + "model": ["test-model"], + "scenario": ["ssp245"], + "region": ["World"], + "variable": ["Emissions|CO2|MAGICC Fossil and Industrial"], + "unit": ["GtC / yr"], + 2050: [10.0], + 2100: [10.0], + })) + + patches = adapter._write_conc_in_files_and_cfg_updates( + scenarios=empty_run, + cfgs=adapter.cfgs, + out_directory=str(tmp_path), + ) + patch = patches[("ssp245", "test-model")] + + # F-gas group: positional array over all 23 names, switch year 10000. + fgas_array = patch["fgas_files_conc"] + assert len(fgas_array) == len(FGAS_NAMES) == 23 + assert patch["fgas_switchfromconc2emis_year"] == 10000 + sf6_slot = fgas_array[FGAS_NAMES.index("SF6")] + hfc134a_slot = fgas_array[FGAS_NAMES.index("HFC134A")] + assert sf6_slot.endswith("_SF6_CONC.IN") and os.path.exists(sf6_slot) + assert hfc134a_slot.endswith("_HFC134A_CONC.IN") + assert os.path.exists(hfc134a_slot) + # Unsupplied F-gas slots are empty (fall back to MAGICC defaults): + assert fgas_array[FGAS_NAMES.index("CF4")] == "" + + # Montreal-halocarbon group: 18 names, CFC12 filled, HALON1202 empty. + mhalo_array = patch["mhalo_files_conc"] + assert len(mhalo_array) == len(MHALO_NAMES) == 18 + assert patch["mhalo_switchfromconc2emis_year"] == 10000 + cfc12_slot = mhalo_array[MHALO_NAMES.index("CFC12")] + assert cfc12_slot.endswith("_CFC12_CONC.IN") and os.path.exists(cfc12_slot) + assert mhalo_array[MHALO_NAMES.index("HALON1202")] == "" + + +@pytest.mark.magicc +def test_fgas_mhalo_name_lists_match_binary_config(): + """The hardcoded FGAS_NAMES / MHALO_NAMES must track the installed + binary's namelist order; the positional arrays bind to it.""" + from openscm_runner.adapters.magicc7._compat import f90nml + from openscm_runner.adapters.magicc7._concentrations_translator import ( + FGAS_NAMES, + MHALO_NAMES, + ) + + run_dir = Path(MAGICC7()._run_dir()) + defaultall = run_dir / "MAGCFG_DEFAULTALL.CFG" + if not defaultall.exists(): + pytest.skip(f"MAGCFG_DEFAULTALL.CFG not found at {defaultall}") + + cfg = f90nml.read(str(defaultall))["nml_allcfgs"] + binary_fgas = tuple(s.strip().upper() for s in cfg["fgas_names"]) + binary_mhalo = tuple(s.strip().upper() for s in cfg["mhalo_names"]) + assert FGAS_NAMES == binary_fgas + assert MHALO_NAMES == binary_mhalo + + +@pytest.mark.magicc +def test_conc_driven_end_to_end_runs_the_binary(test_scenarios): + """Drive the real MAGICC binary in concentration-driven mode. + + The file-writing tests above bypass the binary, so they can't catch + a malformed ``CONC.IN`` (e.g. a sparse / truncated grid that trips + MAGICC's Fortran ``readdata`` end-of-file check). This runs ssp245 + through the binary in both modes and asserts conc-driven produces a + finite GSAT in the same ballpark as emissions-driven. It also checks + that the prescribed CO2 concentration round-trips: requested as an + output, MAGICC should echo what was written into ``CO2_CONC.IN``. + """ + scenarios = test_scenarios.filter(scenario="ssp245") + + gsat = {} + for mode in (RunMode.EMISSIONS_DRIVEN, RunMode.CONCENTRATION_DRIVEN): + adapter = MAGICC7( + cfgs=[ + { + "core_climatesensitivity": 3, + "rcmip3_bundle_path": str(RCMIP3_MINI_BUNDLE), + }, + ], + mode=mode, + output_variables=( + "Surface Air Temperature Change", + "Atmospheric Concentrations|CO2", + ), + ) + res = openscm_runner.run.run([adapter], scenarios=scenarios) + vals = res.filter( + variable="Surface Air Temperature Change", year=2100, + ).values.flatten() + assert len(vals) > 0 + assert all(pd.notna(vals)) + gsat[mode] = float(vals[0]) + + if mode == RunMode.CONCENTRATION_DRIVEN: + # Round-trip: the prescribed CO2 concentration (rcmip3-mini + # ssp245 baseline, ~602.78 ppm at 2100) should come back out. + co2_out = res.filter( + variable="Atmospheric Concentrations|CO2", year=2100, + ).values.flatten() + assert len(co2_out) > 0 + assert co2_out[0] == pytest.approx(602.78, rel=1e-2) + + # Physically-consistent bundle: the two modes should land close. + assert gsat[RunMode.CONCENTRATION_DRIVEN] == pytest.approx( + gsat[RunMode.EMISSIONS_DRIVEN], abs=0.5, + ) + + +@pytest.mark.magicc +def test_conc_driven_requires_rcmip3_bundle_path(): + """No rcmip3_bundle_path -> clear ValueError on conc-driven dispatch.""" + adapter = MAGICC7( + cfgs=[{"core_climatesensitivity": 3}], + mode=RunMode.CONCENTRATION_DRIVEN, + output_variables=("Surface Air Temperature Change",), + ) + with pytest.raises(ValueError, match="rcmip3_bundle_path"): + adapter._resolve_rcmip3_bundle_path(adapter.cfgs) + + @pytest.mark.magicc def test_write_scen_files_and_make_full_cfgs(test_scenarios): adapter = MAGICC7() diff --git a/tests/integration/test_magicc7_ar6_validation.py b/tests/integration/test_magicc7_ar6_validation.py new file mode 100644 index 00000000..1cd9e527 --- /dev/null +++ b/tests/integration/test_magicc7_ar6_validation.py @@ -0,0 +1,125 @@ +"""AR6 Table 7.SM.4 validation for MAGICC7 concentration-driven runs. + +This is a *smoke* harness: it proves the end-to-end pipeline — load the +AR6 probabilistic drawnset, run MAGICC7 concentration-driven, rebase GSAT +to 1995-2014, and reduce to per-period percentiles on the AR6 grid — and +compares the result against the published MAGICC7 column. + +It is gated on a local drawnset (licensed, not vendored) and the MAGICC +binary, so it is skipped by default. To run it:: + + AR6_MAGICC_DRAWNSET=/path/to/...drawnset.json \ + MAGICC_EXECUTABLE_7=/path/to/magicc \ + pytest -m magicc tests/integration/test_magicc7_ar6_validation.py -s + +With only the mini RCMIP3 bundle (the default) the run is +concentration-driven on CO2/CH4/N2O with the remaining species from +ssp245 emissions, so the numbers are a ballpark, not the published +values. Point ``AR6_RCMIP3_BUNDLE`` at the full RCMIP Phase 3 bundle to +scale this up toward a faithful reproduction. + +With the full bundle and the full 600-member drawnset the GSAT *medians* +reproduce the published MAGICC7 column to <=0.01 degC (e.g. SSP2-4.5 +2081-2100: 1.82 vs 1.82). The 5th/95th percentiles run wider — the +drawnset is the AR6 *prior*, whereas the published 7.SM.4 ranges come +from the observationally-constrained (weighted) distribution — so the +test validates the median and only reports the tails. +""" +from __future__ import annotations + +import os +from pathlib import Path + +import pytest + +import openscm_runner.run +from openscm_runner import RunMode + +from _ar6_validation import ( # noqa: E402 (pytest prepends the test dir to sys.path) + GSAT_VARIABLE, + assessed_row, + gsat_percentiles_by_period, + load_ar6_drawnset, + load_table_7sm4, +) + +RCMIP3_MINI_BUNDLE = ( + Path(__file__).parent.parent / "test-data" / "rcmip3-mini" +) +TABLE_7SM4 = ( + Path(__file__).parent.parent / "test-data" / "ar6" / "table_7_SM_4.csv" +) + +_DRAWNSET = os.environ.get("AR6_MAGICC_DRAWNSET") +_SMOKE_N = int(os.environ.get("AR6_SMOKE_N", "10")) + + +@pytest.mark.magicc +@pytest.mark.skipif( + not _DRAWNSET or not Path(_DRAWNSET).exists(), + reason="set AR6_MAGICC_DRAWNSET to a local AR6 probabilistic drawnset JSON", +) +def test_ssp245_gsat_matches_ar6_table_7sm4_ballpark(test_scenarios): + bundle = os.environ.get("AR6_RCMIP3_BUNDLE", str(RCMIP3_MINI_BUNDLE)) + full_bundle = "AR6_RCMIP3_BUNDLE" in os.environ + + cfgs = load_ar6_drawnset(_DRAWNSET, n=_SMOKE_N) + for cfg in cfgs: + cfg["rcmip3_bundle_path"] = bundle + + scenarios = test_scenarios.filter(scenario="ssp245") + res = openscm_runner.run.run( + climate_models_cfgs={"MAGICC7": tuple(cfgs)}, + scenarios=scenarios, + output_variables=(GSAT_VARIABLE,), + mode=RunMode.CONCENTRATION_DRIVEN, + ) + + members = len(res.filter(variable=GSAT_VARIABLE).get_unique_meta("run_id")) + assert members >= max(2, _SMOKE_N // 2), members + + pct = gsat_percentiles_by_period(res) + table = load_table_7sm4(TABLE_7SM4) + + # Report card: our percentiles vs the published MAGICC7 column. + print(f"\nAR6 Table 7.SM.4 — MAGICC7 SSP2-4.5 GSAT (rel 1995-2014), " + f"{members} members, bundle={'FULL' if full_bundle else 'mini/WMGHG-only'}") + print(f"{'period':>10} | {'ours (5/50/95)':>26} | {'AR6 MAGICC7 (5/50/95)':>26}") + for period, (lo, ce, up) in pct.items(): + a_lo, a_ce, a_up = assessed_row(table, "SSP2-4.5", "MAGICC7", period) + print(f"{period:>10} | {lo:7.2f}{ce:8.2f}{up:8.2f} | " + f"{a_lo:7.2f}{a_ce:8.2f}{a_up:8.2f}") + + # --- mechanics (always asserted) --- + for period, (lo, ce, up) in pct.items(): + assert lo <= ce <= up, (period, lo, ce, up) + assert all(map(_finite, (lo, ce, up))), (period, lo, ce, up) + + # --- value check: assert on the MEDIAN, report the tails --- + # The drawnset is the AR6 prior; the published 7.SM.4 ranges come from + # the observationally-*constrained* (weighted) distribution. Running it + # unweighted reproduces the median almost exactly but leaves a heavier + # upper tail (the 95th sits high), so we validate against the median and + # only report the 5th/95th. Reproducing the published range would need + # AR6's constraint weights. + medians = {p: ce for p, (_, ce, _) in pct.items()} + if full_bundle and members >= 100: + # Faithful inputs + a real ensemble: the median should land on the + # published MAGICC7 central across every period. + for period in pct: + _, ar6_central, _ = assessed_row(table, "SSP2-4.5", "MAGICC7", period) + assert medians[period] == pytest.approx(ar6_central, abs=0.1), ( + period, medians[period], ar6_central, + ) + elif full_bundle: + # Full conc inputs but a small ensemble: median is a ballpark. + _, ar6_central, _ = assessed_row(table, "SSP2-4.5", "MAGICC7", "2081-2100") + assert medians["2081-2100"] == pytest.approx(ar6_central, abs=0.3) + else: + # Smoke inputs (WMGHG-only conc): only a loose sanity ballpark. + _, ar6_central, _ = assessed_row(table, "SSP2-4.5", "MAGICC7", "2081-2100") + assert medians["2081-2100"] == pytest.approx(ar6_central, abs=1.0) + + +def _finite(x: float) -> bool: + return x == x and abs(x) != float("inf") diff --git a/tests/test-data/ar6/table_7_SM_4.csv b/tests/test-data/ar6/table_7_SM_4.csv new file mode 100644 index 00000000..74c9d651 --- /dev/null +++ b/tests/test-data/ar6/table_7_SM_4.csv @@ -0,0 +1,111 @@ +category,metric,unit,period,source,lower,central,upper +Key Metrics,ECS,degC,,Assessed,2.0,3.0,5.0 +Key Metrics,ECS,degC,,CICERO-SCM,2.53,3.05,4.09 +Key Metrics,ECS,degC,,FaIRv1.6.2,2.05,2.95,5.07 +Key Metrics,ECS,degC,,MAGICC7,1.93,2.97,4.83 +Key Metrics,ECS,degC,,OSCARv3.1.1,1.84,2.54,3.9 +Key Metrics,TCRE,degC per 1000 GtC,,Assessed,1.0,1.65,2.3 +Key Metrics,TCRE,degC per 1000 GtC,,CICERO-SCM,,, +Key Metrics,TCRE,degC per 1000 GtC,,FaIRv1.6.2,1.29,1.53,1.82 +Key Metrics,TCRE,degC per 1000 GtC,,MAGICC7,1.37,1.73,2.19 +Key Metrics,TCRE,degC per 1000 GtC,,OSCARv3.1.1,1.5,1.52,1.83 +Key Metrics,TCR,degC,,Assessed,1.2,1.8,2.4 +Key Metrics,TCR,degC,,CICERO-SCM,1.38,1.71,2.32 +Key Metrics,TCR,degC,,FaIRv1.6.2,1.36,1.81,2.46 +Key Metrics,TCR,degC,,MAGICC7,1.27,1.88,2.61 +Key Metrics,TCR,degC,,OSCARv3.1.1,1.51,1.82,2.05 +Historical,Ocean Heat Content Change,ZJ,1971-2018,Assessed,329,396,463 +Historical,Ocean Heat Content Change,ZJ,1971-2018,CICERO-SCM,250,288,329 +Historical,Ocean Heat Content Change,ZJ,1971-2018,FaIRv1.6.2,346,381,423 +Historical,Ocean Heat Content Change,ZJ,1971-2018,MAGICC7,325,382,436 +Historical,Ocean Heat Content Change,ZJ,1971-2018,OSCARv3.1.1,174,243,508 +Historical,Total Aerosol ERF,W m-2,2005-2014 rel 1750,Assessed,-2.0,-1.3,-0.6 +Historical,Total Aerosol ERF,W m-2,2005-2014 rel 1750,CICERO-SCM,-1.27,-0.82,-0.54 +Historical,Total Aerosol ERF,W m-2,2005-2014 rel 1750,FaIRv1.6.2,-1.68,-1.15,-0.6 +Historical,Total Aerosol ERF,W m-2,2005-2014 rel 1750,MAGICC7,-1.79,-1.2,-0.55 +Historical,Total Aerosol ERF,W m-2,2005-2014 rel 1750,OSCARv3.1.1,-1.24,-1.11,-0.79 +Historical,WMGHG ERF,W m-2,2019 rel 1750,Assessed,3.03,3.32,3.61 +Historical,WMGHG ERF,W m-2,2019 rel 1750,CICERO-SCM,3.14,3.14,3.14 +Historical,WMGHG ERF,W m-2,2019 rel 1750,FaIRv1.6.2,3.07,3.38,3.66 +Historical,WMGHG ERF,W m-2,2019 rel 1750,MAGICC7,3.1,3.35,3.6 +Historical,WMGHG ERF,W m-2,2019 rel 1750,OSCARv3.1.1,3.06,3.42,3.49 +Historical,Methane ERF,W m-2,2019 rel 1750,Assessed,0.43,0.54,0.65 +Historical,Methane ERF,W m-2,2019 rel 1750,CICERO-SCM,0.56,0.56,0.56 +Historical,Methane ERF,W m-2,2019 rel 1750,FaIRv1.6.2,0.44,0.56,0.67 +Historical,Methane ERF,W m-2,2019 rel 1750,MAGICC7,0.43,0.54,0.67 +Historical,Methane ERF,W m-2,2019 rel 1750,OSCARv3.1.1,0.47,0.54,0.62 +Future Warming GSAT,SSP1-1.9,degC,2021-2040 rel 1995-2014,Assessed,0.38,0.61,0.85 +Future Warming GSAT,SSP1-1.9,degC,2021-2040 rel 1995-2014,CICERO-SCM,0.42,0.58,0.94 +Future Warming GSAT,SSP1-1.9,degC,2021-2040 rel 1995-2014,FaIRv1.6.2,0.39,0.61,0.94 +Future Warming GSAT,SSP1-1.9,degC,2021-2040 rel 1995-2014,MAGICC7,0.39,0.61,0.88 +Future Warming GSAT,SSP1-1.9,degC,2021-2040 rel 1995-2014,OSCARv3.1.1,0.43,0.55,0.64 +Future Warming GSAT,SSP1-1.9,degC,2041-2060 rel 1995-2014,Assessed,0.4,0.71,1.07 +Future Warming GSAT,SSP1-1.9,degC,2041-2060 rel 1995-2014,CICERO-SCM,0.43,0.65,1.15 +Future Warming GSAT,SSP1-1.9,degC,2041-2060 rel 1995-2014,FaIRv1.6.2,0.36,0.66,1.14 +Future Warming GSAT,SSP1-1.9,degC,2041-2060 rel 1995-2014,MAGICC7,0.39,0.71,1.15 +Future Warming GSAT,SSP1-1.9,degC,2041-2060 rel 1995-2014,OSCARv3.1.1,0.45,0.65,0.74 +Future Warming GSAT,SSP1-1.9,degC,2081-2100 rel 1995-2014,Assessed,0.24,0.56,0.96 +Future Warming GSAT,SSP1-1.9,degC,2081-2100 rel 1995-2014,CICERO-SCM,0.21,0.42,0.94 +Future Warming GSAT,SSP1-1.9,degC,2081-2100 rel 1995-2014,FaIRv1.6.2,0.18,0.48,1.0 +Future Warming GSAT,SSP1-1.9,degC,2081-2100 rel 1995-2014,MAGICC7,0.2,0.52,0.99 +Future Warming GSAT,SSP1-1.9,degC,2081-2100 rel 1995-2014,OSCARv3.1.1,0.26,0.51,0.66 +Future Warming GSAT,SSP1-2.6,degC,2021-2040 rel 1995-2014,Assessed,0.41,0.63,0.89 +Future Warming GSAT,SSP1-2.6,degC,2021-2040 rel 1995-2014,CICERO-SCM,0.44,0.6,0.94 +Future Warming GSAT,SSP1-2.6,degC,2021-2040 rel 1995-2014,FaIRv1.6.2,0.42,0.64,0.96 +Future Warming GSAT,SSP1-2.6,degC,2021-2040 rel 1995-2014,MAGICC7,0.4,0.62,0.89 +Future Warming GSAT,SSP1-2.6,degC,2021-2040 rel 1995-2014,OSCARv3.1.1,0.45,0.57,0.64 +Future Warming GSAT,SSP1-2.6,degC,2041-2060 rel 1995-2014,Assessed,0.54,0.88,1.32 +Future Warming GSAT,SSP1-2.6,degC,2041-2060 rel 1995-2014,CICERO-SCM,0.58,0.83,1.34 +Future Warming GSAT,SSP1-2.6,degC,2041-2060 rel 1995-2014,FaIRv1.6.2,0.53,0.86,1.38 +Future Warming GSAT,SSP1-2.6,degC,2041-2060 rel 1995-2014,MAGICC7,0.54,0.89,1.35 +Future Warming GSAT,SSP1-2.6,degC,2041-2060 rel 1995-2014,OSCARv3.1.1,0.62,0.82,0.95 +Future Warming GSAT,SSP1-2.6,degC,2081-2100 rel 1995-2014,Assessed,0.51,0.9,1.48 +Future Warming GSAT,SSP1-2.6,degC,2081-2100 rel 1995-2014,CICERO-SCM,0.5,0.78,1.41 +Future Warming GSAT,SSP1-2.6,degC,2081-2100 rel 1995-2014,FaIRv1.6.2,0.47,0.84,1.49 +Future Warming GSAT,SSP1-2.6,degC,2081-2100 rel 1995-2014,MAGICC7,0.48,0.89,1.49 +Future Warming GSAT,SSP1-2.6,degC,2081-2100 rel 1995-2014,OSCARv3.1.1,0.59,0.82,1.05 +Future Warming GSAT,SSP2-4.5,degC,2021-2040 rel 1995-2014,Assessed,0.44,0.66,0.9 +Future Warming GSAT,SSP2-4.5,degC,2021-2040 rel 1995-2014,CICERO-SCM,0.48,0.63,0.94 +Future Warming GSAT,SSP2-4.5,degC,2021-2040 rel 1995-2014,FaIRv1.6.2,0.47,0.65,0.92 +Future Warming GSAT,SSP2-4.5,degC,2021-2040 rel 1995-2014,MAGICC7,0.45,0.64,0.89 +Future Warming GSAT,SSP2-4.5,degC,2021-2040 rel 1995-2014,OSCARv3.1.1,0.42,0.57,0.63 +Future Warming GSAT,SSP2-4.5,degC,2041-2060 rel 1995-2014,Assessed,0.78,1.12,1.57 +Future Warming GSAT,SSP2-4.5,degC,2041-2060 rel 1995-2014,CICERO-SCM,0.81,1.08,1.62 +Future Warming GSAT,SSP2-4.5,degC,2041-2060 rel 1995-2014,FaIRv1.6.2,0.79,1.11,1.59 +Future Warming GSAT,SSP2-4.5,degC,2041-2060 rel 1995-2014,MAGICC7,0.79,1.13,1.6 +Future Warming GSAT,SSP2-4.5,degC,2041-2060 rel 1995-2014,OSCARv3.1.1,0.84,1.03,1.12 +Future Warming GSAT,SSP2-4.5,degC,2081-2100 rel 1995-2014,Assessed,1.24,1.81,2.59 +Future Warming GSAT,SSP2-4.5,degC,2081-2100 rel 1995-2014,CICERO-SCM,1.22,1.63,2.51 +Future Warming GSAT,SSP2-4.5,degC,2081-2100 rel 1995-2014,FaIRv1.6.2,1.21,1.75,2.63 +Future Warming GSAT,SSP2-4.5,degC,2081-2100 rel 1995-2014,MAGICC7,1.21,1.82,2.67 +Future Warming GSAT,SSP2-4.5,degC,2081-2100 rel 1995-2014,OSCARv3.1.1,1.34,1.74,1.96 +Future Warming GSAT,SSP3-7.0,degC,2021-2040 rel 1995-2014,Assessed,0.45,0.67,0.92 +Future Warming GSAT,SSP3-7.0,degC,2021-2040 rel 1995-2014,CICERO-SCM,0.5,0.64,0.93 +Future Warming GSAT,SSP3-7.0,degC,2021-2040 rel 1995-2014,FaIRv1.6.2,0.51,0.68,0.91 +Future Warming GSAT,SSP3-7.0,degC,2021-2040 rel 1995-2014,MAGICC7,0.49,0.68,0.92 +Future Warming GSAT,SSP3-7.0,degC,2021-2040 rel 1995-2014,OSCARv3.1.1,0.43,0.57,0.65 +Future Warming GSAT,SSP3-7.0,degC,2041-2060 rel 1995-2014,Assessed,0.92,1.28,1.75 +Future Warming GSAT,SSP3-7.0,degC,2041-2060 rel 1995-2014,CICERO-SCM,0.96,1.22,1.74 +Future Warming GSAT,SSP3-7.0,degC,2041-2060 rel 1995-2014,FaIRv1.6.2,0.98,1.28,1.72 +Future Warming GSAT,SSP3-7.0,degC,2041-2060 rel 1995-2014,MAGICC7,0.98,1.33,1.77 +Future Warming GSAT,SSP3-7.0,degC,2041-2060 rel 1995-2014,OSCARv3.1.1,0.99,1.17,1.29 +Future Warming GSAT,SSP3-7.0,degC,2081-2100 rel 1995-2014,Assessed,2.0,2.76,3.75 +Future Warming GSAT,SSP3-7.0,degC,2081-2100 rel 1995-2014,CICERO-SCM,1.99,2.55,3.64 +Future Warming GSAT,SSP3-7.0,degC,2081-2100 rel 1995-2014,FaIRv1.6.2,2.07,2.72,3.72 +Future Warming GSAT,SSP3-7.0,degC,2081-2100 rel 1995-2014,MAGICC7,2.13,2.86,3.97 +Future Warming GSAT,SSP3-7.0,degC,2081-2100 rel 1995-2014,OSCARv3.1.1,2.09,2.59,2.81 +Future Warming GSAT,SSP5-8.5,degC,2021-2040 rel 1995-2014,Assessed,0.51,0.76,1.04 +Future Warming GSAT,SSP5-8.5,degC,2021-2040 rel 1995-2014,CICERO-SCM,0.54,0.71,1.06 +Future Warming GSAT,SSP5-8.5,degC,2021-2040 rel 1995-2014,FaIRv1.6.2,0.56,0.77,1.08 +Future Warming GSAT,SSP5-8.5,degC,2021-2040 rel 1995-2014,MAGICC7,0.55,0.77,1.06 +Future Warming GSAT,SSP5-8.5,degC,2021-2040 rel 1995-2014,OSCARv3.1.1,0.51,0.66,0.73 +Future Warming GSAT,SSP5-8.5,degC,2041-2060 rel 1995-2014,Assessed,1.08,1.54,2.08 +Future Warming GSAT,SSP5-8.5,degC,2041-2060 rel 1995-2014,CICERO-SCM,1.11,1.42,2.07 +Future Warming GSAT,SSP5-8.5,degC,2041-2060 rel 1995-2014,FaIRv1.6.2,1.12,1.55,2.17 +Future Warming GSAT,SSP5-8.5,degC,2041-2060 rel 1995-2014,MAGICC7,1.11,1.57,2.16 +Future Warming GSAT,SSP5-8.5,degC,2041-2060 rel 1995-2014,OSCARv3.1.1,1.19,1.44,1.58 +Future Warming GSAT,SSP5-8.5,degC,2081-2100 rel 1995-2014,Assessed,2.44,3.5,4.82 +Future Warming GSAT,SSP5-8.5,degC,2081-2100 rel 1995-2014,CICERO-SCM,2.54,3.24,4.68 +Future Warming GSAT,SSP5-8.5,degC,2081-2100 rel 1995-2014,FaIRv1.6.2,2.58,3.5,4.89 +Future Warming GSAT,SSP5-8.5,degC,2081-2100 rel 1995-2014,MAGICC7,2.63,3.65,5.16 +Future Warming GSAT,SSP5-8.5,degC,2081-2100 rel 1995-2014,OSCARv3.1.1,2.65,3.35,3.62 diff --git a/tests/test-data/rcmip3-mini-halo/rcmip_phase3_concentrations_v2.0.0.csv b/tests/test-data/rcmip3-mini-halo/rcmip_phase3_concentrations_v2.0.0.csv new file mode 100644 index 00000000..d66c8264 --- /dev/null +++ b/tests/test-data/rcmip3-mini-halo/rcmip_phase3_concentrations_v2.0.0.csv @@ -0,0 +1,10 @@ +Model,Scenario,Region,Variable,Unit,Activity_Id,Type,Priority,Mip_Era,Version,1750,1850,1900,2000,2050,2100 +unspecified,historical,World,Atmospheric Concentrations|CH4,ppb,input4MIPs,non-idealised,1,CMIP6,RCMIP Phase 3 v2.0.0,720.8548,798.79926,907.6626,1778.9182,, +unspecified,historical,World,Atmospheric Concentrations|CO2,ppm,input4MIPs,non-idealised,1,CMIP6,RCMIP Phase 3 v2.0.0,278.00745,284.29727,295.62943,369.2326,, +unspecified,historical,World,Atmospheric Concentrations|N2O,ppb,input4MIPs,non-idealised,1,CMIP6,RCMIP Phase 3 v2.0.0,266.57,271.57,278.32245,315.70068,, +MESSAGE-GLOBIOM,ssp245,World,Atmospheric Concentrations|CH4,ppb,input4MIPs,non-idealised,2,CMIP6,RCMIP Phase 3 v2.0.0,731.405995686849,808.2490285237631,925.5521138509116,1778.0101216634114,2020.2396647135413,1683.1598612467449 +MESSAGE-GLOBIOM,ssp245,World,Atmospheric Concentrations|CO2,ppm,input4MIPs,non-idealised,2,CMIP6,RCMIP Phase 3 v2.0.0,277.1470031738281,284.3169987996419,295.6749954223633,369.1249745686849,506.8749542236328,602.781982421875 +MESSAGE-GLOBIOM,ssp245,World,Atmospheric Concentrations|N2O,ppb,input4MIPs,non-idealised,2,CMIP6,RCMIP Phase 3 v2.0.0,273.86505126953125,273.0210469563802,279.4540481567383,315.75899505615234,356.1659749348958,377.2639795939128 +MESSAGE-GLOBIOM,ssp245,World,Atmospheric Concentrations|F-Gases|HFC|HFC134a,ppt,input4MIPs,non-idealised,2,CMIP6,RCMIP Phase 3 v2.0.0,0.0,0.0,0.0,15.3,118.7,182.4 +MESSAGE-GLOBIOM,ssp245,World,Atmospheric Concentrations|F-Gases|SF6,ppt,input4MIPs,non-idealised,2,CMIP6,RCMIP Phase 3 v2.0.0,0.0,0.0,0.0,4.52,11.9,18.3 +MESSAGE-GLOBIOM,ssp245,World,Atmospheric Concentrations|Montreal Gases|CFC|CFC12,ppt,input4MIPs,non-idealised,2,CMIP6,RCMIP Phase 3 v2.0.0,0.0,0.0,0.0,540.1,402.6,251.8 diff --git a/tests/unit/adapters/test_magicc7_concentrations.py b/tests/unit/adapters/test_magicc7_concentrations.py new file mode 100644 index 00000000..c61ffd1d --- /dev/null +++ b/tests/unit/adapters/test_magicc7_concentrations.py @@ -0,0 +1,186 @@ +"""Unit tests for the MAGICC7 conc-driven translator. + +These tests do NOT require pymagicc or the MAGICC binary; they +exercise the pure-Python overlay / filter / cfg-key logic only. The +end-to-end conc-file writer (which goes through pymagicc) is tested +under ``@pytest.mark.magicc`` in ``tests/integration/test_magicc7.py``. +""" +from __future__ import annotations + +from pathlib import Path + +import pandas as pd +import pytest +from scmdata import ScmRun + +from openscm_runner.adapters.magicc7._concentrations_translator import ( + FGAS_NAMES, + MHALO_NAMES, + RCMIP_TO_MAGICC_SPECIES, + build_concentrations_overlay, + cfg_keys_for_species, + classify_conc_species, + to_magicc_species, +) + +_TEST_DATA = Path(__file__).parent.parent.parent / "test-data" +RCMIP3_MINI_BUNDLE = _TEST_DATA / "rcmip3-mini" +# Variant bundle that additionally ships a few F-gas / Montreal-halocarbon +# ssp245 concentration trajectories (HFC134a, SF6, CFC12) to exercise the +# bundled-array path. Kept separate from the main mini bundle so the +# binary end-to-end GSAT comparison stays on the WMGHGs only. +RCMIP3_MINI_HALO_BUNDLE = _TEST_DATA / "rcmip3-mini-halo" + + +def _user_conc_run( + scenario: str, species_short: str, unit: str, values: dict[int, float], +) -> ScmRun: + """Build a one-row ScmRun with an ``Atmospheric Concentrations|`` trajectory.""" + frame = pd.DataFrame({ + "model": ["unspecified"], + "scenario": [scenario], + "region": ["World"], + "variable": [f"Atmospheric Concentrations|{species_short}"], + "unit": [unit], + **{year: [val] for year, val in values.items()}, + }) + return ScmRun(frame) + + +def test_cfg_keys_for_species_naming_convention(): + """MAGICC7 namelist: ``FILE__CONC`` + ``_SWITCHFROMCONC2EMIS_YEAR``.""" + assert cfg_keys_for_species("CO2") == ( + "file_co2_conc", "co2_switchfromconc2emis_year", + ) + assert cfg_keys_for_species("CH4") == ( + "file_ch4_conc", "ch4_switchfromconc2emis_year", + ) + assert cfg_keys_for_species("N2O") == ( + "file_n2o_conc", "n2o_switchfromconc2emis_year", + ) + + +def test_cfg_keys_for_species_rejects_non_per_gas_species(): + """Per-gas flags are CO2/CH4/N2O only; F-gases / halocarbons raise + (they route through the bundled-array path instead).""" + with pytest.raises(ValueError, match="per-gas species"): + cfg_keys_for_species("HFC134A") + with pytest.raises(ValueError, match="per-gas species"): + cfg_keys_for_species("CFC11") + + +def test_classify_conc_species(): + """Species route to per-gas / fgas / mhalo / unsupported correctly.""" + assert classify_conc_species("CO2") == "per_gas" + assert classify_conc_species("N2O") == "per_gas" + assert classify_conc_species("SF6") == "fgas" + assert classify_conc_species("HFC134A") == "fgas" + assert classify_conc_species("CFC12") == "mhalo" + assert classify_conc_species("HALON1211") == "mhalo" + assert classify_conc_species("NOTAGAS") is None + + +def test_to_magicc_species_normalises_case_and_renames(): + """RCMIP3 leaf names map onto MAGICC's upper-case FGAS/MHALO names.""" + # Pure case differences (uppercase fallback): + assert to_magicc_species("HFC134a") == "HFC134A" + assert to_magicc_species("CCl4") == "CCL4" + assert to_magicc_species("CH3Cl") == "CH3CL" + # Explicit renames win over the fallback: + assert to_magicc_species("H-1211") == "HALON1211" + assert to_magicc_species("HFC4310mee") == "HFC4310" + assert to_magicc_species("cC4F8") == "CC4F8" + # Already-canonical names round-trip: + assert to_magicc_species("CO2") == "CO2" + assert to_magicc_species("SF6") == "SF6" + + +def test_rcmip_to_magicc_species_covers_halons_and_hfc4310mee(): + """Spot-check the species rename table for the cases that bite.""" + assert RCMIP_TO_MAGICC_SPECIES["H-1211"] == "HALON1211" + assert RCMIP_TO_MAGICC_SPECIES["HFC4310mee"] == "HFC4310" + + +def test_fgas_mhalo_name_lists_are_sane(): + """Positional name arrays match the MAGICC namelist cardinality.""" + assert len(FGAS_NAMES) == 23 + assert len(MHALO_NAMES) == 18 + assert len(set(FGAS_NAMES)) == len(FGAS_NAMES) + assert len(set(MHALO_NAMES)) == len(MHALO_NAMES) + assert "HALON1202" == MHALO_NAMES[-1] + + +def test_overlay_uses_rcmip3_baseline_when_user_run_empty(): + """With no user overlay, ssp245 in rcmip3-mini supplies CO2/CH4/N2O.""" + empty_run = None + overlay = build_concentrations_overlay( + empty_run, RCMIP3_MINI_BUNDLE, ["ssp245"], + ) + assert set(overlay["magicc_species"]) == {"CO2", "CH4", "N2O"} + co2 = overlay[overlay["magicc_species"] == "CO2"].iloc[0] + assert co2["scenario"] == "ssp245" + assert co2["unit"] == "ppm" + assert co2["series"].loc[2100] == pytest.approx(602.78, rel=1e-3) + + +def test_overlay_includes_fgas_mhalo_from_bundle(): + """The halo bundle's F-gas / Montreal-halocarbon ssp245 rows survive, + are mapped onto MAGICC's upper-case names, and carry ppt units.""" + empty_run = None + overlay = build_concentrations_overlay( + empty_run, RCMIP3_MINI_HALO_BUNDLE, ["ssp245"], + ) + assert set(overlay["magicc_species"]) == { + "CO2", "CH4", "N2O", "HFC134A", "SF6", "CFC12", + } + sf6 = overlay[overlay["magicc_species"] == "SF6"].iloc[0] + assert sf6["unit"] == "ppt" + assert sf6["series"].loc[2100] == pytest.approx(18.3, rel=1e-3) + # RCMIP3 leaf 'HFC134a' resolved to MAGICC's 'HFC134A': + assert "HFC134a" not in set(overlay["magicc_species"]) + + +def test_overlay_user_row_merges_year_by_year_with_baseline(): + """User CO2 values win per-year; gaps are filled from the bundle.""" + user = _user_conc_run( + "ssp245", "CO2", "ppm", {2050: 999.0, 2100: 999.0}, + ) + overlay = build_concentrations_overlay( + user, RCMIP3_MINI_BUNDLE, ["ssp245"], + ) + co2 = overlay[overlay["magicc_species"] == "CO2"].iloc[0] + # User values win where supplied: + assert co2["series"].loc[2100] == pytest.approx(999.0) + assert co2["series"].loc[2050] == pytest.approx(999.0) + # Baseline fills the gaps the user did NOT supply: + assert co2["series"].loc[1850] == pytest.approx(284.32, rel=1e-3) + # CH4 / N2O still served from the RCMIP3 baseline: + ch4 = overlay[overlay["magicc_species"] == "CH4"].iloc[0] + assert ch4["series"].loc[2100] == pytest.approx(1683.16, rel=1e-3) + + +def test_overlay_is_per_scenario_not_batch_intersection(): + """The conc-driven decision is per-scenario, with no batch-consistency + requirement (unlike FaIR2). rcmip3-mini ships ssp245 conc data but no + ssp126; a batch of both should still drive ssp245's species and simply + leave ssp126 to the emissions path — not blank out ssp245.""" + empty_run = None + overlay = build_concentrations_overlay( + empty_run, RCMIP3_MINI_BUNDLE, ["ssp245", "ssp126"], + ) + assert not overlay.empty + # ssp245 species survive ... + assert set(overlay.loc[overlay["scenario"] == "ssp245", "magicc_species"]) == { + "CO2", "CH4", "N2O", + } + # ... and ssp126 (no baseline data) contributes nothing. + assert (overlay["scenario"] == "ssp126").sum() == 0 + + +def test_overlay_returns_empty_when_no_scenarios(): + """No scenario_names -> empty overlay (no work to do).""" + empty_run = None + overlay = build_concentrations_overlay( + empty_run, RCMIP3_MINI_BUNDLE, [], + ) + assert overlay.empty