<jupyter_start><jupyter_text>ENV / ATM 415: Climate Laboratory Radiative- and Radiative-Convective Equilibrium with `climlab` Thursday March 31, 2016<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc
import climlab<jupyter_output><empty_output><jupyter_text>____________ Using `climlab` to implement the two-layer leaky greenhouse model____________One 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.  ValidationWe 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} $$  Initialize a model in `climlab`The first thing we do is create a new model.The following example code is sparsely commented but will hopefully orient you on the basics of defining and working with a `climlab Process` object.<jupyter_code>#  Test in a 2-layer atmosphere
col = climlab.GreyRadiationModel(num_lev=2)
print col
col.state
col.Ts
#  Set the temperatures to our observed values
col.Ts[:] = 288.
col.Tatm[:] = np.array([230., 275.])
col.state
LW = col.subprocess['LW']
print LW
LW.absorptivity
#  Set the absorptivity to our tuned value
LW.absorptivity = 0.58377
LW.absorptivity
#  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
#  Check OLR against our analytical solution
col.OLR
col.state
# perform a single time step
col.step_forward()
col.state
# integrate out to radiative equilibrium
col.integrate_years(2.)
# Check for equilibrium
col.ASR - col.OLR
#  The temperatures at radiative equilibrium
col.state<jupyter_output><empty_output><jupyter_text>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. ____________ The observed annual, global mean temperature profile____________We want to model the OLR in a column whose temperatures match observations. We'll calculate the global, annual mean air temperature from the NCEP Reanalysis data.<jupyter_code># 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)
#  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
#  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()<jupyter_output><empty_output><jupyter_text>____________ A 30-layer model using the observed temperatures____________<jupyter_code>#  initialize a grey radiation model with 30 levels
col = climlab.GreyRadiationModel()
print col
# 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
# Initialize model with observed temperatures
col.Ts[:] = Tglobal[0]
col.Tatm[:] = Tinterp
#  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
# This should look just like the observations
plot_sounding(col)<jupyter_output><empty_output><jupyter_text>Tune absorptivity to get observed OLR<jupyter_code>col.compute_diagnostics()
col.OLR
# Need to tune absorptivity to get OLR = 239
epsarray = np.linspace(0.01, 0.1, 100)
OLRarray = np.zeros_like(epsarray)
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()<jupyter_output><empty_output><jupyter_text>The necessary value seems to lie near 0.055 or so.A precise numerical search gives $\epsilon = 0.0534$<jupyter_code>eps = 0.0534

col.subprocess.LW.absorptivity = eps
col.subprocess.LW.absorptivity
#  Double check to make sure this gives the right OLR

col.compute_diagnostics()
col.OLR
col.ASR - col.OLR<jupyter_output><empty_output><jupyter_text>____________ Radiative equilibrium in the 30-layer model____________<jupyter_code>#  Make a clone of our first model
re = climlab.process_like(col)
#  To get to equilibrium, we just time-step the model forward long enough
re.integrate_years(2.)
#  Check for energy balance
re.ASR - re.OLR<jupyter_output><empty_output><jupyter_text>Make a plot that compares the **observed temperatures** to the **radiative equilibrium** temperature profile.<jupyter_code>plot_sounding([col, re])<jupyter_output><empty_output><jupyter_text>Some properties of the **radiative equilibrium** temperature profile:- The surface is warmer than observed.- The lower troposphere is colder than observed.- Very cold air is sitting immediately above the warm surface.- There is no tropopause, no stratosphere. ____________ Radiative-Convective Equilibrium in the 30-layer model____________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.<jupyter_code>rce = climlab.RadiativeConvectiveModel(adj_lapse_rate=6.)
print rce<jupyter_output>climlab Process of type <class 'climlab.model.column.RadiativeConvectiveModel'>. 
State variables and domain shapes: 
  Tatm: (30,) 
  Ts: (1,) 
The subprocess tree: 
top: <class 'climlab.model.column.RadiativeConvectiveModel'>
   convective adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>
   LW: <class 'climlab.radiation.greygas.GreyGas'>
   SW: <class 'climlab.radiation.greygas.GreyGasSW'>
   insolation: <class 'climlab.radiation.insolation.FixedInsolation'><jupyter_text>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**.<jupyter_code># Set our tuned absorptivity value
rce.subprocess.LW.absorptivity = eps
#  Run out to equilibrium
rce.integrate_years(2.)
#  Check for energy balance
rce.ASR - rce.OLR
#  Make a plot to compare observations, Radiative Equilibrium, and Radiative-Convective Equilibrium
plot_sounding([col, re, rce])<jupyter_output><empty_output><jupyter_text>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**.  The heating rates due to radiation and convection<jupyter_code>#  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)}
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()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10154964938100597, 'fraction_numerical': 0.019825123365942342, 'mean_word_length': 4.035745422842197, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '214795', 'license': 'MIT', 'n_tokens_mistral': 3878, 'n_tokens_neox': 3555, 'n_words': 1127, 'path': 'notes/RadiativeConvectiveEquilibrium.ipynb', 'repo_name': 'brian-rose/env-415-site', 'size': '11551'}