Climate Modeling Guide

import xarray as xr
import numpy as np

# Open a CMIP6 temperature dataset
ds = xr.open_dataset("tas_Amon_CESM2_ssp585_r1i1p1f1_gn_201501-210012.nc")

print(ds)
# Dimensions:  (time: 1032, lat: 192, lon: 288)
# Variables:   tas (surface air temperature, K)
# Attributes:  CF-1.6 compliant, CMIP6 metadata

# Basic inspection
print(f"Variable: {ds.tas.long_name}")
print(f"Units: {ds.tas.units}")
print(f"Time range: {ds.time.values[0]} to {ds.time.values[-1]}")
print(f"Spatial resolution: {np.diff(ds.lat.values[:2])[0]:.2f} deg")

# Using intake-esm to search the CMIP6 catalog
import intake

# Open the Pangeo CMIP6 catalog (cloud-hosted on Google Cloud)
url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
col = intake.open_esm_datastore(url)

# Search for monthly surface temperature under SSP5-8.5
query = col.search(
    experiment_id="ssp585",
    variable_id="tas",
    table_id="Amon",
    source_id=["CESM2", "GFDL-ESM4", "UKESM1-0-LL", "MPI-ESM1-2-HR"],
    member_id="r1i1p1f1",
)
print(f"Found {len(query)} datasets from {query.nunique()['source_id']} models")

# Load as xarray datasets (lazy, Zarr-backed)
dsets = query.to_dataset_dict(zarr_kwargs={"consolidated": True})

def compute_global_mean_anomaly(ds, baseline_start="1850-01-01",
                                  baseline_end="1900-12-31"):
    """
    Compute area-weighted global mean temperature anomaly
    relative to a baseline period.
    """
    # Area weighting by cosine of latitude
    weights = np.cos(np.deg2rad(ds.lat))
    weights.name = "weights"

    # Weighted global mean time series
    global_mean = ds.tas.weighted(weights).mean(dim=["lat", "lon"])

    # Compute baseline climatology
    baseline = global_mean.sel(time=slice(baseline_start, baseline_end))
    climatology = baseline.groupby("time.month").mean("time")

    # Compute anomalies
    anomaly = global_mean.groupby("time.month") - climatology

    # Annual mean anomaly
    annual_anomaly = anomaly.resample(time="YE").mean()
    return annual_anomaly


def multi_model_ensemble(datasets: dict, baseline_period: tuple):
    """
    Compute multi-model ensemble mean and spread for temperature projections.
    datasets: dict of {model_name: xarray.Dataset}
    Returns ensemble mean and 5th/95th percentile bounds.
    """
    anomalies = []
    for name, ds in datasets.items():
        anom = compute_global_mean_anomaly(ds, *baseline_period)
        anom = anom.assign_coords(model=name)
        anomalies.append(anom)

    ensemble = xr.concat(anomalies, dim="model")
    return {
        "mean": ensemble.mean(dim="model"),
        "p05": ensemble.quantile(0.05, dim="model"),
        "p95": ensemble.quantile(0.95, dim="model"),
    }

def compute_nino34(sst_dataset, baseline="1991-01-01/2020-12-31"):
    """Compute Nino 3.4 index from SST data."""
    # Select Nino 3.4 region
    nino34_region = sst_dataset.tos.sel(
        lat=slice(-5, 5), lon=slice(190, 240)
    )
    # Area-weighted mean
    weights = np.cos(np.deg2rad(nino34_region.lat))
    nino34_ts = nino34_region.weighted(weights).mean(dim=["lat", "lon"])

    # Remove monthly climatology
    clim = nino34_ts.sel(time=slice(*baseline.split("/"))).groupby("time.month").mean()
    nino34_index = nino34_ts.groupby("time.month") - clim

    # 5-month running mean for standard definition
    nino34_smoothed = nino34_index.rolling(time=5, center=True).mean()
    return nino34_smoothed

def quantile_mapping(obs: np.ndarray, model_hist: np.ndarray,
                     model_future: np.ndarray, n_quantiles: int = 100):
    """
    Quantile mapping bias correction.
    Maps model quantiles to observed quantiles for bias correction.
    """
    quantiles = np.linspace(0, 1, n_quantiles + 1)
    obs_q = np.quantile(obs, quantiles)
    hist_q = np.quantile(model_hist, quantiles)

    # For each future value, find its quantile in historical distribution
    # then map to corresponding observed quantile
    corrected = np.interp(model_future, hist_q, obs_q)
    return corrected

def energy_balance_model(S0=1361, albedo=0.30, emissivity=0.612):
    """
    Zero-dimensional energy balance model.
    S0: solar constant (W/m2)
    albedo: planetary albedo
    emissivity: effective atmospheric emissivity
    Returns equilibrium surface temperature (K).
    """
    sigma = 5.67e-8  # Stefan-Boltzmann constant
    # Absorbed solar radiation
    absorbed = S0 * (1 - albedo) / 4
    # Surface temperature with greenhouse effect
    T_surface = (absorbed / (emissivity * sigma)) ** 0.25
    return T_surface

T_eq = energy_balance_model()
print(f"Equilibrium surface temperature: {T_eq:.1f} K ({T_eq - 273.15:.1f} C)")