Hi there good people,
I’m pleased to release neutralocean to the world. It’s a Python package to calculate neutral surfaces in the ocean. Currently, the main job is to compute ω-surfaces, which solve a global optimization problem to reduce the neutrality error by typically 1-2 orders of magnitude compared to potential density surfaces.
So if you ever calculate isopycnals using potential density, specific volume anomaly, Jackett & McDougall’s 1997 γⁿ neutral density, or others, bear in mind that these are often not terribly accurate in terms of being aligned with the directions in which the ocean mixes rapidly (“isopycnal” or “epineutral” mixing), and being orthogonal to the directions in which the ocean mixes very weakly (“diapycnal” or “dianeutral” mixing). Consider using ω-surfaces for such analyses.
Also, if you think you need a 3D density variables (like the above), sometimes these analysis can actually be done using 2D surfaces (like ω-surfaces), perhaps 1 of them, perhaps 2 of them very slightly apart (to get layer thickness and potential vorticity), or perhaps a whole set of them.
And feel free to reach out to me with questions.
Geoff
P.S. There’s also a MATLAB version , which can calculate geostrophic streamfunctions and some other advanced neutral surfaces, but only works on rectilinear grids.