pip install qg-python
The 1.5-layer quasi-geostrophic model is the classic testbed for mesoscale ocean variability, Rossby waves, the inverse energy cascade and ensemble data assimilation. QG-Python is a faithful port of the QG-C model of Sakov (2024) in which the legacy Fortran multi-grid Helmholtz solver is replaced by a direct sparse LU factorization (SuperLU through SciPy). The streamfunction inversion is then exact to machine precision and deterministic, so the only numerical error left is the one introduced by the time integrator, which is exactly what the package is built to measure.
The integrators are exposed under a single API: the classical fourth-order Runge-Kutta scheme, the Dormand-Prince scheme, and the Leapfrog scheme with the Robert-Asselin (RA) or the Robert-Asselin-Williams (RAW) filter. Analysis modules compute the linearized operator, its leading eigenspectrum and the maximum Lyapunov exponent. The default configuration is the 97 × 97 grid of QG-C (mrefin = 6, 9,025 interior degrees of freedom), unchanged, so results are comparable with the body of prior work that uses QG-C as a dynamical core.
One module, qg_model, on top of NumPy and SciPy. A run is a parameter object, a model and one call.
from qg_model import QGModel, QGParams
prm = QGParams(mrefin=6, # 97 x 97
scheme='leapfrog', dt=1.0, tend=110000, dtout=40,
ra_alpha=0.1, raw_filter=True,
outfname='qg_lf_raw.npz', save_q=True)
model = QGModel(prm) # factorizes the Helmholtz operator once
psi, q = model.run('restart', fname='qg_fourier.npz', record=-1)| Call | What it gives you |
|---|---|
QGParams(...) | all parameters as a dataclass: grid (mrefin, 6 gives 97 × 97), physics (rkb, rkh, rkh2, f, r, a, k), time integration (scheme, dt, tend, dtout), Leapfrog filter (ra_alpha, raw_filter) and output (outfname, save_q). Defaults are those of QG-C. |
QGParams.from_file("run.prm") | the same parameters read from a QG-C style parameter file |
QGModel(prm) | builds the model: assembles the Helmholtz operator once and factorizes it (SuperLU) |
model.run(ic, **kw) | integrates from t = 0 to tend, saves a record every dtout to outfname (.npz, or .nc when netCDF4 is installed) and returns (psi, q) |
ic = "zero" | "fourier" | "vortex" | "restart" | initial condition: rest, a seeded Fourier perturbation (amplitude, kmax, seed), a vortex, or a record of a previous run (fname, record) |
scheme = "rk4" | "dp5" | "leapfrog" | "order2" | "order1" | the time integrator; with leapfrog, ra_alpha sets the filter strength and raw_filter=True switches from RA to RAW |
model.set_ic(kind, **kw) | set the state without running, for custom loops |
model.psi, model.q, model.t | the current streamfunction, potential vorticity and time |
qg-python --help | the same model from the command line; this is what run_all.sh calls |
The stability and predictability analysis of the paper lives in the analysis/ scripts of the repository, which import the same model.
| Script | What it computes |
|---|---|
analysis/jacobian_spectrum.py | builds the sparse Jacobian of the right-hand side at an equilibrium state by finite differences and computes its leading eigenvalues with ARPACK; compares them with the von Neumann bounds |
analysis/stability_regions.py | draws the stability regions of RK4, dp5 and Leapfrog in the z-plane and overlays the scaled eigenvalues |
analysis/lyapunov.py | the Benettin algorithm: maximum Lyapunov exponent, its convergence trace and the Lyapunov time |
analysis/uncertainty.py | kinetic energy, enstrophy, energy spectra and KL divergences across schemes |
analysis/plot_alpha_sensitivity.py | the filter sensitivity sweep over the RA strength α |
qg_viewer.py | interactive viewer of the .npz records (needs matplotlib) |
pip install qg-python # the model: numpy and scipy only
pip install "qg-python[all]" # plus matplotlib (viewer) and netCDF4 (.nc output)The package installs the qg_model module and the qg-python command. The experiments of the paper, the analysis scripts and the reproducer live in the repository; the Zenodo archive is the authoritative citation target and also holds the simulation outputs (equilibrium restart fields at both resolutions, eigenvalue clouds, the Benettin trace and the time-resolved diagnostics of every run).
One command regenerates every experiment: the spin-up from rest, the RK4 reference run, the LF + RA and LF + RAW runs, the α sensitivity sweep over {0.05, 0.10, 0.20, 0.40}, the eigenspectra at both resolutions, the Lyapunov estimate and all the figures. About 4 to 5 hours on one laptop core, no GPU, bitwise reproducible given the same NumPy, SciPy and ARPACK versions. Each step can also be run on its own.
git clone https://github.com/enino84/qg-leapfrog-stability
cd qg-leapfrog-stability
bash run_all.shIf you use QG-Python or the dataset, please cite the Data Services Article and the Zenodo record.
Nino-Ruiz, E. D. (2026). QG-Python: A reproducible 1.5-layer quasi-geostrophic modelling service with direct sparse solvers, spectral stability diagnostics and filter sensitivity data. Geoscience Data Journal, 13, e70107. doi:10.1002/gdj3.70107
Nino-Ruiz, E. D. (2026). QG-Python: a reproducible 1.5-layer quasi-geostrophic modelling service with spectral stability and filter sensitivity diagnostics (Version 1.0) [Software and dataset]. Zenodo. doi:10.5281/zenodo.22168803
@article{ninoruiz2026qgpython,
title = {QG-Python: A reproducible 1.5-layer quasi-geostrophic modelling service
with direct sparse solvers, spectral stability diagnostics and filter
sensitivity data},
author = {Nino-Ruiz, Elias D.},
journal = {Geoscience Data Journal},
volume = {13},
pages = {e70107},
year = {2026},
doi = {10.1002/gdj3.70107}
}qg_model