The equation for the model is
$$ C \frac{\partial T}{\partial t} = ASR - OLR + \frac{D}{\cos\phi} \frac{\partial}{\partial \phi} \left( \cos\phi \frac{\partial T}{\partial \phi} \right)$$Where, following our standard notation:
In [ ]:
# We start with the usual import statements
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import climlab
In [ ]:
# create a new model with all default parameters (except the grid size)
mymodel = climlab.EBM_annual(num_lat = 30)
In [ ]:
# What did we just do?
print mymodel
Every model in climlab consists of some state variables, and processes that operate on those state variables.
The state variables are the climate variables that we will step forward in time. For our EBM, it is just the surface temperature:
In [ ]:
mymodel.Ts
We have an array of temperatures in degrees Celsius. Let's see how big this array is:
In [ ]:
mymodel.Ts.shape
Every state variable exists on a spatial grid. In this case, the grid is an array of latitude points:
In [ ]:
mymodel.lat
We can, for example, plot the current temperature versus latitude:
In [ ]:
plt.plot(mymodel.lat, mymodel.Ts)
plt.xlabel('Latitude')
plt.ylabel('Temperature (deg C)')
It is based on a very general concept of a model as a collection of individual, interacting processes.
Every model in climlab is a collection of individual, interacting physical processes.
In this case, the EBM_annual object has a list of 4 sub-processes. Each one is basically one of the terms in our energy budget equation above:
To see this explicitly:
In [ ]:
# The dictionary of sub-processes:
mymodel.subprocess
In [ ]:
# Make a copy of the original temperature array
initial = mymodel.Ts.copy()
# Take a single timestep forward!
mymodel.step_forward()
# Check out the difference
print mymodel.Ts - initial
Looks like the temperature got a bit colder near the equator and warmer near the poles
In [ ]:
# How long is a single timestep?
mymodel.timestep
This value is in seconds. It is actually 1/90th of a year (so, to step forward one year, we need 90 individual steps). This is a default value -- we could change it if we wanted to.
Actually we can see a bunch of important parameters for this model all together in a dictionary called param:
In [ ]:
mymodel.param
In [ ]:
# Each climlab model has a dictionary called diagnostics.
# Let's look at the names of the fields in this dicionary:
mymodel.diagnostics.keys()
In [ ]:
# We can access individual fields either through standard dictionary notation:
mymodel.diagnostics['ASR']
In [ ]:
# Or using the more interactive 'dot' notation:
mymodel.ASR
In [ ]:
# Let's use the diagnostics to make a plot of the current state of the radiation
plt.plot(mymodel.lat, mymodel.ASR, label='ASR')
plt.plot(mymodel.lat, mymodel.OLR, label='OLR')
plt.xlabel('Latitude')
plt.ylabel('W/m2')
plt.legend()
plt.grid()
This plot shows that $ASR > OLR$ (system is gaining extra energy) across the tropics, and $ASR < OLR$ (system is losing energy) near the poles.
It's interesting to ask how close this model is to global energy balance.
For that, we have to take an area-weighted global average.
In [ ]:
climlab.global_mean(mymodel.net_radiation)
In [ ]:
climlab.global_mean(mymodel.Ts)
In [ ]:
# Loop 90 times for 1 year of simulation
for n in range(90):
mymodel.step_forward()
print 'Global mean temperature is %0.1f degrees C.' %climlab.global_mean(mymodel.Ts)
print 'Global mean energy imbalance is %0.1f W/m2.' %climlab.global_mean(mymodel.net_radiation)
Since there is still a significant imbalance, we are not yet at equilibrium. We should step forward again.
This time, let's use a convenient built-in function instead of writing our own loop:
In [ ]:
mymodel.integrate_years(1.)
print 'Global mean temperature is %0.1f degrees C.' %climlab.global_mean(mymodel.Ts)
print 'Global mean energy imbalance is %0.1f W/m2.' %climlab.global_mean(mymodel.net_radiation)
We are now quite close to equilibrium. Let's make some plots
In [ ]:
plt.plot(mymodel.lat, mymodel.Ts)
plt.xlabel('Latitude')
plt.ylabel('Temperature (deg C)')
plt.grid()
In [ ]:
plt.plot(mymodel.lat_bounds, mymodel.heat_transport())
plt.xlabel('Latitude')
plt.ylabel('PW')
plt.grid()
Let's make some nice plots of all the terms in the energy budget.
We'll define a reusable function to make the plots.
In [ ]:
def ebm_plot(e, return_fig=False):
templimits = -20,32
radlimits = -340, 340
htlimits = -6,6
latlimits = -90,90
lat_ticks = np.arange(-90,90,30)
fig = plt.figure(figsize=(8,12))
ax1 = fig.add_subplot(3,1,1)
ax1.plot(e.lat, e.Ts)
ax1.set_ylim(templimits)
ax1.set_ylabel('Temperature (deg C)')
ax2 = fig.add_subplot(3,1,2)
ax2.plot(e.lat, e.ASR, 'k--', label='SW' )
ax2.plot(e.lat, -e.OLR, 'r--', label='LW' )
ax2.plot(e.lat, e.net_radiation, 'c-', label='net rad' )
ax2.plot(e.lat, e.heat_transport_convergence(), 'g--', label='dyn' )
ax2.plot(e.lat, e.net_radiation.squeeze() + e.heat_transport_convergence(), 'b-', label='total' )
ax2.set_ylim(radlimits)
ax2.set_ylabel('Energy budget (W m$^{-2}$)')
ax2.legend()
ax3 = fig.add_subplot(3,1,3)
ax3.plot(e.lat_bounds, e.heat_transport() )
ax3.set_ylim(htlimits)
ax3.set_ylabel('Heat transport (PW)')
for ax in [ax1, ax2, ax3]:
ax.set_xlabel('Latitude')
ax.set_xlim(latlimits)
ax.set_xticks(lat_ticks)
ax.grid()
if return_fig:
return fig
In [ ]:
ebm_plot(mymodel)
In [ ]:
model2 = climlab.EBM_annual(num_lat=30)
print model2
In [ ]:
# The default initial temperature
print model2.Ts
In [ ]:
# Now let's change it to be 15 degrees everywhere
model2.Ts[:] = 15.
print model2.Ts
In [ ]:
model2.compute_diagnostics()
ebm_plot(model2)
Why is the heat transport zero everywhere?
Now let's look at an animation showing how all the terms in the energy budget evolve from this initial condition toward equilibrium.
The code below generates a series of snapshots, that we can stitch together into an animation with a graphics program.
In [ ]:
# Code to generate individual frames for the animation
# You need to specify a valid path for your computer, or else this won't work
# Uncomment the code below to re-create the frames from the animation
#nsteps = 90
#for n in range(nsteps):
# fig = ebm_plot(model2, return_fig=True)
#
# filepath = '/Users/br546577/Desktop/temp/ebm_animation_frames/'
# filename = filepath + 'ebm_animation' + str(n).zfill(4) + '.png'
# fig.savefig(filename)
# plt.close()
#
# model2.step_forward()
We will look at the animation together in class.
In [ ]: