QGS model: MAOOAM with non-linear longwave radiation terms¶

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 and quartic temperature tendencies (Non-linear longwave radiation terms in $T^4$).

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}, T4=True)

# 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 from an initial condition on the attractor, that was obtained after a long transient integration time

In [ ]:

ic = np.array([
    4.64467907e-02, -2.25054437e-02, -3.51762104e-02,  4.65057475e-03,
   -5.08842464e-03,  1.92789965e-02,  3.47849539e-03,  1.86038254e-02,
   -2.99742925e-03,  8.85862363e-03,  1.53577181e+00,  4.29584574e-02,
   -1.88143295e-02, -1.15114614e-02,  1.24843108e-03, -3.61467108e-03,
    2.81459020e-03,  5.59630451e-03,  3.56804517e-03, -1.52654934e-03,
    3.09046486e-03,  1.18185262e-06,  6.76756261e-04,  4.21981618e-06,
    1.06290914e-05, -1.64056994e-05,  3.40932989e-05,  1.90943692e-05,
    1.08189600e-05,  3.16661484e+00, -3.06749318e-04,  1.82036792e-01,
   -1.33924039e-04,  9.77967977e-03, -5.58853689e-03,  2.33218135e-02,
    3.45971396e-05, -5.88894363e-05
])

Now integrate to obtain a trajectory on the attractor (this will take ~10 minutes)

In [ ]:

%%time
integrator.integrate(0., 1000000., 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 [ ]:

model_parameters.oblocks

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})