Here is a quick example of using the climlab.GreyRadiationModel with a latitude dimension and seasonally varying insolation.
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import climlab
In [2]:
np.__version__
Out[2]:
In [3]:
model = climlab.GreyRadiationModel(num_lev=30, num_lat=90)
print model
In [4]:
print model.lat
In [5]:
insolation = climlab.radiation.insolation.DailyInsolation(domains=model.Ts.domain)
In [6]:
model.add_subprocess('insolation', insolation)
model.subprocess.SW.flux_from_space = insolation.insolation
In [7]:
print model
In [8]:
model.step_forward()
In [9]:
plt.plot(model.lat, model.SW_down_TOA)
Out[9]:
In [10]:
model.Tatm.shape
Out[10]:
In [11]:
model.integrate_years(1)
In [12]:
plt.plot(model.lat, model.Ts)
Out[12]:
In [13]:
model.integrate_years(1)
In [14]:
plt.plot(model.lat, model.timeave['Ts'])
Out[14]:
In [15]:
def plot_temp_section(model, timeave=True):
fig = plt.figure()
ax = fig.add_subplot(111)
if timeave:
field = model.timeave['Tatm'].transpose()
else:
field = model.Tatm.transpose()
cax = ax.contourf(model.lat, model.lev, field)
ax.invert_yaxis()
ax.set_xlim(-90,90)
ax.set_xticks([-90, -60, -30, 0, 30, 60, 90])
fig.colorbar(cax)
In [16]:
plot_temp_section(model)
In [17]:
model2 = climlab.RadiativeConvectiveModel(num_lev=30, num_lat=90)
insolation = climlab.radiation.insolation.DailyInsolation(domains=model2.Ts.domain)
model2.add_subprocess('insolation', insolation)
model2.subprocess.SW.flux_from_space = insolation.insolation
In [18]:
model2.step_forward()
In [19]:
model2.integrate_years(1)
In [20]:
model2.integrate_years(1)
In [21]:
plot_temp_section(model2)
In [22]:
import climlab
import numpy as np
# Put in some ozone
import netCDF4 as nc
datapath = "http://ramadda.atmos.albany.edu:8080/repository/opendap/latest/Top/Users/Brian+Rose/CESM+runs/"
endstr = "/entry.das"
topo = nc.Dataset( datapath + 'som_input/USGS-gtopo30_1.9x2.5_remap_c050602.nc' + endstr )
ozone = nc.Dataset( datapath + 'som_input/ozone_1.9x2.5_L26_2000clim_c091112.nc' + endstr )
# Dimensions of the ozone file
lat = ozone.variables['lat'][:]
lon = ozone.variables['lon'][:]
lev = ozone.variables['lev'][:]
# Taking annual, zonal average of the ozone data
O3_zon = np.mean( ozone.variables['O3'],axis=(0,3) )
In [23]:
import climlab
import numpy as np
# make a model on the same grid as the ozone
model3 = climlab.BandRCModel(lev=lev, lat=lat)
insolation = climlab.radiation.insolation.DailyInsolation(domains=model3.Ts.domain)
model3.add_subprocess('insolation', insolation)
model3.subprocess.SW.flux_from_space = insolation.insolation
print model3
In [24]:
# Set the ozone mixing ratio
O3_trans = np.transpose(O3_zon)
# model and ozone data are on the same grid, after the transpose.
print O3_trans.shape
print lev
print model3.lev
In [25]:
# Put in the ozone
model3.absorber_vmr['O3'] = O3_trans
In [26]:
print model3.absorber_vmr['O3'].shape
print model3.Tatm.shape
In [27]:
model3.step_forward()
In [28]:
model3.integrate_years(1.)
In [29]:
model3.integrate_years(1.)
In [30]:
#plt.contour(model3.lat, model3.lev, model3.Tatm.transpose())
plot_temp_section(model3)
This is now working. Will need to do some model tuning.
And start to add dynamics!
In [31]:
testmodel = climlab.BandRCModel(num_lat=90, num_lev=30)
In [32]:
testmodel.step_forward()
testmodel.step_forward()
testmodel.step_forward()
In [33]:
%timeit testmodel.step_forward()
In [34]:
np.__version__
Out[34]:
Definitely get better performance with numpy version 1.9 and above due to vectorization of numpy.tril() (lower triangle operator) for multi-dimensional arrays.
In [35]:
print model2
In [36]:
diffmodel = climlab.process_like(model2)
In [37]:
# thermal diffusivity in W/m**2/degC
D = 0.05
# meridional diffusivity in 1/s
K = D / diffmodel.Tatm.domain.heat_capacity[0]
print K
In [38]:
d = climlab.dynamics.diffusion.MeridionalDiffusion(K=K, state={'Tatm': diffmodel.state['Tatm']}, **diffmodel.param)
In [39]:
diffmodel.add_subprocess('diffusion', d)
In [40]:
print diffmodel
In [41]:
diffmodel.step_forward()
In [42]:
diffmodel.integrate_years(1)
In [43]:
diffmodel.integrate_years(1)
In [44]:
plot_temp_section(model2)
plot_temp_section(diffmodel)
Amazingly... this actually works!
As long as K is a constant.
The diffusion operation is broadcast over all vertical levels without any special code.
In [45]:
def inferred_heat_transport( energy_in, lat_deg ):
'''Returns the inferred heat transport (in PW) by integrating the net energy imbalance from pole to pole.'''
from scipy import integrate
from climlab import constants as const
lat_rad = np.deg2rad( lat_deg )
return ( 1E-15 * 2 * np.math.pi * const.a**2 * integrate.cumtrapz( np.cos(lat_rad)*energy_in,
x=lat_rad, initial=0. ) )
In [46]:
# Plot the northward heat transport in this model
Rtoa = np.squeeze(diffmodel.timeave['ASR'] - diffmodel.timeave['OLR'])
plt.plot(diffmodel.lat, inferred_heat_transport(Rtoa, diffmodel.lat))
Out[46]:
In [47]:
diffband = climlab.process_like(model3)
In [48]:
# thermal diffusivity in W/m**2/degC
D = 0.05
# meridional diffusivity in 1/s
K = D / diffband.Tatm.domain.heat_capacity[0]
print K
In [49]:
d = climlab.dynamics.diffusion.MeridionalDiffusion(K=K, state={'Tatm': diffband.state['Tatm']}, **diffband.param)
diffband.add_subprocess('diffusion', d)
print diffband
In [50]:
diffband.integrate_years(1)
In [51]:
diffband.integrate_years(1)
In [52]:
plot_temp_section(model3)
plot_temp_section(diffband)
In [53]:
plt.plot(diffband.lat, diffband.timeave['ASR'] - diffband.timeave['OLR'])
Out[53]:
In [54]:
# Plot the northward heat transport in this model
Rtoa = np.squeeze(diffband.timeave['ASR'] - diffband.timeave['OLR'])
plt.plot(diffband.lat, inferred_heat_transport(Rtoa, diffband.lat))
Out[54]:
In [ ]: