Recovering the result of Reinhold and Pierrehumbert (1982)

In this example, we will recover the attractor shown in

obtained with a 2-layer channel QG atmosphere truncated at wavenumber 2 on a beta-plane and a simple orography (a mountain and a valley).

Modules import

Loading of some modules…

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

Initializing the random number generator (for reproducibility). – Disable if needed.


Importing the model’s modules

from qgs.params.params import QgParams
from qgs.integrators.integrator import RungeKuttaIntegrator, RungeKuttaTglsIntegrator
from qgs.functions.tendencies import create_tendencies
from qgs.plotting.util import std_plot

and diagnostics

from qgs.diagnostics.streamfunctions import MiddleAtmosphericStreamfunctionDiagnostic
from qgs.diagnostics.variables import VariablesDiagnostic
from qgs.diagnostics.multi import MultiDiagnostic

Systems definition

General parameters

# Time parameters
dt = 0.1
# Saving the model state n steps
write_steps = 5

number_of_trajectories = 1
number_of_perturbed_trajectories = 10

Setting some model parameters

# Model parameters instantiation with some non-default specs
model_parameters = QgParams({'phi0_npi': np.deg2rad(50.)/np.pi, 'hd':0.045})
# Mode truncation at the wavenumber 2 in both x and y spatial coordinate
model_parameters.set_atmospheric_channel_fourier_modes(2, 2)

# Setting the orography depth and the meridional temperature gradient
model_parameters.ground_params.set_orography(0.2, 1)
model_parameters.atemperature_params.set_thetas(0.1, 0)
# Printing the model's parameters
Qgs parameters summary

General Parameters:
'time_unit': days,
'rr': 287.058  [J][kg^-1][K^-1]  (gas constant of dry air),
'sb': 5.67e-08  [J][m^-2][s^-1][K^-4]  (Stefan-Boltzmann constant),

Scale Parameters:
'scale': 5000000.0  [m]  (characteristic space scale (L*pi)),
'f0': 0.0001032  [s^-1]  (Coriolis parameter at the middle of the domain),
'n': 1.3    (aspect ratio (n = 2 L_y / L_x)),
'rra': 6370000.0  [m]  (earth radius),
'phi0_npi': 0.2777777777777778    (latitude expressed in fraction of pi),
'deltap': 50000.0  [Pa]  (pressure difference between the two atmospheric layers),

Atmospheric Parameters:
'kd': 0.1  [nondim]  (atmosphere bottom friction coefficient),
'kdp': 0.01  [nondim]  (atmosphere internal friction coefficient),
'sigma': 0.2  [nondim]  (static stability of the atmosphere),

Atmospheric Temperature Parameters:
'hd': 0.045  [nondim]  (Newtonian cooling coefficient),
'thetas[1]': 0.1    (spectral components 1 of the temperature profile),
'thetas[2]': 0.0    (spectral component 2 of the temperature profile),
'thetas[3]': 0.0    (spectral component 3 of the temperature profile),
'thetas[4]': 0.0    (spectral component 4 of the temperature profile),
'thetas[5]': 0.0    (spectral component 5 of the temperature profile),
'thetas[6]': 0.0    (spectral component 6 of the temperature profile),
'thetas[7]': 0.0    (spectral component 7 of the temperature profile),
'thetas[8]': 0.0    (spectral component 8 of the temperature profile),
'thetas[9]': 0.0    (spectral component 9 of the temperature profile),
'thetas[10]': 0.0    (spectral component 10 of the temperature profile),

Ground Parameters:
'hk[1]': 0.0    (spectral component 1 of the orography),
'hk[2]': 0.2    (spectral components 2 of the orography),
'hk[3]': 0.0    (spectral component 3 of the orography),
'hk[4]': 0.0    (spectral component 4 of the orography),
'hk[5]': 0.0    (spectral component 5 of the orography),
'hk[6]': 0.0    (spectral component 6 of the orography),
'hk[7]': 0.0    (spectral component 7 of the orography),
'hk[8]': 0.0    (spectral component 8 of the orography),
'hk[9]': 0.0    (spectral component 9 of the orography),
'hk[10]': 0.0    (spectral component 10 of the orography),
'orographic_basis': atmospheric,

Creating the tendencies function

f, Df = create_tendencies(model_parameters)

Time integration and plotting a section of the attractor

Defining an integrator

integrator = RungeKuttaIntegrator()

Start on a random initial condition and integrate over a transient time to obtain an initial condition on the attractors

ic = np.random.rand(model_parameters.ndim)*0.1
integrator.integrate(0., 200000., dt, ic=ic, write_steps=0)
time, ic = integrator.get_trajectories()

Now integrate to obtain a trajectory on the attractor

integrator.integrate(0., 100000., dt, ic=ic, write_steps=write_steps)
time, traj = integrator.get_trajectories()

and plot \(\psi_2^{\rm a}\) as a function of \(\psi_3^{\rm a}\)

varx = 2
vary = 1
plt.figure(figsize=(10, 8))

plt.plot(traj[varx], traj[vary], marker='o', ms=0.07, ls='')


to recover the figure 3 of the Reinhold and Pierrehumbert article mentioned above.

Showing the resulting fields

One can use objects called Diagnostic to plot movies of the model fields. Here we show how to use them and plot the simultaneous time evolution of the variables \(\psi_2\) and \(\psi_3\), and the geopotential height field at 500 hPa over the orographic height.

Creating the diagnostics:

  • For the 500hPa geopotential height:

psi = MiddleAtmosphericStreamfunctionDiagnostic(model_parameters, geopotential=True)
  • For the nondimensional variables \(\psi_2\) and \(\psi_3\):

variable_nondim = VariablesDiagnostic([2, 1], model_parameters, False)
# setting also the background
background = VariablesDiagnostic([2, 1], model_parameters, False)
background.set_data(time, traj)

Creating a multi diagnostic with both:

m = MultiDiagnostic(1,2)
m.add_diagnostic(variable_nondim, diagnostic_kwargs={'show_time': False, 'background': background}, plot_kwargs={'ms': 0.2})
m.add_diagnostic(psi, diagnostic_kwargs={'style': 'contour', 'contour_labels': False}, plot_kwargs={'colors': 'k'})
m.set_data(time[:501], traj[:501])

and create the movie:,6), anim_kwargs={'interval': 100, 'frames': 500})