API

Base Model Class

This is the base class from which all other models inherit. All of these initialization arguments are available to all of the other model types. This class is not called directly.

class pyqg.Model(nz=1, nx=64, ny=None, L=1000000.0, W=None, dt=7200.0, twrite=1000.0, tmax=1576800000.0, tavestart=315360000.0, taveint=86400.0, useAB2=False, rek=5.787e-07, filterfac=23.6, f=None, g=9.81, q_parameterization=None, uv_parameterization=None, parameterization=None, diagnostics_list='all', ntd=1, log_level=1, logfile=None)

A generic pseudo-spectral inversion model.

Attributes

nx, nyint

Number of real space grid points in the x, y directions (cython)

nk, nlint

Number of spectral space grid points in the k, l directions (cython)

nzint

Number of vertical levels (cython)

kk, llreal array

Zonal and meridional wavenumbers (nk) (cython)

areal array

inversion matrix (nk, nk, nl, nk) (cython)

qreal array

Potential vorticity in real space (nz, ny, nx) (cython)

qhcomplex array

Potential vorticity in spectral space (nk, nl, nk) (cython)

phcomplex array

Streamfunction in spectral space (nk, nl, nk) (cython)

u, vreal array

Zonal and meridional velocity anomalies in real space (nz, ny, nx) (cython)

Ubgreal array

Background zonal velocity (nk) (cython)

Qyreal array

Background potential vorticity gradient (nk) (cython)

ufull, vfullreal arrays

Zonal and meridional full velocities in real space (nz, ny, nx) (cython)

uh, vhcomplex arrays

Velocity anomaly components in spectral space (nk, nl, nk) (cython)

rekfloat

Linear drag in lower layer (cython)

tfloat

Model time (cython)

tcint

Model timestep (cython)

dtfloat

Numerical timestep (cython)

L, Wfloat

Domain length in x and y directions

filterfacfloat

Amplitdue of the spectral spherical filter

twriteint

Interval for cfl writeout (units: number of timesteps)

tmaxfloat

Total time of integration (units: model time)

tavestartfloat

Start time for averaging (units: model time)

tsnapstartfloat

Start time for snapshot writeout (units: model time)

taveintfloat

Time interval for accumulation of diagnostic averages. (units: model time)

tsnapintfloat

Time interval for snapshots (units: model time)

ntdint

Number of threads to use. Should not exceed the number of cores on your machine.

pmodesreal array

Vertical pressure modes (unitless)

radiireal array

Deformation radii (units: model length)

q_parameterizationfunction or pyqg.Parameterization

Optional Parameterization object or function which takes the model as input and returns a numpy array of shape (nz, ny, nx) to be added to \(\partial_t q\) before stepping forward. This can be used to implement subgrid forcing parameterizations.

uv_parameterizationfunction or pyqg.Parameterization

Optional Parameterization object or function which takes the model as input and returns a tuple of two numpy arrays, each of shape (nz, ny, nx), to be added to the zonal and meridional velocity derivatives (respectively) at each timestep (by adding their curl to \(\partial_t q\)). This can also be used to implemented subgrid forcing parameterizations, but expressed in terms of velocity rather than potential vorticity.

Note

All of the test cases use nx==ny. Expect bugs if you choose these parameters to be different.

Note

All time intervals will be rounded to nearest dt interval.

Parameters

nxint

Number of grid points in the x direction.

nyint

Number of grid points in the y direction (default: nx).

Lnumber

Domain length in x direction. Units: meters.

W

Domain width in y direction. Units: meters (default: L).

reknumber

linear drag in lower layer. Units: seconds ^-1.

filterfacnumber

amplitdue of the spectral spherical filter (originally 18.4, later changed to 23.6).

dtnumber

Numerical timstep. Units: seconds.

twriteint

Interval for cfl writeout. Units: number of timesteps.

tmaxnumber

Total time of integration. Units: seconds.

tavestartnumber

Start time for averaging. Units: seconds.

tsnapstartnumber

Start time for snapshot writeout. Units: seconds.

taveintnumber

Time interval for accumulation of diagnostic averages. Units: seconds. (For performance purposes, averaging does not have to occur every timestep)

tsnapintnumber

Time interval for snapshots. Units: seconds.

ntdint

Number of threads to use. Should not exceed the number of cores on your machine.

q_parameterizationfunction or pyqg.Parameterization

Optional Parameterization object or function which takes the model as input and returns a numpy array of shape (nz, ny, nx) to be added to \(\partial_t q\) before stepping forward. This can be used to implement subgrid forcing parameterizations.

uv_parameterizationfunction or pyqg.Parameterization

Optional Parameterization object or function which takes the model as input and returns a tuple of two numpy arrays, each of shape (nz, ny, nx), to be added to the zonal and meridional velocity derivatives (respectively) at each timestep (by adding their curl to \(\partial_t q\)). This can also be used to implemented subgrid forcing parameterizations, but expressed in terms of velocity rather than potential vorticity.

parameterizationpyqg.Parameterization

An explicit Parameterization object representing either a velocity or potential vorticity parameterization, whose type will be inferred.

describe_diagnostics()

Print a human-readable summary of the available diagnostics.

modal_projection(p, forward=True)

Performs a field p into modal amplitudes pn using the basis [pmodes]. The inverse transform calculates p from pn

property parameterization

Return the model’s parameterization if present (either in terms of PV or velocity, warning if there are both).

Returns

parameterizationpyqg.Parameterization or function

run()

Run the model forward without stopping until the end.

run_with_snapshots(tsnapstart=0.0, tsnapint=432000.0)

Run the model forward, yielding to user code at specified intervals.

Parameters

tsnapstartint

The timestep at which to begin yielding.

tstapintint

The interval at which to yield.

spec_var(ph)

compute variance of p from Fourier coefficients ph

stability_analysis(bottom_friction=False)

Performs the baroclinic linear instability analysis given

given the base state velocity :math: (U, V) and the stretching matrix :math: S:

\[A \Phi = \omega B \Phi,\]

where

\[A = B (U k + V l) + I (k Q_y - l Q_x) + 1j \delta_{N N} r_{ek} I \kappa^2\]

where \(\delta_{N N} = [0,0,\dots,0,1] ,\)

and

\[B = S - I \kappa^2 .\]

The eigenstructure is

\[\Phi\]

and the eigenvalue is

\[`\omega`\]

The growth rate is Im\(\{\omega\}\).

Parameters

bottom_friction: optional inclusion linear bottom drag

in the linear stability calculation (default is False, as if :math: r_{ek} = 0)

Returns

omega: complex array

The eigenvalues with largest complex part (units: inverse model time)

phi: complex array

The eigenvectors associated associated with omega (unitless)

to_dataset()

Convert outputs from model to an xarray dataset

Returns

dsxarray.Dataset

vertical_modes()

Calculate standard vertical modes. Simply the eigenvectors of the stretching matrix S

Specific Model Types

These are the actual models which are run.

class pyqg.QGModel(beta=1.5e-11, rd=15000.0, delta=0.25, H1=500, U1=0.025, U2=0.0, **kwargs)

Two layer quasigeostrophic model.

This model is meant to representflows driven by baroclinic instabilty of a base-state shear \(U_1-U_2\). The upper and lower layer potential vorticity anomalies \(q_1\) and \(q_2\) are

\[\begin{split}q_1 &= \nabla^2\psi_1 + F_1(\psi_2 - \psi_1) \\ q_2 &= \nabla^2\psi_2 + F_2(\psi_1 - \psi_2)\end{split}\]

with

\[\begin{split}F_1 &\equiv \frac{k_d^2}{1 + \delta^2} \\ F_2 &\equiv \delta F_1 \ .\end{split}\]

The layer depth ratio is given by \(\delta = H_1 / H_2\). The total depth is \(H = H_1 + H_2\).

The background potential vorticity gradients are

\[\begin{split}\beta_1 &= \beta + F_1(U_1 - U_2) \\ \beta_2 &= \beta - F_2( U_1 - U_2) \ .\end{split}\]

The evolution equations for \(q_1\) and \(q_2\) are

\[\begin{split}\partial_t\,{q_1} + J(\psi_1\,, q_1) + \beta_1\, {\psi_1}_x &= \text{ssd} \\ \partial_t\,{q_2} + J(\psi_2\,, q_2)+ \beta_2\, {\psi_2}_x &= -r_{ek}\nabla^2 \psi_2 + \text{ssd}\,.\end{split}\]

where ssd represents small-scale dissipation and \(r_{ek}\) is the Ekman friction parameter.

Parameters

betanumber

Gradient of coriolis parameter. Units: meters ^-1 seconds ^-1

reknumber

Linear drag in lower layer. Units: seconds ^-1

rdnumber

Deformation radius. Units: meters.

deltanumber

Layer thickness ratio (H1/H2)

U1number

Upper layer flow. Units: meters seconds ^-1

U2number

Lower layer flow. Units: meters seconds ^-1

set_U1U2(U1, U2)

Set background zonal flow.

Parameters

U1number

Upper layer flow. Units: meters seconds ^-1

U2number

Lower layer flow. Units: meters seconds ^-1

set_q1q2(q1, q2, check=False)

Set upper and lower layer PV anomalies.

Parameters

q1array-like

Upper layer PV anomaly in spatial coordinates.

q1array-like

Lower layer PV anomaly in spatial coordinates.

class pyqg.LayeredModel(beta=1.5e-11, nz=3, rd=15000.0, f=0.0001236812857687059, H=None, U=None, V=None, rho=None, delta=None, **kwargs)

Layered quasigeostrophic model.

This model is meant to represent flows driven by baroclinic instabilty of a base-state shear. The potential vorticity anomalies qi are related to the streamfunction psii through

\[ \begin{align}\begin{aligned}{q_i} = \nabla^2\psi_i + \frac{f_0^2}{H_i} \left(\frac{\psi_{i-1}- \psi_i}{g'_{i-1}}- \frac{\psi_{i}-\psi_{i+1}}{g'_{i}}\right)\,, \qquad i = 2,\textsf{N}-1\,,\\{q_1} = \nabla^2\psi_1 + \frac{f_0^2}{H_1} \left(\frac{\psi_{2}- \psi_1}{g'_{1}}\right)\,, \qquad i =1\,,\\{q_\textsf{N}} = \nabla^2\psi_\textsf{N} + \frac{f_0^2}{H_\textsf{N}} \left(\frac{\psi_{\textsf{N}-1}- \psi_\textsf{N}}{g'_{\textsf{N}}}\right) + \frac{f_0}{H_\textsf{N}}h_b\,, \qquad i =\textsf{N}\,,\end{aligned}\end{align} \]

where the reduced gravity, or buoyancy jump, is

\[g'_i \equiv g \frac{\rho_{i+1}-\rho_i}{\rho_i}\,.\]

The evolution equations are

\[\,{q_{i}}_t + \mathsf{J}\left(\psi_i\,, q_i\right) + \textsf{Q}_y {\psi_i}_x - \textsf{Q}_x {\psi_i}_y = \text{ssd} - r_{ek} \delta_{i\textsf{N}} \nabla^2 \psi_i\,, \qquad i = 1,\textsf{N}\,,\]

where the mean potential vorticy gradients are

\[\textsf{Q}_x = \textsf{S}\textsf{V}\,,\]

and

\[\textsf{Q}_y = \beta\,\textsf{I} - \textsf{S}\textsf{U}\,\,,\]

where S is the stretching matrix, I is the identity matrix, and the background velocity is

\(\vec{\textsf{V}}(z) = \left(\textsf{U},\textsf{V}\right)\).

Parameters

nzinteger number

Number of layers (> 1)

betanumber

Gradient of coriolis parameter. Units: meters ^-1 seconds ^-1

rdnumber

Deformation radius. Units: meters. Only necessary for the two-layer (nz=2) case.

deltanumber

Layer thickness ratio (H1/H2). Only necessary for the two-layer (nz=2) case. Unitless.

Ulist of size nz

Base state zonal velocity. Units: meters seconds ^-1

Varray of size nz

Base state meridional velocity. Units: meters seconds ^-1

Harray of size nz

Layer thickness. Units: meters

rho: array of size nz.

Layer density. Units: kilograms meters ^-3

class pyqg.BTModel(beta=0.0, rd=0.0, H=1.0, U=0.0, **kwargs)

Single-layer (barotropic) quasigeostrophic model. This class can represent both pure two-dimensional flow and also single reduced-gravity layers with deformation radius rd.

The equivalent-barotropic quasigeostrophic evolution equations is

\[\partial_t q + J(\psi, q ) + \beta \psi_x = \text{ssd}\]

The potential vorticity anomaly is

\[q = \nabla^2 \psi - \kappa_d^2 \psi\]

Parameters

betanumber, optional

Gradient of coriolis parameter. Units: meters ^-1 seconds ^-1

rdnumber, optional

Deformation radius. Units: meters.

Unumber, optional

Upper layer flow. Units: meters seconds ^-1.

set_U(U)

Set background zonal flow.

Parameters

Unumber

Upper layer flow. Units: meters seconds ^-1.

class pyqg.SQGModel(beta=0.0, Nb=1.0, f_0=1.0, H=1.0, U=0.0, **kwargs)

Surface quasigeostrophic model.

Parameters

betanumber

Gradient of coriolis parameter. Units: meters ^-1 seconds ^-1

Nbnumber

Buoyancy frequency. Units: seconds ^-1.

Unumber

Background zonal flow. Units: meters seconds ^-1.

set_U(U)

Set background zonal flow

Lagrangian Particles

class pyqg.LagrangianParticleArray2D(x0, y0, periodic_in_x=False, periodic_in_y=False, xmin=- inf, xmax=inf, ymin=- inf, ymax=inf, particle_dtype='f8')

A class for keeping track of a set of lagrangian particles in two-dimensional space. Tries to be fast.

Parameters

x0, y0array-like

Two arrays (same size) representing the particle initial positions.

periodic_in_xbool

Whether the domain wraps in the x direction.

periodic_in_ybool

Whether the domain ‘wraps’ in the y direction.

xmin, xmaxnumbers

Maximum and minimum values of x coordinate

ymin, ymaxnumbers

Maximum and minimum values of y coordinate

particle_dtypedtype

Data type to use for particles

step_forward_with_function(uv0fun, uv1fun, dt)

Advance particles using a function to determine u and v.

Parameters

uv0funfunction

Called like uv0fun(x,y). Should return the velocity field u, v at time t.

uv1fun(x,y)function

Called like uv1fun(x,y). Should return the velocity field u, v at time t + dt.

dtnumber

Timestep.

class pyqg.GriddedLagrangianParticleArray2D(x0, y0, Nx, Ny, grid_type='A', **kwargs)

Lagrangian particles with velocities given on a regular cartesian grid.

Parameters

x0, y0array-like

Two arrays (same size) representing the particle initial positions.

Nx, Ny: int

Number of grid points in the x and y directions

grid_type: {‘A’}

Arakawa grid type specifying velocity positions.

interpolate_gridded_scalar(x, y, c, order=1, pad=1, offset=0)

Interpolate gridded scalar C to points x,y.

Parameters

x, yarray-like

Points at which to interpolate

carray-like

The scalar, assumed to be defined on the grid.

orderint

Order of interpolation

padint

Number of pad cells added

offsetint

???

Returns

ciarray-like

The interpolated scalar

step_forward_with_gridded_uv(U0, V0, U1, V1, dt, order=1)

Advance particles using a gridded velocity field. Because of the Runga-Kutta timestepping, we need two velocity fields at different times.

Parameters

U0, V0array-like

Gridded velocity fields at time t - dt.

U1, V1array-like

Gridded velocity fields at time t.

dtnumber

Timestep.

orderint

Order of interpolation.

Diagnostic Tools

Utility functions for pyqg model data.

pyqg.diagnostic_tools.calc_ispec(model, _var_dens, averaging=True, truncate=True, nd_wavenumber=False, nfactor=1)

Compute isotropic spectrum phr from 2D spectrum of variable signal2d such that signal2d.var() = phr.sum() * (kr[1] - kr[0]).

Parameters

modelpyqg.Model instance

The model object from which var_dens originates

var_denssquared modulus of fourier coefficients like this:

np.abs(signal2d_fft)**2/m.M**2

averaging: If True, spectral density is estimated with averaging over circles,

otherwise summation is used and Parseval identity holds

truncate: If True, maximum wavenumber corresponds to inner circle in Fourier space,

otherwise - outer circle

nd_wavenumber: If True, wavenumber is nondimensional:

minimum wavenumber is 1 and corresponds to domain length/width, otherwise - wavenumber is dimensional [m^-1]

nfactor: width of the bin in sqrt(dk^2+dl^2) units

Returns

krarray

isotropic wavenumber

phrarray

isotropic spectrum

pyqg.diagnostic_tools.diagnostic_differences(m1, m2, reduction='rmse', instantaneous=False)

Compute a dictionary of differences in the diagnostics of two models at possibly different resolutions (e.g. for quantifying the effects of parameterizations). Applies normalization/isotropization to certain diagnostics before comparing them and skips others. Also computes differences for each vertical layer separately.

Parameters

m1pyqg.Model instance

The first model to compare

m2pyqg.Model instance

The second model to compare

reductionstring or function

A function that takes two arrays of diagnostics and computes a distance metric. Defaults to the root mean squared difference (‘rmse’).

instantaneousboolean

If true, compute difference metrics for the instantaneous values of a diagnostic, rather than its time average. Defaults to false.

Returns

diffsdict

A dictionary of diagnostic name => distance. If the diagnostic is defined over multiple layers, separate keys are included with an appended z index.

pyqg.diagnostic_tools.diagnostic_similarities(model, target, baseline, **kw)

Like diagnostic_differences, but returning a dictionary of similarity scores between negative infinity and 1 which quantify how much closer the diagnostics of a given model are to a target with respect to a baseline. Scores approach 1 when the distance between the model and the target is small compared to the baseline and are negative when that distance is greater.

Parameters

modelpyqg.Model instance

The model for which we want to compute similiarity scores (e.g. a parameterized low resolution model)

targetpyqg.Model instance

The target model (e.g. a high resolution model)

baselinepyqg.Model instance

The baseline against which we check for improvement or degradation (e.g. an unparameterized low resolution model)

Returns

simsdict

A dictionary of diagnostic name => similarity. If the diagnostic is defined over multiple layers, separate keys are included with an appended z index.

pyqg.diagnostic_tools.spec_sum(ph2)

Compute total spectral sum of the real spectral quantity``ph^2``.

Parameters

modelpyqg.Model instance

The model object from which ph originates

ph2real array

The field on which to compute the sum

Returns

var_densfloat

The sum of ph2

pyqg.diagnostic_tools.spec_var(model, ph)

Compute variance of p from Fourier coefficients ph.

Parameters

modelpyqg.Model instance

The model object from which ph originates

phcomplex array

The field on which to compute the variance

Returns

var_densfloat

The variance of ph

Parameterizations

class pyqg.Parameterization

A generic class representing a subgrid parameterization. Inherit from this class, \(UVParameterization\), or \(QParameterization\) to define a new parameterization.

abstract __call__(m)

Call the parameterization given a pyqg.Model. Override this function in the subclass when defining a new parameterization.

Parameters

mModel

The model for which we are evaluating the parameterization.

Returns

forcingreal array or tuple

The forcing associated with the model. If the model has been initialized with this parameterization as its q_parameterization, this should be an array of shape (nz, ny, nx). For uv_parameterization, this should be a tuple of two such arrays or a single array of shape (2, nz, ny, nx).

abstract property parameterization_type

Whether the parameterization applies to velocity (in which case this property should return "uv_parameterization") or potential vorticity (in which case this property should return "q_parameterization"). If you inherit from UVParameterization or QParameterization, this will be defined automatically.

Returns

parameterization_typestring

Either "uv_parameterization" or "q_parameterization", depending on how the output should be interpreted.

__add__(other)

Add two parameterizations (returning a new object).

Parameters

otherParameterization

The parameterization to add to this one.

Returns

sumParameterization

The sum of the two parameterizations.

__mul__(constant)

Multiply a parameterization by a constant (returning a new object).

Parameters

constantnumber

Multiplicative factor for scaling the parameterization.

Returns

productParameterization

The parameterization times the constant.

class pyqg.parameterizations.UVParameterization

A generic class representing a subgrid parameterization in terms of velocity. Inherit from this to define a new velocity parameterization.

class pyqg.parameterizations.QParameterization

A generic class representing a subgrid parameterization in terms of potential vorticity. Inherit from this to define a new potential vorticity parameterization.

class pyqg.parameterizations.Smagorinsky(constant=0.1)

Velocity parameterization from Smagorinsky 1963.

This parameterization assumes that due to subgrid stress, there is an effective eddy viscosity

\[\nu = (C_S \Delta)^2 \sqrt{2(S_{x,x}^2 + S_{y,y}^2 + 2S_{x,y}^2)}\]

which leads to updated velocity tendencies \(\Pi_{i}, i \in \{1,2\}\) corresponding to \(x\) and \(y\) respectively (equation is the same in each layer):

\[\Pi_{i} = 2 \partial_i(\nu S_{i,i}) + \partial_{2-i}(\nu S_{i,2-i})\]

where \(C_S\) is a tunable Smagorinsky constant, \(\Delta\) is the grid spacing, and

\[S_{i,j} = \frac{1}{2}(\partial_i \mathbf{u}_j + \partial_j \mathbf{u}_i)\]

Parameters

constantnumber

Smagorinsky constant \(C_S\). Defaults to 0.1.

class pyqg.parameterizations.BackscatterBiharmonic(smag_constant=0.08, back_constant=0.99, eps=1e-32)

PV parameterization based on Jansen and Held 2014 and Jansen et al. 2015 (adapted by Pavel Perezhogin). Assumes that a configurable fraction of Smagorinsky dissipation is scattered back to larger scales in an energetically consistent way.

Parameters

smag_constantnumber

Smagorinsky constant \(C_S\) for the dissipative model. Defaults to 0.08.

back_constantnumber

Backscatter constant \(C_B\) describing the fraction of Smagorinsky-dissipated energy which should be scattered back to larger scales. Defaults to 0.99. Normally should be less than 1, but larger values may still be stable, e.g. due to additional dissipation in the model from numerical filtering.

epsnumber

Small constant to add to the denominator of the backscatter formula to prevent division by zero errors. Defaults to 1e-32.

class pyqg.parameterizations.ZannaBolton2020(constant=- 46761284)

Velocity parameterization derived from equation discovery by Zanna and Bolton 2020 (Eq. 6).

Parameters

constantnumber

Scaling constant \(\kappa_{BC}\). Units: meters ^-2. Defaults to \(\approx -4.68 \times 10^7\), a value obtained by empirically minimizing squared error with respect to the subgrid forcing that results from applying the filtering method of Guan et al. 2022 to a two-layer QGModel with default parameters.