Skip to content

Miyawaki 2022

get_dmse_dt(temp, sphum, height, p_levels, time, zonal_mean=True, spline_smoothing_factor=0)

For a given latitude, this computes the time derivative of the mass weighted vertical integral of the zonal mean moist static energy, \(<[\partial_t m]>\), used in the paper to compute the parameter \(R_1\).

Parameters:

Name Type Description Default
temp DataArray

float [n_time, n_p_levels, n_lon].
Temperature at each coordinate considered. Units: Kelvin.

 required
sphum DataArray

float [n_time, n_p_levels, n_lon].
Specific humidity at each coordinate considered. Units: kg/kg.

 required
height DataArray

float [n_time, n_p_levels, n_lon].
Geopotential height of each level considered.

 required
p_levels DataArray

float [n_p_levels].
Pressure levels of atmosphere, p_levels[0] is top of atmosphere and p_levels[-1] is the surface. Units: Pa.

 required
time DataArray

float [n_time].
Time at which data recorded. Units: second.

 required
zonal_mean bool

bool.
Whether to take zonal average or not. If False, returned arrays will have size [n_time x n_lon]

 True
spline_smoothing_factor float

float.
If positive, a spline will be fit to smooth mse_integ and compute dmse_dt. Typical: 0.001.
If 0, the deriviative will be computed using np.gradient, but in general I recommend using the spline.

 0

Returns:

Type Description
DataArray

mse_integ: float [n_time] or float [n_time x n_lon]
The mass weighted vertical integral of the (zonal mean) moist static energy at each time i.e. \(<[m]>\). Units: \(J/m^2\).
DataArray

dmse_dt: float [n_time] or float [n_time x n_lon]
The time derivative of the mass weighted vertical integral of the (zonal mean) moist static energy at each time i.e. \(<[\partial_t m]>\). Units: \(W/m^2\).
Source code in isca_tools/papers/miyawaki_2022.py 
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
def get_dmse_dt(temp: xr.DataArray, sphum: xr.DataArray, height: xr.DataArray, p_levels: xr.DataArray,
                time: xr.DataArray, zonal_mean: bool = True,
                spline_smoothing_factor: float = 0) -> Tuple[xr.DataArray, xr.DataArray]:
    """
    For a given latitude, this computes the time derivative of the mass weighted vertical integral of the zonal mean
    moist static energy, $<[\\partial_t m]>$, used in the paper to compute the parameter $R_1$.

    Args:
        temp: `float [n_time, n_p_levels, n_lon]`.</br>
            Temperature at each coordinate considered. Units: *Kelvin*.
        sphum: `float [n_time, n_p_levels, n_lon]`.</br>
            Specific humidity at each coordinate considered. Units: *kg/kg*.
        height: `float [n_time, n_p_levels, n_lon]`.</br>
            Geopotential height of each level considered.
        p_levels: `float [n_p_levels]`.</br>
            Pressure levels of atmosphere, `p_levels[0]` is top of atmosphere and `p_levels[-1]` is the surface.
            Units: *Pa*.
        time: `float [n_time]`.</br>
            Time at which data recorded. Units: *second*.
        zonal_mean: `bool`.</br>
            Whether to take zonal average or not. If `False`, returned arrays will have size `[n_time x n_lon]`
        spline_smoothing_factor: `float`.</br>
            If positive, a spline will be fit to smooth `mse_integ` and compute `dmse_dt`. *Typical: 0.001*.</br>
            If 0, the deriviative will be computed using `np.gradient`, but in general I recommend using the spline.

    Returns:
        `mse_integ`: `float [n_time]` or `float [n_time x n_lon]`</br>
            The mass weighted vertical integral of the (zonal mean) moist static energy at each time i.e. $<[m]>$.
            Units: $J/m^2$.
        `dmse_dt`: `float [n_time]` or `float [n_time x n_lon]`</br>
            The time derivative of the mass weighted vertical integral of the (zonal mean) moist static energy at
            each time i.e. $<[\\partial_t m]>$. Units: $W/m^2$.
    """
    # compute zonal mean mse in units of J/kg
    mse = moist_static_energy(temp, sphum, height) * 1000
    if 'lon' in list(temp.coords.keys()) and zonal_mean:
        # Take zonal mean if longitude is a coordinate
        mse = mse.mean(dim='lon')
    mse_integ = integrate.simpson(mse / g, p_levels, axis=1)    # do mass weighted integral over atmosphere
    if spline_smoothing_factor > 0:
        if not zonal_mean:
            raise ValueError('Can only use spline if taking the zonal mean')
        # Divide by mean to fit spline as otherwise numbers too large to fit properly
        spl = UnivariateSpline(time, mse_integ / np.mean(mse_integ), s=spline_smoothing_factor)
        dmse_dt = spl.derivative()(time) * np.mean(mse_integ)      # compute derivative direct from spline fit
        mse_integ = spl(time) * np.mean(mse_integ)    # Set mse_integ to spline version so this is what is returned
    elif spline_smoothing_factor == 0:
        # compute derivative directly from the data
        dmse_dt = np.gradient(mse_integ, time, axis=0)
    else:
        raise ValueError(f'spline_smoothing_factor = {spline_smoothing_factor}, which is negative. It must be >=0.')
    return mse_integ, dmse_dt

get_dvmse_dy(R_a, lh, sh, dmse_dt)

This infers the divergence of moist static energy flux, \(<[\partial_y vm]>\), from equation \((2)\) in the paper:

\(<[\partial_t m]> + <[\partial_y vm]> = [R_a] + [LH] + [SH]\)

It is not computed directly for reasons outlined in section 2b of the paper.

Parameters:

Name Type Description Default
R_a DataArray

float [n_time].
Atmospheric radiative heating rate i.e. difference between top of atmosphere and surface radiative fluxes. Negative means atmosphere is cooling. Units: \(W/m^2\).
This is what is returned by frierson_atmospheric_heating if using the grey radiation with the Frierson scheme.

 required
lh DataArray

float [n_time].
Surface latent heat flux (up is positive). This is saved by Isca if the variable flux_lhe in the mixed_layer module is specified in the diagnostic table. Units: \(W/m^2\).

 required
sh DataArray

float [n_time].
Surface sensible heat flux (up is positive). This is saved by Isca if the variable flux_t in the mixed_layer module is specified in the diagnostic table. Units: \(W/m^2\).

 required
dmse_dt DataArray

float [n_time].
The time derivative of the mass weighted vertical integral of the zonal mean moist static energy at each time, \(<[\partial_t m]>\). Units: \(W/m^2\).

 required

Returns:

Type Description
DataArray

float [n_time].
\(<[\partial_y vm]>\) - Mass weighted vertical integral of the zonal mean meridional divergence of moist static energy flux, \(vm\).
Source code in isca_tools/papers/miyawaki_2022.py 
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
def get_dvmse_dy(R_a: xr.DataArray, lh: xr.DataArray, sh: xr.DataArray, dmse_dt: xr.DataArray) -> xr.DataArray:
    """
    This infers the divergence of moist static energy flux, $<[\\partial_y vm]>$, from equation $(2)$ in the paper:

    $<[\\partial_t m]> + <[\\partial_y vm]> = [R_a] + [LH] + [SH]$

    It is not computed directly for reasons outlined in section 2b of the paper.

    Args:
        R_a: `float [n_time]`.</br>
            Atmospheric radiative heating rate i.e. difference between top of atmosphere and surface radiative fluxes.
            Negative means atmosphere is cooling. Units: $W/m^2$.</br>
            This is what is returned by `frierson_atmospheric_heating` if using the grey radiation with the
            *Frierson* scheme.
        lh: `float [n_time]`.</br>
            Surface latent heat flux (up is positive). This is saved by *Isca* if the variable `flux_lhe` in the
            `mixed_layer` module is specified in the diagnostic table. Units: $W/m^2$.
        sh: `float [n_time]`.</br>
            Surface sensible heat flux (up is positive). This is saved by *Isca* if the variable `flux_t` in the
            `mixed_layer` module is specified in the diagnostic table. Units: $W/m^2$.
        dmse_dt: `float [n_time]`.</br>
            The time derivative of the mass weighted vertical integral of the zonal mean moist static energy at
            each time, $<[\\partial_t m]>$. Units: $W/m^2$.

    Returns:
        `float [n_time]`.</br>
            $<[\\partial_y vm]>$ - Mass weighted vertical integral of the zonal mean meridional divergence of
            moist static energy flux, $vm$.
    """
    return R_a + lh + sh - dmse_dt

get_lapse_dev(temp, pressure, pressure_surf, lev_dim='lev')

Gets the lapse rate deviation from moist adiabat parameters given in the paper.

Parameters:

Name Type Description Default
temp DataArray

float [n_lev].
Temperature on model levels

 required
pressure DataArray

float [n_lev].
Pressure at model levels

 required
pressure_surf DataArray

float.
Surface pressure, used to compute \(\sigma\) coordinate.

 required
lev_dim

Name of model level dimension, over which to vertically average

 'lev'

Returns:

Type Description
DataArray

Aloft lapse rate deviation: \(\langle(\Gamma_m - \Gamma)/\Gamma_m\rangle_{0.7-0.3}\)
Computed between \(\sigma=0.3\) and \(\sigma=0.7\). Given as a %.
DataArray

Boundary layer lapse rate deviation: \(\langle(\Gamma_m - \Gamma)/\Gamma_m \rangle_{1.0-0.9}\)
Computed between \(\sigma=0.9\) and \(\sigma=1.0\). Given as a %.
Source code in isca_tools/papers/miyawaki_2022.py 
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
def get_lapse_dev(temp: xr.DataArray, pressure: xr.DataArray, pressure_surf: xr.DataArray,
                  lev_dim='lev') -> Tuple[xr.DataArray, xr.DataArray]:
    """
    Gets the lapse rate deviation from moist adiabat parameters given in the paper.

    Args:
        temp: `float [n_lev]`.</br>
            Temperature on model levels
        pressure: `float [n_lev]`.</br>
            Pressure at model levels
        pressure_surf: `float`.</br>
            Surface pressure, used to compute $\sigma$ coordinate.
        lev_dim: Name of model level dimension, over which to vertically average

    Returns:
        Aloft lapse rate deviation: $\langle(\Gamma_m - \Gamma)/\Gamma_m\\rangle_{0.7-0.3}$</br>
            Computed between $\sigma=0.3$ and $\sigma=0.7$. Given as a %.
        Boundary layer lapse rate deviation: $\langle(\Gamma_m - \Gamma)/\Gamma_m \\rangle_{1.0-0.9}$</br>
            Computed between $\sigma=0.9$ and $\sigma=1.0$. Given as a %.
    """
    grad_xr = wrap_with_apply_ufunc(np.gradient, input_core_dims=[[lev_dim], [lev_dim]], output_core_dims=[[lev_dim]])
    lapse_env = grad_xr(temp, pressure)
    lapse_moist_xr = wrap_with_apply_ufunc(lapse_moist)
    lapse_adiabat = lapse_moist_xr(temp, pressure, pressure_coords=True)

    # Do mass weighted vertical average
    dp = np.abs(pressure.differentiate(lev_dim))  # Pa
    weights = (dp / g)  # mass weights (kg/m²)

    sigma = pressure / pressure_surf
    use_sigma = {'aloft': np.logical_and(sigma>=0.3, sigma<=0.7),
                 'boundary_layer': np.logical_and(sigma>=0.9, sigma<=1)}
    lapse_dev = {}
    for key in use_sigma:
        lapse_dev[key] = ((lapse_adiabat - lapse_env)*weights).where(use_sigma[key]).sum(dim=lev_dim) / \
                         (lapse_adiabat*weights).where(use_sigma[key]).sum(dim=lev_dim) * 100
    return tuple(lapse_dev[key] for key in lapse_dev)

get_r1(R_a, dmse_dt, dvmse_dy, time=None, spline_smoothing_factor=0)

Returns the non-dimensional number, \(R_1 = \frac{\partial_t m + \partial_y vm}{R_a}\).

Parameters:

Name Type Description Default
R_a DataArray

float [n_time].
Atmospheric radiative heating rate i.e. difference between top of atmosphere and surface radiative fluxes. Negative means atmosphere is cooling. Units: \(W/m^2\).
This is what is returned by frierson_atmospheric_heating if using the grey radiation with the Frierson scheme.

 required
dmse_dt DataArray

float [n_time].
The time derivative of the mass weighted vertical integral of the zonal mean moist static energy at each time, \(<[\partial_t m]>\). Units: \(W/m^2\).

 required
dvmse_dy DataArray

float [n_time].
\(<[\partial_y vm]>\) - Mass weighted vertical integral of the zonal mean meridional divergence of moist static energy flux, \(vm\).

 required
time Optional[DataArray]

None or float [n_time].
Time at which data recorded. Only required if using the spline_smoothing_factor>0. Units: second.

 None
spline_smoothing_factor float

float.
If positive, a spline will be fit to smooth \(R_1\). Typical: 0.2.

 0

Returns:

Type Description
DataArray

float [n_time].
DataArray

The non-dimensional number, \(R_1 = \frac{\partial_t m + \partial_y vm}{R_a}\).
Source code in isca_tools/papers/miyawaki_2022.py 
 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
def get_r1(R_a: xr.DataArray, dmse_dt: xr.DataArray, dvmse_dy: xr.DataArray,
           time: Optional[xr.DataArray] = None, spline_smoothing_factor: float = 0) -> xr.DataArray:
    """
    Returns the non-dimensional number, $R_1 = \\frac{\\partial_t m + \\partial_y vm}{R_a}$.

    Args:
        R_a: `float [n_time]`.</br>
            Atmospheric radiative heating rate i.e. difference between top of atmosphere and surface radiative fluxes.
            Negative means atmosphere is cooling. Units: $W/m^2$.</br>
            This is what is returned by `frierson_atmospheric_heating` if using the grey radiation with the
            *Frierson* scheme.
        dmse_dt: `float [n_time]`.</br>
            The time derivative of the mass weighted vertical integral of the zonal mean moist static energy at
            each time, $<[\\partial_t m]>$. Units: $W/m^2$.
        dvmse_dy: `float [n_time]`.</br>
            $<[\\partial_y vm]>$ - Mass weighted vertical integral of the zonal mean meridional divergence of
            moist static energy flux, $vm$.
        time: `None` or `float [n_time]`.</br>
            Time at which data recorded. Only required if using the `spline_smoothing_factor>0`. Units: *second*.
        spline_smoothing_factor: `float`.</br>
            If positive, a spline will be fit to smooth $R_1$. *Typical: 0.2*.</br>

    Returns:
        `float [n_time]`.</br>
        The non-dimensional number, $R_1 = \\frac{\\partial_t m + \\partial_y vm}{R_a}$.

    """
    r1 = (dmse_dt + dvmse_dy)/R_a
    if spline_smoothing_factor > 0:
        if time is None:
            raise ValueError('Specified spline but no time array given.')
        spl = UnivariateSpline(time, r1, s=spline_smoothing_factor)
        r1 = spl(time)
    elif spline_smoothing_factor < 0:
        raise ValueError(f'spline_smoothing_factor = {spline_smoothing_factor}, which is negative. It must be >=0.')
    return r1