Getting Started
Installation
Simply execute either
(.venv) $ pip install neutralocean
if you use pip, or
(.venv) $ conda install -c conda-forge neutralocean
if you use conda.
Source Code
The open-source neutralocean code is available at https://github.com/geoffstanley/neutralocean
Test Example
The example script run_4casts.py (copied below and on github) creates synthetic salt, temperature, and pressure data for 4 vertical casts, specifies which pairs of these 4 casts are adjacent (see Grids), and then calculates a potential density surface, specific volume anomaly surface, and omega surface. Epsilon neutrality errors are calculated for each link on the final omega surface.
Try running this script yourself, e.g. (changing path as needed)
python /path/to/neutralocean/examples/run_4casts.py
If that runs, your neutralocean package is working correctly.
Input
# Show basic use of neutralocean on an unstructured grid using artificial data.
#
# For more advanced usage of neutralocean, see the run_OCCA.py example.
# In[Imports]
import numpy as np
import neutralocean as no
# In[Create an ocean of 4 water columns]
# Make up simple Salt, Temperature, and Pressure data for 4 casts
nc = 4 # Number of casts (water columns)
nk = 20 # Number of grid points per cast
pbot = 4000.0 # Pressure at bottom grid point
P = np.linspace(0, 1, nk) ** 3 * pbot # pressure increasing cubicly going down
S = np.linspace(34, 36, nk).reshape((1, -1)) # saltier going down
S = S + np.linspace(0, 0.9, nc).reshape((-1, 1)) # some lateral structure
T = np.linspace(14, -2, nk).reshape((1, -1)) # warmer going down
T = T + np.linspace(0, 6, nc).reshape((-1, 1)) # some lateral structure
# Arrange the 4 casts (labelled 0, 1, 2, 3) with connections as follows:
# 0
# | \
# 1--2--3
# That is, cast 0 is connected to casts 1 and 2; cast 1 is connected to casts 0
# and 2; cast 2 is conncted to all casts; cast 3 is connected only to cast 2.
a = np.array([0, 0, 1, 2]) # a[i] and b[i] are a pair of adjacent casts
b = np.array([1, 2, 2, 3])
edges = (a, b)
# Invent distances `dist` between pairs of casts, roughly 100km
dist = np.array([1, 1.4, 1, 1]) * 1e5 # Units: [m]
# Invent distances `distperp` of the interfaces between pairs of casts.
# The product of `dist` and `distperp` gives an area associated to the region
# between casts, which is where the ϵ neutrality errors live. We seek to minimize
# these ϵ neutrality errors, weighted by these areas.
distperp = np.array([1, 1, 1, 1]) * 1e5 # Units: [m]
# Package the grid information into a dict, for neutralocean.
grid = {"edges": edges, "dist": dist, "distperp": distperp}
"""
# Note, more complex examples may have the water column adjacency specified by
# a graph structure, represented as a (sparse) matrix. For example (with a
# dense matrix here given the smallness of this example), suppose we have two
# graphs, each with the same sparsity structure, one storing the distances
# between adjacent water columns and the other storing the length of the
# interfaces between adjacent water columns.
graph_dist = np.zeros((nc, nc))
graph_dist[0, :] = [0.0, 1.0, 1.4, 0.0]
graph_dist[1, :] = [1.0, 0.0, 1.0, 0.0]
graph_dist[2, :] = [1.4, 1.0, 0.0, 2.0]
graph_dist[3, :] = [0.0, 0.0, 2.0, 0.0]
graph_dist *= 1e5
graph_distperp = np.sign(graph_dist) * 1e5
# To convert these graphs into the format for the `grid` argument to
# neutralocean functions, simply call:
from neutralocean.grid.graph import build_grid
grid = build_grid({"dist": graph_dist, "distperp": graph_distperp})
"""
# In[Approx Neutral Surfaces]
# Here, ν = 1/ρ is the TEOS-10 specific volume,
# S is Absolute Salinity,
# Θ is Conservative Temperature,
# S is Absolute Salinity,
# p is pressure
# Potential specific volume surface, with given reference pressure and given isovalue.
# This finds the surface satisfying
# ν(S, Θ, 0 dbar) = (1/1027.5) m3 / kg
s, t, p, d = no.potential_surf(
S,
T,
P,
grid=grid,
ref=0.0,
isoval=1 / 1027.5,
)
print(
f" ** The potential specific volume surface (referenced to {d['ref']}dbar)"
f" with isovalue = {d['isoval']} m3 kg-1"
f" has root-mean-square ϵ neutrality error {d['e_RMS']} m2 kg-1."
)
# In-situ specific volume anomaly, with given reference S and Θ values and given isovalue.
# This finds the surface satisfying
# ν(S, Θ, p) - ν(34.5 g/kg, 4.0°C, p) = 0 m3 / kg
s0, t0 = 34.5, 4.0
s, t, p, d = no.anomaly_surf(
S,
T,
P,
grid=grid,
ref=(s0, t0),
isoval=0.0,
)
print(
f" ** The in-situ specific volume anomaly surface (referenced to {d['ref']})"
f" with isovalue = {d['isoval']} m3 kg-1"
f" has root-mean-square ϵ neutrality error {d['e_RMS']} m2 kg-1."
)
# omega-surface, pinned to be 1500dbar on cast 0, initialized by iteratively
# making Neutral Tangent Plane links from cast 0.
s, t, p, d = no.omega_surf(S, T, P, grid, pin_cast=0, p_init=1500.0)
print(
f" ** The omega-surface"
f" initialized from a potential density surface (referenced to 1500 dbar)"
f" intersecting the cast labelled '0' at pressure 1500 bar"
f" has root-mean-square ϵ neutrality error {d['e_RMS'][-1]} m2 kg-1."
)
# Calculate ϵ neutrality errors on the latest surface, between all pairs of adjacent water columns
e = no.ntp_epsilon_errors(s, t, p, grid)
print("The ϵ neutrality errors on the ω-surface are as follows:")
for i in range(len(a)):
print(f" From cast {a[i]} to cast {b[i]}, ϵ = {e[i]} m2 kg-1")
print(
"Note that the connection between casts 2 and 3 has virtually 0 neutrality "
"error. This is because cast 3 is ONLY connected to cast 2, so this link "
"can be along the (discrete) neutral tangent plane joining cast 2 and 3. "
"The ω-surface finds this."
)
Output
potential done | 4 wet casts | RMS(ϵ) = 2.95366114e-13 | 2.762 sec
** The potential specific volume surface (referenced to 0.0dbar) with isovalue = 0.0009732360097323601 m3 kg-1 has root-mean-square ϵ neutrality error 2.953661135576716e-13 m2 kg-1.
anomaly done | 4 wet casts | RMS(ϵ) = 1.21998943e-13 | 2.327 sec
** The in-situ specific volume anomaly surface (referenced to (34.5, 4.0)) with isovalue = 0.0 m3 kg-1 has root-mean-square ϵ neutrality error 1.2199894286212603e-13 m2 kg-1.
iter | MAV(ϕ) | RMS(Δp) | # wet casts (# new) | RMS(ϵ) | time (s)
0 | | 4 | 7.05305721e-14 | 0.000
1 | 7.62146233e-09 | 1.77640515e+01 | 4 ( 0) | 8.48888373e-15 | 3.417
2 | 6.02796769e-10 | 1.44982488e+00 | 4 ( 0) | 6.12093753e-15 | 0.001
3 | 4.57931777e-11 | 1.11692359e-01 | 4 ( 0) | 6.10611565e-15 | 0.001
4 | 3.38861795e-12 | 8.31849303e-03 | 4 ( 0) | 6.10621916e-15 | 0.001
5 | 2.47457692e-13 | 6.09321842e-04 | 4 ( 0) | 6.10623327e-15 | 0.001
6 | 1.79561251e-14 | 4.42768237e-05 | 4 ( 0) | 6.10623432e-15 | 0.178
** The omega-surface initialized from a potential density surface (referenced to 1500 dbar) intersecting the cast labelled '0' at pressure 1500 bar has root-mean-square ϵ neutrality error 6.1062343244329895e-15 m2 kg-1.
The ϵ neutrality errors on the ω-surface are as follows:
From cast 0 to cast 1, ϵ = -6.946406152330754e-15 m^2 kg^-1
From cast 0 to cast 2, ϵ = 6.9464177373746585e-15 m^2 kg^-1
From cast 1 to cast 2, ϵ = -6.946399353678772e-15 m^2 kg^-1
From cast 2 to cast 3, ϵ = 2.4532753952971413e-20 m^2 kg^-1
Note that the connection between casts 2 and 3 has virtually 0 neutrality error. This is because cast 3 is ONLY connected to cast 2, so this link can be along the (discrete) neutral tangent plane joining cast 2 and 3. The ω-surface finds this.