This model version is a simple 2-layer channel QG atmosphere truncated at wavenumber 2 on a beta-plane with a simple orography (a montain and a valley) , with a radiative and heat exchanges scheme inspired by the MAOOAM model.
More detail can be found in the articles:
or in the documentation.
First, setting the path and loading of some modules
import sys, os
sys.path.extend([os.path.abspath('../')])
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.
np.random.seed(210217)
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
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, 'n':1.5})
# Mode truncation at the wavenumber 2 in both x and y spatial coordinate for the atmosphere
model_parameters.set_atmospheric_channel_fourier_modes(2, 2)
# Same modes for the ground temperature modes
model_parameters.set_ground_channel_fourier_modes()
# Changing (increasing) the orography depth
model_parameters.ground_params.set_orography(0.2, 1)
# Setting the parameters of the heat transfer from the soil
model_parameters.gotemperature_params.set_params({'gamma': 1.6e7})
model_parameters.atemperature_params.set_params({ 'hlambda': 10.})
model_parameters.atemperature_params.set_insolation(300*0.33, 0)
model_parameters.gotemperature_params.set_insolation(300, 0)
# Printing the model's parameters
model_parameters.print_params()
Creating the tendencies function
from qgs.inner_products.symbolic import AtmosphericSymbolicInnerProducts
aip = AtmosphericSymbolicInnerProducts(model_parameters)
%%time
f, Df = create_tendencies(model_parameters)
Defining an integrator
integrator = RungeKuttaIntegrator()
integrator.set_func(f)
Start on a random initial condition and integrate over a transient time to obtain an initial condition on the attractors
%%time
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
%%time
integrator.integrate(0., 200000., dt, ic=ic, write_steps=write_steps)
time, traj = integrator.get_trajectories()
Plotting the result in 3D and 2D
varx = 20
vary = 21
varz = 1
fig = plt.figure(figsize=(10, 8))
axi = fig.gca(projection='3d')
axi.scatter(traj[varx], traj[vary], traj[varz], s=0.2);
axi.set_xlabel('$'+model_parameters.latex_var_string[varx]+'$')
axi.set_ylabel('$'+model_parameters.latex_var_string[vary]+'$')
axi.set_zlabel('$'+model_parameters.latex_var_string[varz]+'$');
varx = 20
vary = 21
plt.figure(figsize=(10, 8))
plt.plot(traj[varx], traj[vary], marker='o', ms=0.07, ls='')
plt.xlabel('$'+model_parameters.latex_var_string[varx]+'$')
plt.ylabel('$'+model_parameters.latex_var_string[vary]+'$');
varx = 2
vary = 1
plt.figure(figsize=(10, 8))
plt.plot(traj[varx], traj[vary], marker='o', ms=0.07, ls='')
plt.xlabel('$'+model_parameters.latex_var_string[varx]+'$')
plt.ylabel('$'+model_parameters.latex_var_string[vary]+'$');
var = 21
plt.figure(figsize=(10, 8))
plt.plot(model_parameters.dimensional_time*time, traj[var])
plt.xlim(0.,1000.)
plt.xlabel('time (days)')
plt.ylabel('$'+model_parameters.latex_var_string[var]+'$');
Instantiating a tangent linear integrator with the model tendencies
tgls_integrator = RungeKuttaTglsIntegrator()
tgls_integrator.set_func(f, Df)
Integrating with slightly perturbed initial conditions
tangent_ic = 0.00001*np.random.randn(number_of_perturbed_trajectories, model_parameters.ndim)
%%time
tgls_integrator.integrate(0., 2000., dt=dt, write_steps=write_steps, ic=ic, tg_ic=tangent_ic)
Obtaining the perturbed trajectories
tt, traj, delta = tgls_integrator.get_trajectories()
pert_traj = traj + delta
Trajectories plot
var = 10
plt.figure(figsize=(10, 8))
plt.plot(model_parameters.dimensional_time*tt, traj[var])
plt.plot(model_parameters.dimensional_time*tt, pert_traj[:,var].T, ls=':')
ax = plt.gca()
plt.ylim(0.,0.1)
plt.xlim(0.,50.)
plt.xlabel('time (days)')
plt.ylabel('$'+model_parameters.latex_var_string[var]+'$');
Mean and standard deviation
var = 10
plt.figure(figsize=(10, 8))
plt.plot(model_parameters.dimensional_time*tt, traj[var])
ax = plt.gca()
std_plot(model_parameters.dimensional_time*tt, np.mean(pert_traj[:,var], axis=0), np.sqrt(np.var(pert_traj[:, var], axis=0)), ax=ax, alpha=0.5)
plt.ylim(0.,0.1)
plt.xlim(0.,50.)
plt.xlabel('time (days)')
plt.ylabel('$'+model_parameters.latex_var_string[var]+'$');