Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

SkewT-LogP Sounding with Python

By: meteoblog on May 15th, 2012  |  syntax: Python  |  size: 3.53 KB  |  hits: 164  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
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
  7. # http://code.google.com/p/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()
  51.     _dry_adiabats()
  52.     _moist_adiabats()
  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.  
  85. def _dry_adiabats():
  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.  
  113. def _moist_adiabats():
  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()
create a new version of this paste RAW Paste Data