Public Function Reference

Approximately Neutral Surfaces

Potential Density Surface

neutralocean.surface.potential_surf(S, T, P, **kw)

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.

eosfunction, Default neutralocean.eos.gsw.specvol

Function taking three inputs corresponding to (S, T, P), and outputting the in-situ density or specific volume. Should be @numba.njit decorated and need not be vectorized.

eos_s_tfunction, neutralocean.eos.gsw.specvol_s_t

Function taking three inputs corresponding to (S, T, P), and outputting a tuple containing the partial derivatives of the equation of state with respect to S and T. Should be @numba.njit decorated and need not be vectorized. Note: This is only used when diags is True.

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.anomaly_surf(S, T, P, **kw)

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, 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_surf(S, T, P, grid, pin_cast, p_init, **kw)

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

pin_castint or tuple of int

Index for cast where surface is kept at fixed pressure or depth.

p_initfloat or ndarray or xarray.DataArray

If array, pressure or depth of the initial approximately neutral surface. Must be the same shape as S less its vertical dimension.

If float, pressure or depth at pin_cast. The initial surface is generated by iteratively wetting the perimeter of the region, beginning with the reference cast, pin_cast.

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, interp, eos, eos_s_t, diags, output, TOL_P_SOLVER

See potential_surf Note: eos_s_t is required.

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 mld with p_ml passed as keyword arguments, enabling control over the parameters in that function. See mld 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.

Grids

Rectilinear

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

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)

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)

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)

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)

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)

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)

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=-1, weights=None)

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

Parameters:

edges2D array of shape (2, E), or length 2 tuple of 1D arrays of length E

The edges in a graph whose nodes are labelled 0, ..., N-1. Denoting a = edges[0] and b = edges[1], - a and b are 1D arrays of type int - 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.

Nint, Default -1

Number of nodes in the graph. If -1, 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)

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 shape (2, E), or length 2 tuple of 1D arrays of length E

The edges in a graph whose nodes are labelled 0, ..., N-1. Denoting a = edges[0] and b = edges[1], - a and b are 1D arrays of type int - 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

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=CPUDispatcher(<function specvol_s_t>), **kw)

Calculate epsilon neutrality errors on an approximately neutral surface

Parameters:

s, t, p, eos_s_t

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=CPUDispatcher(<function specvol_s_t>), **kw)

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, Default neutralocean.eos.gsw.specvol_s_t

Function taking three inputs corresponding to (s, t, p), and outputting a tuple containing the partial derivatives of the equation of state with respect to s and t.

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

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=CPUDispatcher(<function specvol>), **kw)

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.

eosfunction, Default neutralocean.eos.gsw.specvol

Function taking three inputs corresponding to (S, T, P), and outputting the in-situ density or specific volume.

neutralocean.traj.neutral_trajectory(S, T, P, p0, vert_dim=-1, tol_p=0.0001, interp='linear', eos=CPUDispatcher(<function specvol>), **kw)

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

See ntp_bottle_to_cast

Veronis Density

neutralocean.label.veronis(p1, S, T, P, p0=None, p_ref=0.0, b=1.0, dp=1.0, interp='linear', eos=None, eos_s_t=None, **kw)

The surface density plus the integrated vertical gradient of Locally Referenced Potential Density

Determines the Veronis density [1] [2] at vertical position p1 on a cast with hydrographic properties (S, T, P). The Veronis density is the potential density (referenced to p_ref) evaluated at p0 on the cast, plus b times the integral (obtained numerically by trapezoidal integration with step sizes of dP or less) of the vertical (d/dP) derivative of Locally Referenced Potential Density (LRPD) from P = p0 to P = p1. The vertical (d/dP) derivative of LRPD is rho_S dS/dP + rho_T dT/dP where rho_S and rho_T are the partial derivatives of density with respect to S and T, and dS/dP and dT/dP are the derivatives of S and T with respect to P in the water column. If p0 or p1 are outside the range of P, NaN is returned.

In common oceanographic mathematical notation, the above is written as (see Equation (A2) in [2])

ρᵥ(p₁) = ρ(S(p₀), Θ(p₀), p_ref) + b ∫ [f dS/dp + g dΘ/dp] dp

where the integral goes from p₀ to p₁, and where f(p) = (∂ρ/∂S)(S(p),Θ(p),p) and g(p) = (∂ρ/∂Θ)(S(p),Θ(p),p) are the partial derivatives of the equation of state ρ. Note the Conservative (or potential) Temperature Θ is denoted T in the code.

Parameters:

p1float

Pressure or depth at which the Veronis density is evaluated

S, T, P1D ndarray 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:

dfloat

Veronis density

Other Parameters:

p0float, Default P[0]

Pressure or depth at which the potential density is evaluated

p_reffloat, Default 0.0

reference pressure or depth for potential density

bfloat, Default 1.0

Factor by which to multiply the integral.

dpfloat, Default 1.0

Maximum interval of pressure or depth in trapezoidal numerical integration

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.

eosfunction, Default neutralocean.eos.gsw.specvol

Function taking three inputs corresponding to (S, T, P), and outputting the in-situ density or specific volume.

eos_s_tfunction, neutralocean.eos.gsw.specvol_s_t

Function taking three inputs corresponding to (S, T, P), and outputting a tuple containing the partial derivatives of the equation of state with respect to S and T.

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

Notes

The Veronis density at a particular pressure / depth on a reference cast can serve as a neutral density label for an approximately neutral surface (such as an omega surface) intersecting that point, and for temporal neutral trajectories through that point, and for other approximately neutral surfaces (such as omega surfaces) intersecting that temporal neutral trajectory. The Veronis density serves as a useful neutral density label in this way because it perfectly represents the stratification in a single vertical cast, but its horizontal variations are not meaningful.

However, the neutral density label obtained by this Veronis density is NOT the same as a value of the Jackett and McDougall (1997) [3] Neutral Density variable. These two values are different even if you were to provide this function with the same cast that Jackett and McDougall (1997) used to initially label their Neutral Density variable, namely the cast at 188 deg E, 4 deg S, from the Levitus (1982) ocean atlas. Some difference would remain, because of differences in numerics, and because of a subsequent smoothing step in the Jackett and McDougall (1997) algorithm. This function merely allows one to label an approximately neutral surface with a density value that is INTERNALLY consistent within the dataset where one’s surface lives. This function is NOT to compare density values against those from any other dataset, such as 1997 Neutral Density.

See Appendix A of [2] for further information.

Examples

>>> S = np.linspace(34, 35.5, 20)
>>> T = np.linspace(18, 0, 20)
>>> P = np.linspace(0, 4000, 20)
>>> veronis(2000, S, T, P)
0.0009737571001312298

Calculate the Veronis density at 2000 dbar on a water column of linearly varying salinity and potential temperature.

[1]

Veronis, G. (1972). On properties of seawater defined by temperature, salinity, and pressure. Journal of Marine Research, 30(2), 227.

[2] (1,2,3)

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.

[3]

Jackett, D. R., & McDougall, T. J. (1997). A neutral density variable for the world’s oceans. Journal of Physical Oceanography, 27 (2), 237–263.

Mixed Layer

neutralocean.mixed_layer.mld(S, T, P, eos=CPUDispatcher(<function specvol>), pot_dens_diff=0.03, ref_p=100.0, bottle_index=1, interp='linear', vert_dim=-1, **kw)

Calculate the mixed layer pressure or depth

The mixed layer depth (mld) 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).

eosfunction, Default neutralocean.eos.gsw.specvol

Function taking three inputs corresponding to (S, T, P), and outputting a the in-situ density or specific volume.

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:

mldndarray

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

Static Stability

neutralocean.stability.count_unstable(S, T, P, **kw)

Count number of statically unstable bottles.

S, T, P :

As in stabilize_ST

interp_two : function or None, Default None

When None, static stability is calculated as the difference of potential density between pairs of adjacent bottles on each cast, referenced to the average pressure between these two bottles.

When a function, it is used to evaluate the first derivative of an interpolant for each of S and T as a function of P, at each bottle’s P value. The vertical gradient of Locally Referenced Potential Density is then calculated by combining these derivatives with the partial derivatives of the equation of state with respect to S and T at the local P.

eosfunction, Default neutralocean.eos.gsw.specvol

Equation of State; used when interp_two is None.

eos_s_t : function, Default neutralocean.eos.gsw.specvol_s_t

Partial derivatives of the Equation of State (eos above) with respect to S and T; used when interp_two is given.

neutralocean.stability.stabilize_ST(S, T, P, **kw)

Mutate S, T to ensure the vertical gradient of Locally Referenced Potential Density (LRPD) is greater than a given threshold everywhere.

Parameters:

S, Tndarray or xarray.DataArray

practical / Absolute salinity and potential / Conservative temperature.

Pndarray or xarray.DataArray

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

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

eosfunction, Default neutralocean.eos.gsw.rho

Equation of state for density (not specific volume). Takes three inputs corresponding to (S, T, P), and outputs density. Should be @numba.njit decorated and need not be vectorized.

min_dLRPDdpfloat or 1D array of float, Default 1e-6

Minimum vertical gradient of LRPD for the mutated S and T.

If an array, this is the minimum vertical gradient of LRPD between adjacent pairs of bottles in each cast; hence, this should have length one less than the length of the vertical dimension of S and T.

weightfloat, Default 10.0

Multiplicative weighting factor applied to S perturbations, but not to T perturbations. When weight > 1, T perturbations are favoured. When weight < 1, S perturbations are favoured.

vert_dimint or str

Specifies which dimension of S, T (and P if more than 1D) 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.

tolfloat, Default 1e-10

Tolerance for (weighted) S and (unweighted) T perturbations; passed to scipy’s minimize.

methodstr, Default “SLSQP”

Minimization method, passed to scipy’s minimize.

optionsdict, Default {“maxiter”1000}

Options passed to the options argument of minimize.

verbosebool, Default True

Whether to print information any time a cast is mutated.

Notes

This function is similar to the method of Jackett and McDougall (1995) [1]. However, it uses a fixed weighting between S and T, and it solves the optimization problem with non-linear constraints. It also uses a slightly different numerical implementation for the vertical gradient of LRPD.

[1]

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

Equation Of State

Tools

neutralocean.eos.tools.load_eos(eos, derivs='', grav=None, rho_c=None)

Load EOS function from library.

Parameters:

eosstr

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

derivsstr, Default “”

String specifying which partial derivatives of the EOS are desired. Only used when eos is a string. The actual function loaded is named eos + derivs. For example, “” loads the EOS itself, “_p” will load the partial derivative with respect to p, “_s_t” will load the function whose two outputs are the s and t partial derivatives, respectively.

grav, rho_cfloat, Default None

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:

fn: function

Equation of State function accepting three arguments: (salinity, temperature, pressure) when grav or rho_c is None, or (salinity, temperature, depth) otherwise.

Notes

[1] (1,2)

Roquet, F., G. Madec, Trevor J. McDougall, and Paul M. Barker. “Accurate Polynomial Expressions for the Density and Specific Volume of Seawater Using the TEOS-10 Standard.” Ocean Modelling 90 (June 2015): 29-43. https://doi.org/10.1016/j.ocemod.2015.04.002.

[2] (1,2)

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.

[3]

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

[4]

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_bsq(fn, grav, rho_c)

Make a Boussinesq version of a given equation of state (or its partial derivative(s))

Parameters:

fnfunction

Function with (salinity, temperature, pressure, pfac) as inputs. Typically this is the equation of state, returning the density or specific volume. However, it can also be a function for partial derivative(s) of the equation of state with respect to salinity, temperature, or pressure. The 4th argument, pfac, pre-multiplies pressure before the main calculation, and also post-multiplies the output as many times as there are pressure partial derivatives.

gravfloat

Gravitational acceleration [m s-2]

rho_cfloat

Boussinesq reference density [kg m-3]

Returns:

fn_bsqfunction

Boussinesq version of fn. The inputs to fn_bsq are (salinity, temperature, depth).

neutralocean.eos.tools.vectorize_eos(eos)

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

Parameters:

eosfunction

Any function taking three scalar inputs (additional optional arguments are allowed but only take their default value) 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.

(Vertical) Interpolation using Piecewise Polynomials (PP)

Piecewise polynomial interpolation in one dimension.

This package separates 1D interpolation into two steps:

1. Compute coefficients for a piecewise polynomial interpolant,

2. 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)

Make function for piecewise polynomial interpolation.

Parameters:

interpolantstr, Default “linear”

Name of the interpolant to use. Currently, either “linear” or “pchip”.

kindstr, Default “u”

Either “1” (single; make a function that performs one interpolation) or “u” (universal; make a function performs a low level loop over multiple interpolation problems, using numba’s guvectorize decorator).

If “1”, the returned function f operates on independent and dependent data, X and Y, that are 1D arrays, and (if out="interp") the evaluation site x is a scalar.

If “u”, the returned function f operates on independent and dependent data, X and Y, may be N+1 dimensional arrays that are mutually broadcastable, and (if out="interp") the evaluation site x may be a N dimensional array that is broadcastable to the first N dimensions of X or Y.

outstr, Default “coeffs”

Specifies the type of output returned by f.

If "coeffs", then f returns the piecewise polynomial coefficients, and f does not take in the evaluation sites x.

If "interp", then f evaluates the interpolant(s) at the evaluation site x.

nansbool, Default True

If True, the independent and dependent data, X and Y, may have NaNs.

If False, X and Y, must have no NaNs. In this case, a slightly faster function is created that skips the step of checking for the first contiguous block of non-NaN data (performed by valid_range_1_two).

num_dep_varsint, Default 1

If 2, f takes two dependent variables (Y and Z) that are both evaluated at the evaluation site x. This is slightly faster than calling an interpolation function twice, because the location of the evaluation site x within the dependent data X need only be calculated once. Note if out = "coeffs" then num_dep_vars = 1 is required. In that case, create two pairs of piecewise polynomail coefficients Yppc and Zppc then evaluate them at x using ppval_two.

Returns:

ffunction

Interpolating function. The possible inputs to f are

• x the evaluation site(s), (only if out = "interp")

• X the independent data

• Y the dependent data

• Z additional dependent data (only if num_dep_vars = 2)

The outputs of f are

• Yppc the piecewise polynomial coefficients (if out = "coeffs") OR

• y the interpolant evaluated at x (if `out = “interp”), OR

• (y, z) the interpolants evaluated at x (if out = "interp" and `num_dep_vars=2).

neutralocean.ppinterp.ppval(x, X, Yppc, d, y)

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. The size of the last dimension of Yppc is the order of the polynomials.

dint

Order of the derivative to evaluate. When 0, just evaluate the polynomial.

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)

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). Second dimension is the order of the polynomials.

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)

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)

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.

Data

neutralocean.data.load_OCCA()

Load grid, salinity, and potential temperature for OCCA 2004-2006 average

Returns:

gdict

Grid information

Sxarray.DataArray

Salinity

Txarray.DataArray

Potential temperature

neutralocean.data.synthocean(shape, pbot=4000.0, SSSSa=0.3, zonally_uniform=False, wrap=(False, False))

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.