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/neutralocean/examples/run_4casts.py
If that runs, your neutralocean package is working correctly.
Input
# In[Imports]
import numpy as np
from neutralocean.surface import potential_surf, anomaly_surf, omega_surf
from neutralocean.eos import make_eos_s_t
from neutralocean.ntp import ntp_epsilon_errors
# 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 = potential_surf(
S,
T,
P,
grid=grid,
eos="gsw",
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 = anomaly_surf(
S,
T,
P,
grid=grid,
eos="gsw",
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, initialized by a potential density surface, pinning the surface
# to be 1500dbar on cast 0.
s, t, p, d = omega_surf(S, T, P, grid, pin_cast=0, pin_p=1500.0, eos="gsw")
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
eos_s_t = make_eos_s_t("gsw")
e = ntp_epsilon_errors(s, t, p, grid, eos_s_t)
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.