Comparing soundings from NCEP Reanalysis and various models

We are going to plot the global, annual mean sounding (vertical temperature profile) from observations.

Read in the necessary NCEP reanalysis data from the online server.

The catalog is here: http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/catalog.html

In [3]: 
import netCDF4 as nc
ncep_url = "http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/"
ncep_air = nc.Dataset( ncep_url + "pressure/air.mon.1981-2010.ltm.nc" )
level = ncep_air.variables['level'][:]
lat = ncep_air.variables['lat'][:]

Take global averages and time averages.

In [4]: 
Tzon = mean(ncep_air.variables['air'],axis=(0,3))
Tglobal = sum( Tzon * cos(deg2rad(lat)), axis=1 ) / sum( cos(deg2rad(lat) ) )

Here is code to make a nicely labeled sounding plot.

In [5]: 
fig = figure( figsize=(10,8) )
ax = fig.add_subplot(111)
ax.plot( Tglobal + 273.15, log(level/1000) )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_yticks( log(level/1000) )
ax.set_yticklabels( level )
ax.set_title('Global, annual mean sounding from NCEP Reanalysis', fontsize = 24)
ax2 = ax.twinx()
ax2.plot( Tglobal + 273.15, -8*log(level/1000) );
ax2.set_ylabel('Approx. height above surface (km)', fontsize=16 );
ax.grid()

Now compute the Radiative Equilibrium solution for the grey-gas column model

In [6]: 
import ClimateUtils as clim
import ColumnModel
In [7]: 
col = ColumnModel.column()
print col

Instance of column class with surface temperature 288.0 K, 30 atmospheric levels, and parameters: 
{'abs_coeff': 0.0001229, 'num_levels': 30, 'Q': 341.3, 'nul': None, 'adj_lapse_rate': None, 'timestep': 86400.0, 'water_depth': 1.0, 'albedo': 0.299}
In [8]: 
col.Step_Forward(300)
print "After " + str(col.days_elapsed) + " days of integration:"
print "Surface temperature is " + str(col.Ts) + " K."
print "Net energy in to the column is " + str(col.SW_absorbed_sfc - col.OLR) + " W / m2."

After 300.0 days of integration:
Surface temperature is 287.849676171 K.
Net energy in to the column is 0.00654103077261 W / m2.

Plot the radiative equilibrium temperature on the same plot with NCEP reanalysis

In [9]: 
fig = figure( figsize=(10,8) )
ax = fig.add_subplot(111)
ax.plot( Tglobal + 273.15, log(level/1000), 'b-', col.Tatm, log( col.p/clim.ps ), 'r-' )
ax.plot( col.Ts, 0, 'ro', markersize=20 )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_yticks( log(level/1000) )
ax.set_yticklabels( level )
ax.set_title('Temperature profiles: observed (blue) and radiative equilibrium in grey gas model (red)', fontsize = 18)
ax2 = ax.twinx()
ax2.plot( Tglobal + clim.tempCtoK, -8*log(level/1000) );
ax2.set_ylabel('Approx. height above surface (km)', fontsize=16 );
ax.grid()

Now use convective adjustment to compute a Radiative-Convective Equilibrium temperature profile

In [10]: 
dalr_col = ColumnModel.column( params = {'adj_lapse_rate':'DALR'} )
print dalr_col

Instance of column class with surface temperature 288.0 K, 30 atmospheric levels, and parameters: 
{'abs_coeff': 0.0001229, 'num_levels': 30, 'Q': 341.3, 'adj_lapse_rate': 'DALR', 'timestep': 86400.0, 'water_depth': 1.0, 'albedo': 0.299}
In [11]: 
dalr_col.Step_Forward(600)
print "After " + str(dalr_col.days_elapsed) + " days of integration:"
print "Surface temperature is " + str(dalr_col.Ts) + " K."
print "Net energy in to the column is " + str(dalr_col.SW_absorbed_sfc - dalr_col.OLR) + " W / m2."

After 600.0 days of integration:
Surface temperature is 283.044704421 K.
Net energy in to the column is 1.93777987079e-05 W / m2.

Now plot this "Radiative-Convective Equilibrium" on the same graph:

In [12]: 
fig = figure( figsize=(10,8) )
ax = fig.add_subplot(111)
ax.plot( Tglobal + 273.15, log(level/1000), 'b-', col.Tatm, log( col.p/clim.ps ), 'r-' )
ax.plot( col.Ts, 0, 'ro', markersize=16 )
ax.plot( dalr_col.Tatm, log( dalr_col.p / clim.ps ), 'k-' )
ax.plot( dalr_col.Ts, 0, 'ko', markersize=16 )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_yticks( log(level/1000) )
ax.set_yticklabels( level )
ax.set_title('Temperature profiles: observed (blue), RE (red) and dry RCE (black)', fontsize = 18)
ax2 = ax.twinx()
ax2.plot( Tglobal + clim.tempCtoK, -8*log(level/1000) );
ax2.set_ylabel('Approx. height above surface (km)', fontsize=16 );
ax.grid()

The convective adjustment gets rid of the unphysical temperature difference between the surface and the overlying air.

But now the surface is colder! Convection acts to move heat upward, away from the surface.

Also, we note that the observed lapse rate (blue) is always shallower than \(\Gamma_d\) (temperatures decrease more slowly with height).

"Moist" Convective Adjustment

To approximately account for the effects of latent heat release in rising air parcels, we can just adjust to a lapse rate that is a little shallow than \(\Gamma_d\).

We will choose 6 K / km, which gets close to the observed mean lapse rate.

We will also re-tune the longwave absorptivity of the column to get a realistic surface temperature of 288 K:

In [13]: 
rce_col = ColumnModel.column( params = {'adj_lapse_rate':6, 'abs_coeff':1.7E-4})
print rce_col

Instance of column class with surface temperature 288.0 K, 30 atmospheric levels, and parameters: 
{'abs_coeff': 0.00017, 'num_levels': 30, 'Q': 341.3, 'adj_lapse_rate': 6, 'timestep': 86400.0, 'water_depth': 1.0, 'albedo': 0.299}
In [14]: 
rce_col.Step_Forward(600)
print "After " + str(rce_col.days_elapsed) + " days of integration:"
print "Surface temperature is " + str(rce_col.Ts) + " K."
print "Net energy in to the column is " + str(rce_col.SW_absorbed_sfc - rce_col.OLR) + " W / m2."

After 600.0 days of integration:
Surface temperature is 287.909683685 K.
Net energy in to the column is 4.2865378191e-05 W / m2.

Now add this new temperature profile to the graph:

In [15]: 
fig = figure( figsize=(10,8) )
ax = fig.add_subplot(111)
ax.plot( Tglobal + 273.15, log(level/1000), 'b-', col.Tatm, log( col.p/clim.ps ), 'r-' )
ax.plot( col.Ts, 0, 'ro', markersize=16 )
ax.plot( dalr_col.Tatm, log( dalr_col.p / clim.ps ), 'k-' )
ax.plot( dalr_col.Ts, 0, 'ko', markersize=16 )
ax.plot( rce_col.Tatm, log( rce_col.p / clim.ps ), 'm-' )
ax.plot( rce_col.Ts, 0, 'mo', markersize=16 )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_yticks( log(level/1000) )
ax.set_yticklabels( level )
ax.set_title('Temperature profiles: observed (blue), RE (red), dry RCE (black), and moist RCE (magenta)', fontsize = 18)
ax2 = ax.twinx()
ax2.plot( Tglobal + clim.tempCtoK, -8*log(level/1000) );
ax2.set_ylabel('Approx. height above surface (km)', fontsize=16 );
ax.grid()

Adding stratospheric ozone

Our model has no equivalent of the stratosphere, where temperature increases with height. That's because our model has been completely transparent to shortwave radiation up until now.

We can load the observed ozone climatology from the input files for the CESM model:

In [16]: 
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 )
In [17]: 
print ozone.variables['O3']

<type 'netCDF4.Variable'>
float32 O3(time, lev, lat, lon)
    units: mol/mol
    long_name: O3 concentration
    cell_method: time: mean
unlimited dimensions: time
current shape = (12, 26, 96, 144)
In [18]: 
lat = ozone.variables['lat'][:]
lon = ozone.variables['lon'][:]
lev = ozone.variables['lev'][:]

The pressure levels in this dataset are:

In [19]: 
print lev

[   3.544638     7.3888135   13.967214    23.944625    37.23029     53.114605
   70.05915     85.439115   100.514695   118.250335   139.115395   163.66207
  192.539935   226.513265   266.481155   313.501265   368.81798    433.895225
  510.455255   600.5242     696.79629    787.70206    867.16076    929.648875
  970.55483    992.5561   ]

Take the global average of the ozone climatology, and plot it as a function of pressure (or height)

In [20]: 
O3_zon = mean( ozone.variables['O3'],axis=(0,3) )
O3_global = sum( O3_zon * cos(deg2rad(lat)), axis=1 ) / sum( cos(deg2rad(lat) ) )
In [21]: 
O3_global.shape
Out[21]:
(26,)
In [22]: 
ax = figure(figsize=(10,8)).add_subplot(111)
ax.plot( O3_global * 1.E6, log(lev/clim.ps) )
ax.invert_yaxis()
ax.set_xlabel('Ozone (ppm)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
yticks = array([1000., 500., 250., 100., 50., 20., 10., 5.])
ax.set_yticks( log(yticks/1000.) )
ax.set_yticklabels( yticks )
ax.set_title('Global, annual mean ozone concentration', fontsize = 24);

This shows that most of the ozone is indeed in the stratosphere, and peaks near the top of the stratosphere.

Now create a new column model object on the same pressure levels as the ozone data. We are also going set an adjusted lapse rate of 6 K / km, and tune the longwave absorption

In [23]: 
oz_col = ColumnModel.column( p = lev, params={'abs_coeff':1.82E-4, 'adj_lapse_rate':6, 'albedo':0.315} )

Now we will do something new: let the column absorb some shortwave radiation. We will assume that the shortwave absorptivity is proportional to the ozone concentration we plotted above. We need to weight the absorptivity by the pressure (mass) of each layer.

In [24]: 
ozonefactor = 75
sw_abs = flipud(O3_global) * oz_col.dp * ozonefactor
oz_col.set_SW_absorptivity( sw_abs )
oz_col.Longwave_Heating()
oz_col.Shortwave_Heating()
In [25]: 
oz_col.SW_absorbed_atm
Out[25]:
array([ 0.01661879,  0.02963284,  0.05036985,  0.07249398,  0.09443201,
        0.11656803,  0.1290027 ,  0.12549225,  0.11582793,  0.10703984,
        0.10521417,  0.11483168,  0.13617335,  0.16175213,  0.18202097,
        0.19545583,  0.21169751,  0.24088714,  0.31126059,  0.52056369,
        0.98422166,  1.52744332,  2.0163448 ,  2.08143267,  1.49424719,
        1.42232175])

Now run it out to Radiative-Convective Equilibrium, and plot

In [26]: 
oz_col.Step_Forward(600)
print "After " + str(oz_col.days_elapsed) + " days of integration:"
print "Surface temperature is " + str(oz_col.Ts) + " K."
print "Net energy in to the column is " + str(oz_col.SW_absorbed_total - oz_col.OLR) + " W / m2."

After 600.0 days of integration:
Surface temperature is 287.918938151 K.
Net energy in to the column is 1.08835529318e-05 W / m2.
In [27]: 
oz_col.SW_absorbed_total
Out[27]:
239.76495997022022
In [28]: 
fig = figure( figsize=(10,8) )
ax = fig.add_subplot(111)
ax.plot( Tglobal + 273.15, log(level/1000), 'b-', col.Tatm, log( col.p/clim.ps ), 'r-' )
ax.plot( col.Ts, 0, 'ro', markersize=16 )
ax.plot( dalr_col.Tatm, log( dalr_col.p / clim.ps ), 'k-' )
ax.plot( dalr_col.Ts, 0, 'ko', markersize=16 )
ax.plot( rce_col.Tatm, log( rce_col.p / clim.ps ), 'm-' )
ax.plot( rce_col.Ts, 0, 'mo', markersize=16 )
ax.plot( oz_col.Tatm, log( oz_col.p / clim.ps ), 'c-' )
ax.plot( oz_col.Ts, 0, 'co', markersize=16 )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_yticks( log(level/1000) )
ax.set_yticklabels( level )
ax.set_title('Temperature profiles: observed (blue), RE (red), dry RCE (black), moist RCE (magenta), RCE with ozone (cyan)', fontsize = 18)
ax.grid()

And we finally have something that looks looks like the tropopause, with temperature increasing above at about the correct rate. Though the tropopause temperature is off by 15 degrees or so.

Greenhouse warming in the RCE model with ozone

In [29]: 
import copy
oz_col2 = copy.deepcopy( oz_col )
In [30]: 
oz_col2.params['abs_coeff'] = 1.2 * oz_col.params['abs_coeff']
oz_col2.set_LW_emissivity()
In [31]: 
oz_col2.Step_Forward(600)
In [32]: 
fig = figure( figsize=(10,8) )
ax = fig.add_subplot(111)
ax.plot( Tglobal + 273.15, log(level/1000), 'b-' )
ax.plot( oz_col.Tatm, log( oz_col.p / clim.ps ), 'c-' )
ax.plot( oz_col.Ts, 0, 'co', markersize=16 )
ax.plot( oz_col2.Tatm, log( oz_col2.p / clim.ps ), 'c--' )
ax.plot( oz_col2.Ts, 0, 'co', markersize=16 )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_yticks( log(level/1000) )
ax.set_yticklabels( level )
ax.set_title('Temperature profiles: observed (blue), RCE with ozone (cyan)', fontsize = 18)
ax.grid()

And we find that the troposphere warms, while the stratosphere cools!

Vertical structure of greenhouse warming in CESM model

In [33]: 
atmstr = ".cam.h0.clim.nc"
cesm_ctrl = nc.Dataset( datapath + 'som_control/som_control' + atmstr + endstr )
cesm_2xCO2 = nc.Dataset( datapath + 'som_2xCO2/som_2xCO2' + atmstr + endstr )
In [38]: 
T_cesm_ctrl = cesm_ctrl.variables['T'][:]
T_cesm_2xCO2 = cesm_2xCO2.variables['T'][:]
print T_cesm_ctrl.shape

(12, 26, 96, 144)

T_cesm_ctrl.shape

In [40]: 
T_cesm_ctrl_zon = mean( T_cesm_ctrl, axis=(0,3) )
T_cesm_2xCO2_zon = mean( T_cesm_2xCO2, axis=(0,3) )
In [46]: 
T_cesm_2xCO2_glob = empty_like( lev )
T_cesm_ctrl_glob = empty_like( lev )
for n in range( lev.size ):
    T_cesm_ctrl_glob[n] = sum( T_cesm_ctrl_zon[n,:] * cos( deg2rad( lat ) ) ) / sum( cos( deg2rad( lat ) ) )
    T_cesm_2xCO2_glob[n] = sum( T_cesm_2xCO2_zon[n,:] * cos( deg2rad( lat ) ) ) / sum( cos( deg2rad( lat ) ) )
In [49]: 
fig = figure( figsize=(10,8) )
ax = fig.add_subplot(111)
ax.plot( Tglobal + 273.15, log(level/1000), 'b-' )
ax.plot( oz_col.Tatm, log( oz_col.p / clim.ps ), 'c-' )
ax.plot( oz_col.Ts, 0, 'co', markersize=16 )
ax.plot( oz_col2.Tatm, log( oz_col2.p / clim.ps ), 'c--' )
ax.plot( oz_col2.Ts, 0, 'co', markersize=16 )
ax.plot( T_cesm_ctrl_glob, log( lev/clim.ps ), 'r-' )
ax.plot( T_cesm_2xCO2_glob, log( lev/clim.ps ), 'r--' )
ax.invert_yaxis()
ax.set_xlabel('Temperature (K)', fontsize=16)
ax.set_ylabel('Pressure (hPa)', fontsize=16 )
ax.set_yticks( log(level/1000) )
ax.set_yticklabels( level )
ax.set_title('Temperature profiles: observed (blue), RCE with ozone (cyan), CESM (red)', fontsize = 18)
ax.grid()

And we find that CESM has the same tendency for increased CO2: warmer troposphere, colder stratosphere.