In this document an Energy Balance Model (EBM) is set up with the Outgoing Longwave Radiation (OLR) parameterized through the Stefan Boltzmann radiation of a grey body.
In [1]:
# import header
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import climlab
from climlab import constants as const
from climlab.domain.field import global_mean
An EBM model instance is created through
In [2]:
# model creation
ebm_boltz = climlab.EBM()
The model is set up by default with a linearized OLR parameterization (A+BT).
In [3]:
# print model states and suprocesses
print ebm_boltz
The creation of a subprocess needs some information from the model, especially on which model state the subprocess should be defined on.
In [4]:
# create Boltzmann subprocess
LW_boltz = climlab.radiation.Boltzmann(eps=0.69, tau=0.95,
state=ebm_boltz.state,
**ebm_boltz.param)
Note that the model's whole state dictionary is given as input to the subprocess. In case only the temperature field ebm_boltz.state['Ts']
would be given, a new state dictionary would be created which holds the surface temperature with the key 'default'
. That raises an error as the Boltzmann process refers the temperature with key 'Ts'
.
Now the new OLR subprocess has to be merged into the model. Therefore the AplusBT
subprocess has to be removed first.
In [5]:
# remove the old longwave subprocess
ebm_boltz.remove_subprocess('LW')
# add the new longwave subprocess
ebm_boltz.add_subprocess('LW',LW_boltz)
Note that the new OLR subprocess has to have the same key 'LW'
as the old one, as the model refers to this key for radiation balance computation.
That is why the old process has to be removed before the new one is added.
In [6]:
print ebm_boltz
In [7]:
# integrate model for a single timestep
ebm_boltz.step_forward()
In [8]:
# creating plot figure
fig = plt.figure(figsize=(15,10))
# Temperature plot
ax1 = fig.add_subplot(221)
ax1.plot(ebm_boltz.lat,ebm_boltz.Ts)
ax1.set_xticks([-90,-60,-30,0,30,60,90])
ax1.set_xlim([-90,90])
ax1.set_title('Surface Temperature', fontsize=14)
ax1.set_ylabel('(degC)', fontsize=12)
ax1.grid()
# Albedo plot
ax2 = fig.add_subplot(223, sharex = ax1)
ax2.plot(ebm_boltz.lat,ebm_boltz.albedo)
ax2.set_title('Albedo', fontsize=14)
ax2.set_xlabel('latitude', fontsize=10)
ax2.set_ylim([0,1])
ax2.grid()
# Net Radiation plot
ax3 = fig.add_subplot(222, sharex = ax1)
ax3.plot(ebm_boltz.lat, ebm_boltz.OLR, label='OLR',
color='cyan')
ax3.plot(ebm_boltz.lat, ebm_boltz.ASR, label='ASR',
color='magenta')
ax3.plot(ebm_boltz.lat, ebm_boltz.ASR-ebm_boltz.OLR,
label='net radiation',
color='red')
ax3.set_title('Net Radiation', fontsize=14)
ax3.set_ylabel('(W/m$^2$)', fontsize=12)
ax3.legend(loc='best')
ax3.grid()
# Energy Balance plot
net_rad = np.squeeze(ebm_boltz.net_radiation)
transport = ebm_boltz.heat_transport_convergence()
ax4 = fig.add_subplot(224, sharex = ax1)
ax4.plot(ebm_boltz.lat, net_rad, label='net radiation',
color='red')
ax4.plot(ebm_boltz.lat, transport, label='diffusion transport',
color='blue')
ax4.plot(ebm_boltz.lat, net_rad+transport, label='balance',
color='black')
ax4.set_title('Energy', fontsize=14)
ax4.set_xlabel('latitude', fontsize=10)
ax4.set_ylabel('(W/m$^2$)', fontsize=12)
ax4.legend(loc='best')
ax4.grid()
plt.show()
In [9]:
# integrate model until solution converges
ebm_boltz.integrate_converge()
In [10]:
# creating plot figure
fig = plt.figure(figsize=(15,10))
# Temperature plot
ax1 = fig.add_subplot(221)
ax1.plot(ebm_boltz.lat,ebm_boltz.Ts)
ax1.set_xticks([-90,-60,-30,0,30,60,90])
ax1.set_xlim([-90,90])
ax1.set_title('Surface Temperature', fontsize=14)
ax1.set_ylabel('(degC)', fontsize=12)
ax1.grid()
# Albedo plot
ax2 = fig.add_subplot(223, sharex = ax1)
ax2.plot(ebm_boltz.lat,ebm_boltz.albedo)
ax2.set_title('Albedo', fontsize=14)
ax2.set_xlabel('latitude', fontsize=10)
ax2.set_ylim([0,1])
ax2.grid()
# Net Radiation plot
ax3 = fig.add_subplot(222, sharex = ax1)
ax3.plot(ebm_boltz.lat, ebm_boltz.OLR, label='OLR',
color='cyan')
ax3.plot(ebm_boltz.lat, ebm_boltz.ASR, label='ASR',
color='magenta')
ax3.plot(ebm_boltz.lat, ebm_boltz.ASR-ebm_boltz.OLR,
label='net radiation',
color='red')
ax3.set_title('Net Radiation', fontsize=14)
ax3.set_ylabel('(W/m$^2$)', fontsize=12)
ax3.legend(loc='best')
ax3.grid()
# Energy Balance plot
net_rad = np.squeeze(ebm_boltz.net_radiation)
transport = ebm_boltz.heat_transport_convergence()
ax4 = fig.add_subplot(224, sharex = ax1)
ax4.plot(ebm_boltz.lat, net_rad, label='net radiation',
color='red')
ax4.plot(ebm_boltz.lat, transport, label='diffusion transport',
color='blue')
ax4.plot(ebm_boltz.lat, net_rad+transport, label='balance',
color='black')
ax4.set_title('Energy', fontsize=14)
ax4.set_xlabel('latitude', fontsize=10)
ax4.set_ylabel('(W/m$^2$)', fontsize=12)
ax4.legend(loc='best')
ax4.grid()
plt.show()
In [11]:
# creating plot figure
fig = plt.figure(figsize=(15,10))
# Temperature plot
ax1 = fig.add_subplot(221)
ax1.plot(ebm_boltz.lat,ebm_boltz.Ts)
ax1.set_xticks([-90,-60,-30,0,30,60,90])
ax1.set_xlim([-90,90])
ax1.set_xlabel('latitude')
ax1.set_ylabel('surface temperature (degC)', fontsize=12)
ax1.grid()
# Albedo plot
ax2 = fig.add_subplot(223, sharex = ax1)
ax2.plot(ebm_boltz.lat,ebm_boltz.albedo)
ax2.set_ylabel('albedo', fontsize=12)
ax2.set_ylim([0,1])
ax2.grid()
# Net Radiation plot
ax3 = fig.add_subplot(222, sharex = ax1)
ax3.plot(ebm_boltz.lat, ebm_boltz.OLR, label='OLR',
color='cyan')
ax3.plot(ebm_boltz.lat, ebm_boltz.ASR, label='ASR',
color='magenta')
ax3.plot(ebm_boltz.lat, ebm_boltz.ASR-ebm_boltz.OLR,
label='net radiation',
color='red')
ax3.set_ylabel('net radiation (W/m$^2$)', fontsize=12)
ax3.legend(loc='best')
ax3.grid()
# Energy Balance plot
net_rad = np.squeeze(ebm_boltz.net_radiation)
transport = ebm_boltz.heat_transport_convergence()
ax4 = fig.add_subplot(224, sharex = ax1)
ax4.plot(ebm_boltz.lat, net_rad, label='net radiation',
color='red')
ax4.plot(ebm_boltz.lat, transport, label='diffusion transport',
color='blue')
ax4.plot(ebm_boltz.lat, net_rad+transport, label='balance',
color='black')
ax4.set_ylabel('energy (W/m$^2$)', fontsize=12)
ax4.legend(loc='best')
ax4.grid()
plt.show()
In [12]:
print 'The global mean temperature is %s degC with a model
ice edge at %s deg.' \
%(np.round(global_mean(ebm_boltz.Ts),2), np.max(ebm_boltz.icelat))