QGS model: MAOOAM with dynamic reference temperatures¶

Coupled ocean-atmosphere model version¶

This model version is a 2-layer channel QG atmosphere truncated at wavenumber 2 coupled, both by friction and heat exchange, to a shallow water ocean with 8 modes, with dynamic reference temperature evolution in both components.

More details can be found in the articles:

Vannitsem, S., Demaeyer, J., De Cruz, L., & Ghil, M. (2015). Low-frequency variability and heat transport in a low-order nonlinear coupled ocean–atmosphere model. Physica D: Nonlinear Phenomena, 309, 71-85. doi:10.1016/j.physd.2015.07.006
De Cruz, L., Demaeyer, J. and Vannitsem, S. (2016). The Modular Arbitrary-Order Ocean-Atmosphere Model: MAOOAM v1.0, Geosci. Model Dev., 9, 2793-2808. doi:10.5194/gmd-9-2793-2016
TBD

or in the documentation and on readthedocs.

Modules import¶

First, setting the path and loading of some modules

In [ ]:

import sys, os

In [ ]:

sys.path.extend([os.path.abspath('../')])

In [ ]:

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

In [ ]:

from matplotlib import rc
rc('font',**{'family':'serif','sans-serif':['Times'],'size':12})

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

In [ ]:

np.random.seed(210217)

Importing the model's modules

In [ ]:

from qgs.params.params import QgParams
from qgs.integrators.integrator import RungeKuttaIntegrator
from qgs.functions.tendencies import create_tendencies

and diagnostics

In [ ]:

from qgs.diagnostics.streamfunctions import MiddleAtmosphericStreamfunctionDiagnostic, OceanicLayerStreamfunctionDiagnostic
from qgs.diagnostics.temperatures import MiddleAtmosphericTemperatureDiagnostic, OceanicLayerTemperatureDiagnostic
from qgs.diagnostics.variables import VariablesDiagnostic, GeopotentialHeightDifferenceDiagnostic
from qgs.diagnostics.multi import MultiDiagnostic

Systems definition¶

General parameters

In [ ]:

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

number_of_trajectories = 1

Setting some model parameters

In [ ]:

# Model parameters instantiation with some non-default specs
model_parameters = QgParams({'n': 1.5}, dynamic_T=True)

# Modes definition: with dynamic reference temperature, the user has to use the symbolic modes
# Mode truncation at the wavenumber 2 in both x and y spatial
# coordinates for the atmosphere
model_parameters.set_atmospheric_channel_fourier_modes(2, 2, mode="symbolic")
# Mode truncation at the wavenumber 2 in the x and at the 
# wavenumber 4 in the y spatial coordinates for the ocean
model_parameters.set_oceanic_basin_fourier_modes(2, 4, mode="symbolic")

In [ ]:

# Setting MAOOAM parameters according to the publication linked above
model_parameters.set_params({'kd': 0.0290, 'kdp': 0.0290, 'r': 1.e-7,
                             'h': 136.5, 'd': 1.1e-7})
model_parameters.atemperature_params.set_params({'eps': 0.7, 'hlambda': 15.06})
model_parameters.gotemperature_params.set_params({'gamma': 5.6e8})

Setting the short-wave radiation component as in the publication above: $C_{\text{a},1}$ and $C_{\text{o},1}$

In [ ]:

model_parameters.atemperature_params.set_insolation(103., 0)
model_parameters.atemperature_params.set_insolation(103., 1)
model_parameters.gotemperature_params.set_insolation(310., 0)
model_parameters.gotemperature_params.set_insolation(310., 1)

Printing the model's parameters

In [ ]:

model_parameters.print_params()

Creating the tendencies function

In [ ]:

%%time
f, Df = create_tendencies(model_parameters)

Time integration¶

Defining an integrator

In [ ]:

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

In [ ]:

%%time
## Might take several minutes, depending on your cpu computational power.
ic = np.random.rand(model_parameters.ndim)*0.01
ic[29] = 3.  # Setting reasonable initial reference temperature
ic[10] = 1.5
integrator.integrate(0., 2000000.1, dt, ic=ic, write_steps=0)
time, ic = integrator.get_trajectories()

Now integrate to obtain a trajectory on the attractor

In [ ]:

%%time
integrator.integrate(0., 3000000., dt, ic=ic, write_steps=write_steps)
reference_time, reference_traj = integrator.get_trajectories()

Plotting the result in 3D and 2D

In [ ]:

varx = 22
vary = 31
varz = 0

fig = plt.figure(figsize=(10, 8))
axi = fig.add_subplot(111, projection='3d')

axi.scatter(reference_traj[varx], reference_traj[vary], reference_traj[varz], s=0.2);

axi.set_xlabel('$'+model_parameters.latex_var_string[varx]+'$', labelpad=12.)
axi.set_ylabel('$'+model_parameters.latex_var_string[vary]+'$', labelpad=12.)
axi.set_zlabel('$'+model_parameters.latex_var_string[varz]+'$', labelpad=12.);

In [ ]:

varx = 22
vary = 31
plt.figure(figsize=(10, 8))

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

plt.xlabel('$'+model_parameters.latex_var_string[varx]+'$')
plt.ylabel('$'+model_parameters.latex_var_string[vary]+'$');

In [ ]:

var = 10
plt.figure(figsize=(10, 8))

plt.plot(model_parameters.dimensional_time*reference_time, 2* model_parameters.temperature_scaling * reference_traj[var], label='$'+model_parameters.latex_var_string[var]+'$ [K]')
plt.plot(model_parameters.dimensional_time*reference_time, model_parameters.temperature_scaling * reference_traj[29], label='$'+model_parameters.latex_var_string[29]+'$ [K]')


plt.xlabel('time (days)')
#plt.ylabel();
plt.legend();

In [ ]:

var = 22
plt.figure(figsize=(10, 8))

plt.plot(model_parameters.dimensional_time*reference_time, reference_traj[var])


plt.xlabel('time (days)')
plt.ylabel('$'+model_parameters.latex_var_string[var]+'$');

Showing the resulting fields (animation)¶

This is an advanced feature showing the time evolution of diagnostic of the model. It shows simultaneously a scatter plot of the variables $\psi_{{\rm a}, 1}$, $\psi_{{\rm o}, 2}$ and $\delta T_{{\rm o}, 2}$, with the corresponding atmospheric and oceanic streamfunctions and temperature at 500 hPa. Please read the documentation for more information.

Creating the diagnostics:

For the 500hPa geopotential height:

In [ ]:

psi_a = MiddleAtmosphericStreamfunctionDiagnostic(model_parameters, geopotential=True, delta_x=0.2, delta_y=0.2)

For the 500hPa atmospheric temperature:

In [ ]:

theta_a = MiddleAtmosphericTemperatureDiagnostic(model_parameters, delta_x=0.2, delta_y=0.2)

For the ocean streamfunction:

In [ ]:

psi_o = OceanicLayerStreamfunctionDiagnostic(model_parameters, delta_x=0.2, delta_y=0.2)

For the ocean temperature:

In [ ]:

theta_o = OceanicLayerTemperatureDiagnostic(model_parameters, delta_x=0.2, delta_y=0.2)

For the nondimensional variables $\psi_{{\rm a}, 1}$, $\psi_{{\rm o}, 2}$ and $\delta T_{{\rm o}, 2}$:

In [ ]:

variable_nondim = VariablesDiagnostic([22, 31, 0], model_parameters, False)

For the geopotential height difference between North and South:

In [ ]:

geopot_dim = GeopotentialHeightDifferenceDiagnostic([[[np.pi/model_parameters.scale_params.n, np.pi/4], [np.pi/model_parameters.scale_params.n, 3*np.pi/4]]],
                                                    model_parameters, True)

In [ ]:

# setting also the background
background = VariablesDiagnostic([22, 31, 0], model_parameters, False)
background.set_data(reference_time, reference_traj)

Selecting a subset of the data to plot:

In [ ]:

stride = 10
time = reference_time[10000:10000+5200*stride:stride]
traj = reference_traj[:, 10000:10000+5200*stride:stride]

Creating a multi diagnostic with all the diagnostics:

In [ ]:

m = MultiDiagnostic(2,3)
m.add_diagnostic(geopot_dim, diagnostic_kwargs={'style':'moving-timeserie'})
m.add_diagnostic(theta_a, diagnostic_kwargs={'show_time':False})
m.add_diagnostic(theta_o, diagnostic_kwargs={'show_time':False})
m.add_diagnostic(variable_nondim, diagnostic_kwargs={'show_time':False, 'background': background, 'style':'3Dscatter'}, plot_kwargs={'ms': 0.2})
m.add_diagnostic(psi_a, diagnostic_kwargs={'show_time':False})
m.add_diagnostic(psi_o, diagnostic_kwargs={'show_time':False})

m.set_data(time, traj)

and show an interactive animation:

In [ ]:

rc('font',**{'family':'serif','sans-serif':['Times'],'size':12})
m.animate(figsize=(23,12))

or a movie (may takes some minutes to compute):

In [ ]:

%%time
rc('font',**{'family':'serif','sans-serif':['Times'],'size':12})
m.movie(figsize=(23.5,12), anim_kwargs={'interval': 100, 'frames':2000})

In [ ]: