Check out the Pastebin Gadgets Shop. We have thousands of fun, geeky & affordable gadgets on sale :-)Want more features on Pastebin? Sign Up, it's FREE!
tweet

SkewT-LogP Sounding with Python

By: meteoblog on May 15th, 2012  |  syntax: Python  |  size: 3.53 KB  |  views: 478  |  expires: Never
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
1. import math
2. import numpy as np
3. import pylab as pl
4.
5.
6. # most of the code adapted from pywrfplot
8.
9. # purpose: plot raw data of Temperature and Dewpoint in skewT-log(p)-diagram
10.
11. skewness = 37.5
12. # Defines the ranges of the plot, do not confuse with P_bot and P_top
13. P_b = 105000.
14. P_t = 10000.
15. dp = 100.
16. plevs = np.arange(P_b,P_t-1,-dp)
17.
18. T_zero = 273.15
19.
20. # P_top must be the same as what is used in the WRF simulation
21. P_top = 10**4
22. P_bot = 10**5
23.
24.
25. L = 2.501e6 # latent heat of vaporization
26. R = 287.04  # gas constant air
27. Rv = 461.5  # gas constant vapor
28. eps = R/Rv
29.
30. cp = 1005.
31. cv = 718.
32. kappa = (cp-cv)/cp
33. g = 9.81
34.
35.
36. # constants used to calculate moist adiabatic lapse rate
37. # See formula 3.16 in Rogers&Yau
38. a = 2./7.
39. b = eps*L*L/(R*cp)
40. c = a*L/R
41.
42.
43.
44.
45.
46. def SkewTPlot():
47.
48.     pl.figure()
49.     _isotherms()
50.     _isobars()
53.
54.     _plot_data()
55.
56.
57.     pl.axis([-40,50,P_b,P_t])
58.     pl.xlabel('Temperatur (\$^{\circ}\! C\$)')
59.     xticks = np.arange(-40,51,5)
60.     pl.xticks(xticks,['' if tick%10!=0 else str(tick) for tick in xticks])
61.     pl.ylabel('Pressure (hPa)')
62.     yticks = np.arange(P_bot,P_t-1,-10**4)
63.     pl.yticks(yticks,yticks/100)
64.
65.     pl.show()
66.     pl.close()
67.     pl.clf()
68.
69.
70. def _skewnessTerm(P):
71.     return skewness * np.log(P_bot/P)
72.
73.
74. def _isotherms():
75.     for temp in np.arange(-100,50,10):
76.         pl.semilogy(temp + _skewnessTerm(plevs), plevs,  basey=math.e, \
77.                      color = ('blue'), \
78.                      linestyle=('solid' if temp == 0 else 'dashed'), linewidth = .5)
79.
80. def _isobars():
81.     for n in np.arange(P_bot,P_t-1,-10**4):
82.         pl.plot([-40,50], [n,n], color = 'black', linewidth = .5)
83.
84.
86.     for tk in T_zero+np.arange(-30,210,10):
87.         dry_adiabat = tk * (plevs/P_bot)**kappa - T_zero + _skewnessTerm(plevs)
88.         pl.semilogy(dry_adiabat, plevs, basey=math.e, color = 'brown', \
89.                      linestyle='dashed', linewidth = .5)
90.
91.
92.
93. def es(T):
94.     """Returns saturation vapor pressure (Pascal) at temperature T (Celsius)
95.    Formula 2.17 in Rogers&Yau"""
96.     return 611.2*np.exp(17.67*T/(T+243.5))
97.
98.
99. def gamma_s(T,p):
100.     """Calculates moist adiabatic lapse rate for T (Celsius) and p (Pa)
101.    Note: We calculate dT/dp, not dT/dz
102.    See formula 3.16 in Rogers&Yau for dT/dz, but this must be combined with
103.    the dry adiabatic lapse rate (gamma = g/cp) and the
104.    inverse of the hydrostatic equation (dz/dp = -RT/pg)"""
105.     esat = es(T)
106.     wsat = eps*esat/(p-esat) # Rogers&Yau 2.18
107.     numer = a*(T+T_zero) + c*wsat
108.     denom = p * (1 + b*wsat/((T+T_zero)**2))
109.     return numer/denom # Rogers&Yau 3.16
110.
111.
112.
114.     ps = [p for p in plevs if p<=P_bot]
115.     for temp in np.concatenate((np.arange(-40.,10.1,5.),np.arange(12.5,45.1,2.5))):
116.         moist_adiabat = []
117.         for p in ps:
118.             temp -= dp*gamma_s(temp,p)
119.             moist_adiabat.append(temp + _skewnessTerm(p))
120.         pl.semilogy(moist_adiabat, ps, basey=math.e, color = 'green', \
121.                      linestyle = 'dashed', linewidth = .5)
122.
123.
124.
125. def _plot_data():
126.     p,h,T,Td = np.loadtxt('/Users/daniel/tmp/soundingdata.txt',usecols=range(0,4),unpack=True)
127.
128.     pl.semilogy(T+ _skewnessTerm(p*100),p*100,basey=math.e, color=('black'),linestyle=('solid'),linewidth= 1.5)
129.
130.     pl.semilogy(Td+ _skewnessTerm(p*100),p*100,basey=math.e, color=('black'),linestyle=('solid'),linewidth= 1.5)
131.
132.     return 1
133.
134.
135. SkewTPlot()
clone this paste RAW Paste Data
Top