diff --git a/examples/swm.ipynb b/examples/swm.ipynb
new file mode 100644
index 0000000000..6f851ea438
--- /dev/null
+++ b/examples/swm.ipynb
@@ -0,0 +1,453 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "GT4Py - GridTools Framework\n",
+    "\n",
+    "Copyright (c) 2014-2024, ETH Zurich\n",
+    "All rights reserved.\n",
+    "\n",
+    "Please, refer to the LICENSE file in the root directory.\n",
+    "SPDX-License-Identifier: BSD-3-Clause"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Shallow Water Model\n",
+    "\n",
+    "A doubly periodic shallow water model on an Arakawa C-grid, integrated with a leapfrog\n",
+    "scheme and a Robert-Asselin time filter. It is the NCAR\n",
+    "[SWM](https://github.com/NCAR/SWM) mini-app, a long-standing benchmark for comparing\n",
+    "programming models, written here as a single `gt4py.next` field operator.\n",
+    "\n",
+    "The prognostic equations are\n",
+    "\n",
+    "$$\n",
+    "\\frac{\\partial \\mathbf{V}}{\\partial t} + (\\zeta + f)\\,\\mathbf{k} \\times \\mathbf{V}\n",
+    "  + \\nabla\\left(P + \\tfrac{1}{2}\\mathbf{V}\\cdot\\mathbf{V}\\right) = 0,\n",
+    "\\qquad\n",
+    "\\frac{\\partial P}{\\partial t} + \\nabla\\cdot(P\\mathbf{V}) = 0 ,\n",
+    "$$\n",
+    "\n",
+    "discretised with the potential-enstrophy-conserving finite differences of Sadourny (1975,\n",
+    "his Eq. 4). The four intermediates below are his mass fluxes and diagnostics,\n",
+    "\n",
+    "$$\n",
+    "U = \\overline{P}^{\\,x} u, \\qquad V = \\overline{P}^{\\,y} v, \\qquad\n",
+    "H = P + \\tfrac{1}{2}\\left(\\overline{u^2}^{\\,x} + \\overline{v^2}^{\\,y}\\right), \\qquad\n",
+    "\\eta = \\frac{\\delta_x v - \\delta_y u}{\\overline{P}^{\\,xy}} ,\n",
+    "$$\n",
+    "\n",
+    "named `cu`, `cv`, `h` and `z` here.\n",
+    "\n",
+    "**Provenance.** The *scheme* is Sadourny's. The *configuration* --- grid size, spacing,\n",
+    "time step, filter coefficient and initial condition --- comes from the\n",
+    "[NCAR/SWM](https://github.com/NCAR/SWM) benchmark rather than from the paper; Sadourny\n",
+    "damps the leapfrog scheme by averaging odd and even time levels every $N$ steps, not with\n",
+    "the per-step Robert-Asselin filter used here."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "import gt4py.next as gtx\n",
+    "\n",
+    "import swm_numpy"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Grid and backend\n",
+    "\n",
+    "Fields carry **one halo cell on every side**, so the domain runs from `-1` to `M+1`. All\n",
+    "three prognostic fields share this one index space; the C-grid staggering is expressed by\n",
+    "the shifts inside the operators rather than by different array shapes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "I = gtx.Dimension(\"I\")\n",
+    "J = gtx.Dimension(\"J\")\n",
+    "\n",
+    "IJField = gtx.Field[gtx.Dims[I, J], gtx.float64]\n",
+    "\n",
+    "M = swm_numpy.M  # 16\n",
+    "N = swm_numpy.N  # 16\n",
+    "\n",
+    "DX = swm_numpy.DX  # 100 km\n",
+    "DY = swm_numpy.DY\n",
+    "DT = swm_numpy.DT  # 90 s\n",
+    "ALPHA = swm_numpy.ALPHA  # Robert-Asselin filter coefficient\n",
+    "\n",
+    "backend = None\n",
+    "# backend = gtx.gtfn_cpu\n",
+    "# backend = gtx.gtfn_gpu"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Differential operators\n",
+    "\n",
+    "Each operator comes in two variants that differ only in shift direction. The forward form\n",
+    "($+1$) maps an unstaggered quantity onto a staggered location, the backward form ($-1$)\n",
+    "maps back. Which one a term needs is fixed by where its operands live: `p` and `h` sit at\n",
+    "cell centres, `u` and `cu` on $x$-faces, `v` and `cv` on $y$-faces, and `z` on corners.\n",
+    "\n",
+    "Keeping the two variants apart by hand is exactly the bookkeeping that\n",
+    "[staggered dimensions](../docs/development/ADRs/next/0026-Staggered_Dimensions.md) are\n",
+    "designed to remove; a follow-up rewrites these eight operators as four."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@gtx.field_operator\n",
+    "def avg_x(f: IJField) -> IJField:\n",
+    "    \"\"\"Average onto the next x-face.\"\"\"\n",
+    "    return 0.5 * (f(I + 1) + f)\n",
+    "\n",
+    "\n",
+    "@gtx.field_operator\n",
+    "def avg_y(f: IJField) -> IJField:\n",
+    "    \"\"\"Average onto the next y-face.\"\"\"\n",
+    "    return 0.5 * (f(J + 1) + f)\n",
+    "\n",
+    "\n",
+    "@gtx.field_operator\n",
+    "def avg_x_staggered(f: IJField) -> IJField:\n",
+    "    \"\"\"Average an x-staggered field back onto cell centres.\"\"\"\n",
+    "    return 0.5 * (f(I - 1) + f)\n",
+    "\n",
+    "\n",
+    "@gtx.field_operator\n",
+    "def avg_y_staggered(f: IJField) -> IJField:\n",
+    "    \"\"\"Average a y-staggered field back onto cell centres.\"\"\"\n",
+    "    return 0.5 * (f(J - 1) + f)\n",
+    "\n",
+    "\n",
+    "@gtx.field_operator\n",
+    "def delta_x(dx: gtx.float64, f: IJField) -> IJField:\n",
+    "    \"\"\"Forward difference in x.\"\"\"\n",
+    "    return (1.0 / dx) * (f(I + 1) - f)\n",
+    "\n",
+    "\n",
+    "@gtx.field_operator\n",
+    "def delta_y(dy: gtx.float64, f: IJField) -> IJField:\n",
+    "    \"\"\"Forward difference in y.\"\"\"\n",
+    "    return (1.0 / dy) * (f(J + 1) - f)\n",
+    "\n",
+    "\n",
+    "@gtx.field_operator\n",
+    "def delta_x_staggered(dx: gtx.float64, f: IJField) -> IJField:\n",
+    "    \"\"\"Backward difference in x, for an x-staggered field.\"\"\"\n",
+    "    return (1.0 / dx) * (f - f(I - 1))\n",
+    "\n",
+    "\n",
+    "@gtx.field_operator\n",
+    "def delta_y_staggered(dy: gtx.float64, f: IJField) -> IJField:\n",
+    "    \"\"\"Backward difference in y, for a y-staggered field.\"\"\"\n",
+    "    return (1.0 / dy) * (f - f(J - 1))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The timestep\n",
+    "\n",
+    "One fused field operator produces all six updated fields. The first three are the leapfrog\n",
+    "update of the prognostic variables; the last three are the time-filtered previous level,\n",
+    "\n",
+    "$$ u^{\\text{old}} \\leftarrow u + \\alpha\\,(u^{\\text{new}} - 2u + u^{\\text{old}}) , $$\n",
+    "\n",
+    "which suppresses the computational mode of the leapfrog scheme."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@gtx.field_operator\n",
+    "def timestep(\n",
+    "    u: IJField,\n",
+    "    v: IJField,\n",
+    "    p: IJField,\n",
+    "    uold: IJField,\n",
+    "    vold: IJField,\n",
+    "    pold: IJField,\n",
+    "    dx: gtx.float64,\n",
+    "    dy: gtx.float64,\n",
+    "    tdt: gtx.float64,\n",
+    "    alpha: gtx.float64,\n",
+    ") -> tuple[IJField, IJField, IJField, IJField, IJField, IJField]:\n",
+    "    \"\"\"Advance one leapfrog step and return (unew, vnew, pnew, uold, vold, pold).\"\"\"\n",
+    "    cu = avg_x(p) * u\n",
+    "    cv = avg_y(p) * v\n",
+    "    z = (delta_x(dx, v) - delta_y(dy, u)) / avg_x(avg_y(p))\n",
+    "    h = p + 0.5 * (avg_x_staggered(u * u) + avg_y_staggered(v * v))\n",
+    "\n",
+    "    unew = uold + avg_y_staggered(z) * avg_y_staggered(avg_x(cv)) * tdt - delta_x(dx, h) * tdt\n",
+    "    vnew = vold - avg_x_staggered(z) * avg_x_staggered(avg_y(cu)) * tdt - delta_y(dy, h) * tdt\n",
+    "    pnew = pold - delta_x_staggered(dx, cu) * tdt - delta_y_staggered(dy, cv) * tdt\n",
+    "\n",
+    "    uold_new = u + alpha * (unew - 2.0 * u + uold)\n",
+    "    vold_new = v + alpha * (vnew - 2.0 * v + vold)\n",
+    "    pold_new = p + alpha * (pnew - 2.0 * p + pold)\n",
+    "\n",
+    "    return unew, vnew, pnew, uold_new, vold_new, pold_new"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Periodicity\n",
+    "\n",
+    "The domain wraps in both directions. Here the halo is refreshed **outside** GT4Py, in\n",
+    "NumPy, after each step: the operator computes the interior only, and the halo is then\n",
+    "overwritten with the opposite edge.\n",
+    "\n",
+    "Rows are copied before columns, which is what fills the corners correctly.\n",
+    "\n",
+    "Expressing the wrap *inside* the field operator instead --- with `concat_where` over a\n",
+    "domain boundary --- is the subject of a follow-up; it needs shift offsets that are only\n",
+    "known at compile time rather than when the operator is traced."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def apply_periodic(field: gtx.Field, m: int = M, n: int = N) -> None:\n",
+    "    \"\"\"Refresh the one-cell halo of `field` in place from the opposite edge.\"\"\"\n",
+    "    a = field.ndarray\n",
+    "    a[0, :] = a[m, :]\n",
+    "    a[m + 1, :] = a[1, :]\n",
+    "    a[:, 0] = a[:, n]\n",
+    "    a[:, n + 1] = a[:, 1]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Time loop\n",
+    "\n",
+    "`uold`, `vold` and `pold` are updated in place and never swapped; only the prognostic\n",
+    "fields ping-pong with their `new` counterparts. The first cycle is a forward Euler step\n",
+    "with the filter switched off, which is how the leapfrog scheme is started."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def run(itmax: int, m: int = M, n: int = N) -> tuple[gtx.Field, gtx.Field, gtx.Field]:\n",
+    "    \"\"\"Integrate `itmax` steps and return the final `u`, `v`, `p` fields.\"\"\"\n",
+    "    domain = gtx.domain({I: (-1, m + 1), J: (-1, n + 1)})\n",
+    "    interior = gtx.domain({I: (0, m), J: (0, n)})\n",
+    "\n",
+    "    u0, v0, p0 = swm_numpy.initial_conditions(m, n, DX, DY)\n",
+    "\n",
+    "    def as_halo_field(a: np.ndarray) -> gtx.Field:\n",
+    "        return gtx.as_field(domain, np.pad(a, 1, mode=\"wrap\"), allocator=backend)\n",
+    "\n",
+    "    u, v, p = as_halo_field(u0), as_halo_field(v0), as_halo_field(p0)\n",
+    "    uold, vold, pold = as_halo_field(u0), as_halo_field(v0), as_halo_field(p0)\n",
+    "    unew, vnew, pnew = (gtx.zeros(domain, dtype=gtx.float64, allocator=backend) for _ in range(3))\n",
+    "\n",
+    "    for cycle in range(itmax):\n",
+    "        tdt = DT if cycle == 0 else 2.0 * DT\n",
+    "        alpha = 0.0 if cycle == 0 else ALPHA\n",
+    "\n",
+    "        timestep(\n",
+    "            u,\n",
+    "            v,\n",
+    "            p,\n",
+    "            uold,\n",
+    "            vold,\n",
+    "            pold,\n",
+    "            DX,\n",
+    "            DY,\n",
+    "            tdt,\n",
+    "            alpha,\n",
+    "            out=(unew, vnew, pnew, uold, vold, pold),\n",
+    "            domain=interior,\n",
+    "            offset_provider={},\n",
+    "        )\n",
+    "\n",
+    "        for field in (unew, vnew, pnew):\n",
+    "            apply_periodic(field, m, n)\n",
+    "\n",
+    "        u, unew = unew, u\n",
+    "        v, vnew = vnew, v\n",
+    "        p, pnew = pnew, p\n",
+    "\n",
+    "    return u, v, p\n",
+    "\n",
+    "\n",
+    "def interior_of(field: gtx.Field, m: int = M, n: int = N) -> np.ndarray:\n",
+    "    return field.asnumpy()[1 : m + 1, 1 : n + 1]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Validation against NumPy\n",
+    "\n",
+    "The companion module `swm_numpy.py` implements the same scheme with `numpy.roll`, which\n",
+    "gives periodicity for free and needs no halo at all. Agreement between the two is a check\n",
+    "that the halo bookkeeping here is right."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "STEPS = 100\n",
+    "\n",
+    "# The two implementations sum the same terms in a different order, so they agree to\n",
+    "# rounding rather than bit-for-bit. Tolerances are absolute because u and v cross zero.\n",
+    "ATOL_UV = 1e-11  # peak |u|, |v| ~ 4\n",
+    "ATOL_P = 1e-8  # peak |p| ~ 5e4\n",
+    "\n",
+    "\n",
+    "def test_matches_numpy():\n",
+    "    u, v, p = run(STEPS)\n",
+    "    u_ref, v_ref, p_ref = swm_numpy.run(STEPS)\n",
+    "\n",
+    "    np.testing.assert_allclose(interior_of(u), u_ref, atol=ATOL_UV)\n",
+    "    np.testing.assert_allclose(interior_of(v), v_ref, atol=ATOL_UV)\n",
+    "    np.testing.assert_allclose(interior_of(p), p_ref, atol=ATOL_P)\n",
+    "\n",
+    "\n",
+    "test_matches_numpy()\n",
+    "print(\"GT4Py and NumPy agree after\", STEPS, \"steps\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Validation against the NCAR reference\n",
+    "\n",
+    "`swm_reference.npz` holds `u`, `v`, `p` after the full 4000-step benchmark, taken from\n",
+    "[NCAR/SWM](https://github.com/NCAR/SWM) `ref/16x16` (Apache-2.0). Reproducing it end to\n",
+    "end is the real correctness check, but it takes far longer than a notebook should in CI,\n",
+    "so it is opt-in."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "RUN_FULL_BENCHMARK = False\n",
+    "\n",
+    "\n",
+    "def test_matches_reference():\n",
+    "    u, v, p = run(swm_numpy.ITMAX)\n",
+    "    u_out, v_out, p_out = swm_numpy.to_reference_layout(\n",
+    "        interior_of(u), interior_of(v), interior_of(p)\n",
+    "    )\n",
+    "\n",
+    "    reference = np.load(\"swm_reference.npz\")\n",
+    "    np.testing.assert_allclose(u_out, reference[\"u\"], atol=ATOL_UV)\n",
+    "    np.testing.assert_allclose(v_out, reference[\"v\"], atol=ATOL_UV)\n",
+    "    np.testing.assert_allclose(p_out, reference[\"p\"], atol=ATOL_P)\n",
+    "\n",
+    "\n",
+    "if RUN_FULL_BENCHMARK:\n",
+    "    test_matches_reference()\n",
+    "    print(\"GT4Py reproduces the NCAR reference after\", swm_numpy.ITMAX, \"steps\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The flow\n",
+    "\n",
+    "The initial stream function sets up a pair of counter-rotating vortices; after a few\n",
+    "hundred steps the pressure field has developed the finer structure the scheme is meant to\n",
+    "handle without spurious accumulation of energy at the grid scale."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "u, v, p = run(STEPS)\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 3, figsize=(13, 3.6))\n",
+    "for ax, field, name in zip(axes, (u, v, p), (\"u\", \"v\", \"p\")):\n",
+    "    image = ax.contourf(interior_of(field).T, levels=20)\n",
+    "    ax.set_title(f\"{name} after {STEPS} steps\")\n",
+    "    ax.set_aspect(\"equal\")\n",
+    "    fig.colorbar(image, ax=ax)\n",
+    "fig.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## References\n",
+    "\n",
+    "Sadourny, R. (1975). The dynamics of finite-difference models of the shallow-water\n",
+    "equations. *Journal of the Atmospheric Sciences*, 32(4), 680-689.\n",
+    "[doi:10.1175/1520-0469(1975)032<0680:TDOFDM>2.0.CO;2](https://doi.org/10.1175/1520-0469(1975)032%3C0680:TDOFDM%3E2.0.CO;2)\n",
+    "\n",
+    "The benchmark itself --- its configuration, initial condition and reference data --- is\n",
+    "[NCAR/SWM](https://github.com/NCAR/SWM), which implements the same model across a range of\n",
+    "languages and programming models."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/examples/swm_numpy.py b/examples/swm_numpy.py
new file mode 100644
index 0000000000..7db27aa980
--- /dev/null
+++ b/examples/swm_numpy.py
@@ -0,0 +1,164 @@
+# GT4Py - GridTools Framework
+#
+# Copyright (c) 2014-2024, ETH Zurich
+# All rights reserved.
+#
+# Please, refer to the LICENSE file in the root directory.
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""
+NumPy reference for the shallow water model example.
+
+Companion to `swm.ipynb`. The model is doubly periodic, so periodicity is expressed here
+with `numpy.roll` on interior-only arrays of shape `(M, N)` -- there is no halo. The
+notebook's GT4Py version instead carries one halo cell per side and updates it explicitly,
+which is what makes the two implementations worth comparing.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+
+# Model constants of the NCAR/SWM benchmark configuration.
+M = 16
+N = 16
+DX = 100000.0
+DY = 100000.0
+DT = 90.0
+ALPHA = 0.001
+A = 1000000.0
+P_MEAN = 50000.0
+ITMAX = 4000
+
+
+def avg_x(f: np.ndarray) -> np.ndarray:
+    return 0.5 * (np.roll(f, -1, axis=0) + f)
+
+
+def avg_y(f: np.ndarray) -> np.ndarray:
+    return 0.5 * (np.roll(f, -1, axis=1) + f)
+
+
+def avg_x_staggered(f: np.ndarray) -> np.ndarray:
+    return 0.5 * (np.roll(f, 1, axis=0) + f)
+
+
+def avg_y_staggered(f: np.ndarray) -> np.ndarray:
+    return 0.5 * (np.roll(f, 1, axis=1) + f)
+
+
+def delta_x(dx: float, f: np.ndarray) -> np.ndarray:
+    return (np.roll(f, -1, axis=0) - f) / dx
+
+
+def delta_y(dy: float, f: np.ndarray) -> np.ndarray:
+    return (np.roll(f, -1, axis=1) - f) / dy
+
+
+def delta_x_staggered(dx: float, f: np.ndarray) -> np.ndarray:
+    return (f - np.roll(f, 1, axis=0)) / dx
+
+
+def delta_y_staggered(dy: float, f: np.ndarray) -> np.ndarray:
+    return (f - np.roll(f, 1, axis=1)) / dy
+
+
+def initial_conditions(
+    m: int = M, n: int = N, dx: float = DX, dy: float = DY, a: float = A, p_mean: float = P_MEAN
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Return the initial `u`, `v`, `p` on the interior grid, each of shape `(m, n)`.
+
+    The velocity field is derived from a doubly periodic stream function
+    `psi = a * sin(2*pi*(i+0.5)/m) * sin(2*pi*(j+0.5)/n)`.
+    """
+    d_i = 2.0 * np.pi / m
+    d_j = 2.0 * np.pi / n
+    el = n * dx
+    pcf = (np.pi * np.pi * a * a) / (el * el)
+
+    psi = (
+        a
+        * np.sin((np.arange(0, m + 1).reshape(-1, 1) + 0.5) * d_i)
+        * np.sin((np.arange(0, n + 1) + 0.5) * d_j)
+    )
+    p = (
+        pcf
+        * (np.cos(2.0 * np.arange(0, m).reshape(-1, 1) * d_i) + np.cos(2.0 * np.arange(0, n) * d_j))
+        + p_mean
+    )
+    u = -(psi[1:, 1:] - psi[1:, :-1]) / dy
+    v = (psi[1:, 1:] - psi[:-1, 1:]) / dx
+
+    return u, v, p
+
+
+def timestep(
+    u: np.ndarray,
+    v: np.ndarray,
+    p: np.ndarray,
+    uold: np.ndarray,
+    vold: np.ndarray,
+    pold: np.ndarray,
+    dx: float,
+    dy: float,
+    tdt: float,
+    alpha: float,
+) -> tuple[np.ndarray, ...]:
+    """One leapfrog step with a Robert-Asselin time filter. Returns the six updated fields."""
+    cu = avg_x(p) * u
+    cv = avg_y(p) * v
+    z = (delta_x(dx, v) - delta_y(dy, u)) / avg_x(avg_y(p))
+    h = p + 0.5 * (avg_x_staggered(u * u) + avg_y_staggered(v * v))
+
+    unew = uold + avg_y_staggered(z) * avg_y_staggered(avg_x(cv)) * tdt - delta_x(dx, h) * tdt
+    vnew = vold - avg_x_staggered(z) * avg_x_staggered(avg_y(cu)) * tdt - delta_y(dy, h) * tdt
+    pnew = pold - delta_x_staggered(dx, cu) * tdt - delta_y_staggered(dy, cv) * tdt
+
+    uold_new = u + alpha * (unew - 2.0 * u + uold)
+    vold_new = v + alpha * (vnew - 2.0 * v + vold)
+    pold_new = p + alpha * (pnew - 2.0 * p + pold)
+
+    return unew, vnew, pnew, uold_new, vold_new, pold_new
+
+
+def run(
+    itmax: int,
+    m: int = M,
+    n: int = N,
+    dx: float = DX,
+    dy: float = DY,
+    dt: float = DT,
+    alpha: float = ALPHA,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Integrate `itmax` steps and return the final interior `u`, `v`, `p`."""
+    u, v, p = initial_conditions(m, n, dx, dy)
+    uold, vold, pold = u.copy(), v.copy(), p.copy()
+
+    for cycle in range(itmax):
+        # The first step is forward Euler and unfiltered; afterwards leapfrog with tdt = 2*dt.
+        tdt = dt if cycle == 0 else 2.0 * dt
+        step_alpha = 0.0 if cycle == 0 else alpha
+        unew, vnew, pnew, uold, vold, pold = timestep(
+            u, v, p, uold, vold, pold, dx, dy, tdt, step_alpha
+        )
+        u, v, p = unew, vnew, pnew
+
+    return u, v, p
+
+
+def to_reference_layout(
+    u: np.ndarray, v: np.ndarray, p: np.ndarray
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Convert interior `(M, N)` fields to the `(M+1, N+1)` layout of the NCAR/SWM dumps.
+
+    The extra row and column are the periodic images the original Fortran carried
+    explicitly; which side they sit on differs per field because of the C-grid staggering.
+    """
+    return (
+        np.pad(u, ((1, 0), (0, 1)), mode="wrap"),
+        np.pad(v, ((0, 1), (1, 0)), mode="wrap"),
+        np.pad(p, ((0, 1), (0, 1)), mode="wrap"),
+    )
diff --git a/examples/swm_reference.npz b/examples/swm_reference.npz
new file mode 100644
index 0000000000..6a9ca4af2f
Binary files /dev/null and b/examples/swm_reference.npz differ