Rad_calc solution


In [9]:
from IPython.display import Image
import sys

Weighting functions in the $CO_2$ 15 $\mu m$ absorption band

Below is a plot of radiance (or intensity) (left axis) and brightness temperature (right axis) vs. wavenumber near the main $CO_2$ absorption band. Wavenumber is defined as $1/\lambda$; the center of the absorption band is at $\lambda = 15\ \mu m$ which is a wavenumber of 1/0.0015 = 666 $cm^{-1}$. The VTPR (vertical temperature profiling radiometer) has six channels: Channel 1 is at the center of the band -- it has the lowest transmissivity and is measuring photons coming from 45 km at the top of the stratosphere (See Stull Chapter 1, Figure 1.10). As the channel number increases from 2-6 the transmissivity also increases, and the photons originate from increasing lower levels of the atmosphere with increasing kinetic temperatures. Note that the different heights for the peaks in each of the weighting functions.


In [10]:
Image('figures/wallace4_33.png',width=500)


Out[10]:

Assignment 9

This notebook uses the five standard atmospheres from hydrostatic.ipynb to show how to use Stull eq. 8.4 to calculate radiance at the top of the atmosphere for wavelengths similar to the 6 sounder channels shown in the above figure. I define a new function find_tau to calculate the optical thickness for a $CO_2$-like absorbing gas, and get the weighting functions at 7 wavelengths from the transmission. I then ask you to plot the brightness temperatures at each wavelength for each of the five standard atmospheres.


In [11]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticks
from a301utils.a301_readfile import download
import h5py
#
# from the hydrostatic noteboo
#
from a301lib.thermo import calcDensHeight

In [12]:
filename='std_soundings.h5'
download(filename)


std_soundings.h5 already exists, need to either delete or rename for new download

In [13]:
from pandas import DataFrame
with h5py.File(filename) as infile:
    sound_dict={}
    print('soundings: ',list(infile.keys()))
    #
    # names are separated by commas, so split them up
    # and strip leading blanks
    #
    column_names=infile.attrs['variable_names'].split(',')
    column_names = [item.strip() for item in column_names]
    column_units = infile.attrs['units'].split(',')
    column_units = [item.strip() for item in column_units]
    for name in infile.keys():
        data = infile[name][...]
        sound_dict[name]=DataFrame(data,columns=column_names)


soundings:  ['midsummer', 'midwinter', 'subsummer', 'subwinter', 'tropics']

mass absorption coefficient for fake $CO_2$

To keep things simple I'm going to make up a set of 7 absorption coefficients that will give weighting functions that look something like the VPTR.

We have been working with the volume absorption coefficient:

$$\beta_\lambda = n b\ (m^{-1})$$

where $n$ is the absorber number concentration in $\#\,m^{-3}$ and $b$ is the absorption cross section of a molecule in $m^2$. For absorbing gasses like $CO_2$ that obey the ideal gas law $n$ depends inversely on temperature -- which changes rapidly with height.

A better route is to use the mass absorption coefficient $k_\lambda$:

$$k_\lambda = \frac{n b}{\rho_{air}}$$

units: $m^2/kg$. With this definition the optical depth is:

$$\tau_\lambda = \int_0^z \rho_{air} k_\lambda dz$$

Why is this an improvement? We now have an absorption coefficient that can is roughly constant with height and can be taken out of the integral.


In [14]:
Rd = 287.
def find_tau(r_gas,k_lambda,df):
    """
    given a data frame df with a standard sounding, return the optical depth assuming an
    absorbing gas with a constant mixing ration r_gas and absorption coefficient k_lambda
    
    Parameters
    ----------
    
    r_gas: float
       absorber mixing ratio, kg/kg
       
    k_lambda: float
       mass absorption coefficient, m^2/kg
       
    df:  dataframe
       sounding with height levels as rows, columns 'temp': temperature in K
       'press': pressure in Pa
       
    Returns
    -------
    tau: vector (float)
      optical depth measured from the surface at each height in df
    """
    #
    # density scale height
    #
    Hdens = calcDensHeight(df)
    #
    # surface density
    #
    rho0 = df['press']/(Rd*df['temp'])
    rho = rho0*np.exp(-df['z']/Hdens)
    height = df['z'].values
    tau=np.empty_like(rho)
    #
    # start from the surface
    #
    tau[0]=0
    num_levels=len(rho)
    num_layers=num_levels-1
    for index in range(num_layers):
        delta_z=height[index+1] - height[index]
        delta_tau=r_gas*rho[index]*k_lambda*delta_z
        tau[index+1]=tau[index] + delta_tau
    return tau

In [15]:
%matplotlib inline
r_gas=0.01  #kg/kg
k_lambda=0.01  #m^2/kg
df=sound_dict['tropics']
top = 20.e3
df = df.loc[df['z']<top]
height = df['z']
press = df['press']
tau=find_tau(r_gas,k_lambda,df)
fig1,axis1=plt.subplots(1,1)
axis1.plot(tau,height*1.e-3)
axis1.set_title('vertical optical depth vs. height')
axis1.set_ylabel('height (km)')
axis1.set_xlabel('optical depth (no units)')
fig2,axis2=plt.subplots(1,1)
axis2.plot(tau,press*1.e-3)
axis2.invert_yaxis()
axis2.set_title('vertical optical depth vs. pressure')
axis2.set_ylabel('pressure (kPa)')
axis2.set_xlabel('optical depth (no units)')


Out[15]:
<matplotlib.text.Text at 0x1177b2390>

In [16]:
r_gas=0.01  #kg/kg
#
# assign the 7 k_lambdas to 7 CO2 absorption band wavelengths
# (see Wallace and Hobbs figure 4.33)
#
wavenums=np.linspace(666,766,7)  #wavenumbers in cm^{-1}
wavelengths=1/wavenums*1.e-2  #wavelength in m
wavelengths_um = wavelengths*1.e6  # in microns
print('channel wavelengths (microns) ',wavelengths_um)  #microns
df=sound_dict['tropics']
top = 20.e3  #stop at 20 km
df = df.loc[df['z']< top]
height = df['z'].values
mid_height = (height[1:] - height[:-1])/2.
#
# here are the mass absorption coefficients for each of the 7 wavelengths
# in m^2/kg
#
k_lambda_list=np.array([ 0.175 , 0.15 ,  0.125 , 0.1  ,  0.075,  0.05 ,  0.025])
legend_string=["{:5.3f}".format(item) for item in k_lambda_list]
#
# make a list of tuples of k_lambda and its label
# using zip
#
k_vals=zip(k_lambda_list,legend_string)
#
#  find the height at mid-layer
#
mid_height=(height[1:] + height[:-1])/2.
fig1,axis_list=plt.subplots(2,2,figsize=(14,14))
#
# turn axis_list from 2 x 2 to 4 x 1
#
axis_list=np.array(axis_list).ravel()
#
# throw away the 4th plot using delaxes -- only plotting 3 graphs
#
fig1.delaxes(axis_list[3])
axis_list=axis_list[:3] 
heightkm=height*1.e-3
mid_heightkm=mid_height*1.e-3
for k_lambda,k_label in k_vals:
    tau=find_tau(r_gas,k_lambda,df)
    tau_tot=tau[-1]
    axis_list[0].plot(tau,heightkm,label=k_label)
    trans=np.exp(-(tau_tot - tau))
    axis_list[1].plot(trans,heightkm,label=k_label)
    del_trans=np.diff(trans)
    axis_list[2].plot(del_trans,mid_heightkm,label=k_label)
titles = ['optical depth for 7 values of $k_\lambda$',
        'transmittance for 7 values of $k_\lambda$',
        'weighting function for 7 values of $k_\lambda$']
xlabels = ['optical depth',
          'transmittance',
          'weighting function']
for axis,title,xlabel in zip(axis_list,titles,xlabels):
    axis.set(title=title,xlabel=xlabel)
[axis.set_ylabel('height (km)') for axis in axis_list]
_=[axis.legend(loc='best') for axis in axis_list]
fig1.savefig('trans_plots.png')


channel wavelengths (microns)  [ 15.01501502  14.6484375   14.2993327   13.96648045  13.64877161
  13.34519573  13.05483029]

In [17]:
from collections import defaultdict
from a301lib.radiation import Blambda,planckInvert
Tbright=defaultdict(list)
for the_sound,the_df in sound_dict.items():
    for the_wave_m,the_k in zip(wavelengths,k_lambda_list):
        tau = find_tau(r_gas,the_k,the_df)
        tau_tot=tau[-1]
        temps=the_df['temp'].values
        Bsfc=Blambda(the_wave_m,temps[0])
        sfc_contrib = Bsfc*np.exp(-tau_tot)
        mid_temps=(temps[1:] + temps[:-1])/2.
        Batm=Blambda(the_wave_m,mid_temps)
        trans = np.exp(-(tau_tot - tau))
        weights = np.diff(trans)
        atm_contrib = np.sum(Batm*weights)
        Lsat = sfc_contrib + atm_contrib
        the_temp = planckInvert(the_wave_m,Lsat)
        Tbright[the_sound].append(the_temp)
    Tbright[the_sound]=np.array(Tbright[the_sound])

fig, ax = plt.subplots(1,1)
the_wave_um = wavelengths*1.e6
for the_sound,the_temps in Tbright.items():
    ax.plot(the_wave_um,the_temps,label=the_sound)
ax.legend(loc='best')
ax.set(ylabel="Tbright (K)",xlabel='wavelength (microns)')
print(the_wave_um)


[ 15.01501502  14.6484375   14.2993327   13.96648045  13.64877161
  13.34519573  13.05483029]

In [ ]: