<jupyter_start><jupyter_text>Rossby Wave Packet Diagnostic (Ghinassi et al. 2018)A reproduction of the results of section 4, Large-amplitude Rossby wave packets in an idealized barotropic model simulation, by [Ghinassi et al. (2018)](https://doi.org/10.1175/MWR-D-18-0068.1).<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import barotropic as bt<jupyter_output><empty_output><jupyter_text>Model setup and integration<jupyter_code># 2° resolution of latitude-longitude grid
grid = bt.Grid(resolution=2.)
# No forcing
rhs = None
# Model with hyperdiffusion
model = bt.BarotropicModel(rhs, diffusion_order=2, diffusion_coeff=1.0e15)
# Start from Held and Phillips (1987) zonal profile
background_u = bt.init.held_1985(grid)

# Add perturbation according to eqs. (24) to (26)
# λ_2 = 3π/2 here instead of 2π/3 as in the paper, which seems to be a typo
def perturbation_vorticity(grid, s1, s2, zeta0, phi0, lam1, lam2, sigma_phi, sigma_lam):
    L1 = np.exp(-((grid.lam - lam1)/sigma_lam)**2) * np.sin(s1 * grid.lam)
    L2 = np.exp(-((grid.lam - lam2)/sigma_lam)**2) * np.cos(s2 * grid.lam)
    return zeta0 * np.cos(grid.phi) * np.exp(-((grid.phi - phi0)/sigma_phi)**2) * (L1 + L2)

perturbation = perturbation_vorticity(grid, s1=6, s2=9, zeta0=8.0e-5, phi0=np.pi/4,
        lam1=np.pi/3, lam2=3*np.pi/2, sigma_phi=np.deg2rad(10), sigma_lam=np.deg2rad(30))

# Add vorticity perturbation to background state to obtain initial state
init = background_u.add_vorticity(perturbation)

# Run model for 5 days with a timestep of 900 seconds and save state every 24 hours
last, states = model.run(init, 900, 5*bt.DAY, save_every=1*bt.DAY)<jupyter_output><empty_output><jupyter_text>Figure 5Absolute vorticity in the barotropic model simulation (a) at initial time and (b) after 5 days of integration.<jupyter_code>def style_lonlat(ax):
    ax.set_xticks([0, 60, 120, 180, 240, 300])
    ax.set_xticklabels(["0°", "60°E", "120°E", "180°", "120°W", "60°W"])
    ax.set_ylim([0, 90])
    ax.set_yticks([0, 10, 20, 30, 40, 50, 60, 70, 80])
    ax.set_yticklabels(["0°", "10°N", "20°N", "30°N", "40°N", "50°N", "60°N", "70°N", "80°N"])
    ax.grid(linestyle="dotted", color="k")
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(6, 5), gridspec_kw=dict(height_ratios=(8, 8, 1)))

pv_levels = [0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6]
pv_cmap = "viridis"

ax1.set_title("(a)", loc="left")
ax1.set_title("Absolute vort (t = {:.0f} h)".format(init.time / bt.HOUR), loc="center")
ax1.contourf(grid.lon, grid.lat, init.pv * 10000,
             levels=pv_levels, cmap=pv_cmap, extend="both")
style_lonlat(ax1)

ax2.set_title("(b)", loc="left")
ax2.set_title("Absolute vort (t = {:.0f} h)".format(last.time / bt.HOUR), loc="center")
pv2 = ax2.contourf(grid.lon, grid.lat, last.pv * 10000,
                   levels=pv_levels, cmap=pv_cmap, extend="both")
plt.colorbar(pv2, cax=ax3, orientation="horizontal")
ax3.set_title("$10^{-4}\;\mathrm{s}^{-1}$", loc="right")
style_lonlat(ax2)

fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Figure 6Evolution of (left) filtered LWA and (right) the envelope of the meridional wind in the barotropic model simulation, shown as snapshots every 24h.<jupyter_code>fig, axs = plt.subplots(7, 2, figsize=(10, 12), gridspec_kw=dict(height_ratios=(7, 7, 7, 7, 7, 7, 1)))

panels = list("lfkejdichbga")

for state, (axl, axr) in zip(states, axs):
    axl.set_title("(" + panels.pop() + ")", loc="left")
    axl.set_title("Filtered LWA (t = {:.0f}h)".format(state.time / bt.HOUR), loc="center")
    ml = axl.contourf(grid.lon, grid.lat, state.falwa_filtered,
                      levels=[0, 4, 8, 12, 16, 20], extend="max", cmap="YlOrRd")
    axr.set_title("(" + panels.pop() + ")", loc="left")
    axr.set_title("Envelope of v (t = {:.0f}h)".format(state.time / bt.HOUR), loc="center")
    mr = axr.contourf(grid.lon, grid.lat, state.envelope_hilbert,
                      levels=[4, 8, 12, 16, 20], extend="both", cmap="Purples")

plt.colorbar(ml, cax=axs[-1,0], orientation="horizontal")
plt.colorbar(mr, cax=axs[-1,1], orientation="horizontal")

axs[-1,0].set_title("$\mathrm{m}\;\mathrm{s}^{-1}$", loc="right")
axs[-1,1].set_title("$\mathrm{m}\;\mathrm{s}^{-1}$", loc="right")

for ax in axs[:-1,:].flatten():
    style_lonlat(ax)

fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Figure A1Initial state of the barotropic model simulation: (a) meridional wind $v$ and (b) local dominant wavenumber from wavelet analysis.<jupyter_code>fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(6, 6), gridspec_kw=dict(height_ratios=(8, 1, 8, 1)))

ax1.set_title("(a)", loc="left")
ax1.set_title("Meridional wind (t = {:.0f} h)".format(init.time / bt.HOUR), loc="center")
m1 = ax1.contourf(grid.lon, grid.lat, init.v,
                  levels=[-30, -20, -10, -5, 5, 10, 20, 30], cmap="RdBu_r", extend="both")
plt.colorbar(m1, cax=ax2, orientation="horizontal")
ax2.set_title("$\mathrm{m}\;\mathrm{s}^{-1}$", loc="right")
style_lonlat(ax1)

ax3.set_title("(b)", loc="left")
ax3.set_title("Local dominant wavenumber (t = {:.0f} h)".format(init.time / bt.HOUR), loc="center")
m3 = ax3.contourf(grid.lon, grid.lat, init.dominant_wavenumber,
                  levels=[1, 2, 3, 4, 5, 6, 7, 8, 9], cmap="Spectral_r", extend="both")
plt.colorbar(m3, cax=ax4, orientation="horizontal")
ax4.set_title("dimensionless", loc="right")
style_lonlat(ax3)

fig.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18798235942668137, 'fraction_numerical': 0.05880191106210952, 'mean_word_length': 4.884324324324324, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '301322', 'license': 'Apache-2.0', 'n_tokens_mistral': 2406, 'n_tokens_neox': 2175, 'n_words': 478, 'path': 'examples/rwp-diagnostic-ghinassi-et-al-2018.ipynb', 'repo_name': 'chpolste/barotropic', 'size': '5442'}