In [9]:
from IPython.display import Image
import sys
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]:
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)
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)
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]:
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')
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)
In [ ]: