爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
123
返回列表 发新帖
楼主: shuihu88888

[分享资料] 如何利用GRADS和1*1 Ncep再分析资料画CAPE对流有效位能

[复制链接]
发表于 2018-3-8 11:04:01 | 显示全部楼层
happyouc 发表于 2018-3-8 09:21
我的qq289619047 谢谢

我没有专门算CAPE的程序,传个skewplot的gs程序,可能会输出CAPE,前半截只要改一下路径和变量名就好了,然后直接套用到plotskew()函数里面,后半截是skewplot.gs程序本身:
  1. 'reinit'
  2. 'open D:\gradsdata\ncepfnl_1_6h\fnl_20120719_00_00_c.ctl'
  3. 'set grid off'
  4. 'set grads off'
  5. 'set z  1 21'
  6. 'set lat 37.75'
  7. 'set lon 103.92'
  8. tt=1
  9. while(tt<=13)
  10. 'c'
  11. 'set grads off'
  12. 'set t 'tt' '
  13. 'define tmp=tmpprs-273.16'
  14. 'define RH=rhprs'
  15. 'td=tmp-((14.55+0.114*tmp)*(1-0.01*RH) + pow((2.5+0.007*tmp)*(1-0.01*RH),3) + (15.9+0.37*tmp)*pow((1-0.01*RH),14))'
  16. 'define uv=mag(vgrdprs,ugrdprs)'
  17. 'define pi=3.14159'
  18. if (ugrdprs>=0.)
  19. if (vgrdprs>=0.)
  20. 'define dir=180+atan2(ugrdprs,vgrdprs)*180/pi'
  21. endif
  22. endif
  23. if (vgrdprs<=0.)
  24. if (ugrdprs>=0.)
  25. 'define dir=270+atan2(ugrdprs,vgrdprs)*180/pi'
  26. endif
  27. endif
  28. if (ugrdprs<=0.)
  29. if (vgrdprs>=0.)
  30. 'define dir=180-atan2(ugrdprs,vgrdprs)*180/pi'
  31. endif
  32. endif
  33. if (vgrdprs<0.)
  34. if (ugrdprs<=0.)
  35. 'define dir=270-atan2(ugrdprs,vgrdprs)*180/pi'
  36. endif
  37. endif
  38. *if(ugrdprs > 0.);'define dir=90-180/pi*atan2(vgrdprs,ugrdprs)';endif
  39. *if(ugrdprs < 0.);'define dir=270-180/pi*atan2(vgrdprs,ugrdprs)';endif
  40. rc=plotskew(tmp,td,uv,dir)

  41. 'printim D:\gradsdata\ncepfnl_1_6h\0719\站点垂直剖面\观测点'%tt%'.gif  white'
  42. tt=tt+1
  43. endwhile
  44. ;

  45. function plotskew(sndtemp,snddewp,sndspd,snddir)
  46. *************************************************************************
  47. *
  48. * GrADS Script to Plot a SkewT/LogP Diagram      
  49. *
  50. * Bob Hart
  51. * Penn State University / Dept of Meteorology
  52. * Last Update:  January 23, 2001
  53. *
  54. * Recent Changes:
  55. *
  56. * 01/23/01 - Fixed a small bug in the theta-e calculation.  
  57. *            Errors averaged 0.5-3K.  Thank you George Bryan.
  58. *
  59. * 11/10/99 - Change in calculation method for CAPE/CIN.  Trapezoid
  60. *            integration method is now used.  Speeds up execution
  61. *            by 25%, and increases accuracy by 5-10%.
  62. *
  63. * 10/18/99 - Minor glitch fixed that occasionally caused crash.
  64. *
  65. *  8/26/99 - Datasets with missing data can now be used.
  66. *
  67. * Features:
  68. *   - All features of standard skewt/logp plot
  69. *   - RH sounding
  70. *   - LCL location
  71. *   - Parcel trajectory for both sfc based convection and elevated from
  72. *     most unstable level (highest theta-e level reported)
  73. *   - Stability indices and precipitable water calculations
  74. *   - CAPE & CIN Calculations
  75. *   - Wind Profile
  76. *   - Hodograph / Hodograph scaling
  77. *   - Helicity and SR Helicity Calculations and Display
  78. *   - Color aspects of output
  79. *   - Line Thickness, style aspects of output
  80. *   - Can be run in either PORTRAIT or LANDSCAPE mode.
  81. *
  82. * There are numerous tunable parameters below to change the structure
  83. * and output for the diagram.
  84. *
  85. * Function Arguments:
  86. *    sndtemp - temperature data (Celsius) as a function of pressure
  87. *    snddewp - dewpoint data (Celsius) as a function of pressure      
  88. *    sndspd  - wind speed data (knots) as a function of pressure
  89. *    snddir  - wind direction data as a function of pressure
  90. *
  91. * Use '-1' for any of the above 4 arguments to indicate that you
  92. * are not passing that variable.  The appropriate options will
  93. * be ignored based on your specifying '-1' for that variable.
  94. *
  95. * NOTE:  Make sure to set the vertical range of the plot before running.
  96. *        I.e., "SET LEV 1050 150", for example.   This does not have to
  97. *        be limited to the pressure range of your data.
  98. *
  99. * Labelling:  Pressure/Height is labelled along left side.  Temperature is
  100. *             labelled along bottom.  Mixing ratio is labelled along right
  101. *             side/top.
  102. *
  103. *
  104. * PROBLEMS:  First check out the web page for the script (which also
  105. *            has a link to a FAQ with answers to many common questions
  106. *            about using the script):
  107. *            http://www.ems.psu.edu/~hart/skew.html
  108. *
  109. * Please send any further problems, comments, or suggestions to
  110. * <hart@ems.psu.edu>
  111. *
  112. * ACKNOWLEDGMENTS:  Thanks go to the innumerable users who have helped
  113. * fine tune the script from the horrible mess from which it began.
  114. * In particular, thanks go out to Steve Lord (NCEP), Mike Fiorino (ECMWF),
  115. * George Bryan (PSU), Davide Sacchetti (CMIRL), and Enrico Minguzzi (CMIRL).
  116. *
  117. **************************************************************************
  118. *           !!!!!   BEGINNING OF USER-SPECIFIED OPTIONS  !!!!!!
  119. **************************************************************************
  120. *
  121. * --------------------- Initialization options  ----------------------
  122. *
  123. * ClrScrn = Whether to clear the screen before drawing diagram
  124. *           [1 = yes, 0 = no]

  125. ClrScrn = 1

  126. *
  127. * ------------------- Define Skew-T Diagram Shape/Slope-----------------
  128. *
  129. * (P1,T1) = Pres, Temp of some point on left-most side
  130. * (P2,T2) = Pres, Temp of some point on right-most side
  131. * (P3,T3) = Pres, Temp of some point in diagram which is mid-point
  132. *           in the horizontal between 1 and 2.
  133. *
  134. * P1, P2, P3 are in hPa ; T1, T2, T3 are in Celsius
  135. *
  136. * These define the SLOPE and WIDTH of the diagram as you see it but DO NOT
  137. * DEFINE THE HEIGHT of the diagram as you see it.  In other words,
  138. * 1 and 2 do NOT necessarily need to be at the bottom of the diagram and
  139. * 3 does NOT necessarily need to be at the top.  THE VERTICAL PRESSURE
  140. * RANGE OF THE SKEWT AS YOU SEE IT IS DETERMINED BY YOUR 'SET Z ...'  
  141. * COMMAND OR THE 'SET LEV ...' COMMAND BEFORE RUNNING THIS SCRIPT.
  142. *
  143. *    _______________________
  144. *   |                       |
  145. *   |                       |
  146. *   |           3           |
  147. *   |                       |
  148. *   |                       |
  149. *   |                       |
  150. *   |                       |
  151. *   |                       |
  152. *   |                       |
  153. *   |                       |
  154. *   |                       |
  155. *   |1                     2|
  156. *   |                       |
  157. *   |_______________________|
  158. *   
  159. *
  160. * A good set of defining points are given below.   Feel free
  161. * to experiment with variations.


  162. P1 = 1000
  163. T1 = -40

  164. P2 = 1000
  165. T2 = 40

  166. P3 = 200
  167. T3 = -50

  168. * Another good set of defining points suggested by George Bryan (PSU)
  169. * are:
  170. *
  171. * P1 = 1000
  172. * T1 = -30
  173. *
  174. * P2 = 1000
  175. * T2 = 40
  176. *
  177. * P3 = 500
  178. * T3 = -18

  179. * ------------------- Contour Intervals / Levels --------------------------
  180. *
  181. * All variables below are contour intervals/levels for diagram
  182. *
  183. * Thetaint = interval for potential temperature lines
  184. * Thetwint = interval for moist pseudo adiabats
  185. * tempint  = interval for temperature lines
  186. * wsclevs  = contour LEVELS for mixing ratio lines
  187. *
  188. *
  189. thetaint= 10
  190. thetwint= 5
  191. tempint = 10
  192. wsclevs = "1 2 3 4 6 8 10 15 20 25 30 35 40"
  193. *
  194. *
  195. * ------------------------ Output Options --------------------------------
  196. *
  197. * All variables below are logical .. 1=yes, 0=no, unless otherwise
  198. * specified.
  199. *
  200. * DrawBarb = Draw wind barbs along right side of plot
  201. * DrawThet = Draw dry adiabats
  202. * DrawThtw = Draw moist pseudo-adiabats
  203. * DrawTemp = Draw temperature lines
  204. * DrawMix  = Draw mixing ratio lines
  205. * DrawTSnd = Draw temperature sounding
  206. * DrawDSnd = Draw dewpoint sounding
  207. * DrawRH   = Draw relative humidity sounding
  208. * DrawPrcl = Draw parcel path from surface upward
  209. * DrawPMax = Draw parcel path from most unstable level upward
  210. * DrawIndx = Display stability indices & CAPE
  211. * DrawHeli = Calculate and display absolute and storm-relative helicity
  212. * DrawHodo = Draw hodograph
  213. * DrawPLev = Draw Pressure Levels
  214. * DrawZLev = Draw height levels and lines
  215. *            0 = no lines
  216. *            1 = above ground level (AGL)
  217. *            2 = above sea level (ASL)
  218. * DrawZSTD = Draw Height levels using standard atm lapse rate           
  219. * LblAxes  = Label the x,y axes (temperature, pressure,mixing ratio)
  220. *
  221. * ThtwStop = Pressure level at which to stop drawing Theta-w lines
  222. * MixStop  = Pressure level at which to stop drawing Mixratio lines

  223. DrawBarb= 1
  224. DrawThet= 1
  225. DrawThtw= 0
  226. DrawTemp= 1
  227. DrawMix = 1
  228. DrawTSnd= 1
  229. DrawDSnd= 1
  230. DrawRH  = 0
  231. DrawPrcl= 1
  232. DrawPMax= 0
  233. DrawIndx= 1
  234. DrawHeli= 1
  235. DrawHodo= 1
  236. DrawPLev= 1
  237. DrawZLev= 0
  238. DrawZSTD= 0
  239. LblAxes = 1

  240. ThtwStop = 200
  241. MixStop  = 600


  242. *
  243. * -----------------  Sounding Geography options ------------------------
  244. *
  245. * SfcElev = Elevation above sea-level (meters) of lowest level reported
  246. *           in sounding.  Used only if DrawZLev = 2

  247. SfcElev = 0


  248. *
  249. * ------------------ Thermodynamic Index Options --------------------
  250. *
  251. * All variables here are in inches.  Use -1 for the default values.
  252. *
  253. *  Text1XC = X-location of midpoint of K,TT,PW output box
  254. *  Text1YC = Y-location of midpoint of K,TT,PW output box
  255. *  Text2XC = X-Location of midpoint of surface indices output box
  256. *  Text2YC = Y-location of midpoint of surface indices output box
  257. *  Text3XC = X-Location of midpoint of most unstable level-based indices
  258. *            output box
  259. *  Text3YC = Y-location of midpoint of most unstable level-based indices
  260. *            output box

  261. Text1XC = -1
  262. Text1YC = -1
  263. Text2XC = -1
  264. Text2YC = -1
  265. Text3XC = -1
  266. Text3YC = -1

  267. *
  268. * ----------------- Wind Barb Profile Options ----------------------------
  269. *
  270. * All variables here are in units of inches, unless otherwise specified
  271. *
  272. *  barbint = Interval for plotting barbs (in units of levels)
  273. *  poleloc = X-Location of profile.  Choose -1 for the default.
  274. *  polelen = Length of wind-barb pole
  275. *  Len05   = Length of each 5-knot barb
  276. *  Len10   = Length of each 10-knot barb
  277. *  Len50   = Length of each 50-knot flag
  278. *  Wid50   = Width of base of 50-knot flag
  279. *  Spac50  = Spacing between 50-knot flag and next barb/flag
  280. *  Spac10  = Spacing between 10-knot flag and next flag
  281. *  Spac05  = Spacing between 5-knot flag and next flag
  282. *  Flagbase= Draw flagbase (filled circle) for each windbarb [1=yes, 0 =no]
  283. *  Fill50  = Solid-fill 50-knot flag [1=yes, 0=no]
  284. *  barbline= Draw a vertical line connecting all the wind barbs [1=yes, 0=no]
  285. *
  286. barbint = 1
  287. poleloc = -1
  288. polelen = 0.35
  289. len05   = 0.07
  290. len10   = 0.15
  291. len50   = 0.15
  292. wid50   = 0.06
  293. spac50  = 0.07
  294. spac10  = 0.05
  295. spac05  = 0.05
  296. Fill50  = 1
  297. flagbase= 1
  298. barbline= 1

  299. *
  300. *
  301. *---------------- Hodograph Options -------------------------------------
  302. *
  303. * All variables here are in units of inches, unless otherwise specified
  304. *
  305. * HodXcent= x-location of hodograph center.  Use -1 for default location.
  306. * HodYcent= y-location of hodograph center.  Use -1 for default location.
  307. * HodSize = Size of hodograph in inches
  308. * NumRing = Number of rings to place in hodograph (must be at least 1)
  309. * HodRing = Wind speed increment of each hodograph ring
  310. * HodoDep = Depth (above lowest level in hPa) of end of hodograph trace
  311. * TickInt = Interval (in kts) at which tick marks are drawn along the axes
  312. *           Use 0 for no tick marks.
  313. * TickSize= Size of tick mark in inches
  314. * Text4XC = X-location of midpoint of hodograph text output. Use -1 for default.
  315. * Text4YC = Y-location of midpoint of hodograph text output. Use -1 for default.

  316. HodXcent= -1
  317. HodYcent= -1
  318. HodSize = 2
  319. NumRing = 3
  320. HodRing = 15
  321. HodoDep = 300
  322. TickInt = 5
  323. TickSize= 0.05
  324. Text4XC = -1
  325. Text4YC = -1

  326. *--------------- Helicity Options ---------------------------------------
  327. *
  328. * MeanVTop = Top pressure level (hPa) of mean-wind calculation
  329. * MeanVBot = Bottom pressure level (hPa) of mean-wind calculation
  330. * HelicDep = Depth in hPa (above ground) of helicity integration
  331. * StormMot = Type of storm motion estimation scheme.  Use following:  
  332. *            0 = No departure from mean wind.
  333. *            1 = Davies-Jones (1990) approach
  334. * FillArrw = Whether to fill the arrowhead of the storm motion vector
  335. *            [1 = yes, 0 = no]

  336. MeanVTop= 300
  337. MeanVBot= 850
  338. HelicDep= 300
  339. StormMot= 1
  340. FillArrw= 1

  341. *
  342. *---------------- Color Options ------------------------------------------
  343. *
  344. * ThetCol = Color of dry adiabats
  345. * TempCol = Color of temperature lines
  346. * MixCol  = Color of mixing ratio lines
  347. * ThtwCol = Color of moist adiabats
  348. * TSndCol = Color of Temperature Sounding
  349. * DSndCol = Color of Dewpoint Sounding
  350. * RHCol   = Color of RH Sounding
  351. * PrclCol = Color of parcel trace
  352. * BarbCol = Color of wind barbs (choose -1 for color according to speed)
  353. * HodoCol = Color of hodograph trace

  354. ThetCol = 2
  355. TempCol = 4
  356. MixCol  = 7
  357. ThtwCol = 3
  358. TSndCol = 1
  359. DSndCol = 1
  360. RHCol   = 3
  361. PrclCol = 5
  362. BarbCol = -1
  363. HodoCol = 1

  364. *
  365. *-------------------- Line Style Options ------------------------------------
  366. *
  367. * GrADS Styles: 1=solid;2=long dash;3=short dash;4=long,short dashed;
  368. *               5=dotted;6=dot dash;7=dot dot dash
  369. *
  370. * ThetLine = Line Style of dry adiabats
  371. * TempLine = Line Style of temperature lines
  372. * MixLine  = Line Style of mixing ratio lines
  373. * ThtwLine = Line Style of moist adiabats
  374. * TSndLine = Line Style of Temperature Sounding
  375. * DSndLine = Line Style of Dewpoint Sounding
  376. * RHLine   = Line Style of RH sounding
  377. * PrclLine = Line Style of parcel trace
  378. * HodoLine = Line Style of hodograph trace
  379. *

  380. ThetLine = 1
  381. TempLine = 1
  382. MixLine  = 5
  383. ThtwLine = 3
  384. TSndLine = 1
  385. DSndLine = 2
  386. RHLine   = 3
  387. PrclLine = 1
  388. HodoLine = 1

  389. *
  390. *------------------- Line Thickness Options---------------------------------
  391. * GrADS Line Thickness: increases with increasing number. Influences
  392. *                       hardcopy output more strongly than screen output.
  393. *
  394. *
  395. * ThetThk = Line Thickness of dry adiabats
  396. * TempThk = Line Thickness of temperature lines
  397. * MixThk  = Line Thickness of mixing ratio lines
  398. * ThtwThk = Line Thickness of moist adiabats
  399. * TSndThk = Line Thickness of temperature sounding
  400. * DSndThk = Line thickness of dewpoint sounding
  401. * RHThk   = Line thickness of RH sounding
  402. * PrclThk = Line thickness of parcel trace
  403. * HodoThk = Line thickness of hodograph trace
  404. * BarbThk = Line thickness of wind barbs

  405. ThetThk = 3
  406. TempThk = 3
  407. MixThk  = 3
  408. ThtwThk = 3
  409. TSndThk = 8
  410. DSndThk = 8
  411. RHThk   = 8
  412. PrclThk = 6
  413. HodoThk = 3
  414. BarbThk = 2

  415. *
  416. *------------------- Data Point Marker Options -----------------------------
  417. * GrADS Marker Types: 0 = none ; 1 = cross ; 2 = open circle ;
  418. *                     3 = closed circle ; 4 = open square ; 5 = closed square
  419. *                     6 = X ; 7 = diamond ; 8 = triangle ; 9 = none
  420. *                    10 = open circle with vertical line ; 11 = open oval
  421. *
  422. * TSndMrk = Mark type of data point marker for temperature sounding
  423. * DSndMrk = Mark type of data point marker for dewpoint sounding
  424. * RHMrk   = Mark type of data point marker for relative humidity sounding
  425. * MrkSize = Mark size (inches) of each data marker

  426. TSndMrk = 0
  427. DSndMrk = 0
  428. RHMrk   = 0
  429. MrkSize = 0.1


  430. * !!!!! YOU SHOULD NOT NEED TO CHANGE ANYTHING BELOW HERE !!!!!
  431. ****************************************************************************

  432. *-------------------------------------------
  433. * grab user-specified environment dimensions
  434. *-------------------------------------------

  435. "q dims"
  436. rec=sublin(result,2)
  437. _xtype=subwrd(rec,3)
  438. _xval=subwrd(rec,9)
  439. rec=sublin(result,3)
  440. _yval=subwrd(rec,9)
  441. _ytype=subwrd(rec,3)
  442. rec=sublin(result,4)
  443. _ptype=subwrd(rec,3)
  444. _pmax=subwrd(rec,6)
  445. _pmin=subwrd(rec,8)
  446. _zmin=subwrd(rec,11)
  447. _zmax=subwrd(rec,13)
  448. rec=sublin(result,5)
  449. _ttype=subwrd(rec,3)
  450. _tval=subwrd(rec,9)

  451. "q file"
  452. rec=sublin(result,5)
  453. _zmaxfile=subwrd(rec,9)

  454. *-------------------------------------------------------------
  455. * Check to ensure that dimensions are valid.  Warn & exit if not.
  456. *--------------------------------------------------------------

  457. dimrc=0
  458. If (_xtype != "fixed")
  459.   say "X-Dims Error:  Not fixed.  Use 'set lon' or 'set x' to specify a value."
  460.   dimrc=-1
  461. Endif

  462. If (_ytype != "fixed")
  463.   say "Y-Dims Error:  Not fixed.  Use 'set lat' or 'set y' to specify a value"
  464.   dimrc=-1
  465. Endif

  466. If (_ptype != "varying")
  467.    say "Z-Dims Error:  Not varying.  Use 'set lev' or 'set z' to specify a range."
  468.    dimrc=-1
  469. Endif

  470. If (_ttype != "fixed")
  471.   say "Time Error:     Not fixed.  Use 'set time' or 'set t' to specify a value"
  472.   dimrc=-1
  473. Endif


  474. If (dimrc < 0)
  475.   Return(-1)
  476. Endif


  477. *
  478. * A few global variables used in units conversion
  479. *

  480. _pi=3.14159265
  481. _dtr=_pi/180
  482. _rtd=1/_dtr
  483. _ktm=0.514444
  484. _mtk=1/_ktm

  485. * A few global constants used in thermo calcs

  486. _C0=0.99999683
  487. _C1=-0.90826951/100
  488. _C2= 0.78736169/10000
  489. _C3=-0.61117958/1000000
  490. _C4= 0.43884187/pow(10,8)
  491. _C5=-0.29883885/pow(10,10)
  492. _C6= 0.21874425/pow(10,12)
  493. _C7=-0.17892321/pow(10,14)
  494. _C8= 0.11112018/pow(10,16)         
  495. _C9=-0.30994571/pow(10,19)

  496. * A pressure array of power calculations which should be performed
  497. * only once to reduce execution time.

  498. zz=1100
  499. while (zz > 10)
  500.     subscr=zz/10
  501.     _powpres.subscr=pow(zz,0.286)
  502.     zz=zz-10
  503. endwhile

  504. *
  505. * Turn off options not available due to user data limitations
  506. *

  507. If (ClrScrn = 1)
  508.   "clear"
  509. Endif

  510. If (sndspd = -1 | snddir = -1)
  511.   DrawBarb = 0
  512.   DrawHodo = 0
  513.   DrawHeli = 0
  514. Endif

  515. If (snddewp = -1)
  516.   DrawDSnd = 0
  517.   DrawRH   = 0
  518.   DrawPrcl = 1
  519.   DrawPMax = 0
  520.   DrawIndx = 0
  521. Endif

  522. If (sndtemp = -1)
  523.   DrawTSnd = 0
  524.   DrawRH   = 0
  525.   DrawPrcl = 1
  526.   DrawPMax = 0
  527.   DrawIndx = 0
  528.   DrawZLev = 0
  529. Endif

  530. If (NumRing < 1)
  531.   DrawHodo = 0
  532. Endif
  533.   
  534. "q gxinfo"
  535. rec=sublin(result,2)
  536. xsize=subwrd(rec,4)

  537. If (xsize = 11)
  538.    PageType = "Landscape"
  539. Else
  540.    PageType = "Portrait"
  541. Endif
  542.   
  543. *------------------------------------------------------
  544. * calculate constants determining slope/shape of diagram
  545. * based on temp/pressure values given by user
  546. *-------------------------------------------------------

  547. "set x 1"
  548. "set y 1"
  549. "set z 1"
  550. "set t 1"
  551. _m1=(T1+T2-2*T3)/(2*log10(P2/P3))
  552. _m2=(T2-T3-_m1*log10(P2/P3))/50
  553. _m3=(T1-_m1*log10(P1))

  554. "set z "_zmin" "_zmax           
  555. "set zlog on"
  556. "set xlab off"

  557. *-------------------------------------------------
  558. * perform coordinate transformation to Skew-T/LogP
  559. *-------------------------------------------------

  560. "set gxout stat"
  561. "set x "_xval
  562. "set y "_yval
  563. "set t "_tval
  564. "define tempx=("sndtemp"-"_m1"*log10(lev)-"_m3")/"_m2
  565. "define dewpx=("snddewp"-"_m1"*log10(lev)-"_m3")/"_m2

  566. If (PageType = "Portrait")
  567.    "set parea 0.7 7 0.75 10.5"
  568. Else
  569.    "set parea 0.7 6.5 0.5 8"
  570. Endif
  571.   
  572. "set axlim 0 100"
  573. "set lon 0 100"
  574. "set grid on 1 1"

  575. "set z "_zmin " " _zmax
  576. "set lon 0 100"
  577. "set clevs -900"
  578. "set gxout contour"

  579. *-------------------------------------
  580. * Draw pressure lines
  581. *-------------------------------------

  582. If (DrawPLev = 0)
  583.    "set ylab off"
  584. Else
  585.    "set ylab on"
  586.    "set ylopts 1 3 0.10"
  587.    "set xlopts 1 3 0.125"
  588. Endif

  589. "d lon"

  590. *--------------------------------------
  591. * Determine corners of skewt/logp frame
  592. *--------------------------------------

  593. "q w2xy 100 "_pmin
  594. rxloc=subwrd(result,3)
  595. tyloc=subwrd(result,6)
  596. "q w2xy 0 "_pmax
  597. lxloc=subwrd(result,3)
  598. byloc=subwrd(result,6)

  599. If (DrawPLev = 1 & LblAxes = 1)
  600.    "set strsiz 0.10"
  601.    "set string 1 c 3 0"
  602.    If (PageType = "Portrait")
  603.       "draw string 0.5 10.85 hPa"
  604.    Else
  605.       "draw string 0.5 8.35 hPa"
  606.    Endif
  607. Endif

  608. *---------------------------------------------------
  609. * Calculate & draw actual height lines using temp data
  610. *---------------------------------------------------

  611. If (DrawZLev > 0)
  612.    say "Calculating observed height levels from temp/pressure data."
  613.    zz=1
  614.    "set gxout stat"
  615.    "set x "_xval
  616.    "set y "_yval
  617.    "set t "_tval
  618.    count=0
  619.    while (zz < _zmax)
  620.       "set z "zz
  621.       pp.zz=subwrd(result,4)
  622.       lpp.zz=log(pp.zz)
  623.       "d "sndtemp
  624.       rec=sublin(result,8)
  625.       tt=subwrd(rec,4)
  626.       if (tt > -900)
  627.          tk=tt+273.15
  628.          count=count+1
  629.          zzm=zz-1
  630.          If (count = 1)
  631.             If (DrawZLev = 2)
  632.                htlb="ASL"
  633.                height.zz=SfcElev
  634.             Else
  635.                htlb="AGL"
  636.                height.zz=0
  637.             Endif
  638.             sfcz=height.zz
  639.          Else
  640.             DZ=29.2857*(lpp.zzm-lpp.zz)*(lpp.zz*tk+lpp.zzm*tkold)/(lpp.zz+lpp.zzm)
  641.             height.zz=height.zzm+DZ
  642.             highz=height.zz
  643.          Endif
  644.       else
  645.          height.zz = -9999
  646.       endif
  647.       tkold=tk
  648.       zz=zz+1
  649.    endwhile

  650.    maxht=int(highz/1000)
  651.    if (int(sfcz/1000) = sfcz/1000)
  652.       minht=int(sfcz/1000)
  653.    else
  654.       minht=1+int(sfcz/1000)
  655.    endif

  656.    ht=minht
  657.    "set line 1 3 1"
  658.    "set strsiz 0.10"
  659.    "set string 1 l 3 0"
  660.    while (ht <= maxht)
  661.        zz=1
  662.        while (height.zz/1000 <= ht)
  663.           zz=zz+1
  664.        endwhile
  665.        zzm=zz-1
  666.        PBelow=pp.zzm
  667.        PAbove=pp.zz
  668.        HBelow=height.zzm
  669.        HAbove=height.zz
  670.        DZ=HAbove-HBelow
  671.        DP=PAbove-PBelow
  672.        Del=ht*1000-HBelow
  673.        Est=PBelow+Del*DP/DZ
  674.        If (Est >= _pmin & Est <= _pmax)
  675.           "q w2xy 1 " Est
  676.           yloc=subwrd(result,6)
  677.           "draw line " lxloc " " yloc " " rxloc " " yloc
  678.           "draw string 0.22 "yloc-0.05" "ht
  679.        Endif
  680.        ht=ht+1
  681.    endwhile
  682.    "set strsiz 0.10"
  683.    "set string 1"
  684.    If (LblAxes = 1)
  685.       If (PageType = "Portrait")
  686.          "draw string 0.25 10.85 km"
  687.          "draw string 0.25 10.75 "htlb
  688.          "draw string 0.25 10.65 OBS"
  689.       Else
  690.          "draw string 0.25 8.35 km"
  691.          "draw string 0.25 8.25 "htlb
  692.          "draw string 0.25 8.15 OBS"
  693.       Endif
  694.    Endif
  695. Endif


  696. *---------------------------------------------------
  697. * Draw height levels (height above MSL using Std Atm)
  698. *---------------------------------------------------

  699. If (DrawZSTD = 1)
  700.    "set strsiz 0.10"
  701.    minht=30.735*(1-pow(_pmax/1013.26,0.287))
  702.    minht=int(minht+0.5)
  703.    maxht=30.735*(1-pow(_pmin/1013.26,0.287))
  704.    maxht=int(maxht)
  705.    "set gxout stat"
  706.    zcount=minht        
  707.    while (zcount <= maxht)
  708.       plev=1013.26*pow((1-zcount/30.735),3.4843)
  709.       "q w2xy 0 "plev
  710.       yloc=subwrd(result,6)
  711.       "draw string 0 "yloc-0.05" "zcount
  712.       zcount=zcount+1
  713.    endwhile
  714.    "set strsiz 0.10"
  715.    If (LblAxes = 1)
  716.       If (PageType = "Portrait")
  717.          "draw string 0 10.85 km"
  718.          "draw string 0 10.75 ASL"
  719.          "draw string 0 10.65 STD"
  720.       Else
  721.          "draw string 0 8.35 km"
  722.          "draw string 0 8.25 ASL"
  723.          "draw string 0 8.15 STD"
  724.       Endif
  725.   Endif
  726. Endif


  727. *-----------------------
  728. * Plot temperature lines
  729. *-----------------------

  730. If (DrawTemp = 1)
  731.    "set strsiz 0.1"
  732.    "set z "_zmin " " _zmax
  733.    "set line "TempCol " " TempLine " "TempThk
  734.    "set string 1 c 3 0"
  735.    "set gxout stat"
  736.    maxtline=GetTemp(100,_pmax)
  737.    mintline=GetTemp(0,_pmin)

  738.    maxtline=tempint*int(maxtline/tempint)
  739.    mintline=tempint*int(mintline/tempint)

  740.    tloop=mintline
  741.    While (tloop <= maxtline)
  742.        Botxtemp=GetXLoc(tloop,_pmax)
  743.        "q w2xy "Botxtemp " " _pmax
  744.        Botxloc=subwrd(result,3)
  745.        Botyloc=byloc           
  746.        Topxtemp=GetXLoc(tloop,_pmin)
  747.         "q w2xy "Topxtemp " " _pmin
  748.        Topxloc=subwrd(result,3)
  749.        Topyloc=tyloc            
  750.        If (Botxtemp <= 100 | Topxtemp <= 100)
  751.           If (Topxtemp > 100)
  752.              Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)
  753.              b=Topyloc-Slope*Topxtemp
  754.              Topyloc=Slope*100+b
  755.              Topxloc=rxloc         
  756.           Endif
  757.           If (Botxtemp < 0)
  758.              Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)
  759.              b=Botyloc-Slope*Botxtemp
  760.              Botyloc=b
  761.              Botxloc=lxloc
  762.           Else
  763.              "draw string " Botxloc-0.05 " " Botyloc-0.15 " " tloop
  764.           Endif
  765.           "draw line "Botxloc " " Botyloc " " Topxloc " " Topyloc
  766.        Endif
  767.        tloop=tloop+tempint
  768.    EndWhile
  769.    If (LblAxes = 1)
  770.       "set strsiz 0.15"
  771.       "set string 1 c"
  772.       If (PageType = "Portrait")
  773.          "draw string 4.0 0.35 Temperature (`3.`0C)"
  774.       Else
  775.          "draw string 3.5 0.15 Temperature (`3.`0C)"
  776.       Endif
  777.    Endif
  778. Endif


  779. *------------------
  780. * Plot dry adiabats
  781. *------------------

  782. If (DrawThet = 1)
  783.    temp=GetTemp(100,_pmin)
  784.    maxtheta=GetThet2(temp,-100,_pmin)
  785.    maxtheta=thetaint*int(maxtheta/thetaint)
  786.    temp=GetTemp(0,_pmax)
  787.    mintheta=GetThet2(temp,-100,_pmax)
  788.    mintheta=thetaint*int(mintheta/thetaint)
  789.    
  790.    "set lon 0 100"
  791.    "set y 1"
  792.    "set z 1"
  793.    tloop=mintheta
  794.    "set line "ThetCol" "ThetLine " "ThetThk
  795.    While (tloop <= maxtheta)
  796.      PTemp=LiftDry(tloop,1000,_pmin,1,_pmin,_pmax)     
  797.      tloop=tloop+thetaint
  798.    Endwhile
  799. Endif

  800. *------------------------
  801. * Plot mixing ratio lines
  802. *------------------------

  803. If (DrawMix = 1)
  804.    If (MixStop < _pmin)
  805.       MixStop = _pmin
  806.    Endif
  807.    "set string 1 l"
  808.    "set z "_zmin " " _zmax
  809.    "set cint 1"
  810.    "set line "MixCol" " MixLine " "MixThk
  811.    cont = 1
  812.    mloop=subwrd(wsclevs,1)
  813.    count = 1
  814.    While (cont = 1)
  815.        BotCoef=log(mloop*_pmax/3801.66)
  816.        BotTval=-245.5*BotCoef/(BotCoef-17.67)
  817.        Botxtemp=GetXLoc(BotTval,_pmax)
  818.        "q w2xy "Botxtemp " " _pmax
  819.        Botxloc=subwrd(result,3)
  820.        Botyloc=byloc            
  821.        TopCoef=log(mloop*MixStop/3801.66)
  822.        TopTval=-245.5*TopCoef/(TopCoef-17.67)
  823.        Topxtemp=GetXLoc(TopTval,MixStop)
  824.        "q w2xy "Topxtemp " " MixStop
  825.        Topxloc=subwrd(result,3)
  826.        Topyloc=subwrd(result,6)
  827.        "set string "MixCol" l 3"
  828.        "set strsiz 0.09"
  829.        If (Botxtemp <= 100 | Topxtemp <= 100)
  830.           If (Topxtemp > 100)
  831.              Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)
  832.              b=Topyloc-Slope*Topxtemp
  833.              Topyloc=Slope*100+b
  834.              Topxloc=rxloc
  835.              "draw string " Topxloc+0.05 " " Topyloc  " " mloop
  836.           Else
  837.              "draw string " Topxloc " " Topyloc+0.1 " " mloop
  838.           Endif
  839.           If (Botxtemp < 0)
  840.              Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)
  841.              b=Botyloc-Slope*Botxtemp
  842.              Botyloc=b
  843.              Botxloc=lxloc        
  844.           Endif
  845.           "draw line "Botxloc " " Botyloc " " Topxloc " " Topyloc
  846.        Endif
  847.        count=count+1
  848.        mloop=subwrd(wsclevs,count)
  849.        If (mloop = "" | count > 50)
  850.           cont = 0
  851.        Endif
  852.    EndWhile
  853.    If (LblAxes = 1)
  854.       "set strsiz 0.15"
  855.       "set string 1 c 3 90"
  856.       If (PageType = "Portrait")
  857.          "draw string 7.40 4.75 Mixing Ratio (g/kg)"
  858.       Else
  859.          "draw string 6.90 4.25 Mixing Ratio (g/kg)"
  860.       Endif
  861.       "set string 1 c 3 0"
  862.    Endif
  863. Endif

  864. *-----------------------------
  865. * Plot moist (pseudo) adiabats
  866. *-----------------------------

  867. If (DrawThtw = 1)
  868.    "set lon 0 100"
  869.    "set y 1"
  870.    "set z 1"
  871.    "set gxout stat"
  872.    tloop=80
  873.    "set line "ThtwCol" "ThtwLine " "ThtwThk
  874.    While (tloop > -80)
  875.      PTemp=LiftWet(tloop,1000,ThtwStop,1,_pmin,_pmax)     
  876.      tloop=tloop-thetwint
  877.    Endwhile
  878. Endif


  879. *-----------------------------------------------------
  880. * Plot transformed user-specified temperature sounding
  881. *-----------------------------------------------------

  882. If (DrawTSnd = 1)
  883.    say "Drawing temperature sounding."
  884.    "set gxout line"
  885.    "set x "_xval
  886.    "set y "_yval
  887.    "set z "_zmin" "_zmax     
  888.    "set ccolor "TSndCol
  889.    "set cstyle "TSndLine
  890.    "set cmark "TSndMrk
  891.    "set digsiz "MrkSize
  892.    "set cthick "TSndThk
  893.    "set missconn on"
  894.    "d tempx"
  895. Endif

  896. *---------------------------------------------------
  897. * Plot transformed user-specified dewpoint sounding
  898. *---------------------------------------------------

  899. If (DrawDSnd = 1)
  900.    say "Drawing dewpoint sounding."
  901.    "set gxout line"
  902.    "set x "_xval
  903.    "set y "_yval
  904.    "set z "_zmin" "_zmax
  905.    "set cmark "DSndMrk
  906.    "set digsiz "MrkSize
  907.    "set ccolor "DSndCol
  908.    "set cstyle "DSndLine
  909.    "set cthick "DSndThk
  910.    "set missconn on"
  911.    "d dewpx"
  912. Endif

  913. *----------------------------------------
  914. * Determine lowest level of reported  data
  915. *----------------------------------------

  916. If (DrawTSnd = 1)
  917.    field=sndtemp
  918. Else
  919.    field=sndspd
  920. Endif

  921. "set gxout stat"
  922. "set x "_xval
  923. "set y "_yval
  924. "set t "_tval
  925. "set lev " _pmax " " _pmin
  926. "d maskout(lev,"field"+300)"
  927. rec=sublin(result,1)
  928. check=substr(rec,1,6)
  929. If (check = "Notice")
  930.     rec=sublin(result,9)
  931. Else
  932.     rec=sublin(result,8)
  933. Endif
  934. SfcPlev=subwrd(rec,5)

  935. If (DrawTSnd = 1 & DrawDSnd = 1)
  936.    "set lev "SfcPlev
  937.    "d "sndtemp
  938.    rec=sublin(result,8)
  939.    Sfctemp=subwrd(rec,4)
  940.    "d "snddewp
  941.    rec=sublin(result,8)
  942.    Sfcdewp=subwrd(rec,4)
  943.    SfcThee=Thetae(Sfctemp,Sfcdewp,SfcPlev)

  944. *------------------------------------------
  945. * Calculate temperature and pressure of LCL
  946. *------------------------------------------

  947.    TLcl=Templcl(Sfctemp,Sfcdewp)
  948.    PLcl=Preslcl(Sfctemp,Sfcdewp,SfcPlev)
  949. Endif

  950. *----------------------------------------------------------
  951. * Plot parcel path from surface to LCL and up moist adiabat
  952. *----------------------------------------------------------

  953. If (DrawPrcl = 1)
  954.    say "Drawing parcel path from surface upward."
  955.    If (PageType = "Portrait")
  956.       xloc=7.15
  957.    Else
  958.       xloc=6.65
  959.    Endif
  960.    "q w2xy 1 "PLcl
  961.    rec=sublin(result,1)
  962.    yloc=subwrd(rec,6)
  963.    "set strsiz 0.1"
  964.    If (PLcl < _pmax)
  965.       "set string 1 l"
  966.       "draw string "xloc" "yloc" LCL"
  967.       "set line 1 1 1"
  968.       "draw line "xloc-0.15" "yloc" "xloc-0.05" "yloc
  969.    Endif
  970.    "set lon 0 100"
  971.    "set gxout stat"
  972.    "set line "PrclCol" "PrclLine " " PrclThk
  973.    PTemp=LiftDry(Sfctemp,SfcPlev,PLcl,1,_pmin,_pmax)
  974.    Ptemp=LiftWet(TLcl,PLcl,_pmin,1,_pmin,_pmax)
  975. Endif

  976. *-------------------------------------------------------
  977. * Determine level within lowest 250hPa having highest
  978. * theta-e value
  979. *-------------------------------------------------------

  980. If (DrawTSnd = 1 & DrawDSnd = 1)
  981.   "set x "_xval
  982.   "set y "_yval
  983.   "set t "_tval
  984.    zz=1
  985.    MaxThee=-999
  986.    "set gxout stat"
  987.    while (zz <= _zmax & pp > _pmax-250)
  988.        "set z "zz
  989.        pp=subwrd(result,4)
  990.        "d "sndtemp
  991.        rec=sublin(result,8)
  992.        tt=subwrd(rec,4)
  993.        "d "snddewp
  994.        rec=sublin(result,8)
  995.        dd=subwrd(rec,4)
  996.        If (abs(tt) < 130 & abs(dd) < 130)
  997.           Thee=Thetae(tt,dd,pp)
  998.           If (Thee > MaxThee)
  999.              MaxThee=Thee
  1000.              TMaxThee=tt
  1001.              DMaxThee=dd
  1002.              PMaxThee=pp
  1003.           Endif
  1004.        endif
  1005.        zz=zz+1
  1006.    Endwhile
  1007.    If (PMaxThee = SfcPlev-250)
  1008.       PMaxThee = SfcPlev
  1009.    Endif
  1010. *------------------------------------------------------
  1011. * Calculate temperature and pressure of LCL from highest   
  1012. * theta-e level
  1013. *------------------------------------------------------
  1014.    If (SfcPlev != PMaxThee)
  1015.       TLclMax=Templcl(TMaxThee,DMaxThee)
  1016.       PLclMax=Preslcl(TMaxThee,DMaxThee,PMaxThee)
  1017.    Endif
  1018. Endif

  1019. *----------------------------------------------------------
  1020. * Plot parcel path from highest theta-e level to LCL and up
  1021. * moist adiabat
  1022. *----------------------------------------------------------

  1023. If (DrawPMax = 1 & SfcPlev != PMaxThee)
  1024.    say "Drawing parcel path from most unstable level upward."
  1025.    If (PageType = "Portrait")
  1026.       xloc=7.15
  1027.    Else
  1028.       xloc=6.65
  1029.    Endif
  1030.    "q w2xy 1 "PLclMax
  1031.    rec=sublin(result,1)
  1032.    yloc=subwrd(rec,6)
  1033.    "set strsiz 0.1"
  1034.    If (PLclMax < _pmax)
  1035.       "set string 1 l"
  1036.       "draw string "xloc" "yloc" LCL"
  1037.       "set line 1 1 1"
  1038.       "draw line "xloc-0.15" "yloc" "xloc-0.05" "yloc
  1039.    Endif
  1040.    "set lon 0 100"
  1041.    "set gxout stat"
  1042.    "set line "PrclCol" "PrclLine " " PrclThk
  1043.    PTemp=LiftDry(TMaxThee,PMaxThee,PLclMax,1,_pmin,_pmax)
  1044.    Ptemp=LiftWet(TLclMax,PLclMax,_pmin,1,_pmin,_pmax)
  1045. Endif

  1046. *--------------------------------
  1047. * Draw thermodynamic indices
  1048. *--------------------------------

  1049. If (DrawIndx = 1)
  1050.    "set string 1 l"
  1051.    "set strsiz 0.10"
  1052.    "set x "_xval
  1053.    "set y "_yval
  1054.    "set t "_tval
  1055.    say "Calculating precipitable water."
  1056.    pw=precipw(sndtemp,snddewp,_pmax,_pmin)
  1057.    say "Calculating thermodynamic indices."
  1058.    Temp850=interp(sndtemp,850)
  1059.    Temp700=interp(sndtemp,700)
  1060.    Temp500=interp(sndtemp,500)
  1061.    Dewp850=interp(snddewp,850)
  1062.    Dewp700=interp(snddewp,700)
  1063.    Dewp500=interp(snddewp,500)
  1064.    If (Temp850>-900 & Dewp850>-900 & Dewp700>-900 & Temp700>-900 & Temp500>-900)
  1065.       K=Temp850+Dewp850+Dewp700-Temp700-Temp500
  1066.    Else
  1067.       K=-999
  1068.    Endif
  1069.    If (Temp850 > -900 & Dewp850 > -900 & Temp500 > -900)
  1070.       tt=Temp850+Dewp850-2*Temp500
  1071.    Else
  1072.       tt=-999
  1073.    Endif
  1074.    Temp500V=virtual2(Temp500+273.15,Dewp500+273.15,500)-273.15
  1075.    PclTemp=LiftWet(TLcl,PLcl,500,0)
  1076.    PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,500)-273.15
  1077.    SLI=Temp500V-PclTempV
  1078.    rec=CAPE(TLcl,PLcl,100,sndtemp,snddewp)
  1079.    Pos=subwrd(rec,1)
  1080.    CIN=subwrd(rec,2)

  1081.    If (SfcPlev != PMaxThee)
  1082.       PclTemp=LiftWet(TLclMax,PLclMax,500,0)
  1083.       PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,500)-273.15
  1084.       LIMax=Temp500V-PclTempV
  1085.       rec=CAPE(TLclMax,PLclMax,100,sndtemp,snddewp)
  1086.       PosMax=subwrd(rec,1)
  1087.       CINMax=subwrd(rec,2)
  1088.    Else
  1089.       LIMax=SLI
  1090.       PosMax=Pos
  1091.       CINMax=CIN
  1092.       MaxThee=SfcThee
  1093.    Endif

  1094.    If (PageType = "Portrait")
  1095.       If (Text1XC = -1)
  1096.          Text1XC=rxloc-0.75
  1097.       Endif
  1098.       If (Text1YC = -1)
  1099.          Text1YC=tyloc-2.25
  1100.       Endif
  1101.       If (Text2XC = -1)
  1102.          Text2XC=rxloc-0.75
  1103.       Endif
  1104.       If (Text2YC = -1)
  1105.          Text2YC=tyloc-3.25
  1106.       Endif
  1107.       If (Text3XC = -1)
  1108.           Text3XC=rxloc-0.75
  1109.       Endif
  1110.       If (Text3YC = -1)
  1111.          Text3YC=tyloc-4.40
  1112.       Endif
  1113.    Else
  1114.       If (Text1XC = -1)
  1115.          Text1XC=rxloc+2.50
  1116.       Endif
  1117.       If (Text1YC = -1)
  1118.          Text1YC=tyloc-3.00
  1119.       Endif
  1120.       If (Text2XC = -1)
  1121.          Text2XC=rxloc+2.50
  1122.       Endif
  1123.       If (Text2YC = -1)
  1124.          Text2YC=tyloc-4.00
  1125.       Endif
  1126.       If (Text3XC = -1)
  1127.          Text3XC=rxloc+2.50
  1128.       Endif
  1129.       If (Text3YC = -1)
  1130.          Text3YC=tyloc-5.10
  1131.       Endif
  1132.    Endif
  1133.    "set string 1 l 3"
  1134.    "set line 0 1 3"
  1135.    "draw recf  "Text1XC-0.75 " " Text1YC-0.40 " " Text1XC+0.75 " " Text1YC+0.25
  1136.    "set line 1 1 3"
  1137.    "draw rec  "Text1XC-0.75 " " Text1YC-0.40 " " Text1XC+0.75 " " Text1YC+0.25
  1138.    "draw string "Text1XC-0.70 " " Text1YC+0.10"  K"
  1139.    "draw string "Text1XC+0.25 " " Text1YC+0.10" " int(K)      
  1140.    "draw string "Text1XC-0.70 " " Text1YC-0.10 "  TT"
  1141.    "draw string "Text1XC+0.25 " " Text1YC-0.10 " " int(tt)
  1142.    "draw string "Text1XC-0.70 " " Text1YC-0.25 "  PW(cm)"
  1143.    "draw string "Text1XC+0.25 " " Text1YC-0.25 " " int(pw*100)/100
  1144.    "set line 0 1 3"
  1145.    "draw recf  "Text2XC-0.75 " " Text2YC-0.60 " " Text2XC+0.75 " " Text2YC+0.60
  1146.    "set line 1 1 3"
  1147.    "draw rec  "Text2XC-0.75 " " Text2YC-0.60 " " Text2XC+0.75 " " Text2YC+0.60
  1148.    "draw string "Text2XC-0.35 " " Text2YC+0.50 " Surface"
  1149.    "draw string "Text2XC-0.70 " " Text2YC+0.30 "  Temp(`3.`0C)"
  1150.    "draw string "Text2XC+0.25 " " Text2YC+0.30 " " int(Sfctemp*10)/10
  1151.    "draw string "Text2XC-0.70 " " Text2YC+0.15 "  Dewp(`3.`0C)"
  1152.    "draw string "Text2XC+0.25 " " Text2YC+0.15 " " int(Sfcdewp*10)/10
  1153.    "draw string "Text2XC-0.70 " " Text2YC "   `3z`0`bE`n(K)"
  1154.    "draw string "Text2XC+0.25 " " Text2YC " " int(SfcThee)
  1155.    "draw string "Text2XC-0.70 " " Text2YC-0.15 "  LI"
  1156.    "draw string "Text2XC+0.25 " " Text2YC-0.15 " " round(SLI)
  1157.    "draw string "Text2XC-0.70 " " Text2YC-0.30 "  CAPE(J)"
  1158.    "draw string "Text2XC+0.25 " " Text2YC-0.30 " " int(Pos)   
  1159.    "draw string "Text2XC-0.70 " " Text2YC-0.45 "  CIN(J)"
  1160.    "draw string "Text2XC+0.25 " " Text2YC-0.45 " " int(CIN)      
  1161.    "set line 0 1 3"
  1162.    "draw recf  "Text3XC-0.75 " " Text3YC-0.55 " "  Text3XC+0.75 " " Text3YC+0.55
  1163.    "set line 1 1 3"
  1164.    "draw rec  "Text3XC-0.75 " " Text3YC-0.55 " "  Text3XC+0.75 " " Text3YC+0.55
  1165.    "draw string "Text3XC-0.60 " " Text3YC+0.45 " Most Unstable"
  1166.    "draw string "Text3XC-0.70 " " Text3YC+0.20 "  Press(hPa)"
  1167.    "draw string "Text3XC+0.25 " " Text3YC+0.20 " " int(PMaxThee)
  1168.    "draw string "Text3XC-0.70 " " Text3YC+0.05 " `3z`0`bE`n(K)"
  1169.    "draw string "Text3XC+0.25 " " Text3YC+0.05 " " int(MaxThee)
  1170.    "draw string "Text3XC-0.70 " " Text3YC-0.10 " LI"
  1171.    "draw string "Text3XC+0.25 " " Text3YC-0.10 " "round(LIMax)
  1172.    "draw string "Text3XC-0.70 " " Text3YC-0.25 " CAPE(J)"
  1173.    "draw string "Text3XC+0.25 " " Text3YC-0.25 " "int(PosMax)
  1174.    "draw string "Text3XC-0.70 " " Text3YC-0.40 " CIN(J)"
  1175.    "draw string "Text3XC+0.25 " " Text3YC-0.40 " " int(CINMax)
  1176. Endif

  1177. *-----------------------------
  1178. * Draw wind profile along side
  1179. *-----------------------------

  1180. If (DrawBarb = 1)
  1181.    say "Drawing Wind Profile."
  1182.    If (poleloc = -1)
  1183.       If (PageType = "Portrait")
  1184.          poleloc = 8.0
  1185.       Else
  1186.          poleloc = 7.5
  1187.       Endif
  1188.    Endif
  1189.    If (barbline = 1)
  1190.       "set line 1 1 3"
  1191.       "draw line "poleloc " " byloc " " poleloc " " tyloc
  1192.    Endif
  1193.    If (BarbCol = -1)
  1194.       'set rgb 41 255 0 132'
  1195.       'set rgb 42 255 0 168'
  1196.       'set rgb 43 255 0 204'
  1197.       'set rgb 44 255 0 240'
  1198.       'set rgb 45 255 0 255'
  1199.       'set rgb 46 204 0 255'
  1200.       'set rgb 47 174 0 255'
  1201.       'set rgb 48 138 0 255'
  1202.       'set rgb 49 108 0 255'
  1203.       'set rgb 50 84 0 255'
  1204.       'set rgb 51 40 0 255'
  1205.       'set rgb 52 0 0  255'
  1206.       'set rgb 53 0 42 255'
  1207.       'set rgb 54 0 84 255'
  1208.       'set rgb 55 0 120 255'
  1209.       'set rgb 56 0 150 255'
  1210.       'set rgb 57 0 192 255'
  1211.       'set rgb 58 0 240 255'
  1212.       'set rgb 59 0 255 210'
  1213.       'set rgb 60 0 255 160'
  1214.       'set rgb 61 0 255 126'
  1215.       'set rgb 62 0 255 78'
  1216.       'set rgb 63 84 255 0'
  1217.       'set rgb 64 138 255 0'
  1218.       'set rgb 65 188 255 0'
  1219.       'set rgb 66 236 255 0'
  1220.       'set rgb 67 255 255 0'
  1221.       'set rgb 68 255 222 0'
  1222.       'set rgb 69 255 192 0'
  1223.       'set rgb 70 255 162 0'
  1224.       'set rgb 71 255 138 0'
  1225.       'set rgb 72 255 108 0'
  1226.       'set rgb 73 255 84 0'
  1227.       'set rgb 74 255 54 0'
  1228.       'set rgb 75 255 12 0'
  1229.       'set rgb 76 255 0 34'
  1230.       'set rgb 77 255 0 70'
  1231.       'set rgb 78 255 0 105'
  1232.       'set rgb 79 255 0 140'
  1233.       'set rgb 80 255 0 175'
  1234.       'set rgb 81 255 0 215'
  1235.       'set rgb 82 255 0 255'
  1236.       'set rgb 83 255 255 255'

  1237.       col1='83 83 83 83 83 83 83 83 83 83 82 81 80 79 78'
  1238.       col2='77 76 75 74 73 72 71 70 69 68 67 66 65 64 63'
  1239.       col3='62 61 60 59 58 57 56 55 54 53 52 51 50 49 48'
  1240.       'set rbcols 'col1' 'col2' 'col3
  1241.    Endif
  1242.    "set z "_zmin" "_zmax
  1243.    "set gxout stat"
  1244.    zz=1
  1245.    wspd=-999
  1246.    cont=1
  1247.    While (cont = 1 & zz < _zmax)
  1248.       "set z "zz
  1249.       pres=subwrd(result,4)
  1250.       "d "sndspd
  1251.       rec=sublin(result,8)
  1252.       wspd=subwrd(rec,4)
  1253.       if (wspd < 0 | pres > _pmax)
  1254.           zz=zz+1
  1255.       else
  1256.           cont=0
  1257.       Endif
  1258.    Endwhile
  1259.    While (zz <= _zmax)
  1260.       "d "sndspd"(z="zz")"
  1261.       rec=sublin(result,8)
  1262.       wspd=subwrd(rec,4)
  1263.       If (BarbCol >= 0)
  1264.          "set line "BarbCol " 1 "BarbThk
  1265.       Else
  1266.          tempbcol=55+wspd/5     
  1267.          If (tempbcol > 83)
  1268.             tempbcol = 83
  1269.          Endif
  1270.          "set line "tempbcol " 1 "BarbThk
  1271.       Endif
  1272.       "d "snddir"(z="zz")"
  1273.       rec=sublin(result,8)
  1274.       wdir=subwrd(rec,4)
  1275.       xwind=GetUWnd(wspd,wdir)
  1276.       ywind=GetVWnd(wspd,wdir)
  1277.       "query gr2xy 5 "zz
  1278.       y1=subwrd(result,6)
  1279.       if (wspd > 0)
  1280.          cc=polelen/wspd
  1281.          xendpole=poleloc-xwind*cc
  1282.          yendpole=y1-ywind*cc
  1283.       endif
  1284.       if (xendpole>0 & wspd >= 0.5)
  1285.         if (flagbase = 1)
  1286.            "draw mark 3 "poleloc " " y1 " 0.05"
  1287.         endif
  1288.         "draw line " poleloc " " y1 "  " xendpole " " yendpole
  1289.         flagloop=wspd/10
  1290.         windcount=wspd
  1291.         flagcount=0
  1292.         xflagstart=xendpole
  1293.         yflagstart=yendpole
  1294.         dx=cos((180-wdir)*_dtr)
  1295.         dy=sin((180-wdir)*_dtr)
  1296.         while (windcount > 47.5)
  1297.            flagcount=flagcount+1
  1298.            dxflag=-len50*dx
  1299.            dyflag=-len50*dy
  1300.            xflagend=xflagstart+dxflag
  1301.            yflagend=yflagstart+dyflag
  1302.            windcount=windcount-50
  1303.            x1=xflagstart+0.5*wid50*xwind/wspd
  1304.            y1=yflagstart+0.5*wid50*ywind/wspd
  1305.            x2=xflagstart-0.5*wid50*xwind/wspd
  1306.            y2=yflagstart-0.5*wid50*ywind/wspd
  1307.            If (Fill50 = 1)
  1308.               "draw polyf "x1" "y1" "x2" "y2" "xflagend" "yflagend" "x1" "y1
  1309.            Else
  1310.               "draw line "x1 " "y1 " " xflagend " " yflagend " "  
  1311.               "draw line "x2 " "y2 " " xflagend " " yflagend
  1312.               "draw line "x1 " "y1 " " x2 " " y2
  1313.            Endif
  1314.            xflagstart=xflagstart+spac50*xwind/wspd
  1315.            yflagstart=yflagstart+spac50*ywind/wspd
  1316.         endwhile
  1317.         while (windcount > 7.5 )
  1318.            flagcount=flagcount+1
  1319.            dxflag=-len10*dx
  1320.            dyflag=-len10*dy
  1321.            xflagend=xflagstart+dxflag
  1322.            yflagend=yflagstart+dyflag
  1323.            windcount=windcount-10
  1324.            "draw line " xflagstart " " yflagstart " " xflagend " " yflagend
  1325.            xflagstart=xflagstart+spac10*xwind/wspd
  1326.            yflagstart=yflagstart+spac10*ywind/wspd
  1327.         endwhile
  1328.         if (windcount > 2.5)
  1329.            flagcount=flagcount+1
  1330.            if (flagcount = 1)
  1331.               xflagstart=xflagstart+spac05*xwind/wspd
  1332.               yflagstart=yflagstart+spac05*ywind/wspd
  1333.            endif
  1334.            dxflag=-len05*dx
  1335.            dyflag=-len05*dy
  1336.            xflagend=xflagstart+dxflag
  1337.            yflagend=yflagstart+dyflag
  1338.            windcount=windcount-5
  1339.            "draw line " xflagstart " " yflagstart " " xflagend " " yflagend
  1340.         endif
  1341.       else
  1342.         if (wspd < 0.5 & wspd >= 0)
  1343.            "draw mark 2 " poleloc " " y1 " 0.08"
  1344.         endif
  1345.       endif
  1346.       zz=zz+barbint
  1347.    endwhile
  1348. Endif

  1349. *----------------
  1350. * Draw Hodograph
  1351. *----------------

  1352. If (DrawHodo = 1)
  1353.    say "Drawing Hodograph."

  1354.    If (HodXcent = -1 | HodYcent = -1)
  1355.       If (PageType = "Portrait")
  1356.          HodXcent=6
  1357.          HodYcent=9.5
  1358.       Else
  1359.          HodXcent=9
  1360.          HodYcent=7.0
  1361.       Endif
  1362.    Endif
  1363.    HodL=HodXcent-HodSize/2.0
  1364.    HodR=HodXcent+HodSize/2.0
  1365.    HodT=HodYcent+HodSize/2.0
  1366.    HodB=HodYcent-HodSize/2.0
  1367.    RingSpac=HodSize/(NumRing*2)
  1368.    "set line 0"
  1369.    "draw recf "HodL" "HodB" "HodR" "HodT
  1370.    "set line "HodoCol" 1 6"
  1371.    "draw rec "HodL" "HodB" "HodR" "HodT
  1372.    "set line 1 1 3"
  1373.    "set string 1 c"
  1374.    "draw mark 1 "HodXcent " "HodYcent " " HodSize
  1375.    i=1
  1376.    While (i <= NumRing)
  1377.      "set strsiz 0.10"
  1378.      "set string 1 c 3 45"
  1379.      uwnd=-i*HodRing*cos(45*_dtr)
  1380.      xloc=HodXcent+uwnd*RingSpac/HodRing
  1381.      yloc=HodYcent+uwnd*RingSpac/HodRing
  1382.   
  1383.      "draw mark 2 "HodXcent " " HodYcent " " i*HodSize/NumRing
  1384.      "draw string "xloc " " yloc " " HodRing*i
  1385.      i=i+1
  1386.    Endwhile
  1387.    "set string 1 l 3 0"
  1388.    If (TickInt > 0)
  1389.       i=0
  1390.       while (i < HodRing*NumRing)
  1391.          dist=i*RingSpac/HodRing
  1392.          hrxloc=HodXcent+dist               
  1393.          hlxloc=HodXcent-dist                     
  1394.          htyloc=HodYcent+dist
  1395.          hbyloc=HodYcent-dist
  1396.          "set line 1 1 3"
  1397.          "draw line "hrxloc " " HodYcent-TickSize/2 " " hrxloc " " HodYcent+TickSize/2
  1398.          "draw line "hlxloc " " HodYcent-TickSize/2 " " hlxloc " " HodYcent+TickSize/2
  1399.          "draw line "HodXcent+TickSize/2 " " htyloc " " HodXcent-TickSize/2 " " htyloc
  1400.          "draw line "HodXcent+TickSize/2 " " hbyloc " " HodXcent-TickSize/2 " " hbyloc
  1401.          i=i+TickInt
  1402.       endwhile
  1403.    Endif
  1404.    "set line "HodoCol " " HodoLine " "HodoThk
  1405.    "draw string "HodL+0.05 " " HodT-0.1 " knots"
  1406.    zloop=_zmin
  1407.    xold=-999
  1408.    yold=-999
  1409.    count=0
  1410.    Depth=0
  1411.    While (zloop < _zmax & Depth < HodoDep)
  1412.       "set z "zloop
  1413.       pres=subwrd(result,4)
  1414.       "d "sndspd
  1415.       rec=sublin(result,8)
  1416.       wspd=subwrd(rec,4)
  1417.       "d "snddir
  1418.       rec=sublin(result,8)
  1419.       wdir=subwrd(rec,4)
  1420.       uwnd=GetUWnd(wspd,wdir)
  1421.       vwnd=GetVWnd(wspd,wdir)
  1422.       If (wspd >= 0)
  1423.          xloc=HodXcent+uwnd*RingSpac/HodRing
  1424.          yloc=HodYcent+vwnd*RingSpac/HodRing
  1425.          If (xloc > 0 & yloc > 0 & xold > 0 & yold > 0)
  1426.             Depth=Depth+pold-pres
  1427.             count=count+1
  1428.             If (count = 1)
  1429.                "draw mark 3 "xold " " yold " 0.05"
  1430.             Endif
  1431.             "draw line "xold" "yold" "xloc" "yloc
  1432.          Endif
  1433.          xold=xloc
  1434.          yold=yloc
  1435.       Endif
  1436.       zloop=zloop+1
  1437.       pold=pres
  1438.    EndWhile

  1439.    If (count > 0)
  1440.       "draw mark 3 "xold " " yold " 0.05"
  1441.    Endif
  1442. Endif

  1443. *----------------------------------------------
  1444. * Calculate and Display Absolute & S-R Helicity
  1445. *----------------------------------------------

  1446. If (DrawHeli = 1)
  1447.    say "Calculating Helicity & SR Helicity."
  1448.    delp=10
  1449.    UTotal=0
  1450.    VTotal=0

  1451. * First, calculate mass-weighted mean wind
  1452. * Since delp is a constant, and mass is proportional to
  1453. * delp, this is a simple sum.

  1454.    "set lev "_pmax " " _pmin
  1455.    "define uwndarr="sndspd"*cos((270-"snddir")*"_dtr")"
  1456.    "define vwndarr="sndspd"*sin((270-"snddir")*"_dtr")"

  1457.    pres=MeanVBot
  1458.    While (pres >= MeanVTop)
  1459.       uwnd=interp(uwndarr,pres)*_ktm
  1460.       vwnd=interp(vwndarr,pres)*_ktm
  1461.       If (uwnd > -900 & vwnd > -900)
  1462.          UTotal=UTotal+uwnd
  1463.          VTotal=VTotal+vwnd
  1464.       Endif
  1465.       pres=pres-delp
  1466.    EndWhile
  1467.    vcount=1+(MeanVBot-MeanVTop)/delp
  1468.    Umean=UTotal/vcount
  1469.    Vmean=VTotal/vcount
  1470.    Spdmean=GetWSpd(Umean,Vmean)
  1471.    MeanDir=GetWDir(Umean,Vmean)

  1472. * Now, rotate and reduce mean wind to get storm motion

  1473.    If (StormMot = 1)
  1474.       If (Spdmean < 15)
  1475.          Reduct=0.25
  1476.          Rotate=30
  1477.       Else
  1478.          Reduct=0.20
  1479.          Rotate=20
  1480.       Endif
  1481.    Else
  1482.       Reduct=0.0
  1483.       Rotate=0.0
  1484.    Endif

  1485.    UReduce=(1-Reduct)*Umean
  1486.    VReduce=(1-Reduct)*Vmean
  1487.    StormSpd=GetWSpd(UReduce,VReduce)

  1488.    StormDir=GetWDir(UReduce,VReduce)+Rotate
  1489.    If (StormDir >= 360)
  1490.       StormDir=StormDir-360
  1491.    Endif

  1492.    StormU=GetUWnd(StormSpd,StormDir)
  1493.    StormV=GetVWnd(StormSpd,StormDir)

  1494. * Draw Storm Motion Vector

  1495.    xloc=HodXcent+_mtk*StormU*RingSpac/HodRing
  1496.    yloc=HodYcent+_mtk*StormV*RingSpac/HodRing

  1497.    "set line 1 1 4"
  1498.    "draw line "HodXcent " " HodYcent " " xloc " " yloc
  1499.    Arr1U=GetUWnd(HodRing/10,StormDir+30)
  1500.    Arr1V=GetVWnd(HodRing/10,StormDir+30)
  1501.    Arr2U=GetUWnd(HodRing/10,StormDir-30)
  1502.    Arr2V=GetVWnd(HodRing/10,StormDir-30)

  1503.    xloc2=xloc-Arr1U/HodRing
  1504.    xloc3=xloc-Arr2U/HodRing
  1505.    yloc2=yloc-Arr1V/HodRing
  1506.    yloc3=yloc-Arr2V/HodRing

  1507.    "set line 1 1 3"

  1508.    If (FillArrw = 0)
  1509.       "draw line "xloc" "yloc" "xloc2" "yloc2
  1510.       "draw line "xloc" "yloc" "xloc3" "yloc3
  1511.    Else
  1512.       "draw polyf "xloc" "yloc" "xloc2" "yloc2" "xloc3" "yloc3" "xloc" "yloc
  1513.    Endif


  1514. * Now, calculate SR and Environmental Helicity

  1515.    helic=0
  1516.    SRhelic=0
  1517.    MinP=SfcPlev-HelicDep
  1518.    pres=SfcPlev
  1519.    uwndold=-999
  1520.    vwndold=-999
  1521.    While (pres >= MinP)
  1522.       uwnd=interp(uwndarr,pres)*_ktm
  1523.       vwnd=interp(vwndarr,pres)*_ktm
  1524.       If (uwnd > -900 & uwndold > -900)
  1525.           du=uwnd-uwndold
  1526.           dv=vwnd-vwndold
  1527.           ubar=0.5*(uwnd+uwndold)
  1528.           vbar=0.5*(vwnd+vwndold)
  1529.           uhelic=-dv*ubar                  
  1530.           vhelic=du*vbar                  
  1531.           SRuhelic=-dv*(ubar-StormU)
  1532.           SRvhelic=du*(vbar-StormV)
  1533.           SRhelic=SRhelic+SRuhelic+SRvhelic
  1534.           helic=helic+uhelic+vhelic
  1535.       Endif
  1536.       uwndold=uwnd
  1537.       vwndold=vwnd
  1538.       pres=pres-delp
  1539.    EndWhile

  1540.    "set strsiz 0.1"
  1541.    "set string 1 l 3"
  1542.    If (PageType = "Portrait")
  1543.       If (Text4XC = -1)
  1544.          Text4XC=rxloc-0.75
  1545.       Endif
  1546.       If (Text4YC = -1)
  1547.          Text4YC=tyloc-5.45
  1548.       Endif
  1549.    Else
  1550.       If (Text4XC = -1)
  1551.          Text4XC=rxloc+2.50
  1552.       Endif
  1553.       If (Text4YC = -1)
  1554.          Text4YC=tyloc-6.10
  1555.       Endif
  1556.    Endif
  1557.    "set line 0 1 3"
  1558.    "draw recf  "Text4XC-0.75 " "Text4YC-0.5 " " Text4XC+0.75 " " Text4YC+0.5
  1559.    "set line 1 1 3"
  1560.    "draw rec  "Text4XC-0.75 " "Text4YC-0.5 " " Text4XC+0.75 " " Text4YC+0.5
  1561.    "draw string "Text4XC-0.45 " " Text4YC+0.40 " Hodograph"
  1562.    "draw string "Text4XC-0.70 " " Text4YC+0.20 " EH"
  1563.    "draw string "Text4XC+0.25 " " Text4YC+0.20 " "int(helic)
  1564.    "draw string "Text4XC-0.70 " " Text4YC+0.05 " SREH"
  1565.    "draw string "Text4XC+0.25 " " Text4YC+0.05 " " int(SRhelic)
  1566.    "draw string "Text4XC-0.70 " " Text4YC-0.20 " StmDir"
  1567.    "draw string "Text4XC+0.25 " " Text4YC-0.20 " " int(StormDir)"`3.`0"
  1568.    "draw string "Text4XC-0.70 " " Text4YC-0.35 " StmSpd(kt)"
  1569.    "draw string "Text4XC+0.25 " " Text4YC-0.35 " " int(_mtk*StormSpd)
  1570. Endif

  1571. *---------------------------------------------------
  1572. * Plot RH profile.
  1573. *---------------------------------------------------

  1574. If (DrawRH = 1)
  1575.   "set z "_zmin" "_zmax
  1576.   "set x "_xval
  1577.   "set y "_yval
  1578.   "set t "_tval
  1579.   "set zlog on"
  1580.   "set gxout line"
  1581.   "set ccolor "RHCol
  1582.   "set cstyle "RHLine
  1583.   "set cmark "RHMrk
  1584.   "set digsiz "MrkSize
  1585.   "set missconn on"
  1586.   "set xlab on"
  1587.   "set frame off"
  1588.   "set vrange 0 350"
  1589.   "set xlpos 0 t"
  1590.   "set xlevs 25 50 75 100"
  1591.   "set grid vertical 5"
  1592.   "define rh=100*exp((17.2694*"snddewp")/("snddewp"+237.3)-(17.2694*"sndtemp")/("sndtemp"+237.3))"
  1593.   "d rh"
  1594.    If (LblAxes = 1)
  1595.      "set string 1 c 3 0"
  1596.      "set strsiz 0.125"
  1597.      If (PageType = "Portrait")
  1598.        "draw string 1.5 10.85 RH (%)"
  1599.      Else
  1600.        "draw string 1.75 8.35 RH (%)"
  1601.      Endif
  1602.    Endif
  1603. Endif

  1604. *------------------------------------------
  1605. * Reset environment to original dimensions
  1606. *------------------------------------------

  1607. "set t "_tval
  1608. "set x "_xval
  1609. "set y "_yval
  1610. "set z "_zmin " "_zmax

  1611. say "Done."

  1612. Return(0)

  1613. *************************************************************************
  1614. function Templcl(temp,dewp)

  1615. *------------------------------------------------------
  1616. * Calculate the temp at the LCL given temp & dewp in C
  1617. *------------------------------------------------------

  1618. tempk=temp+273.15
  1619. dewpk=dewp+273.15
  1620. Parta=1/(dewpk-56)
  1621. Partb=log(tempk/dewpk)/800
  1622. Tlcl=1/(Parta+Partb)+56
  1623. return(Tlcl-273.15)

  1624. **************************************************************************

  1625. function Preslcl(temp,dewp,pres)

  1626. *-------------------------------------------------------
  1627. * Calculate press of LCL given temp & dewp in C and pressure
  1628. *-------------------------------------------------------

  1629. Tlclk=Templcl(temp,dewp)+273.15
  1630. tempk=temp+273.15
  1631. theta=tempk*pow(1000/pres,0.286)
  1632. plcl=1000*pow(Tlclk/theta,3.48)
  1633. return(plcl)

  1634. **************************************************************************
  1635. function LiftWet(startt,startp,endp,display,Pmin,Pmax)

  1636. *--------------------------------------------------------------------
  1637. * Lift a parcel moist adiabatically from startp to endp.
  1638. * Init temp is startt in C.  If you wish to see the parcel's
  1639. * path plotted, display should be 1.  Returns temp of parcel at endp.
  1640. *--------------------------------------------------------------------

  1641. temp=startt
  1642. pres=startp
  1643. cont = 1
  1644. delp=10
  1645. While (pres >= endp & cont = 1)
  1646.     If (display = 1)
  1647.        xtemp=GetXLoc(temp,pres)
  1648.        "q w2xy "xtemp" "pres
  1649.        xloc=subwrd(result,3)
  1650.        yloc=subwrd(result,6)
  1651.        If (xtemp < 0 | xtemp > 100)
  1652.           cont=0
  1653.        Else
  1654.           If (pres >= Pmin & pres < Pmax & pres < startp)  
  1655.              "draw line "xold" "yold" "xloc" "yloc
  1656.           Endif
  1657.        Endif
  1658.     Endif
  1659.     xold=xloc
  1660.     yold=yloc
  1661.     temp=temp-100*delp*gammaw(temp,pres-delp/2,100)
  1662.     pres=pres-delp
  1663. EndWhile
  1664. return(temp)


  1665. **************************************************************************
  1666. function LiftDry(startt,startp,endp,display,Pmin,Pmax)

  1667. *--------------------------------------------------------------------
  1668. * Lift a parcel dry adiabatically from startp to endp.
  1669. * Init temp is startt in C.  If you wish to see the parcel's
  1670. * path plotted, display should be 1.  Returns temp of parcel at endp.
  1671. *--------------------------------------------------------------------

  1672. starttk=startt+273.15
  1673. cont = 1
  1674. delp=10
  1675. round=int(startp/10)*10
  1676. subscr=0.1*round
  1677. powstart=pow(startp,-0.286)
  1678. temp=starttk*_powpres.subscr*powstart-273.15
  1679. pres=round-10
  1680. While (pres >= endp & cont = 1)
  1681.     subscr=0.1*pres
  1682.     temp=starttk*_powpres.subscr*powstart-273.15
  1683.     If (display = 1)
  1684.        xtemp=GetXLoc(temp,pres)
  1685.        "q w2xy "xtemp" "pres
  1686.        xloc=subwrd(result,3)
  1687.        yloc=subwrd(result,6)
  1688.        If (xtempold > 0 & xtempold < 100 & xtemp > 0 & xtemp < 100)
  1689.           If (pres >= Pmin & pres < Pmax & pres < startp)  
  1690.              "draw line "xold" "yold" "xloc" "yloc
  1691.           Endif
  1692.        Endif
  1693.     Endif
  1694.     xold=xloc
  1695.     xtempold=xtemp
  1696.     yold=yloc
  1697.     pres=pres-delp
  1698. EndWhile
  1699. return(temp)

  1700. **************************************************************************
  1701. function CAPE(startt,startp,endp,sndtemp,snddewp)

  1702. *---------------------------------------------------------------------
  1703. * Returns all postive area and convective inhibition above LCL.
  1704. * Parcel is lifted from LCL at startt,startp and is halted
  1705. * at endp.   Integration method used is trapezoid method.
  1706. *---------------------------------------------------------------------

  1707. pres=startp
  1708. PclTemp=startt
  1709. PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15
  1710. delp=10
  1711. Pos=0
  1712. Neg=0
  1713. Neg2=0

  1714. count=0
  1715. While (pres >= endp)
  1716.    EnvTemp=interp(sndtemp,pres)
  1717.    EnvDewp=interp(snddewp,pres)
  1718.    EnvTempV=virtual2(EnvTemp+273.15,EnvDewp+273.15,pres)-273.15
  1719.    DelT=PclTempV-EnvTempV
  1720.    If (abs(EnvTempV) < 130 & abs(PclTempV) < 130)
  1721.      count=count+1
  1722.      If (count > 1)
  1723.        Val=DelT/pres+Prev
  1724.        If (Val > 0)
  1725.           Pos=Pos+Val
  1726.           Neg2=0
  1727.        Else
  1728.           Neg=Neg+abs(Val)
  1729.           Neg2=Neg2+abs(Val)
  1730.        Endif
  1731.      Endif
  1732.      Prev=DelT/pres
  1733.    Endif
  1734.    pres=pres-delp
  1735.    PclTemp=PclTemp-100*delp*gammaw(PclTemp,pres,100)
  1736.    PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15
  1737. Endwhile

  1738. Pos=0.5*Pos*287*delp
  1739. CIN=0.5*(Neg-Neg2)*287*delp

  1740. return(Pos" "CIN)

  1741. ***************************************************************************
  1742. function gammaw(tempc,pres,rh)

  1743. *-----------------------------------------------------------------------
  1744. * Function to calculate the moist adiabatic lapse rate (deg C/Pa) based
  1745. * on the temperature, pressure, and rh of the environment.
  1746. *----------------------------------------------------------------------

  1747. tempk=tempc+273.15
  1748. es=satvap2(tempc)
  1749. ws=mixratio(es,pres)
  1750. w=rh*ws/100
  1751. tempv=virtual(tempk,w)
  1752. latent=latentc(tempc)

  1753. A=1.0+latent*ws/(287*tempk)
  1754. B=1.0+0.622*latent*latent*ws/(1005*287*tempk*tempk)
  1755. Density=100*pres/(287*tempv)
  1756. lapse=(A/B)/(1005*Density)
  1757. return(lapse)

  1758. *************************************************************************
  1759. function latentc(tempc)

  1760. *-----------------------------------------------------------------------
  1761. * Function to return the latent heat of condensation in J/kg given
  1762. * temperature in degrees Celsius.
  1763. *-----------------------------------------------------------------------

  1764. val=2502.2-2.43089*tempc

  1765. return(val*1000)

  1766. *************************************************************************
  1767. function precipw(sndtemp,snddewp,startp,endp)

  1768. *-----------------------------------------------------------------------
  1769. * Function to calculate the precipitable water (cm) in a sounding
  1770. * starting at pressure level startp and ending at pressure level endp.
  1771. *-----------------------------------------------------------------------

  1772. ppold=-999
  1773. ttold=-999
  1774. ddold=-999
  1775. delp=10
  1776. Int=0
  1777. mix=0
  1778. pres=startp
  1779. logpp=log(pres)
  1780. logppm=log(pres-delp)
  1781. while (pres >= endp)
  1782.    tt=interp(sndtemp,pres)
  1783.    dd=interp(snddewp,pres)
  1784.    if (tt>-900 & ttold>-900 & dd>-900 & ddold>-900)
  1785.       e=satvap2(dd)
  1786.       mix=mixratio(e,pres)
  1787.       mixavg=(logpp*mix+logppm*mixold)/(logpp+logppm)
  1788.       Int=Int+1.020408*mixavg*delp
  1789.    endif
  1790.    ttold=tt
  1791.    ddold=dd
  1792.    ppold=pp
  1793.    mixold=mix
  1794.    pres=pres-delp
  1795.    logpp=logppm
  1796.    logppm=log(pres-delp)
  1797. endwhile

  1798. return(Int)

  1799. *************************************************************************

  1800. function virtual(temp,mix)

  1801. *------------------------------------------------------------
  1802. * Function to return virtual temperature given temperature in
  1803. * kelvin and mixing ratio in g/g.
  1804. *-------------------------------------------------------------

  1805. tempv=temp*(1.0+0.6*mix)

  1806. return (tempv)

  1807. ************************************************************************

  1808. function virtual2(temp,dewp,pres)
  1809.   
  1810. *------------------------------------------------------------
  1811. * Function to return virtual temperature in kelvin given temperature in
  1812. * kelvin and dewpoint in kelvin and pressure in hPa
  1813. *-------------------------------------------------------------

  1814. if (temp > 0 & dewp > 0)
  1815.   vap=satvap2(dewp-273.15)
  1816.   mix=mixratio(vap,pres)
  1817.   tempv=virtual(temp,mix)
  1818. else
  1819.   tempv=-9999
  1820. endif

  1821. return (tempv)
  1822.   
  1823. ************************************************************************

  1824. function satvapor(temp)

  1825. *---------------------------------------------------------------
  1826. * Given temp in Celsius, returns saturation vapor pressure in hPa
  1827. *---------------------------------------------------------------

  1828. pol=_C0+temp*(_C1+temp*(_C2+temp*(_C3+temp*(_C4+temp*(_C5+temp*(_C6+temp*(_C7+temp*(_C8+temp*(_C9)))))))))

  1829. return(6.1078/pow(pol,8))

  1830. ************************************************************************

  1831. function satvap2(temp)

  1832. *---------------------------------------------------------------
  1833. * Given temp in Celsius, returns saturation vapor pressure in hPa*---------------------------------------------------------------

  1834. es=6.112*exp(17.67*temp/(temp+243.5))

  1835. return(es)

  1836. *************************************************************************

  1837. function mixratio(e,p)

  1838. *------------------------------------------------------
  1839. * Given vapor pressure and pressure, return mixing ratio
  1840. *-------------------------------------------------------

  1841. mix=0.622*e/(p-e)

  1842. return(mix)

  1843. ************************************************************************

  1844. function getrh(temp,dewp,pres)

  1845. tempk=temp+273.15
  1846. dewpk=dewp+273.15

  1847. es=satvap2(temp)

  1848. If (temp > 0)
  1849.    A=2.53*pow(10,9)
  1850.    B=5420
  1851. Else
  1852.    A=3.41*pow(10,10)
  1853.    B=6130
  1854. Endif

  1855. w=A*0.622*exp(-B/dewpk)/pres
  1856. ws=mixratio(es,pres)

  1857. return(100*w/ws)

  1858. ************************************************************************
  1859. function interp(array,pres)

  1860. *------------------------------------------------------------------------
  1861. * Interpolate inside array for pressure level pres.
  1862. * Returns estimated value of array at pressure pres.
  1863. *------------------------------------------------------------------------

  1864. "set gxout stat"
  1865. "set lev "pres
  1866. altpres=subwrd(result,4)
  1867. "q dims"
  1868. rec=sublin(result,4)
  1869. zlev=subwrd(rec,9)

  1870. If (zlev < 2 | zlev > _zmaxfile)
  1871.   Vest = -9999.0
  1872. Else
  1873.   If (altpres > pres)
  1874.     zlev=zlev+1
  1875.   Endif
  1876.   "set z "zlev
  1877.   PAbove=subwrd(result,4)
  1878.   "d "array"(lev="PAbove")"
  1879.   rec=sublin(result,8)
  1880.   VAbove=subwrd(rec,4)
  1881.   "set z "zlev-1
  1882.   PBelow=subwrd(result,4)
  1883.   "d "array"(lev="PBelow")"
  1884.   rec=sublin(result,8)
  1885.   VBelow=subwrd(rec,4)

  1886. * Now if we are in a region of missing data, find next good level.

  1887.   If (abs(VAbove) > 130 & zlev > 1 & zlev < _zmaxfile)
  1888.      zz=zlev
  1889.      While (abs(VAbove) > 130 & zz < _zmaxfile)
  1890.        zz=zz+1
  1891.        "set z "zz
  1892.        PAbove=subwrd(result,4)
  1893.        "d "array"(lev="PAbove")"
  1894.        rec=sublin(result,8)
  1895.        VAbove=subwrd(rec,4)
  1896.      EndWhile
  1897.   Endif

  1898.   If (abs(VBelow) > 130 & zlev > 1 & zlev < _zmaxfile)
  1899.       zz=zlev-1
  1900.       While (abs(VBelow) > 130 & zz > 1)
  1901.         zz=zz-1
  1902.         "set z "zz
  1903.         PBelow=subwrd(result,4)
  1904.         "d "array"(lev="PBelow")"
  1905.         rec=sublin(result,8)
  1906.         VBelow=subwrd(rec,4)
  1907.       EndWhile
  1908.   Endif

  1909.   If (abs(VAbove) < 130 & abs(VBelow) < 130)
  1910.      Vest=VBelow+log(PBelow/pres)*(VAbove-VBelow)/log(PBelow/PAbove)
  1911.   Else
  1912.      Vest=-9999.0
  1913.   Endif

  1914. Endif

  1915. Return(Vest)

  1916. ****************************************************************************

  1917. function GetUWnd(wspd,wdir)

  1918. *------------------------
  1919. * Get x-component of wind.
  1920. *------------------------


  1921. If (wspd >= 0)
  1922.    xwind=wspd*cos((270-wdir)*_dtr)
  1923. Else
  1924.    xwind = -9999.0
  1925. Endif
  1926. return(xwind)

  1927. **************************************************************************

  1928. function GetVWnd(wspd,wdir)

  1929. *-----------------------
  1930. * Get y-component of wind
  1931. *------------------------

  1932. If (wspd >= 0)
  1933.    ywind=wspd*sin((270-wdir)*_dtr)
  1934. Else
  1935.    ywind = -9999.0
  1936. Endif
  1937. return(ywind)


  1938. *************************************************************************

  1939. function GetWSpd(xwind,ywind)


  1940. "set gxout stat"
  1941. "d mag("xwind","ywind")"
  1942. rec=sublin(result,8)
  1943. val=subwrd(rec,4)

  1944. return (val)

  1945. *************************************************************************

  1946. function GetWDir(xwind,ywind)

  1947. * Return wind direction given x and y components

  1948. "set gxout stat"
  1949. "define theta=270-"_rtd"*atan2("ywind","xwind")"
  1950. "d theta"
  1951. rec=sublin(result,8)
  1952. Dir=subwrd(rec,4)

  1953. If (Dir > 360)
  1954.    Dir=Dir-360
  1955. Endif

  1956. If (Dir < 0)
  1957.    Dir=360+Dir
  1958. Endif

  1959. return(Dir)

  1960. *************************************************************************

  1961. function GetXLoc(temp,pres)

  1962. *-------------------------------------------------
  1963. * Get x-location on skew-t based on temp, pressure
  1964. *-------------------------------------------------

  1965. xloc=(temp-_m1*log10(pres)-_m3)/_m2
  1966. return(xloc)

  1967. *************************************************************************

  1968. function GetTemp(xloc,pres)

  1969. *-------------------------------------------------
  1970. * Return temperature at location given by xloc,pres
  1971. *-------------------------------------------------

  1972. tempval=_m1*log10(pres)+_m2*xloc+_m3
  1973. return(tempval)

  1974. **************************************************************************

  1975. function GetTheta(temp,pres)         

  1976. *---------------------------------------------------
  1977. * Calculate theta for a given temperature and pressure
  1978. *---------------------------------------------------

  1979. theta=(temp+273.15)*pow(1000/pres,0.286)-273.15
  1980. return(theta)


  1981. *************************************************************************

  1982. function GetThet2(temp,dewp,pres)         

  1983. *---------------------------------------------------
  1984. * Calculate theta for a given temperature,dewp, and pressure
  1985. *---------------------------------------------------

  1986. tempk=273.15+temp
  1987. dewpk=273.15+dewp

  1988. es=satvap2(temp)
  1989. ws=mixratio(es,pres)

  1990. mix=10*getrh(temp,dewp,pres)*ws

  1991. exponent=0.2854*(1.0-0.00028*mix)
  1992. theta=(temp+273.15)*pow(1000/pres,exponent)-273.15
  1993. return(theta)

  1994. *************************************************************************

  1995. function Thetae(temp,dewp,pres)

  1996. *--------------------------------------------------------------
  1997. * Return equiv. pot. temp in Kelvin given temp, dewp in celsius
  1998. * From Bolton (1980) Mon Wea Rev
  1999. *--------------------------------------------------------------

  2000. es=satvap2(temp)
  2001. ws=mixratio(es,pres)
  2002. mix=10*getrh(temp,dewp,pres)*ws
  2003. theta=GetThet2(temp,dewp,pres)+273.15
  2004. TLcl=Templcl(temp,dewp)+273.15
  2005. thetae=theta*exp((3.376/TLcl-0.00254)*mix*(1.0+0.00081*mix))

  2006. return(thetae)

  2007. **************************************************************************


  2008. function int(i0)

  2009. *--------------------------
  2010. * Return integer of i0
  2011. *--------------------------
  2012.   i=0
  2013.   while(i<12)
  2014.     i=i+1
  2015.     if(substr(i0,i,1)='.')
  2016.       i0=substr(i0,1,i-1)
  2017.       break
  2018.     endif
  2019.   endwhile
  2020. return(i0)

  2021. *************************************************************************

  2022. function abs(i)

  2023. *----------------------------
  2024. * return absolute value of i
  2025. *----------------------------

  2026.   if (i < 0)
  2027.      absval=-i
  2028.   else
  2029.      absval=i
  2030.   endif

  2031. return(absval)

  2032. *************************************************************************

  2033. function log(i)

  2034. *---------------------------
  2035. * return natural log of i
  2036. *---------------------------

  2037. "set gxout stat"
  2038. "d log("i")"
  2039. rec=sublin(result,8)
  2040. val=subwrd(rec,4)
  2041. return(val)

  2042. *************************************************************************

  2043. function log10(i)

  2044. *--------------------------
  2045. * return log base 10 of i
  2046. *--------------------------

  2047. "set gxout stat"
  2048. "d log10("i")"
  2049. rec=sublin(result,8)
  2050. val=subwrd(rec,4)
  2051. return(val)

  2052. *************************************************************************

  2053. function pow(i,j)

  2054. *-------------------------------
  2055. * return power of i raised to j
  2056. *-------------------------------

  2057. "set gxout stat"
  2058. "d pow("i","j")"
  2059. rec=sublin(result,8)
  2060. val=subwrd(rec,4)
  2061. return(val)

  2062. ************************************************************************

  2063. function cos(i)

  2064. *-----------------------------------------
  2065. * return cosine of i, where i is in radians
  2066. *------------------------------------------

  2067. "set gxout stat"
  2068. "d cos("i")"
  2069. rec=sublin(result,8)
  2070. val=subwrd(rec,4)
  2071. return(val)

  2072. ************************************************************************

  2073. function sin(i)

  2074. *------------------------------------------
  2075. * return sine of i, where i is in radians
  2076. *------------------------------------------

  2077. "set gxout stat"
  2078. "d sin("i")"
  2079. rec=sublin(result,8)
  2080. val=subwrd(rec,4)
  2081. return(val)

  2082. ************************************************************************

  2083. function exp(i)

  2084. *------------------------------------------
  2085. * return exponential of i
  2086. *------------------------------------------

  2087. "set gxout stat"
  2088. "d exp("i")"
  2089. rec=sublin(result,8)
  2090. val=subwrd(rec,4)
  2091. return(val)

  2092. ***********************************************************************
  2093. function round(i)

  2094. rr=abs(1.0*i)
  2095. rr=int(rr+0.5)
  2096. if (i < 0)
  2097.    rr=-1*rr      
  2098. endif
  2099. return(rr)



复制代码


密码修改失败请联系微信:mofangbao
发表于 2018-3-27 09:30:44 | 显示全部楼层
schliezer 发表于 2014-3-10 15:53
这个变量我画过,也拿EC和NCEP的相关数据用skewplot.gs算过,一对比三者量级相差很大,直接画出来的是最 ...

能把您计算cape的脚本发给我吗 我也是想利于基础数据计算cape 565653979@qq 万分感谢
密码修改失败请联系微信:mofangbao
发表于 2018-4-2 10:58:18 | 显示全部楼层
schliezer 发表于 2014-3-10 15:53
这个变量我画过,也拿EC和NCEP的相关数据用skewplot.gs算过,一对比三者量级相差很大,直接画出来的是最 ...

能把cape的计算方法发给我一份吗,我希望通过微波辐射计计算cape,微波辐射计也是只有基础数据
密码修改失败请联系微信:mofangbao
发表于 2019-10-17 09:23:02 | 显示全部楼层
happyouc 发表于 2018-3-8 09:20
请教 能不能把画cape的程序分享一下 我想拿ec再分析资料画 谢谢

你好 你用EC的计算过cape么 为什么我计算得到的值和EC提供的cape值相差很大
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

发表于 2020-3-24 17:04:48 | 显示全部楼层
求问1°×1°的怎么下载的!我只找到2.5°×2.5°的再分析资料
QQ截图.png
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表