/usr/lib/python3/dist-packages/gsw/stability.py is in python3-gsw 3.2.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 | """
Vertical stability functions.
These work with ndarrays of profiles; use the `axis` keyword
argument to specify the axis along which pressure varies.
For example, the default, following the Matlab versions, is
`axis=0`, meaning the pressure varies along the first dimension.
Use `axis=-1` if pressure varies along the last dimension--that
is, along a row, as the column index increases, in the 2-D case.
Docstrings will be added later, either manually or via
an automated mechanism.
"""
import numpy as np
from ._utilities import match_args_return, axis_slicer
from ._gsw_ufuncs import grav, specvol_alpha_beta
__all__ = ['Nsquared',
'Turner_Rsubrho',
'IPV_vs_fNsquared_ratio',
]
# In the following, axis=0 matches the Matlab behavior.
@match_args_return
def Nsquared(SA, CT, p, lat=None, axis=0):
"""
Calculate the square of the buoyancy frequency.
Parameters
----------
SA : array-like
Absolute Salinity, g/kg
CT : array-like
Conservative Temperature (ITS-90), degrees C
p : array-like
Sea pressure (absolute pressure minus 10.1325 dbar), dbar
lat : array-like, 1-D, optional
Latitude, degrees.
axis : int, optional
The dimension along which pressure increases.
Returns
-------
N2 : array
Buoyancy frequency-squared at pressure midpoints, 1/s.
The shape along the pressure axis dimension is one
less than that of the inputs.
p_mid : array
Pressure at midpoints of p, dbar.
The array shape matches N2.
"""
if lat is not None:
if np.any((lat < -90) | (lat > 90)):
raise ValueError('lat is out of range')
SA, CT, p, lat = np.broadcast_arrays(SA, CT, p, lat)
g = grav(lat, p)
else:
SA, CT, p = np.broadcast_arrays(SA, CT, p)
g = 9.7963 # (Griffies, 2004)
db_to_pa = 1e4
shallow = axis_slicer(SA.ndim, slice(-1), axis)
deep = axis_slicer(SA.ndim, slice(1, None), axis)
if lat is not None:
g_local = 0.5 * (g[shallow] + g[deep])
else:
g_local = g
dSA = SA[deep] - SA[shallow]
dCT = CT[deep] - CT[shallow]
dp = p[deep] - p[shallow]
SA_mid = 0.5 * (SA[shallow] + SA[deep])
CT_mid = 0.5 * (CT[shallow] + CT[deep])
p_mid = 0.5 * (p[shallow] + p[deep])
specvol_mid, alpha_mid, beta_mid = specvol_alpha_beta(SA_mid,
CT_mid, p_mid)
N2 = ((g_local**2) / (specvol_mid * db_to_pa * dp))
N2 *= (beta_mid*dSA - alpha_mid*dCT)
return N2, p_mid
@match_args_return
def Turner_Rsubrho(SA, CT, p, axis=0):
"""
Calculate the Turner Angle and the Stability Ratio.
Parameters
----------
SA : array-like
Absolute Salinity, g/kg
CT : array-like
Conservative Temperature (ITS-90), degrees C
p : array-like
Sea pressure (absolute pressure minus 10.1325 dbar), dbar
axis : int, optional
The dimension along which pressure increases.
Returns
-------
Tu : array
Turner Angle at pressure midpoints, degrees.
The shape along the pressure axis dimension is one
less than that of the inputs.
Rsubrho : array
Stability Ratio, dimensionless.
The shape matches Tu.
p_mid : array
Pressure at midpoints of p, dbar.
The array shape matches Tu.
"""
SA = np.clip(SA, 0, 50)
SA, CT, p = np.broadcast_arrays(SA, CT, p)
shallow = axis_slicer(SA.ndim, slice(-1), axis)
deep = axis_slicer(SA.ndim, slice(1, None), axis)
dSA = -SA[deep] + SA[shallow]
dCT = -CT[deep] + CT[shallow]
SA_mid = 0.5 * (SA[shallow] + SA[deep])
CT_mid = 0.5 * (CT[shallow] + CT[deep])
p_mid = 0.5 * (p[shallow] + p[deep])
_, alpha, beta = specvol_alpha_beta(SA_mid, CT_mid, p_mid)
Tu = np.arctan2((alpha*dCT + beta*dSA), (alpha*dCT - beta*dSA))
Tu = np.degrees(Tu)
igood = (dSA != 0)
Rsubrho = np.zeros_like(dSA)
Rsubrho.fill(np.nan)
Rsubrho[igood] = (alpha[igood]*dCT[igood])/(beta[igood]*dSA[igood])
return Tu, Rsubrho, p_mid
@match_args_return
def IPV_vs_fNsquared_ratio(SA, CT, p, p_ref=0, axis=0):
"""
Calculates the ratio of the vertical gradient of potential density to
the vertical gradient of locally-referenced potential density. This
is also the ratio of the planetary Isopycnal Potential Vorticity
(IPV) to f times N^2, hence the name for this variable,
IPV_vs_fNsquared_ratio (see Eqn. (3.20.17) of IOC et al. (2010)).
Parameters
----------
SA : array-like
Absolute Salinity, g/kg
t : array-like
In-situ temperature (ITS-90), degrees C
p : array-like
Sea pressure (absolute pressure minus 10.1325 dbar), dbar
p_ref : float
Reference pressure, dbar
Returns
-------
IPV_vs_fNsquared_ratio : array
The ratio of the vertical gradient of
potential density referenced to p_ref, to the vertical
gradient of locally-referenced potential density, dimensionless.
p_mid : array
Pressure at midpoints of p, dbar.
The array shape matches IPV_vs_fNsquared_ratio.
"""
SA = np.clip(SA, 0, 50)
SA, CT, p = np.broadcast_arrays(SA, CT, p)
shallow = axis_slicer(SA.ndim, slice(-1), axis)
deep = axis_slicer(SA.ndim, slice(1, None), axis)
dSA = -SA[deep] + SA[shallow]
dCT = -CT[deep] + CT[shallow]
SA_mid = 0.5 * (SA[shallow] + SA[deep])
CT_mid = 0.5 * (CT[shallow] + CT[deep])
p_mid = 0.5 * (p[shallow] + p[deep])
_, alpha, beta = specvol_alpha_beta(SA_mid, CT_mid, p_mid)
_, alpha_pref, beta_pref = specvol_alpha_beta(SA_mid, CT_mid, p_ref)
num = dCT*alpha_pref - dSA*beta_pref
den = dCT*alpha - dSA*beta
igood = (den != 0)
IPV_vs_fNsquared_ratio = np.zeros_like(dSA)
IPV_vs_fNsquared_ratio.fill(np.nan)
IPV_vs_fNsquared_ratio[igood] = num[igood] / den[igood]
return IPV_vs_fNsquared_ratio, p_mid
|