Thursday March 31, 2016
In [39]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
import climlab
climlab to implement the two-layer leaky greenhouse modelOne of the things that climlab is set up to do is the grey-radiation modeling we have already been discussing.
Since we already derived a complete analytical solution to the two-layer leaky greenhouse model, we will use this to validate the climlab code.
We want to verify that the model reproduces the observed OLR given observed temperatures, and the absorptivity that we tuned in the analytical model. The target numbers are:
\begin{align} T_s &= 288 \text{ K} \\ T_0 &= 275 \text{ K} \\ T_1 &= 230 \text{ K} \\ \end{align}$$ \epsilon = 0.58377 $$$$ OLR = 239 \text{ W m}^{-2} $$
In [50]:
# Test in a 2-layer atmosphere
col = climlab.GreyRadiationModel(num_lev=2)
print col
In [51]:
col.state
Out[51]:
In [52]:
col.Ts
Out[52]:
In [53]:
# Set the temperatures to our observed values
col.Ts[:] = 288.
col.Tatm[:] = np.array([230., 275.])
col.state
Out[53]:
In [54]:
LW = col.subprocess['LW']
print LW
In [55]:
LW.absorptivity
Out[55]:
In [56]:
# Set the absorptivity to our tuned value
LW.absorptivity = 0.58377
LW.absorptivity
Out[56]:
In [57]:
# This does all the calculations that would be performed at each time step,
# but doesn't actually update the temperatures
col.compute_diagnostics()
# Let's see what's in the diagnostics dictionary
col.diagnostics
Out[57]:
In [58]:
# Check OLR against our analytical solution
col.OLR
Out[58]:
In [59]:
col.state
Out[59]:
In [60]:
# perform a single time step
col.step_forward()
In [61]:
col.state
Out[61]:
In [62]:
# integrate out to radiative equilibrium
col.integrate_years(2.)
In [63]:
# Check for equilibrium
col.ASR - col.OLR
Out[63]:
In [64]:
# The temperatures at radiative equilibrium
col.state
Out[64]:
Compare these to the analytical solutions for radiative equilibrium with $\epsilon = 0.58$:
\begin{align} T_1 &= 234 \text{ K} \\ T_0 &= 262 \text{ K} \\ T_s &= 296 \text{ K} \\ \end{align}So it looks like climlab agrees with our analytical results. That's good.
In [65]:
# This will try to read the data over the internet.
ncep_filename = 'air.mon.1981-2010.ltm.nc'
#ncep_url = "http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/"
#ncep_air = nc.Dataset( ncep_url + 'pressure/' + ncep_filename )
# Or to read from local disk
ncep_air = nc.Dataset( ncep_filename )
level = ncep_air.variables['level'][:]
lat = ncep_air.variables['lat'][:]
# A log-pressure height coordinate
zstar = -np.log(level/1000)
In [66]:
# Take averages of the temperature data
Tzon = np.mean(ncep_air.variables['air'][:],axis=(0,3))
Tglobal = np.average( Tzon , weights=np.cos(np.deg2rad(lat)), axis=1) + climlab.constants.tempCtoK
# Note the useful conversion factor. climlab.constants has lots of commonly used constant pre-defined
In [67]:
# Here we are plotting with respect to log(pressure) but labeling the axis in pressure units
fig = plt.figure( figsize=(8,6) )
ax = fig.add_subplot(111)
ax.plot( Tglobal , zstar )
yticks = np.array([1000., 750., 500., 250., 100., 50., 20., 10.])
ax.set_yticks(-np.log(yticks/1000.))
ax.set_yticklabels(yticks)
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_title('Global, annual mean sounding from NCEP Reanalysis', fontsize = 24)
ax.grid()
In [68]:
# initialize a grey radiation model with 30 levels
col = climlab.GreyRadiationModel()
print col
In [69]:
# interpolate to 30 evenly spaced pressure levels
lev = col.lev
Tinterp = np.interp(lev, np.flipud(level), np.flipud(Tglobal))
print Tinterp
# Need to 'flipud' because the interpolation routine needs the pressure data to be in increasing order
In [70]:
# Initialize model with observed temperatures
col.Ts[:] = Tglobal[0]
col.Tatm[:] = Tinterp
In [71]:
# A handy re-usable routine for making a plot of the temperature profiles
# We will plot temperatures with respect to log(pressure) to get a height-like coordinate
def plot_sounding(collist):
color_cycle=['r', 'g', 'b', 'y']
# col is either a column model object or a list of column model objects
if isinstance(collist, climlab.Process):
# make a list with a single item
collist = [collist]
fig = plt.figure()
ax = fig.add_subplot(111)
for i, col in enumerate(collist):
zstar = -np.log(col.lev/climlab.constants.ps)
ax.plot(col.Tatm, zstar, color=color_cycle[i])
ax.plot(col.Ts, 0, 'o', markersize=12, color=color_cycle[i])
#ax.invert_yaxis()
yticks = np.array([1000., 750., 500., 250., 100., 50., 20., 10.])
ax.set_yticks(-np.log(yticks/1000.))
ax.set_yticklabels(yticks)
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('Pressure (hPa)')
ax.grid()
return ax
In [72]:
# This should look just like the observations
plot_sounding(col)
Out[72]:
In [73]:
col.compute_diagnostics()
col.OLR
Out[73]:
In [74]:
# Need to tune absorptivity to get OLR = 239
epsarray = np.linspace(0.01, 0.1, 100)
OLRarray = np.zeros_like(epsarray)
In [75]:
for i in range(epsarray.size):
col.subprocess['LW'].absorptivity = epsarray[i]
col.compute_diagnostics()
OLRarray[i] = col.OLR
plt.plot(epsarray, OLRarray)
plt.grid()
The necessary value seems to lie near 0.055 or so.
A precise numerical search gives $\epsilon = 0.0534$
In [76]:
eps = 0.0534
col.subprocess.LW.absorptivity = eps
col.subprocess.LW.absorptivity
Out[76]:
In [77]:
# Double check to make sure this gives the right OLR
col.compute_diagnostics()
col.OLR
Out[77]:
In [81]:
col.ASR - col.OLR
Out[81]:
In [82]:
# Make a clone of our first model
re = climlab.process_like(col)
In [83]:
# To get to equilibrium, we just time-step the model forward long enough
re.integrate_years(2.)
In [84]:
# Check for energy balance
re.ASR - re.OLR
Out[84]:
Make a plot that compares the observed temperatures to the radiative equilibrium temperature profile.
In [85]:
plot_sounding([col, re])
Out[85]:
Some properties of the radiative equilibrium temperature profile:
We recognize that the large drop in temperature just above the surface is unphysical. Parcels of air in direct contact with the ground will be warmed by mechansisms other than radiative transfer.
These warm air parcels will then become buoyant, and will convect upward, mixing their heat content with the environment.
We parameterize the statistical effects of this mixing through a convective adjustment.
At each timestep, our model checks for any locations at which the lapse rate exceeds some threshold. Unstable layers are removed through an energy-conserving mixing formula.
This process is assumed to be fast relative to radiative heating. In the model, it is instantaneous.
In [86]:
rce = climlab.RadiativeConvectiveModel(adj_lapse_rate=6.)
print rce
This model is exactly like our previous models, except for one additional subprocess called convective adjustment.
We passed a parameter adj_lapse_rate (in K / km) that sets the neutrally stable lapse rate -- in this case, 6 K / km.
This number is chosed to very loosely represent the net effect of moist convection.
In [87]:
# Set our tuned absorptivity value
rce.subprocess.LW.absorptivity = eps
In [88]:
# Run out to equilibrium
rce.integrate_years(2.)
In [89]:
# Check for energy balance
rce.ASR - rce.OLR
Out[89]:
In [90]:
# Make a plot to compare observations, Radiative Equilibrium, and Radiative-Convective Equilibrium
plot_sounding([col, re, rce])
Out[90]:
Introducing convective adjustment into the model cools the surface quite a bit (compared to Radiative Equilibrium, in green here) -- and warms the lower troposphere. It gives us a MUCH better fit to observations.
But of course we still have no stratosphere.
The missing ingredient is absorption of shortwave UV radiation by ozone.
In [91]:
# Here we plotting the heating rates in K / day
degrees_per_day_atm = {'LW': (rce.subprocess['LW'].heating_rate['Tatm'] /
rce.Tatm.domain.heat_capacity *
climlab.constants.seconds_per_day),
'SW': (rce.subprocess['SW'].heating_rate['Tatm'] /
rce.Tatm.domain.heat_capacity *
climlab.constants.seconds_per_day),
'Convection': (rce.subprocess['convective adjustment'].adjustment['Tatm'] /
rce.timestep *
climlab.constants.seconds_per_day)}
degrees_per_day_sfc = {'LW': (rce.subprocess['LW'].heating_rate['Ts'] /
rce.Ts.domain.heat_capacity *
climlab.constants.seconds_per_day),
'SW': (rce.subprocess['SW'].heating_rate['Ts'] /
rce.Ts.domain.heat_capacity *
climlab.constants.seconds_per_day),
'Convection': (rce.subprocess['convective adjustment'].adjustment['Ts'] /
rce.timestep *
climlab.constants.seconds_per_day)}
In [92]:
color_cycle=['r', 'g', 'b', 'y']
fig = plt.figure( figsize=(8,6) )
ax = fig.add_subplot(111)
zstar = -np.log(rce.lev/climlab.constants.ps)
for i, item in enumerate(degrees_per_day_atm):
ax.plot(degrees_per_day_atm[item], zstar, color=color_cycle[i], label=item)
for i, item in enumerate(degrees_per_day_sfc):
ax.plot(degrees_per_day_sfc[item], 0, 'o', markersize=12, color=color_cycle[i])
yticks = np.array([1000., 750., 500., 250., 100., 50., 20., 10.])
ax.set_yticks(-np.log(yticks/1000.))
ax.set_yticklabels(yticks)
ax.set_xlabel('Heating rates (degrees per day)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.legend()
ax.grid()
What do we see here?
Collectively the LW and SW radiation are trying to push the temperatures towards radiative equilibrium while the convection is moving heat from the surface to the troposphere.
In [ ]: