Skip to content

Circulation

get_stream(v, p, lat, ax_p=0)

Computes the streamfunction \(\psi\) at latitude \(\phi\) and pressure \(p\) according to:

\(\psi(\phi, p) = \int_0^z v\rho dx dz = 2\pi a \cos \phi\int_0^z v\rho dz = -\frac{2\pi a \cos \phi}{g} \int_{p_{surf}}^p vdp\)

Parameters:

Name Type Description Default
v ndarray

float [n_p_levels, n_lat].
Meridional velocity at each pressure and latitude. Units: m/s.

 required
p ndarray

float [n_p_levels].
Pressure levels corresponding to meridional velocity.
Must be descending so first value is closest to surface. Units: Pa.

 required
lat ndarray

float [n_lat].
Latitude corresponding to meridional velocity. Units: degrees.

 required
ax_p int

Axis corresponding to pressure levels in \(v\).

 0

Returns:

Type Description
ndarray

float [n_lat].
Streamfunction at pressure level given by p[-1]. Units: kg/s.
Source code in isca_tools/utils/circulation.py 
 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
def get_stream(v: np.ndarray, p: np.ndarray, lat: np.ndarray, ax_p: int = 0) -> np.ndarray:
    """
    Computes the streamfunction $\psi$ at latitude $\phi$ and pressure $p$
    [according to](https://sandrolubis.wordpress.com/2012/06/17/mass-streamfunction/):

    $\psi(\phi, p) = \int_0^z v\\rho dx dz = 2\pi a \cos \phi\int_0^z v\\rho dz =
    -\\frac{2\pi a \cos \phi}{g} \int_{p_{surf}}^p vdp$

    Args:
        v: `float [n_p_levels, n_lat]`. </br>
            Meridional velocity at each pressure and latitude. Units: *m/s*.
        p: `float [n_p_levels]`.</br>
            Pressure levels corresponding to meridional velocity.</br>
            Must be descending so first value is closest to surface. Units: *Pa*.
        lat: `float [n_lat]`.</br>
            Latitude corresponding to meridional velocity. Units: *degrees*.
        ax_p: Axis corresponding to pressure levels in $v$.

    Returns:
        `float [n_lat]`.</br>
            Streamfunction at pressure level given by `p[-1]`. Units: *kg/s*.

    """
    if len(p) > 1:
        if p[1] > p[0]:
            raise ValueError(f'Pressure is not in correct order, expect p[0] to be surface value but it is {p[0]}Pa')
    cos_lat = np.cos(np.deg2rad(lat))
    stream = -2 * np.pi * radius_earth * cos_lat / g * integrate.simpson(v, p, axis=ax_p)
    return stream

get_u_thermal(temp, p, lat, ax_p=0, ax_lat=1)

Computes thermal wind at pressure \(p\) and latitude \(\phi\) according to Equation 1 in shaw_2023 paper:

\(u_T(p, \phi) = \int_{p_s}^{p}\frac{R}{fap'}\frac{\partial T}{\partial \phi} dp'\)

where \(p_s\) is near-surface pressure, \(f\) is the coriolis parameter, \(R\) is the gas constant for dry air, \(a\) is the radius of the Earth and \(T\) is temperature.

Parameters:

Name Type Description Default
temp ndarray

float [n_p_levels, n_lat, ...].
Temperature at each pressure and latitude. Units: K.

 required
p ndarray

float [n_p_levels].
Pressure levels corresponding to meridional velocity.
Must be descending so first value is closest to surface. Units: Pa.

 required
lat ndarray

float [n_lat].
Latitude corresponding to meridional velocity. Units: degrees.

 required
ax_p int

Axis corresponding to pressure levels in temp.
Must be \(0\) or \(1\).

 0
ax_lat int

Axis corresponding to latitude in temp.
Must be \(0\) or \(1\).

 1

Returns:

Type Description
ndarray

float [n_lat, ...].
Thermal wind at pressure level given by p[-1]. Units: m/s.
Source code in isca_tools/utils/circulation.py 
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
def get_u_thermal(temp: np.ndarray, p: np.ndarray, lat: np.ndarray, ax_p: int = 0, ax_lat: int = 1) -> np.ndarray:
    """
    Computes thermal wind at pressure $p$ and latitude $\phi$ according to Equation 1
    in [shaw_2023](https://www.nature.com/articles/s41558-023-01884-1) paper:

    $u_T(p, \phi) = \int_{p_s}^{p}\\frac{R}{fap'}\\frac{\partial T}{\partial \phi} dp'$

    where $p_s$ is near-surface pressure, $f$ is the coriolis parameter, $R$ is the gas constant for dry air,
    $a$ is the radius of the Earth and $T$ is temperature.

    Args:
        temp: `float [n_p_levels, n_lat, ...]`. </br>
            Temperature at each pressure and latitude. Units: *K*.
        p: `float [n_p_levels]`.</br>
            Pressure levels corresponding to meridional velocity.</br>
            Must be descending so first value is closest to surface. Units: *Pa*.
        lat: `float [n_lat]`.</br>
            Latitude corresponding to meridional velocity. Units: *degrees*.
        ax_p: Axis corresponding to pressure levels in `temp`.</br>
            Must be $0$ or $1$.
        ax_lat: Axis corresponding to latitude in `temp`.</br>
            Must be $0$ or $1$.

    Returns:
        `float [n_lat, ...]`.</br>
            Thermal wind at pressure level given by `p[-1]`. Units: *m/s*.
    """
    if len(p) > 1:
        if p[1] > p[0]:
            raise ValueError(f'Pressure is not in correct order, expect p[0] to be surface value but it is {p[0]}Pa')
    n_ax = len(temp.shape)  # all axis 2 or higher are not p or lat
    ax_expand_dims = list(np.append(np.asarray([ax_lat]), np.arange(2, n_ax)))
    integrand = np.gradient(temp, np.deg2rad(lat), axis=ax_lat) / np.expand_dims(p, axis=ax_expand_dims)
    f_coriolis = 2 * rot_earth * np.sin(np.deg2rad(lat).to_numpy())
    if n_ax > 2:
        f_coriolis = np.expand_dims(f_coriolis, list(np.arange(2, n_ax) - 1))
    return integrate.simpson(integrand, p, axis=ax_p) * R / (radius_earth * f_coriolis)