In [1]:
from IPython.display import Image
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 [2]:
Image('figures/wallace4_33.png',width=500)
Out[2]:
In this assignmet we'll work with 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 below to calculate the optical thickness for a $CO_2$-like absorbing gas. I ask your to find the transmissivities $t=\exp(-\tau_{tot} - \tau)$ and the weighting functions $\Delta t$, at 7 wavelengths, and use those plus $B_\lambda$ to calculate the radiance for the 7 channels for each of the 5 atmospheres.
In [3]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticks
from a301utils.a301_readfile import download
import h5py
from a301lib.radiation import Blambda,planckInvert
#
# from the hydrostatic noteboo
#
from a301lib.thermo import calcDensHeight
In [4]:
filename='std_soundings.h5'
download(filename)
In [5]:
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)
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 [6]:
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 [7]:
%matplotlib inline
r_gas=0.01 #kg/kg
k_lambda=0.01 #m^2/kg
df=sound_dict['tropics']
#
# set the top of the atmosphere at 20 km
#
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[7]:
In [8]:
#
# 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]
#
# find the height at mid-layer
#
mid_height=(height[1:] + height[:-1])/2.
#
# make a list of tuples of k_lambda and its label
# using zip
#
k_vals=zip(k_lambda_list,legend_string)
fig1,ax=plt.subplots(1,1,figsize=(10,10))
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)
ax.plot(tau,heightkm,label=k_label)
ax.legend()
1) make two more plots for the tropical atmosphere:
a) a plot of the transmissivity $t_\lambda = exp(-(\tau_{tot} - \tau)$ for each wavelength vs. height in km, with a legend, axes labels and a title
b) a plot of the weighting functions $\Delta t_\lambda$ vs. height for each wavelength, with a legend, axes labels and a title
2) Use Stull 8.4 with Blambda and planckInvert imported above to calculate the brighntess temperature for each of the 5 atmospheres at each of the 7 wavelengths. Make a plot of wavelength vs. brightness temperature with 5 lines on it, one for each atmosphere. The plot should have a title, axes labels and a legend for the 5 lines
Hint: I wrote two nested loops, looping through the dataframes in sound_dict and then looping over the 7 wavelengths with 7 k_lambdas giving me 7 taus, transmissivities and weighting functions. I initialized a dictionary with each sounding name as a key containing empty lists to hold the 7 brightness temperatures. Then in my wavelength loop I used Stull 8.4 to calculate the radiance at 20 km (our top of atmosphere), inverted to get brightness temperature and appended that to the list for that sounding. The brightness temperatures should range between about 220 K and 290 K
In [ ]: