Public Function Reference

Approximately Neutral Surfaces

Potential Density Surface

neutralocean.surface.isopycnal.potential_surf(S, T, P, **kwargs)[source]

Calculate a potential density (or specific volume) surface.

Given practical / Absolute salinity S, potential / Conservative temperature T, and pressure (when non-Boussinesq) or depth (when Boussinesq) P, and given a reference pressure / depth ref, calculate an isosurface of eos(S, T, ref) where eos is the equation of state.

Parameters:

S, Tndarray or xarray.DataArray

3D practical / Absolute salinity and potential / Conservative temperature — that is, a 2D array of 1D water columns or “casts”

Pndarray or xarray.DataArray

In the non-Boussinesq case, P is the 3D pressure, sharing the same dimensions as S and T.

In the Boussinesq case, P is the depth and can be 3D with the same structure as S and T, or can be 1D with as many elements as there are in the vertical dimension of S and T.

reffloat

The reference pressure or depth, in the same units as P.

If ref is None, the reference value is taken as pin_p.

See Examples section.

isovalfloat

Isovalue of the potential density surface. Units are same as returned by eos.

See Examples section.

pin_pfloat

Pressure or depth at which the surface intersects the cast pin_cast.

See Examples section.

pin_casttuple or list of int of length 2

Index for cast where surface is at pressure or depth pin_p.

See Examples section.

Returns:

s, t, pndarray or xarray.DataArray

practical / Absolute salinity, potential / Conservative temperature, and pressure / depth on the surface

ddict

Diagnostics.

"e_MAV" : float

Mean Absolute Value of the ϵ neutrality error on the surface, area-weighted. Units are those of eos return values divided by those of dist* inputs.

"e_RMS" : float

As "e_MAV" but for the Root Mean Square.

"n_wet": float

Number of wet casts (surface points).

"timer" : float

Time spent on the whole algorithm, excluding set-up and diagnostics.

"ref" : float

Reference pressure / depth for surface (matching input ref if given, otherwise this is calculated internally).

"isoval" : float

Isovalue of potential density for surface (matching input isoval if given, otherwise this is calculated internally).

Other Parameters:

griddict

Containing the following:

edgestuple of length 2: Each element is an array of int of length E, where E is the number of edges in the grid’s graph, i.e. the number of pairs of adjacent water columns (including land) in the grid. If edges = (a, b), the nodes (water columns) whose linear indices are a[i] and b[i] are adjacent. Required if diags is True
dist1d array: Horizontal distance between adjacent water columns (nodes). dist[i] is the distance between nodes whose linear indices are edges[0][i] and edges[1][i]. If absent, a value of 1.0 is assumed for all edges.
distperp1d array: Horizontal distance of the face between adjacent water columns (nodes). distperp[i] is the distance of the interface between nodes whose linear indices are edges[0][i] and edges[1][i]. If absent, a value of 1.0 is assumed for all edges.
For a rectilinear grid (e.g. latitude-longitude), use: neutralocean.grid.rectilinear.build_grid
For a tiled rectilinear grid, such as works with XGCM, use: neutralocean.grid.xgcm.build_grid
For a general grid given as a graph, use: neutralocean.grid.graph.build_grid

Also see the examples in neutralocean.examples.

vert_dimint or str, Default -1

Specifies which dimension of S, T (and P if 3D) is vertical.

If S and T are ndarray, then vert_dim is the int indexing the vertical dimension of S and T (e.g. -1 indexes the last dimension).

If S and T are xarray.DataArray, then vert_dim is a str naming the vertical dimension of S and T.

Ideally, vert_dim is -1. See Notes.

eosstr or function or tuple of functions, Default ‘gsw’

Specification for the equation of state.

If a str, can be any of the strings accepted by neutralocean.eos.tools.make_eos, e.g. 'jmd95', 'jmdfwg06', 'gsw'.

If a function, must take three inputs corresponding to S, T, and P, and output the density (or specific volume). This form is not allowed when diags is True. This can be made as, e.g., eos = neutralocean.eos.make_eos('gsw') for a non-Boussinesq ocean, or as eos = neutralocean.eos.make_eos('gsw', grav, rho_c) for a Boussinesq ocean with grav and rho_c (see inputs below).

If a tuple of functions, the first element must be a function for the equation of state as above, and the second element must be a function taking the same three inputs as above and returning two outputs, namely the partial derivatives of the equation of state with respect to S and T. The second element can be made as, e.g., eos_s_t = neutralocean.eos.make_eos_s_t('gsw', grav, rho_c)

The function (or the first element of the tuple of functions) should be @numba.njit decorated and need not be vectorized – it will be called many times with scalar inputs.

gravfloat, Default None

Gravitational acceleration [m s-2]. When non-Boussinesq, pass None.

rho_cfloat, Default None

Boussinesq reference density [kg m-3]. When non-Boussinesq, pass None.

interpstr, Default ‘linear’

Method for vertical interpolation. Use 'linear' for linear interpolation, and 'pchip' for Piecewise Cubic Hermite Interpolating Polynomials. Other interpolants can be added through the subpackage, ppinterp.

diagsbool, Default True

If True, calculate diagnostics (4th output). If False, 4th output is an empty dict.

outputbool, Default True

If True, prints diagnostic output during computation. diags must be True for this to have any effect. To redirect this output to a file, do the following >>> import sys >>> tmp = sys.stdout >>> sys.stdout = file_id = open(‘myfile.txt’, ‘w’) >>> # Now call this function … >>> sys.stdout = tmp >>> file_id.close()

TOL_P_SOLVERfloat, Default 1e-4

Error tolerance when root-finding to update the pressure or depth of the surface in each water column. Units are the same as P.

Notes

This code will internally re-arrange S, T, P to have the vertical dimension last, so that the data for an individual water column is contiguous in memory. If you call this function many times, consider using lib._process_casts to pre-process your S, T, P inputs to have the vertical dimension last.

Examples

The output surface must be specified by some combination of reference pressure / depth ref, isovalue isoval, pinning cast pin_cast and pinning pressure pin_p. The following methods are valid, and listed in order of precedence (e.g. pin_cast and pin_p are not used if both ref and isoval are given).

>>> potential_surf(S, T, P, ref=___, isoval=___, ...)

This finds the surface with given reference pressure / depth and the given isovalue.

>>> potential_surf(S, T, P, ref=___, pin_p=___, pin_cast=___, ...)

This finds the surface with given reference pressure / depth that intersects the given cast at the given pressure or depth.

>>> potential_surf(S, T, P, pin_p=___, pin_cast=___, ...)

This is as for the previous method, but selects the reference pressure / depth as pin_p (i.e. the local P value at the given cast’s given pressure or depth).

Specific Volume Anomaly Surface

neutralocean.surface.isopycnal.anomaly_surf(S, T, P, **kwargs)[source]

Calculate a specific volume (or in-situ density) anomaly surface.

Given practical / Absolute salinity S, potential / Conservative temperature T, and pressure / depth P, and given reference values S0 and T0, calculate an isosurface of eos(S, T, P) - eos(S0, T0, P) where eos is the equation of state.

In a non-Boussinesq ocean, P is pressure. Also, if one is computing geostrophic streamfunctions, it is most convenient if eos provides the specific volume.

In a Boussinesq ocean, P is depth. Also, if one is computing geostrophic streamfunctions, it is most convenient if eos provides the in-situ density.

Parameters:

S, T, Pndarray or xarray.DataArray

See potential_surf

reftuple of float of length 2

The reference S and T values.

If ref is None or has a None element, the reference values are taken from the local S and T at the pressure or depth pin_p on the pinning cast pin_cast.

isoval, pin_p, pin_cast

See potential_surf

Returns:

s, t, p, d: See potential_surf. Note d["ref"] returns a 2 element tuple, namely ref as here.

Other Parameters:

grid, vert_dim, eos, grav, rho_c, interp, diags, output, TOL_P_SOLVER: See potential_surf

Notes

See potential_surf.

Examples

The output surface must be specified by some combination of reference salinity and temperature ref, isovalue isoval, pinning cast pin_cast and pinning pressure pin_p. The following methods are valid, and listed in order of precedence (e.g. pin_cast and pin_p are not used if both ref and isoval are given).

>>> anomaly_surf(S, T, P, ref=___, isoval=___, ...)

This finds the surface with given reference salinity and temperature and the given isovalue.

>>> anomaly_surf(S, T, P, ref=___, pin_p=___, pin_cast=___, ...)

This finds the surface with given reference salinity and temperature that intersects the given cast at the given pressure or depth.

>>> anomaly_surf(S, T, P, pin_p=___, pin_cast=___, ...)

This is as for the previous method, but selects the reference salinity and temperature from the local S and T values at the given cast’s given pressure or depth.

Omega Surface

neutralocean.surface.omega.omega_surf(S, T, P, grid, **kw)[source]

Calculate an omega surface from structured ocean data.

Given 3D salinity, temperature, and pressure or depth data arranged on a rectilinear grid, calculate a 2D omega surface [1] [2], which is a highly accurate approximately neutral surface.

Parameters:

S, T, Pndarray or xarray.DataArray

See potential_surf

p_initndarray, Default None

Pressure or depth on the initial approximately neutral surface.

See Examples section.

reffloat, or tuple of float of length 2

If p_init is None, the reference value(s) for the initial potential density surface or in-situ density (specific volume) anomaly surface that initializes the omega surface algorithm. If ref is a scalar, a potential density urface is used, and if ref is None, the reference P is pin_p (i.e. taken local to the pinning location). If ref is a tuple of length two, a in-situ density anomaly surface is used, and if ref is (None, None), then the reference S and T values are taken local to the pinning location (pressure or depth pin_p on the pinning cast pin_cast).

See Examples section.

isovalfloat

Isovalue for the initial potential density or in-situ density anomaly surface when p_init is not given. Units are same as returned by the function specified by eos.

See Examples section.

pin_p, pin_cast

See potential_surf

Returns:

s, t, pndarray or xarray.DataArray

practical / Absolute salinity, potential / Conservative temperature, and pressure / depth on surface

ddict

Diagnostics. The first four listed below give information going into the i’th iteration (e.g. the 0’th element is about the initial surface). The rest give information about what the i’th iteration did (and hence their 0’th elements are irrelevant).

"e_MAV" : array of float

Mean Absolute Value of the ϵ neutrality error on the surface, area-weighted. Units are those of eos return values divided by those of dist* inputs.

"e_RMS" : array of float

As "e_MAV" but for the Root Mean Square.

"n_wet": array of float

Number of wet casts (surface points).

"timer" : array of float

Time spent on each iteration, excluding set-up (approximately) and diagnostics.

"ϕ_MAV" : array of float

Mean Absolute Value of the Locally Referenced Potential Density perturbation, per iteration

"Δp_MAV" : array of float

Mean Absolute Value of the pressure or depth change from one iteration to the next

"Δp_RMS" : array of float

Root Mean Square of the pressure or depth change from one iteration to the next

"Δp_Linf" : array of float

Maximum absolute value (infinity norm) of the pressure or depth change from one iteration to the next

"n_newly_wet" : array of int

Number of casts that are newly wet, per iteration

"timer_bfs" : array of float

Time spent in Breadth-First Search including wetting, per iteration.

"timer_mat" : array of float

Time spent building and solving the matrix problem, per iteration.

"timer_update" : array of float

Time spent vertically updating the surface.

Other Parameters:

grid, vert_dim, grav, rho_c, interp, diags, output, TOL_P_SOLVER

See potential_surf

eosstr or tuple of functions, Default ‘gsw’

As in potential_surf, excluding the option to pass a single function. The omega surface algorithm relies on knowing the partial derivatives of the equation of state with respect to salinity and temperature, so the eos_s_t function is also required if given as a tuple of functions.

ITER_MINint, Default 1

Minimum number of iterations.

ITER_MAXint, Default 10

Maximum number of iterations.

ITER_START_WETTINGint, Default 1

Iteration on which wetting begins. Set to np.inf (ITER_MAX + 1 would also do) to deactivate.

ITER_STOP_WETTINGint, Default 5

The last iteration on which to perform wetting. This can be useful to avoid pesky water columns that repeatedly wet then dry.

TOL_LRPD_MAVfloat, Default 1e-7

Exit iterations when the mean absolute value of the Locally Referenced Potential Density perturbation that updates the surface from one iteration to the next is less than this value. Units are [kg m-3], even if eos returns a specific volume. Set to 0 to deactivate.

TOL_P_CHANGE_RMSfloat, Default 0.0

Exit iterations when the root mean square of the pressure or depth change on the surface from one iteration to the next is less than this value. Set to 0 to deactivate. Units are the same as P [dbar or m].

p_mlndarray or dict, Default None

If a dict, the pressure or depth at the base of the mixed layer is computed using mixed_layer with p_ml passed as keyword arguments, enabling control over the parameters in that function. See mixed_layer for details.

If an ndarray (of the same shape as the lateral dimensions of S), the pressure or depth at the base of the mixed layer in each water column.

When the surface’s pressure is shallower than p_ml in any water column, it is set to NaN (a “dry” water column). This is not applied to the initial surface, but only to the surface after the first iteration, as the initial surface could be very far from neutral.

If None, the mixed layer is not removed.

OMEGA_FORMULATIONstr, Default ‘poisson’

Specify how the matrix problem is set up and solved. Options are

'poisson', to solve the Poisson problem as in [1] with Cholesky, or
'gradient', to solve the overdetermined gradient equations as in [2] using LSQR.

Notes

See potential_surf Notes.

[1] (1,2)

Stanley, G. J., McDougall, T. J., & Barker, P. M. (2021). Algorithmic Improvements to Finding Approximately Neutral Surfaces. Journal of Advances in Modeling Earth Systems, 13(5), e2020MS002436.

[2] (1,2)

Klocker, A., McDougall, T. J., & Jackett, D. R. (2009). A new method for forming approximately neutral surfaces. Ocean Science, 5 (2), 155-172.

Examples

omega surfaces require a pinning cast and initial surface. The surface is iteratively updated while remaining fixed at the pinning cast. The initial surface can be provided directly, as the surface with pressure or depth given by p_init, in the following method:

>>> omega_surf(S, T, P, pin_cast, p_init, ...)

Alternatively, a potential density surface or an in-situ density (specific volume) anomaly surface can be used as the initial surface. To do this, use one of the following two methods:

>>> omega_surf(S, T, P, ref, isoval, pin_cast, ...)

>>> omega_surf(S, T, P, ref, pin_p, pin_cast, ...)

For more info on these methods, see the Examples section of potential_surf. Note that pin_cast is always a required input. Note that ref is needed to distinguish which of the two types of isopycnals will be used as the initial surface.

Grids

Rectilinear

neutralocean.grid.rectilinear.build_grid(dims, periodic, dxC=1.0, dyC=1.0, dyG=1.0, dxG=1.0)[source]

Build pairs of adjacent grid points and associated geometry for a rectilinear grid

Parameters:

dimstuple of int: Dimensions of the grid. dims[n] is the number of grid cells in the n’th direction, for n = 1,2.
periodictuple of bool: Specifies periodicity. The n’th dimension is periodic when periodic[n] is True.
dxCfloat or ndarray, Default 1.0: Distance between adjacent grid points in the 1st (‘x’) dimension. Specifically, dxC[i,j] is the distance between the centers of cells (i,j) and (i-1,j).
dyCfloat or ndarray, Default 1.0: Distance between adjacent grid points in the 2nd (‘y’) dimension. Specifically, dyC[i,j] is the distance between the centers of cells (i,j) and (i,j-1).
dyGfloat or ndarray, Default 1.0: Distance (in the 2nd, ‘y’ dimension) of the face between grid points that are adjacent in the 1st (‘x’) dimension. Lives at same location as dxC. Specifically, dyG[i,j] is the distance of the interface between cells (i,j) and (i-1,j).
dxGfloat or ndarray, Default 1.0: Distance (in the 1st, ‘x’ dimension) of the face between grid points that are adjacent in the 2nd (‘y’) dimension. Lives at same location as dyC. Specifically, dxG[i,j] is the distance of the interface between cells (i,j) and (i,j-1).

Returns:

griddict

Containing the following:

edgestuple of length 2: Each element is an array of int of length E, where E is the number of edges in the grid’s graph, i.e. the number of pairs of adjacent water columns (including land) in the grid. If edges = (a, b), the nodes (water columns) whose linear indices are a[i] and b[i] are adjacent.
dist1d array: Horizontal distance between adjacent water columns (nodes). dist[i] is the distance between nodes whose linear indices are edges[0][i] and edges[1][i].
distperp1d array: Horizontal distance of the face between adjacent water columns (nodes). distperp[i] is the distance of the interface between nodes whose linear indices are edges[0][i] and edges[1][i].

Notes

The above notes for dxC, dyC, dyG, dxG are valid for 0 <= i <= ni-1 and 0 <= j <= nj-1, where dims == (ni, nj), using Python’s -1 indexing notation:

dxC[0,j] is the distance between the centers of cells (0,j) and (ni-1,j)
dyG[0,j] is the distance of the face between cells (0,j) and (ni-1,j)
dyC[0,j] is the distance between the centers of cells (i,0) and (i,nj-1)
dxG[0,j] is the distance of the face between cells (i,0) and (i,nj-1)

If periodic[0] is False, dxC[0,j] and dyG[0,j] are not used.

If periodic[1] is False, dyC[i,0] and dxG[i,0] are not used.

Broadcasting: The shape of dxC, dyC, dyG, dxG will be broadcast to shape dims == (ni, nj).

Example on a latitude-longitude grid: the meridional grid distances are constant, so dyC and dyG can be passed as scalars rather than 2D arrays; likewise, the zonal grid distances only depend on latitude, so dxC and dxG can be passed as 1D arrays of length nj.

neutralocean.grid.rectilinear.edgedata_to_maps(edgedata, dims, periodic)[source]

Convert 1D array of data living on edges into two 2D arrays, one for each spatial dimension, for a rectilinear grid

Parameters:

edgedataarray: 1D array of data that lives on edges of the grid’s graph, given in the same order as the edges cosntructed from _build_edges. If a, b = _build_edges(dims, periodic), then edgedata[i] is the data that lives at the interface between the grid cells indexed by a[i] and b[i].
dims, periodic: See build_grid

Returns:

Fi, Fjndarray: 2D arrays of shape dims. Fi contains the values from edgedata that correspond to data living between grid points in the first dimension, and similarly for Fj but in the second dimension.

xgcm (tiled rectilinear)

neutralocean.grid.xgcm.build_grid(n, face_connections, dims, xsh, ysh, dxC=1.0, dyC=1.0, dyG=1.0, dxG=1.0)[source]

Make a list of edges between adjacent nodes in a graph consisting of rectilinear tiles, with adjacency between tiles specified as for xgcm; also generate lists of distances associated with these edges.

Here, the grid is thought of as a graph, where each node represents a water column, and the graph has an edge between two nodes when the corresponding water columns are adjacent in physical space.

Parameters:

nint

Number of grid points on each side of a square tile. (Note xgcm only handles square tiles.)

face_connectionsdict

xgcm-style specification for how tiles (faces) are connected. See https://xgcm.readthedocs.io/en/latest/grid_topology.html

dimstuple of str

A tuple of length 3, with one entry being the same as the only key in the face_connections dict, and the other two being ‘i’ and ‘j’.

This names, in order, the dimensions of a DataArray whose data lives on each (tracer) grid point in the horizontal, and has no vertical dimension.

If D is a DataArray giving the Sea Surface Temperature, for example, then this argument should be D.dims.

Example: in ECCOv4r4, dims == ('tile', 'j', 'i').

xshstr

The direction of shifting in the ‘i’ dimension. Can be either “left” or “right”.

If “left”, then dxC is the distance between the local grid point (i, j) and the previous grid point in the ‘i’ dimension, at (i-1, j).

If “right”, then dxC is the distance between the local grid point (i, j) and the next grid point in the ‘i’ dimension, at (i+1, j).

yshstr

The direction of shifting in the ‘j’ dimension. Can be either “left” or “right”.

If “left”, then dyC is the distance between the local grid point (i, j) and the previous grid point in the ‘j’ dimension, at (i, j+1).

If “right”, then dyC is the distance between the local grid point (i, j) and the next grid point in the ‘j’ dimension, at (i, j+1).

dxCfloat or ndarray, Default 1.0

Distance between adjacent grid points in the ‘i’ dimension.

dyCfloat or ndarray, Default 1.0

Distance between adjacent grid points in the ‘j’ dimension

dyGfloat or ndarray, Default 1.0

Distance (in the ‘j’ dimension) of the face between grid points that are adjacent in the ‘i’ dimension. Lives at the same location as dxC.

dxGfloat or ndarray, Default 1.0

Distance (in the ‘i’ dimension) of the face between grid points that are adjacent in the ‘j’ dimension. Lives at the same location as dyC.

Returns:

griddict

Containing the following:

edgestuple of length 2: Each element is an array of int of length E, where E is the number of edges in the grid’s graph, i.e. the number of pairs of adjacent water columns (including land) in the grid. If edges = (a, b), the nodes (water columns) whose linear indices are a[i] and b[i] are adjacent.
dist1d array: Horizontal distance between adjacent water columns (nodes). dist[i] is the distance between nodes whose linear indices are edges[0][i] and edges[1][i].
distperp1d array: Horizontal distance of the face between adjacent water columns (nodes). distperp[i] is the distance of the interface between nodes whose linear indices are edges[0][i] and edges[1][i].

neutralocean.grid.xgcm.edgedata_to_maps(edgedata, n, face_connections, dims, xsh, ysh)[source]

Convert 1D array of data living on edges into two nD arrays, one for each spatial dimension, for a rectilinear grid

Parameters:

edgedataarray: 1D array of data that lives on edges of the grid’s graph, given in the same order as the edges cosntructed from build_grid. If a, b, ... = build_grid(...), then edgedata[i] is the data that lives at the interface between the grid cells indexed by a[i] and b[i].
n, face_connections, dims, xsh, ysh: See build_grid

Returns:

Fi, Fjndarray: Fi contains the values from edgedata that correspond to data living between grid points in the first dimension, and similarly for Fj but in the second dimension.

tripolar

neutralocean.grid.tripolar.build_grid(dims, e1u=1.0, e2v=1.0, e2u=1.0, e1v=1.0)[source]

Build pairs of adjacent grid points and associated geometry for a tripolar grid

Parameters:

dimstuple of int: Dimensions of the grid. Letting dims == (nj, ni), then nj is the number of grid cells in the y (latitude-like) direction, & ni is the number of grid cells in the x (longitude-like) direction.
e1ufloat or 2D array, Default 1.0: Distance in the x dimension that lives on the U cell.
e2vfloat or 2D array, Default 1.0: Distance in the y dimension that lives on the V cell.
e2ufloat or 2D array, Default 1.0: Distance in the y dimension that lives on the U cell.
e1vfloat or 2D array, Default 1.0: Distance in the x dimension that lives on the V cell.

Returns:

griddict

Containing the following:

edgestuple of length 2: Each element is an array of int of length E, where E is the number of edges in the grid’s graph, i.e. the number of pairs of adjacent water columns (including land) in the grid. If edges = (a, b), the nodes (water columns) whose linear indices are a[i] and b[i] are adjacent.
dist1d array: Horizontal distance between adjacent water columns (nodes). dist[i] is the distance between nodes whose linear indices are edges[0][i] and edges[1][i].
distperp1d array: Horizontal distance of the face between adjacent water columns (nodes). distperp[i] is the distance of the interface between nodes whose linear indices are edges[0][i] and edges[1][i].

Notes

e1u, e2v, e2u, e1v can have shape (nj,ni) or (nj+1, ni+2). The latter form has 1 cell of padding in the y dimension to handle the tripolar grid, and 1 cell of padding in the x dimension to handle the longitudinal periodicity. If given the latter form, this extra padding is removed internally. When of shape (nj, ni),

e1u[j,i] is the x distance on the U grid at [j,i+1/2],
e2v[j,i] is the y distance on the V grid at [j+1/2,i],
e2u[j,i] is the y distance on the U grid at [j,i+1/2],
e1v[j,i] is the x distance on the V grid at [j+1/2,i].

neutralocean.grid.tripolar.edgedata_to_maps(edgedata, dims)[source]

Convert 1D array of data living on edges into two 2D arrays, one for each spatial dimension, for a tripolar grid

Parameters:

edgedataarray: 1D array of data that lives on edges of the grid’s graph, given in the same order as the edges cosntructed from _build_edges. If a, b = _build_edges(dims, periodic), then edgedata[i] is the data that lives at the interface between the grid cells indexed by a[i] and b[i].
dims: See build_grid

Returns:

v, u2D array: 2D arrays of shape dims. v and u contains the values from edgedata that correspond to data living on the V and U grids, respectively.

Graph

neutralocean.grid.graph.build_grid(G)[source]

Convert a dict of (undirected) graphs into a grid dict used by neutralocean.

Parameters:

Gmatrix or dict

If a (sparse) matrix, this specifies the edges in the graph:: an edge exists between nodes m and n iff G[m,n] is nonzero or G[n,m] is nonzero. If G is symmetric, only its upper triangle is used. If G is not symmetric, it must have only one of G[m,n] or G[n,m] be non-zero; this condition is not checked for. The geometric distances are both taken to be 1.0, regardless of the data in G.
If a dict, must have 'dist' and 'distperp' entries, which are both: (sparse) matrices with the same sparsity structure. The edges in the graph are determined from G['dist'] as above. The distance between adjacent nodes m and n is given by G['dist'][m,n] or G['dist'][n,m]. The distance of the interface between adjacent nodes m and n is given by G['distperp'][m,n] or G['distperp'][n,m].

Returns:

griddict

grid[‘edges’]tuple of length 2: Each element is an array of int of length E, where E is the number of edges in the grid’s graph, i.e. the number of pairs of adjacent water columns (including land) in the grid. If grid['edges'] = (a, b), the nodes (water columns) whose linear indices are a[i] and b[i] are adjacent.
grid[‘dist’]1d array: grid['dist'][i] is the distance between nodes a[i] and b[i], where grid['edges'] = (a, b).
grid[‘distperp’]1d array: grid['distperp'][i] is the distance of the interface between nodes a[i] and b[i], where grid['edges'] = (a, b).

neutralocean.grid.graph.edges_to_graph(edges, N=None, weights=None)[source]

Convert list of pairs of adjacent nodes into a sparse matrix graph representation

Parameters:

edgestuple: Each element is an array of int of length E, where E is the number of edges in the grid’s graph, i.e. the number of pairs of adjacent water columns (including land) in the grid. If edges = (a, b), the nodes (water columns) whose linear indices are a[i] and b[i] are adjacent.
Nint, Default None: Number of nodes in the graph. If None, the number of nodes is assumed to be the largest value in edges, plus one since elements in edges are zero-indexed, i.e. nodes are laballed 0, 1, …, N-1.
weightsarray, Default None: Weights of the edges in the graph. If None, a weight of 1 is given for each edge.

Returns:

Gcsr_matrix: Symmetric, sparse, N by N matrix representation of the undirected graph. G[m,n] = G[n,m] = weights[i] where i is such that edges[0][i] = m and edges[1][i] = n.

neutralocean.grid.graph.graph_binary_fcn(edges, nodevals, binary_fcn)[source]

Apply a binary function to all pairs of data for all pairs of adjacent nodes in a graph specified by a sparse COO matrix.

Parameters:

edges2D array of int of shape (2, E): The edges in a graph whose nodes are labelled 0, ..., N-1. Denoting a = edges[0] and b = edges[1], - the nodes labelled a[i] and b[i] are adjacent, for 0 <= i <= E-1 - must have 0 <= a[i], b[i] <= N-1 for 0 <= i <= E-1, so edges is referring to valid nodes.
nodevalsarray-like: values at the nodes in a graph
binary_fcnfunction: A numba.njit function accepting two numbers and returning one number.

Returns:

ev1d array: The value of binary_fcn evaluated on each edge. That is, ev[e] = binary_fcn(data[i], data[j]) where i = row[e] and j = col[e]. If either i or j is < 0, ev[e] is nan.

Notes

edges can instead be a tuple of length 2, with each element a tuple of int of length E. If edges is a 2D array of shape (E, 2) or is a tuple of length E with each element a tuple of int of length 2, then one must convert it as follows:

>>> edges = numpy.array(edges).transpose()

Neutrality Error

neutralocean.ntp.ntp_epsilon_errors(s, t, p, grid, eos_s_t='gsw', grav=None, rho_c=None)[source]

Calculate epsilon neutrality errors on an approximately neutral surface

Parameters:

s, t, p, eos_s_t, grav, rho_c: See ntp_epsilon_errors_norms.
griddict: See ntp_epsilon_errors_norms. Need not have the distperp element, as this is not used. If the dist element is missing, a value of 1.0 will be used. Can alternatively pass a 2 element tuple that is just grid['edges'], in which case dist will be taken as 1.0.

Returns:

earray: The epsilon neutrality errors on the surface. e[i] is the neutrality error from node a[i] to node b[i], where a, b = grid['edges'].

neutralocean.ntp.ntp_epsilon_errors_norms(s, t, p, grid, eos_s_t='gsw', grav=None, rho_c=None)[source]

Calculate norms of the epsilon neutrality errors on an approximately neutral surface

Parameters:

s, t, pndarray

practical / Absolute Salinity, potential / Conservative temperature, and pressure / depth on the surface

griddict

Containing the following:

edgestuple of length 2, Required: Each element is an array of int of length E, where E is the number of edges in the grid’s graph, i.e. the number of pairs of adjacent water columns (including land) in the grid. If edges = (a, b), the nodes (water columns) whose linear indices are a[i] and b[i] are adjacent.
dist1d array, Default 1.0: Horizontal distance between adjacent water columns (nodes). dist[i] is the distance between nodes whose linear indices are edges[0][i] and edges[1][i].
distperp1d array, Default 1.0: Horizontal distance of the face between adjacent water columns (nodes). distperp[i] is the distance of the interface between nodes whose linear indices are edges[0][i] and edges[1][i].

eos_s_tfunction

Equation of state for the partial derivatives of density or specific volume with respect to S and T as a function of S, T, and P inputs.

gravfloat, Default None

Gravitational acceleration [m s-2]. When non-Boussinesq, pass None.

rho_cfloat, Default None

Boussinesq reference density [kg m-3]. When non-Boussinesq, pass None.

Returns:

e_RMSfloat: Area-weighted root-mean-square of the epsilon neutrality error on the surface
e_MAVTYPE: Area-weighted mean-absolute-value of the epsilon neutrality error on the surface

Neutral Trajectory

neutralocean.traj.ntp_bottle_to_cast(sB, tB, pB, S, T, P, tol_p=0.0001, interp='linear', eos='gsw', grav=None, rho_c=None)[source]

Find the neutral tangent plane from a bottle to a cast

Finds a point on the cast salinity, temperature, and pressure (S, T, P) where the salinity, temperature, and pressure (s, t, p) is neutrally related to a bottle of salinity, temperature and pressure (sB, tB, pB). That is, the density of (s, t, p_avg) very nearly equals the density of (sB, tB, p_avg), where p_avg = (p + pB) / 2. Within tol_p of this point on the cast, there is a point where these two densities are exactly equal.

Parameters:

sB, tB, pBfloat: practical / Absolute salinity, potential / Conservative temperature, and pressure or depth of the bottle
S, T, P1D array of float: practical / Absolute salinity, potential / Conservative temperature, and pressure or depth of data points on the cast. P must increase monotonically along its last dimension.

Returns:

s, t, pfloat: practical / Absolute salinity, potential / Conservative temperature, and pressure or depth at a point on the cast that is very nearly neutrally related to the bottle.

Other Parameters:

tol_pfloat, Default 1e-4

Error tolerance in terms of pressure or depth when searching for a root of the nonlinear equation. Units are the same as P.

interpstr, Default ‘linear’

Method for vertical interpolation. Use 'linear' for linear interpolation, and 'pchip' for Piecewise Cubic Hermite Interpolating Polynomials. Other interpolants can be added through the subpackage, ppinterp.

eosstr or function, Default ‘gsw’

The equation of state for the density or specific volume as a function of S, T, and pressure (if non-Boussinesq) or depth(if Boussinesq).

If a str, can be any of the strings accepted by neutralocean.eos.tools.make_eos, e.g. 'jmd95', 'jmdfwg06', 'gsw'.

If a function, this should be @numba.njit decorated and need not be vectorized, as it will be called many times with scalar inputs.

gravfloat, Default None

Gravitational acceleration [m s-2]. When non-Boussinesq, pass None.

rho_cfloat, Default None

Boussinesq reference density [kg m-3]. When non-Boussinesq, pass None.

neutralocean.traj.neutral_trajectory(S, T, P, p0, vert_dim=-1, tol_p=0.0001, interp='linear', eos='gsw', grav=None, rho_c=None)[source]

Calculate a neutral trajectory through a sequence of casts.

Given a sequence of casts with hydrographic properties (S, T, P), calculate a neutral trajectory starting from the first cast at pressure p0, or starting from a bottle prior to the first cast with hydrographic properties (s0, t0, p0).

Parameters:

S, T, P2D ndarray or xarray: 1D data specifying the practical / Absolute salinity, and potential / Conservative temperature, and pressure / depth down a 1D sequence of casts. The first dimension specifies the cast number, while the second provides data on that cast; e.g. S[i, :] is the salinity down cast i.
p0float: The pressure / depth at which to begin the neutral trajectory on the first cast

Returns:

s, t, p1D ndarray: practical / Absolute Salinity, potential / Conservative Temperature, and pressure / depth along the neutral trajectory.

Other Parameters:

vert_dimint or str, Default -1

Specifies which dimension of S, T (and P if 3D) is vertical.

If S and T are ndarray, then vert_dim is the int indexing the vertical dimension of S and T (e.g. -1 indexes the last dimension).

If S and T are xarray.DataArray, then vert_dim is a str naming the vertical dimension of S and T.

Ideally, vert_dim is -1. See Notes in potential_surf.

tol_p, interp, eos, grav, rho_c

See ntp_bottle_to_cast

Veronis Density

Mixed Layer

neutralocean.mixed_layer.mixed_layer(S, T, P, eos='gsw', grav=None, rho_c=None, pot_dens_diff=0.03, ref_p=100.0, bottle_index=1, interp='linear', vert_dim=-1)[source]

Calculate the pressure or depth at the bottom of the mixed layer

The mixed layer is the pressure or depth at which the potential density (referenced to ref_p) exceeds the potential density near the surface (specifically, the bottle indexed by bottle_index on each cast) by an amount of pot_dens_diff.

Parameters:

S, Tndarray or xarray.DataArray

3D practical / Absolute salinity and potential / Conservative temperature — that is, a 2D array of 1D water columns or “casts”

Pndarray or xarray.DataArray

In the non-Boussinesq case, P is the 3D pressure, sharing the same dimensions as S and T.

In the Boussinesq case, P is the depth and can be 3D with the same structure as S and T, or can be 1D with as many elements as there are in the vertical dimension of S and T.

Must increase monotonically along the first dimension (i.e. downwards).

eosstr or function, Default ‘gsw’

Specification for the equation of state.

If a str, can be any of the strings accepted by neutralocean.eos.tools.make_eos, e.g. 'jmd95', 'jmdfwg06', 'gsw'.

If a function, must take three inputs corresponding to S, T, and P, and output the density (or specific volume). This form is not allowed when diags is True. This can be made as, e.g., eos = neutralocean.eos.make_eos('gsw') for a non-Boussinesq ocean, or as eos = neutralocean.eos.make_eos('gsw', grav, rho_c) for a Boussinesq ocean with grav and rho_c (see inputs below).

The function should be @numba.njit decorated and need not be vectorized – it will be called many times with scalar inputs.

gravfloat, Default None

Gravitational acceleration [m s-2]. When non-Boussinesq, pass None.

rho_cfloat, Default None

Boussinesq reference density [kg m-3]. When non-Boussinesq, pass None.

pot_dens_difffloat, Default 0.03

Difference in potential density [kg m-3] between the near-surface and the base of the mixed layer.

ref_pfloat, Default 100.0

Reference pressure or depth for potential density [dbar or m].

bottle_indexint, Default 1

The index for the bottle on each cast where the “near surface” potential density is calculated. Note this index is 0-based. The Default of 1 therefore indexes the second bottle in each cast.

interpstr, Default ‘linear’

Method for vertical interpolation. Use 'linear' for linear interpolation, and 'pchip' for Piecewise Cubic Hermite Interpolating Polynomials. Other interpolants can be added through the subpackage, interp1d.

vert_dimint or str, Default -1

Specifies which dimension of S, T (and P if 3D) is vertical.

If S and T are ndarray, then vert_dim is the int indexing the vertical dimension of S and T (e.g. -1 indexes the last dimension).

If S and T are xarray.DataArray, then vert_dim is a str naming the vertical dimension of S and T.

Ideally, vert_dim is -1. See Notes section of potential_surf.

Returns:

MLndarray: A 2D array giving the pressure [dbar] or depth [m, positive] of the mixed layer

Equation Of State

Tools

neutralocean.eos.tools.make_eos(eos, grav=None, rho_c=None)[source]

Make an equation of state function, possibly modified for Boussinesq

Parameters:

eosstr or function

If a str, can be 'gsw' to generate the TEOS-10 specific volume [1], 'jmd95' to generate the Jackett and McDougall (1995) in-situ density [2], or 'jmdfwg06' to generate the Jackett et al (2006) in-situ density [3].

If a function, should be an equation of state as a function of practical / Absolute salinity, potential / Conservative temperature, and pressure (for non-Boussesinq) / depth (for Boussinesq).

grav, rho_cfloat

Gravitational acceleration [m s-2] and Boussinesq reference density [kg m-3]. If both are provided, the equation of state is modified as appropriate for the Boussinesq approximation, in which the third argument is depth, not pressure. Specifically, a depth z is converted to 1e-4 * grav * rho_c * z, which is the hydrostatic pressure [dbar] at depth z [m] caused by a water column of density rho_c under gravity grav.

Returns:

eosfunction: The desired equation of state.

[1]

McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127, SBN 978-0-646-55621-5.

[2]

Jackett and McDougall, 1995, JAOT 12(4), pp. 381-388 ..

[3]

Jackett, D. R., McDougall, T. J., Feistel, R., Wright, D. G., & Griffies, S. M. (2006). Algorithms for Density, Potential Temperature, Conservative Temperature, and the Freezing Temperature of Seawater. Journal of Atmospheric and Oceanic Technology, 23(12), 1709–1728. https://doi.org/10.1175/JTECH1946.1

neutralocean.eos.tools.make_eos_s_t(eos, grav=None, rho_c=None)[source]

Make a function for the partial S and T derivatives of an equation of state

Parameters:

eos, grav, rho_c: See make_eos

Returns:

eos_s_tfunction: Function returning two outputs, namely the partial derivatives with respect to its first two arguments (practical / Absolute salinity and potential / Conservative temperature) of the desired equation of state.

neutralocean.eos.tools.make_eos_p(eos, grav=None, rho_c=None)[source]

Make a function for the partial P derivative of an equation of state

Parameters:

eos, grav, rho_c: See make_eos

Returns:

eos_pfunction: Function returning the partial derivative with respect to the third argument (pressure) of the desired equation of state.

neutralocean.eos.tools.vectorize_eos(eos)[source]

Convert an eos function that takes scalar inputs into one taking arrays.

Parameters:

eosfunction: Any function taking three scalar inputs and returning one scalar output, such as the equation of state. Note this does not work for functions returning multiple outputs (i.e. a tuple), such as a function returning the partial derivatives of the equation of state.

Returns:

eos_vecfunction: A @numba.vectorize’d version of eos, which can take array inputs and returns one array output. The array inputs’ shape need not match exactly, but must be broadcastable to each other.

JMD95

Density of Sea Water using the Jackett and McDougall 1995 [1] function

Functions:

rho :: computes in-situ density from salinity, potential temperature and: pressure
rho_s_t :: compute the partial derivatives of in-situ density with: respect to salinity and potential temperature
rho_p :: compute the partial derivative of in-situ density with: respect to pressure

Notes: To make Boussinesq versions of these functions, see neutralocean.eos.tools.make_eos_bsq.

To make vectorized versions of these functions, see neutralocean.eos.tools.vectorize_eos.

[1]

Jackett and McDougall, 1995, JAOT 12[4], pp. 381-388

neutralocean.eos.jmd95.rho(s, t, p)[source]

Fast JMD95 [1] in-situ density.

Parameters:

sfloat: Practical salinity [PSS-78]
tfloat: Potential temperature [IPTS-68]
pfloat: Pressure [dbar]

Returns:

rhofloat: JMD95 in-situ density [kg m-3]

Notes

This function is derived from densjmd95.m, documented below. Input checks and expansion of variables have been removed. That code used arrays to store coefficients, and also began by multiplying the input pressure by 10 (dbar -> bar). The coefficients, modified to not need this multiplication by 10, have been taken out of arrays and simply hardcoded into the later expressions. The polynomial calculations have also been optimized to favour faster, nested multiplications.

As such, the output of this function differs from the output of the original densjmd95.m function, though the difference is at the level of machine precision.

% DENSJMD95    Density of sea water
%=========================================================================
%
% USAGE:  dens = densjmd95(S,Theta,P)
%
% DESCRIPTION:
%    Density of Sea Water using Jackett and McDougall 1995 (JAOT 12)
%    polynomial (modified UNESCO polynomial).
%
% INPUT:  (all must have same dimensions)
%   S     = salinity    [psu      (PSS-78)]
%   Theta = potential temperature [degree C (IPTS-68)]
%   P     = pressure    [dbar]
%       (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) )
%
% OUTPUT:
%   dens = density  [kg/m^3]
%
% AUTHOR:  Martin Losch 2002-08-09  (mlosch@mit.edu)

% Jackett and McDougall, 1995, JAOT 12(4), pp. 381-388

% created by mlosch on 2002-08-09
% $Header: /u/gcmpack/MITgcm/utils/matlab/densjmd95.m,v 1.2 2007/02/17 23:49:43 jmc Exp $
% $Name:  $

[1]

Jackett and McDougall, 1995, JAOT 12[4], pp. 381-388

neutralocean.eos.jmd95.rho_s_t(s, t, p)[source]

Fast salinity and potential temperature partial derivatives of JMD95 in-situ density

Parameters:

s, t, pfloat: See rho

Returns:

rho_sfloat: Partial derivative of JMD95 in-situ density with respect to practical salinity s [kg m-3 psu-1]
rho_tfloat: Partial derivative of JMD95 in-situ density with respect to potential temperature t [kg m-3 degC-1]

Notes

This function is derived from rho.

neutralocean.eos.jmd95.rho_p(s, t, p)[source]

Fast pressure derivative of JMD95 in-situ density.

Parameters:

s, t, pfloat: See rho

Returns:

rho_pfloat: Partial derivative of JMD95 in-situ density with respect to pressure p [kg m-3 dbar-1]

Notes

This function is derived from rho.

JMDFWG06

Density of Sea Water using the Jackett et al. (2006) [1] function

Functions:

rho :: computes in-situ density from salinity, potential temperature and: pressure
rho_s_t :: compute the partial derivatives of in-situ density with: respect to salinity and potential temperature
rho_p :: compute the partial derivative of in-situ density with: respect to pressure

Notes: To make Boussinesq versions of these functions, see neutralocean.eos.tools.make_eos_bsq.

To make vectorized versions of these functions, see neutralocean.eos.tools.vectorize_eos.

[1]

Jackett, D. R., McDougall, T. J., Feistel, R., Wright, D. G., & Griffies, S. M. (2006). Algorithms for Density, Potential Temperature, Conservative Temperature, and the Freezing Temperature of Seawater. Journal of Atmospheric and Oceanic Technology, 23(12), 1709–1728. https://doi.org/10.1175/JTECH1946.1

History

Code adapted from MOM5.1 (Griffies et al) originally written in Fortran 2022 January 14 - David Hutchinson - Translated into python 2022 January 14 - Geoff Stanley - code optimization, partial derivatives, check vals

neutralocean.eos.jmdfwg06.rho(s, t, p)[source]

Parameters:

sfloat: Practical salinity [PSS-78]
tfloat: Potential temperature [ITS-90]
pfloat: Pressure [dbar]

Returns:

rhofloat: In-situ density [kg m-3]

neutralocean.eos.jmdfwg06.rho_s_t(s, t, p)[source]

Parameters:

sfloat: Practical salinity [PSS-78]
tfloat: Potential temperature [ITS-90]
pfloat: Pressure [dbar]

Returns:

rho_sfloat: Partial derivative of in-situ density with respect to salinity [kg m-3 psu-1]
rho_tfloat: Partial derivative of in-situ density with respect to temperature [kg m-3 degc-1]

neutralocean.eos.jmdfwg06.rho_p(s, t, p)[source]

Parameters:

sfloat: Practical salinity [PSS-78]
tfloat: Potential temperature [ITS-90]
pfloat: Pressure [dbar]

Returns:

rho_pfloat: Partial derivative of in-situ density with respect to pressure [kg m-3 dbar-1]

TEOS-10 GSW

Density of Sea Water using TEOS-10 Gibbs Sea Water [1] in pure Python

Functions:

specvol :: compute specific volume from Absolute Salinity, Conservative Temperature and: pressure

specvol_first_derivs :: compute first order partial derivatives of specvol

specvol_second_derivs :: compute second order partial derivatives of specvol

specvol_third_derivs :: compute third order partial derivatives of specvol

specvol_s_t :: compute the partial derivatives of specvol with respect to: Absolute Salinity and Conservative Temperature

specvol_p :: compute the partial derivative of specvol with respect to pressure

specvol_s_t_ss_st_tt_sp_tp :: compute various partial derivatives of specvol,: namely with respect to s, t, ss, st, tt, sp, tp.
specvol_s_t_ss_st_tt_sp_tp_sss_sst_stt_ttt_ssp_stp_ttp_spp_tpp :: compute various: partial derivatives of specvol, namely with respect to s, t, ss, st, tt, sp, tp, sss, sst, stt, ttt, ssp, stp, ttp, spp, tpp.

Notes: To make Boussinesq versions of these functions, see neutralocean.eos.tools.make_eos_bsq.

To make vectorized versions of these functions, see neutralocean.eos.tools.vectorize_eos.

All functions here are derived from gsw_specvol, documented below.

# gsw_specvol                            specific volume (75-term equation)
#==========================================================================
# 
# USAGE:  
#  specvol = gsw_specvol(SA,CT,p)
# 
# DESCRIPTION:
#  Calculates specific volume from Absolute Salinity, Conservative 
#  Temperature and pressure, using the computationally-efficient 75-term 
#  polynomial expression for specific volume (Roquet et al., 2015).
#
#  Note that the 75-term equation has been fitted in a restricted range of 
#  parameter space, and is most accurate inside the "oceanographic funnel" 
#  described in McDougall et al. (2003).  The GSW library function 
#  "gsw_infunnel(SA,CT,p)" is available to be used if one wants to test if 
#  some of one's data lies outside this "funnel".  
#
# INPUT:
#  SA  =  Absolute Salinity                                        [ g/kg ]
#  CT  =  Conservative Temperature (ITS-90)                       [ deg C ]
#  p   =  sea pressure                                             [ dbar ]
#         ( i.e. absolute pressure - 10.1325 dbar )
#
#  SA & CT need to have the same dimensions.
#  p may have dimensions 1x1 or Mx1 or 1xN or MxN, where SA & CT are MxN.
#
# OUTPUT:
#  specvol  =  specific volume                                   [ m^3/kg ]
#
# AUTHOR: 
#  Fabien Roquet, Gurvan Madec, Trevor McDougall & Paul Barker
#                                                      [ help@teos-10.org ]
#
# VERSION NUMBER: 3.06 (30th August, 2018)
# 
# REFERENCES:
#  IOC, SCOR and IAPSO, 2010: The international thermodynamic equation of 
#   seawater - 2010: Calculation and use of thermodynamic properties.  
#   Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
#   UNESCO (English), 196 pp.  Available from http://www.TEOS-10.org.  
# 
#  McDougall, T.J., D.R. Jackett, D.G. Wright and R. Feistel, 2003: 
#   Accurate and computationally efficient algorithms for potential 
#   temperature and density of seawater.  J. Atmos. Ocean. Tech., 20, 
#   730-741.
#
#  Roquet, F., G. Madec, T.J. McDougall, P.M. Barker, 2015: Accurate
#   polynomial expressions for the density and specifc volume of seawater
#   using the TEOS-10 standard. Ocean Modelling., 90, pp. 29-43.
#
#  This software is available from http://www.TEOS-10.org
# 
#==========================================================================

[1]

McDougall, T.J. and P.M. Barker, 2011: Getting started with TEOS-10 and the Gibbs Seawater (GSW) Oceanographic Toolbox, 28pp., SCOR/IAPSO WG127, ISBN 978-0-646-55621-5.

neutralocean.eos.gsw.specvol(SA, CT, p)[source]

GSW specific volume.

Parameters:

SAfloat: Absolute Salinity [g/kg]
CTfloat: Conservative Temperature [deg C]
pfloat: sea pressure (i.e. absolute pressure - 10.1325 dbar) [dbar]

Returns:

specvolfloat: Specific volume [m3 kg-1]

neutralocean.eos.gsw.specvol_s_t(SA, CT, p)[source]

Partial derivatives of GSW specific volume with respect to salinity & temperature

Parameters:

SAfloat: Absolute Salinity [g/kg]
CTfloat: Conservative Temperature [deg C]
pfloat: sea pressure (i.e. absolute pressure - 10.1325 dbar) [dbar]

Returns:

sfloat: Partial deriv of specific volume w.r.t. SA [m3 kg-1 / (g/kg)]
tfloat: Partial deriv of specific volume w.r.t. CT [m3 kg-1 / (deg C)]

neutralocean.eos.gsw.specvol_p(SA, CT, p)[source]

Partial derivative of GSW specific volume with respect to pressure

Parameters:

SAfloat: Absolute Salinity [g/kg]
CTfloat: Conservative Temperature [deg C]
pfloat: sea pressure (i.e. absolute pressure - 10.1325 dbar) [dbar]

Returns:

pfloat: Partial deriv of specific volume w.r.t. p [m3 kg-1 / (dbar)]

neutralocean.eos.gsw.specvol_s_t_ss_st_tt_sp_tp(SA, CT, p)[source]

Select partial derivatives of GSW specific volume up to second order

Parameters:

SAfloat: Absolute Salinity [g/kg]
CTfloat: Conservative Temperature [deg C]
pfloat: sea pressure (i.e. absolute pressure - 10.1325 dbar) [dbar]

Returns:

s, t, ss, st, tt, sp, tpfloat: Partial derivatives of specific volume w.r.t. SA, CT, SA*SA, SA*CT, CT*CT, SA*p, CT*p.

neutralocean.eos.gsw.specvol_s_t_ss_st_tt_sp_tp_sss_sst_stt_ttt_ssp_stp_ttp_spp_tpp(SA, CT, p)[source]

Select partial derivatives of GSW specific volume up to third order

Parameters:

SAfloat: Absolute Salinity [g/kg]
CTfloat: Conservative Temperature [deg C]
pfloat: sea pressure (i.e. absolute pressure - 10.1325 dbar) [dbar]

Returns:

s, t, ss, st, tt, sp, tp, sss, sst, stt, ttt, ssp, stp, ttp, spp, tppfloat: Partial derivatives of specific volume.

(Vertical) Interpolation using Piecewise Polynomials (PP)

Piecewise polynomial interpolation in one dimension.

This package separates 1D interpolation into two steps:

Compute coefficients for a piecewise polynomial interpolant,
Evaluate the interpolant.

If the interpolant will be evaluated many times, then it is fastest to do Step 1 once, then repeat Step 2 as many times as needed. To do this, build the coefficients using *coeffs* methods, such as linear_coeffs and pchip_coeffs. Then evaluate the interpolants using ppval methods.

If the interpolant will only be evaluated once or a few times, then it will be faster to combine Steps 1 and 2, building the polynomial (not a piecewise polynomial) for just the interval containing the evaluation site, using only as much data around the evaluation site as is required. To do this, use the *interp* methods, such as linear_interp and pchip_interp.

The core numerical methods handle just one interpolation problem; that is, the evaluation site (x) is a scalar and the input independent and dependent data (X and Y) are 1D arrays. These functions are named with "_1".

Functions that lack this "_1" are “universal”, accepting an N-dimensional array for x and (N+1)-dimensional arrays for X, and Y, so long as they are mutually broadcastable from the basic 1D problem, thus allowing multiple 1D interpolation problems to be solved with one function call. Effectively, a low-level for loop is used to wrap the single 1D problem. Operationally, this is achieved using numba’s guvectorize decorator.

If the interpolant must be evaluated many times but the dataset is large and memory is limited, then the best strategy is to compute the piecewise polynomial coefficients for a single 1D interpolation problem (Step 1), then evaluate it (Step 2) as many times as needed. This can be achieved by a numba.njit accelerated for loop over each 1D problem, using the "_1" function variants. For example, neutralocean.surface._vertsolve does this to solve many individual 1D non-linear root finding problems.

If performing a single 1D interpolation problem and all input data is finite (no NaN’s), the "_nonan" variants provide a small speed advantage.

If there are two dependent variables that share the same dependent variable, the "_two" variants offer a speed advantage over calling the basic functions twice, since the "_two" variants do the work of locating where the evaluation site x lies within the dependent data X only once. This could easily be extended to handle more dependent variables, of course.

The make_pp factory function will help select the desired variant of the *coeffs* and *interp* functions.

neutralocean.ppinterp.tools.make_pp(interpolant='linear', kind='u', out='coeffs', nans=True, num_dep_vars=1)[source]: Make function for piecewise polynomial interpolation.

neutralocean.ppinterp.ppval(x, X, Yppc, d, y)[source]

Evaluate a given Piecewise Polynomial

Parameters:

xndarray

Sites at which to evaluate the independent variable.

Must be broadcastable to the shape of X or Y less their final dimension.

Xndarray

Independent variable. Must be broadcastable to the shape of Y.

X must be monotonically increasing in its last dimension. That is, X[*i,:] must be monotonically increasing for any i tuple indexing all but the last dimension.

Yndarray

Dependent variable. Must be broadcastable to the shape of X.

Yppcndarray

Piecewise Polynomial Coefficients for the interpolant. The shape of Yppc less its final dimension must be broadcastable to the shape of X and/or Y.

Returns:

yndarray: Value of the piecewise polynomial interpolant (or its d’th derivative) at x.

Notes

If X and Y have NaN’s, only the first contiguous block of non-NaN data between both X and Y is used.

A binary search is performed to find the indices of X that x lies between. As such, nonsense results will arise if X is not sorted along its last dimension.

neutralocean.ppinterp.ppval_1(x, X, Yppc, d=0)[source]

Evaluate a single Piecewise Polynomial (PP).

Parameters:

xfloat: Evaluation site
Xndarray, 1d: Independent variable. Must be monotonically increasing.
Yppcndarray, 2d: Piecewise Polynomial Coefficients. First dimension must match len(X).
dint, Default 0: Number of derivatives to take. If 0, simply evaluate the PP.

Returns:

yfloat: The value of the PP (or its d’th derivative) at X = x.

Notes

If X and Y have NaN’s, only the first contiguous block of non-NaN data between both X and Y is used.

neutralocean.ppinterp.ppval_two(x, X, Yppc, Zppc, d, y, z)[source]

Evaluate two given Piecewise Polynomials

The inputs and outputs are as for ppval, but another input, Zppc, gives the second piecewise polynomial, for which a second output z is given.

Zppc must have the same NaN-structure as Yppc.

Notes

Calling ppval_two is faster than calling ppval twice, since the former calls np.searchsorted half as many times.

neutralocean.ppinterp.ppval_1_two(x, X, Yppc, Zppc, d=0)[source]

Evaluate two piecewise polynomials.

As ppval_1 but a second input, Zppc, provides a second set of Piecewise Polynomial Coefficients. Correspondingly, a second output, z, is returned.

Zppc must have the same NaN-structure as Yppc.

neutralocean.ppinterp.linear.linear_interp(x, X, Y, d, y)[source]

neutralocean.ppinterp.linear.linear_interp_1(x, X, Y, d=0)[source]: Build and evaluate a piecewise polynomial, building only what’s needed.

neutralocean.ppinterp.linear.linear_coeffs(X, Y)[source]

Piecewise Polynomial Coefficients for a linear interpolant

Parameters:

Xndarray

Independent variable. Must be same dimensions as Y, with any dimension possibly a singleton: that is, must be broadcastable to the size of Y.

X must be monotonically increasing in its last dimension. That is, X[*i,:] should be monotonically increasing for any i tuple indexing all but the last dimension.

Yndarray

Dependent variable. Must be same dimensions as X, with any dimension possibly a singleton: that is, must be broadcastable to the size of X.

Returns:

Yppcndarray: Coefficients for a linear interpolant of Y to X.

Notes

If X and Y have NaN’s, only the first contiguous block of non-NaN data between both X and Y is used.

Evaluate the piecewise polynomial at x as

>>> y = ppval(x, X, Yppc)

neutralocean.ppinterp.linear.linear_coeffs_1(X, Y)[source]

Coefficients for a single linear interpolant

Inputs and outputs are as for linear_coeffs, but X and Y both 1D arrays of length n, and so Yppc is a 2D array of size (n, 2).

neutralocean.ppinterp.linear.linear_interp_two(x, X, Y, Z, d, y, z)[source]

neutralocean.ppinterp.linear.linear_interp_1_two(x, X, Y, Z, d=0)[source]

neutralocean.ppinterp.pchip.pchip_interp(x, X, Y, d=0)[source]

neutralocean.ppinterp.pchip.pchip_interp_1(x, X, Y, d=0)[source]

neutralocean.ppinterp.pchip.pchip_coeffs(X, Y)[source]

Piecewise Polynomial Coefficients for a PCHIP (Piecewise Cubic Hermite Interpolating Polynomial)

Parameters:

Xndarray

Independent variable. Must be same dimensions as Y, with any dimension possibly a singleton: that is, must be broadcastable to the size of Y.

X must be monotonically increasing in its last dimension. That is, X[*i,:] should be monotonically increasing for any i tuple indexing all but the last dimension.

Yndarray

Dependent variable. Must be same dimensions as X, with any dimension possibly a singleton: that is, must be broadcastable to the size of X.

Returns:

Yppcndarray: Coefficients for a PCHIP interpolant of Y to X.

Notes

If X and Y have NaN’s, only the first contiguous block of non-NaN data between both X and Y is used.

Evaluate the piecewise polynomial at x as

>>> y = ppval(x, X, Yppc)

neutralocean.ppinterp.pchip.pchip_coeffs_1(X, Y)[source]

Coefficients for a single PCHIP when the data may contain NaNs

Inputs and outputs are as for pchip_coeffs, but X and Y both 1D arrays of length n, and so Yppc is a 2D array of size (n, 4).

neutralocean.ppinterp.pchip.pchip_coeffs_1_nonan(X, Y)[source]

Coefficients for a single PCHIP when the data has no NaNs

Inputs and outputs are as for pchip_coeffs, but X and Y both 1D arrays of length n, and so Yppc is a 2D array of size (n, 4).

neutralocean.ppinterp.pchip.pchip_interp_two(x, X, Y, Z, d=0)[source]

neutralocean.ppinterp.pchip.pchip_interp_1_two(x, X, Y, Z, d=0)[source]

Synthetic Ocean

neutralocean.synthocean.synthocean(shape, pbot=4000.0, SSSSa=0.3, zonally_uniform=False, wrap=(False, False))[source]

Synthetic idealization of the Pacific and Southern Ocean with tuneable parameters

Parameters:

shapetuple of int: A three element tuple giving the dimensions of the output. The elements specify the number of points in longitude, latitude, and depth space, respectively.
pbotfloat, Default 4000.0: Value of P at the bottommost data point, i.e. P[-1]
SSSSafloat, Default 0.3: Southern Sea Surface Salinity anomaly: the sea surface salinity is increased by SSSS in the southernmost grid cells, and by 0 in the northernmost grid cells, and linearly in between. The degree to which the southern casts are statically unstable can be controlled by this parameter. Increase to make more unstable.
zonally_uniformbool, Default False: If True, the synthetic ocean is uniform in longitude, and the neutral helicity is zero. If False, some zonal structure is created which results in non-zero neutral helicity.
wraptuple of bool, Default (False, False): Specify periodicity of the lateral dimensions. When wrap[0] is False, S[0, :, :] = T[0, :, :] = nan. When wrap[1] is False, S[:, 0, :] = T[:, 0, :] = nan.

Returns:

S, Tndarray: The practical / Absolute Salinity and potential / Conservative Temperature as a 3D array.
Pndarray: The pressure / depth as a 1D array.
gdict: The geometry of the grid, including latitude, longitude, cell widths and areas. See code for details.