There are five different average profiles for the tropics, subarctic summer, subarctic winter, midlatitude summer, midlatitude winter. These are called the US Standard Atmospheres. This notebook shows how to read and plot the soundings, and calculate the pressure and density scale heights.
In [1]:
from matplotlib import pyplot as plt
import matplotlib.ticker as ticks
import urllib
import numpy as np
from a301utils.a301_readfile import download
import h5py
In [2]:
filename='std_soundings.h5'
download(filename)
Use HDFViewer to check the layout of std_soundings.h5. There are five soundings. The soundings have six columns and 33 rows (i.e. 33 height levels). The variables are z, press, temp, rmix, den, o3den -- where rmix is the mixing ratio of water vapor, den is the dry air density and o3den is the ozone density. The units are m, pa, K, kg/kg, kg/m^3, kg/m^3
I will read the 6 column soundings into a pandas (panel data) DataFrame, which is like a matrix except the columns can be accessed by column name in addition to column number. The main advantage for us is that it's easier to keep track of which variables we're plotting
In [3]:
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)
We use these keys to get a dataframe with 6 columns, and 33 levels. Here's an example for the midsummer sounding
In [4]:
midsummer=sound_dict['midsummer']
print(midsummer.head())
list(midsummer.columns)
Out[4]:
In [5]:
%matplotlib inline
plt.style.use('ggplot')
meters2km=1.e-3
plt.close('all')
fig,(ax1,ax2)=plt.subplots(1,2,figsize=(11,8))
for a_name,df in sound_dict.items():
ax1.plot(df['temp'],df['z']*meters2km,label=a_name)
ax1.set(ylim=(0,40),title='Temp soundings',ylabel='Height (km)',
xlabel='Temperature (K)')
ax2.plot(df['rmix']*1.e3,df['z']*meters2km,label=a_name)
ax2.set(ylim=(0,8),title='Vapor soundings',ylabel='Height (km)',
xlabel='vapor mixing ratio (g/kg)')
ax1.legend()
_=ax2.legend()
Here is equation 1 of the hydrostatic balance notes
where
$$H=R_d T/g$$and here is the Python code to do that integral:
In [6]:
g=9.8 #don't worry about g(z) for this exercise
Rd=287. #kg/m^3
def calcScaleHeight(T,p,z):
"""
Calculate the pressure scale height H_p
Parameters
----------
T: vector (float)
temperature (K)
p: vector (float) of len(T)
pressure (pa)
z: vector (float) of len(T
height (m)
Returns
-------
Hbar: vector (float) of len(T)
pressure scale height (m)
"""
dz=np.diff(z)
TLayer=(T[1:] + T[0:-1])/2.
oneOverH=g/(Rd*TLayer)
Zthick=z[-1] - z[0]
oneOverHbar=np.sum(oneOverH*dz)/Zthick
Hbar = 1/oneOverHbar
return Hbar
Similarly, equation (5) of the hydrostatic balance notes is:
$$\frac{d\rho }{\rho} = - \left ( \frac{1 }{H} + \frac{1 }{T} \frac{dT }{dz} \right ) dz \equiv - \frac{dz }{H_\rho} $$Which leads to
$$\frac{ 1}{\overline{H_\rho}} = \frac{\int_{0 }^{z}\!\left [ \frac{1}{H} + \frac{1 }{T} \frac{dT }{dz} \right ] dz^\prime }{z-0} $$and the following python function:
In [7]:
def calcDensHeight(T,p,z):
"""
Calculate the density scale height H_rho
Parameters
----------
T: vector (float)
temperature (K)
p: vector (float) of len(T)
pressure (pa)
z: vector (float) of len(T
height (m)
Returns
-------
Hbar: vector (float) of len(T)
density scale height (m)
"""
dz=np.diff(z)
TLayer=(T[1:] + T[0:-1])/2.
dTdz=np.diff(T)/np.diff(z)
oneOverH=g/(Rd*TLayer) + (1/TLayer*dTdz)
Zthick=z[-1] - z[0]
oneOverHbar=np.sum(oneOverH*dz)/Zthick
Hbar = 1/oneOverHbar
return Hbar
In [8]:
sounding='tropics'
#
# grab the dataframe and get the sounding columns
#
df=sound_dict[sounding]
z=df['z'].values
Temp=df['temp'].values
press=df['press'].values
In [9]:
#
# limit calculation to bottom 10 km
#
hit=z<10000.
zL,pressL,TempL=(z[hit],press[hit],Temp[hit])
rhoL=pressL/(Rd*TempL)
Hbar= calcScaleHeight(TempL,pressL,zL)
Hrho= calcDensHeight(TempL,pressL,zL)
print("pressure scale height for the {} sounding is {:5.2f} km".format(sounding,Hbar*1.e-3))
print("density scale height for the {} is {:5.2f} km".format(sounding,Hrho*1.e-3))
In [10]:
theFig,theAx=plt.subplots(1,1)
theAx.semilogy(Temp,press/100.)
#
# need to flip the y axis since pressure decreases with height
#
theAx.invert_yaxis()
tickvals=[1000,800, 600, 400, 200, 100, 50,1]
theAx.set_yticks(tickvals)
majorFormatter = ticks.FormatStrFormatter('%d')
theAx.yaxis.set_major_formatter(majorFormatter)
theAx.set_yticklabels(tickvals)
theAx.set_ylim([1000.,50.])
theAx.set_title('{} temperature profile'.format(sounding))
theAx.set_xlabel('Temperature (K)')
_=theAx.set_ylabel('pressure (hPa)')
Now check the hydrostatic approximation by plotting the pressure column against
$$p(z) = p_0 \exp \left (-z/\overline{H_p} \right )$$vs. the actual sounding p(T):
In [11]:
fig,theAx=plt.subplots(1,1)
hydroPress=pressL[0]*np.exp(-zL/Hbar)
theAx.plot(pressL/100.,zL/1000.,label='sounding')
theAx.plot(hydroPress/100.,zL/1000.,label='hydrostat approx')
theAx.set_title('height vs. pressure for tropics')
theAx.set_xlabel('pressure (hPa)')
theAx.set_ylabel('height (km)')
theAx.set_xlim([500,1000])
theAx.set_ylim([0,5])
tickVals=[500, 600, 700, 800, 900, 1000]
theAx.set_xticks(tickVals)
theAx.set_xticklabels(tickVals)
_=theAx.legend(loc='best')
Again plot the hydrostatic approximation
$$\rho(z) = \rho_0 \exp \left (-z/\overline{H_\rho} \right )$$vs. the actual sounding $\rho(z)$:
In [12]:
fig,theAx=plt.subplots(1,1)
hydroDens=rhoL[0]*np.exp(-zL/Hrho)
theAx.plot(rhoL,zL/1000.,label='sounding')
theAx.plot(hydroDens,zL/1000.,label='hydrostat approx')
theAx.set_title('height vs. density for the tropics')
theAx.set_xlabel('density ($kg\,m^{-3}$)')
theAx.set_ylabel('height (km)')
theAx.set_ylim([0,5])
_=theAx.legend(loc='best')
Add cells to this notebook to:
1. Print out the density and pressure scale heights for each of the five soundings
2. Define a function that takes a sounding dataframe and returns the "total precipitable water", which is defined as:
$$W = \int_0^{z_{top}} \rho_v dz $$Do a change of units to convert $kg\,m^{-2}$ to $cm\,m^{-2}$ using the density of liquid water (1000 $kg\,m^{-3}$) -- that is, turn the kg of water in the 1 square meter column into cubic meters and turn that into $cm/m^{-2}$
3. Use your function to print out W for all five soundings
In [ ]: