爱气象，爱气象家园！ 
·设为首页 开启辅助访问 充值积分新帖子新回复热门贴精华帖

气象家园

 找回密码
 立即注册

QQ登录

只需一步，快速开始

搜索
气象家园»论坛 › 编程作图 › GrADS › 如何利用GRADS和1*1 Ncep再分析资料画CAPE对流有效位能
123
返回列表 发新帖
楼主: shuihu88888

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

[复制链接]
schliezer
发表于 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)

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

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

  128. ClrScrn = 1

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


  168. P1 = 1000
  169. T1 = -40

  171. P2 = 1000
  172. T2 = 40

  174. P3 = 200
  175. T3 = -50

  177. * Another good set of defining points suggested by George Bryan (PSU)
  178. * are:
  179. *
  180. * P1 = 1000
  181. * T1 = -30
  182. *
  183. * P2 = 1000
  184. * T2 = 40
  185. *
  186. * P3 = 500
  187. * T3 = -18

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

  234. DrawBarb= 1
  235. DrawThet= 1
  236. DrawThtw= 0
  237. DrawTemp= 1
  238. DrawMix = 1
  239. DrawTSnd= 1
  240. DrawDSnd= 1
  241. DrawRH  = 0
  242. DrawPrcl= 1
  243. DrawPMax= 0
  244. DrawIndx= 1
  245. DrawHeli= 1
  246. DrawHodo= 1
  247. DrawPLev= 1
  248. DrawZLev= 0
  249. DrawZSTD= 0
  250. LblAxes = 1

  252. ThtwStop = 200
  253. MixStop  = 600


  256. *
  257. * -----------------  Sounding Geography options ------------------------
  258. *
  259. * SfcElev = Elevation above sea-level (meters) of lowest level reported
  260. *           in sounding.  Used only if DrawZLev = 2

  262. SfcElev = 0


  265. *
  266. * ------------------ Thermodynamic Index Options --------------------
  267. *
  268. * All variables here are in inches.  Use -1 for the default values.
  269. *
  270. *  Text1XC = X-location of midpoint of K,TT,PW output box
  271. *  Text1YC = Y-location of midpoint of K,TT,PW output box
  272. *  Text2XC = X-Location of midpoint of surface indices output box
  273. *  Text2YC = Y-location of midpoint of surface indices output box
  274. *  Text3XC = X-Location of midpoint of most unstable level-based indices
  275. *            output box
  276. *  Text3YC = Y-location of midpoint of most unstable level-based indices
  277. *            output box

  279. Text1XC = -1
  280. Text1YC = -1
  281. Text2XC = -1
  282. Text2YC = -1
  283. Text3XC = -1
  284. Text3YC = -1

  286. *
  287. * ----------------- Wind Barb Profile Options ----------------------------
  288. *
  289. * All variables here are in units of inches, unless otherwise specified
  290. *
  291. *  barbint = Interval for plotting barbs (in units of levels)
  292. *  poleloc = X-Location of profile.  Choose -1 for the default.
  293. *  polelen = Length of wind-barb pole
  294. *  Len05   = Length of each 5-knot barb
  295. *  Len10   = Length of each 10-knot barb
  296. *  Len50   = Length of each 50-knot flag
  297. *  Wid50   = Width of base of 50-knot flag
  298. *  Spac50  = Spacing between 50-knot flag and next barb/flag
  299. *  Spac10  = Spacing between 10-knot flag and next flag
  300. *  Spac05  = Spacing between 5-knot flag and next flag
  301. *  Flagbase= Draw flagbase (filled circle) for each windbarb [1=yes, 0 =no]
  302. *  Fill50  = Solid-fill 50-knot flag [1=yes, 0=no]
  303. *  barbline= Draw a vertical line connecting all the wind barbs [1=yes, 0=no]
  304. *
  305. barbint = 1
  306. poleloc = -1
  307. polelen = 0.35
  308. len05   = 0.07
  309. len10   = 0.15
  310. len50   = 0.15
  311. wid50   = 0.06
  312. spac50  = 0.07
  313. spac10  = 0.05
  314. spac05  = 0.05
  315. Fill50  = 1
  316. flagbase= 1
  317. barbline= 1

  319. *
  320. *
  321. *---------------- Hodograph Options -------------------------------------
  322. *
  323. * All variables here are in units of inches, unless otherwise specified
  324. *
  325. * HodXcent= x-location of hodograph center.  Use -1 for default location.
  326. * HodYcent= y-location of hodograph center.  Use -1 for default location.
  327. * HodSize = Size of hodograph in inches
  328. * NumRing = Number of rings to place in hodograph (must be at least 1)
  329. * HodRing = Wind speed increment of each hodograph ring
  330. * HodoDep = Depth (above lowest level in hPa) of end of hodograph trace
  331. * TickInt = Interval (in kts) at which tick marks are drawn along the axes
  332. *           Use 0 for no tick marks.
  333. * TickSize= Size of tick mark in inches
  334. * Text4XC = X-location of midpoint of hodograph text output. Use -1 for default.
  335. * Text4YC = Y-location of midpoint of hodograph text output. Use -1 for default.

  337. HodXcent= -1
  338. HodYcent= -1
  339. HodSize = 2
  340. NumRing = 3
  341. HodRing = 15
  342. HodoDep = 300
  343. TickInt = 5
  344. TickSize= 0.05
  345. Text4XC = -1
  346. Text4YC = -1

  348. *--------------- Helicity Options ---------------------------------------
  349. *
  350. * MeanVTop = Top pressure level (hPa) of mean-wind calculation
  351. * MeanVBot = Bottom pressure level (hPa) of mean-wind calculation
  352. * HelicDep = Depth in hPa (above ground) of helicity integration
  353. * StormMot = Type of storm motion estimation scheme.  Use following:  
  354. *            0 = No departure from mean wind.
  355. *            1 = Davies-Jones (1990) approach
  356. * FillArrw = Whether to fill the arrowhead of the storm motion vector
  357. *            [1 = yes, 0 = no]

  359. MeanVTop= 300
  360. MeanVBot= 850
  361. HelicDep= 300
  362. StormMot= 1
  363. FillArrw= 1

  365. *
  366. *---------------- Color Options ------------------------------------------
  367. *
  368. * ThetCol = Color of dry adiabats
  369. * TempCol = Color of temperature lines
  370. * MixCol  = Color of mixing ratio lines
  371. * ThtwCol = Color of moist adiabats
  372. * TSndCol = Color of Temperature Sounding
  373. * DSndCol = Color of Dewpoint Sounding
  374. * RHCol   = Color of RH Sounding
  375. * PrclCol = Color of parcel trace
  376. * BarbCol = Color of wind barbs (choose -1 for color according to speed)
  377. * HodoCol = Color of hodograph trace

  379. ThetCol = 2
  380. TempCol = 4
  381. MixCol  = 7
  382. ThtwCol = 3
  383. TSndCol = 1
  384. DSndCol = 1
  385. RHCol   = 3
  386. PrclCol = 5
  387. BarbCol = -1
  388. HodoCol = 1

  390. *
  391. *-------------------- Line Style Options ------------------------------------
  392. *
  393. * GrADS Styles: 1=solid;2=long dash;3=short dash;4=long,short dashed;
  394. *               5=dotted;6=dot dash;7=dot dot dash
  395. *
  396. * ThetLine = Line Style of dry adiabats
  397. * TempLine = Line Style of temperature lines
  398. * MixLine  = Line Style of mixing ratio lines
  399. * ThtwLine = Line Style of moist adiabats
  400. * TSndLine = Line Style of Temperature Sounding
  401. * DSndLine = Line Style of Dewpoint Sounding
  402. * RHLine   = Line Style of RH sounding
  403. * PrclLine = Line Style of parcel trace
  404. * HodoLine = Line Style of hodograph trace
  405. *

  407. ThetLine = 1
  408. TempLine = 1
  409. MixLine  = 5
  410. ThtwLine = 3
  411. TSndLine = 1
  412. DSndLine = 2
  413. RHLine   = 3
  414. PrclLine = 1
  415. HodoLine = 1

  417. *
  418. *------------------- Line Thickness Options---------------------------------
  419. * GrADS Line Thickness: increases with increasing number. Influences
  420. *                       hardcopy output more strongly than screen output.
  421. *
  422. *
  423. * ThetThk = Line Thickness of dry adiabats
  424. * TempThk = Line Thickness of temperature lines
  425. * MixThk  = Line Thickness of mixing ratio lines
  426. * ThtwThk = Line Thickness of moist adiabats
  427. * TSndThk = Line Thickness of temperature sounding
  428. * DSndThk = Line thickness of dewpoint sounding
  429. * RHThk   = Line thickness of RH sounding
  430. * PrclThk = Line thickness of parcel trace
  431. * HodoThk = Line thickness of hodograph trace
  432. * BarbThk = Line thickness of wind barbs

  434. ThetThk = 3
  435. TempThk = 3
  436. MixThk  = 3
  437. ThtwThk = 3
  438. TSndThk = 8
  439. DSndThk = 8
  440. RHThk   = 8
  441. PrclThk = 6
  442. HodoThk = 3
  443. BarbThk = 2

  445. *
  446. *------------------- Data Point Marker Options -----------------------------
  447. * GrADS Marker Types: 0 = none ; 1 = cross ; 2 = open circle ;
  448. *                     3 = closed circle ; 4 = open square ; 5 = closed square
  449. *                     6 = X ; 7 = diamond ; 8 = triangle ; 9 = none
  450. *                    10 = open circle with vertical line ; 11 = open oval
  451. *
  452. * TSndMrk = Mark type of data point marker for temperature sounding
  453. * DSndMrk = Mark type of data point marker for dewpoint sounding
  454. * RHMrk   = Mark type of data point marker for relative humidity sounding
  455. * MrkSize = Mark size (inches) of each data marker

  457. TSndMrk = 0
  458. DSndMrk = 0
  459. RHMrk   = 0
  460. MrkSize = 0.1


  463. * !!!!! YOU SHOULD NOT NEED TO CHANGE ANYTHING BELOW HERE !!!!!
  464. ****************************************************************************

  466. *-------------------------------------------
  467. * grab user-specified environment dimensions
  468. *-------------------------------------------

  470. "q dims"
  471. rec=sublin(result,2)
  472. _xtype=subwrd(rec,3)
  473. _xval=subwrd(rec,9)
  474. rec=sublin(result,3)
  475. _yval=subwrd(rec,9)
  476. _ytype=subwrd(rec,3)
  477. rec=sublin(result,4)
  478. _ptype=subwrd(rec,3)
  479. _pmax=subwrd(rec,6)
  480. _pmin=subwrd(rec,8)
  481. _zmin=subwrd(rec,11)
  482. _zmax=subwrd(rec,13)
  483. rec=sublin(result,5)
  484. _ttype=subwrd(rec,3)
  485. _tval=subwrd(rec,9)

  487. "q file"
  488. rec=sublin(result,5)
  489. _zmaxfile=subwrd(rec,9)

  491. *-------------------------------------------------------------
  492. * Check to ensure that dimensions are valid.  Warn & exit if not.
  493. *--------------------------------------------------------------

  495. dimrc=0
  496. If (_xtype != "fixed")
  497.   say "X-Dims Error:  Not fixed.  Use 'set lon' or 'set x' to specify a value."
  498.   dimrc=-1
  499. Endif

  501. If (_ytype != "fixed")
  502.   say "Y-Dims Error:  Not fixed.  Use 'set lat' or 'set y' to specify a value"
  503.   dimrc=-1
  504. Endif

  506. If (_ptype != "varying")
  507.    say "Z-Dims Error:  Not varying.  Use 'set lev' or 'set z' to specify a range."
  508.    dimrc=-1
  509. Endif

  511. If (_ttype != "fixed")
  512.   say "Time Error:     Not fixed.  Use 'set time' or 'set t' to specify a value"
  513.   dimrc=-1
  514. Endif


  517. If (dimrc < 0)
  518.   Return(-1)
  519. Endif


  522. *
  523. * A few global variables used in units conversion
  524. *

  526. _pi=3.14159265
  527. _dtr=_pi/180
  528. _rtd=1/_dtr
  529. _ktm=0.514444
  530. _mtk=1/_ktm

  532. * A few global constants used in thermo calcs

  534. _C0=0.99999683
  535. _C1=-0.90826951/100
  536. _C2= 0.78736169/10000
  537. _C3=-0.61117958/1000000
  538. _C4= 0.43884187/pow(10,8)
  539. _C5=-0.29883885/pow(10,10)
  540. _C6= 0.21874425/pow(10,12)
  541. _C7=-0.17892321/pow(10,14)
  542. _C8= 0.11112018/pow(10,16)         
  543. _C9=-0.30994571/pow(10,19)

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

  548. zz=1100
  549. while (zz > 10)
  550.     subscr=zz/10
  551.     _powpres.subscr=pow(zz,0.286)
  552.     zz=zz-10
  553. endwhile

  555. *
  556. * Turn off options not available due to user data limitations
  557. *

  559. If (ClrScrn = 1)
  560.   "clear"
  561. Endif

  563. If (sndspd = -1 | snddir = -1)
  564.   DrawBarb = 0
  565.   DrawHodo = 0
  566.   DrawHeli = 0
  567. Endif

  569. If (snddewp = -1)
  570.   DrawDSnd = 0
  571.   DrawRH   = 0
  572.   DrawPrcl = 1
  573.   DrawPMax = 0
  574.   DrawIndx = 0
  575. Endif

  577. If (sndtemp = -1)
  578.   DrawTSnd = 0
  579.   DrawRH   = 0
  580.   DrawPrcl = 1
  581.   DrawPMax = 0
  582.   DrawIndx = 0
  583.   DrawZLev = 0
  584. Endif

  586. If (NumRing < 1)
  587.   DrawHodo = 0
  588. Endif
  589.   
  590. "q gxinfo"
  591. rec=sublin(result,2)
  592. xsize=subwrd(rec,4)

  594. If (xsize = 11)
  595.    PageType = "Landscape"
  596. Else
  597.    PageType = "Portrait"
  598. Endif
  599.   
  600. *------------------------------------------------------
  601. * calculate constants determining slope/shape of diagram
  602. * based on temp/pressure values given by user
  603. *-------------------------------------------------------

  605. "set x 1"
  606. "set y 1"
  607. "set z 1"
  608. "set t 1"
  609. _m1=(T1+T2-2*T3)/(2*log10(P2/P3))
  610. _m2=(T2-T3-_m1*log10(P2/P3))/50
  611. _m3=(T1-_m1*log10(P1))

  613. "set z "_zmin" "_zmax           
  614. "set zlog on"
  615. "set xlab off"

  617. *-------------------------------------------------
  618. * perform coordinate transformation to Skew-T/LogP
  619. *-------------------------------------------------

  621. "set gxout stat"
  622. "set x "_xval
  623. "set y "_yval
  624. "set t "_tval
  625. "define tempx=("sndtemp"-"_m1"*log10(lev)-"_m3")/"_m2
  626. "define dewpx=("snddewp"-"_m1"*log10(lev)-"_m3")/"_m2

  628. If (PageType = "Portrait")
  629.    "set parea 0.7 7 0.75 10.5"
  630. Else
  631.    "set parea 0.7 6.5 0.5 8"
  632. Endif
  633.   
  634. "set axlim 0 100"
  635. "set lon 0 100"
  636. "set grid on 1 1"

  638. "set z "_zmin " " _zmax
  639. "set lon 0 100"
  640. "set clevs -900"
  641. "set gxout contour"

  643. *-------------------------------------
  644. * Draw pressure lines
  645. *-------------------------------------

  647. If (DrawPLev = 0)
  648.    "set ylab off"
  649. Else
  650.    "set ylab on"
  651.    "set ylopts 1 3 0.10"
  652.    "set xlopts 1 3 0.125"
  653. Endif

  655. "d lon"

  657. *--------------------------------------
  658. * Determine corners of skewt/logp frame
  659. *--------------------------------------

  661. "q w2xy 100 "_pmin
  662. rxloc=subwrd(result,3)
  663. tyloc=subwrd(result,6)
  664. "q w2xy 0 "_pmax
  665. lxloc=subwrd(result,3)
  666. byloc=subwrd(result,6)

  668. If (DrawPLev = 1 & LblAxes = 1)
  669.    "set strsiz 0.10"
  670.    "set string 1 c 3 0"
  671.    If (PageType = "Portrait")
  672.       "draw string 0.5 10.85 hPa"
  673.    Else
  674.       "draw string 0.5 8.35 hPa"
  675.    Endif
  676. Endif

  678. *---------------------------------------------------
  679. * Calculate & draw actual height lines using temp data
  680. *---------------------------------------------------

  682. If (DrawZLev > 0)
  683.    say "Calculating observed height levels from temp/pressure data."
  684.    zz=1
  685.    "set gxout stat"
  686.    "set x "_xval
  687.    "set y "_yval
  688.    "set t "_tval
  689.    count=0
  690.    while (zz < _zmax)
  691.       "set z "zz
  692.       pp.zz=subwrd(result,4)
  693.       lpp.zz=log(pp.zz)
  694.       "d "sndtemp
  695.       rec=sublin(result,8)
  696.       tt=subwrd(rec,4)
  697.       if (tt > -900)
  698.          tk=tt+273.15
  699.          count=count+1
  700.          zzm=zz-1
  701.          If (count = 1)
  702.             If (DrawZLev = 2)
  703.                htlb="ASL"
  704.                height.zz=SfcElev
  705.             Else
  706.                htlb="AGL"
  707.                height.zz=0
  708.             Endif
  709.             sfcz=height.zz
  710.          Else
  711.             DZ=29.2857*(lpp.zzm-lpp.zz)*(lpp.zz*tk+lpp.zzm*tkold)/(lpp.zz+lpp.zzm)
  712.             height.zz=height.zzm+DZ
  713.             highz=height.zz
  714.          Endif
  715.       else
  716.          height.zz = -9999
  717.       endif
  718.       tkold=tk
  719.       zz=zz+1
  720.    endwhile

  722.    maxht=int(highz/1000)
  723.    if (int(sfcz/1000) = sfcz/1000)
  724.       minht=int(sfcz/1000)
  725.    else
  726.       minht=1+int(sfcz/1000)
  727.    endif

  729.    ht=minht
  730.    "set line 1 3 1"
  731.    "set strsiz 0.10"
  732.    "set string 1 l 3 0"
  733.    while (ht <= maxht)
  734.        zz=1
  735.        while (height.zz/1000 <= ht)
  736.           zz=zz+1
  737.        endwhile
  738.        zzm=zz-1
  739.        PBelow=pp.zzm
  740.        PAbove=pp.zz
  741.        HBelow=height.zzm
  742.        HAbove=height.zz
  743.        DZ=HAbove-HBelow
  744.        DP=PAbove-PBelow
  745.        Del=ht*1000-HBelow
  746.        Est=PBelow+Del*DP/DZ
  747.        If (Est >= _pmin & Est <= _pmax)
  748.           "q w2xy 1 " Est
  749.           yloc=subwrd(result,6)
  750.           "draw line " lxloc " " yloc " " rxloc " " yloc
  751.           "draw string 0.22 "yloc-0.05" "ht
  752.        Endif
  753.        ht=ht+1
  754.    endwhile
  755.    "set strsiz 0.10"
  756.    "set string 1"
  757.    If (LblAxes = 1)
  758.       If (PageType = "Portrait")
  759.          "draw string 0.25 10.85 km"
  760.          "draw string 0.25 10.75 "htlb
  761.          "draw string 0.25 10.65 OBS"
  762.       Else
  763.          "draw string 0.25 8.35 km"
  764.          "draw string 0.25 8.25 "htlb
  765.          "draw string 0.25 8.15 OBS"
  766.       Endif
  767.    Endif
  768. Endif


  771. *---------------------------------------------------
  772. * Draw height levels (height above MSL using Std Atm)
  773. *---------------------------------------------------

  775. If (DrawZSTD = 1)
  776.    "set strsiz 0.10"
  777.    minht=30.735*(1-pow(_pmax/1013.26,0.287))
  778.    minht=int(minht+0.5)
  779.    maxht=30.735*(1-pow(_pmin/1013.26,0.287))
  780.    maxht=int(maxht)
  781.    "set gxout stat"
  782.    zcount=minht        
  783.    while (zcount <= maxht)
  784.       plev=1013.26*pow((1-zcount/30.735),3.4843)
  785.       "q w2xy 0 "plev
  786.       yloc=subwrd(result,6)
  787.       "draw string 0 "yloc-0.05" "zcount
  788.       zcount=zcount+1
  789.    endwhile
  790.    "set strsiz 0.10"
  791.    If (LblAxes = 1)
  792.       If (PageType = "Portrait")
  793.          "draw string 0 10.85 km"
  794.          "draw string 0 10.75 ASL"
  795.          "draw string 0 10.65 STD"
  796.       Else
  797.          "draw string 0 8.35 km"
  798.          "draw string 0 8.25 ASL"
  799.          "draw string 0 8.15 STD"
  800.       Endif
  801.   Endif
  802. Endif


  805. *-----------------------
  806. * Plot temperature lines
  807. *-----------------------

  809. If (DrawTemp = 1)
  810.    "set strsiz 0.1"
  811.    "set z "_zmin " " _zmax
  812.    "set line "TempCol " " TempLine " "TempThk
  813.    "set string 1 c 3 0"
  814.    "set gxout stat"
  815.    maxtline=GetTemp(100,_pmax)
  816.    mintline=GetTemp(0,_pmin)

  818.    maxtline=tempint*int(maxtline/tempint)
  819.    mintline=tempint*int(mintline/tempint)

  821.    tloop=mintline
  822.    While (tloop <= maxtline)
  823.        Botxtemp=GetXLoc(tloop,_pmax)
  824.        "q w2xy "Botxtemp " " _pmax
  825.        Botxloc=subwrd(result,3)
  826.        Botyloc=byloc           
  827.        Topxtemp=GetXLoc(tloop,_pmin)
  828.         "q w2xy "Topxtemp " " _pmin
  829.        Topxloc=subwrd(result,3)
  830.        Topyloc=tyloc            
  831.        If (Botxtemp <= 100 | Topxtemp <= 100)
  832.           If (Topxtemp > 100)
  833.              Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)
  834.              b=Topyloc-Slope*Topxtemp
  835.              Topyloc=Slope*100+b
  836.              Topxloc=rxloc         
  837.           Endif
  838.           If (Botxtemp < 0)
  839.              Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)
  840.              b=Botyloc-Slope*Botxtemp
  841.              Botyloc=b
  842.              Botxloc=lxloc
  843.           Else
  844.              "draw string " Botxloc-0.05 " " Botyloc-0.15 " " tloop
  845.           Endif
  846.           "draw line "Botxloc " " Botyloc " " Topxloc " " Topyloc
  847.        Endif
  848.        tloop=tloop+tempint
  849.    EndWhile
  850.    If (LblAxes = 1)
  851.       "set strsiz 0.15"
  852.       "set string 1 c"
  853.       If (PageType = "Portrait")
  854.          "draw string 4.0 0.35 Temperature (`3.`0C)"
  855.       Else
  856.          "draw string 3.5 0.15 Temperature (`3.`0C)"
  857.       Endif
  858.    Endif
  859. Endif


  862. *------------------
  863. * Plot dry adiabats
  864. *------------------

  866. If (DrawThet = 1)
  867.    temp=GetTemp(100,_pmin)
  868.    maxtheta=GetThet2(temp,-100,_pmin)
  869.    maxtheta=thetaint*int(maxtheta/thetaint)
  870.    temp=GetTemp(0,_pmax)
  871.    mintheta=GetThet2(temp,-100,_pmax)
  872.    mintheta=thetaint*int(mintheta/thetaint)
  873.    
  874.    "set lon 0 100"
  875.    "set y 1"
  876.    "set z 1"
  877.    tloop=mintheta
  878.    "set line "ThetCol" "ThetLine " "ThetThk
  879.    While (tloop <= maxtheta)
  880.      PTemp=LiftDry(tloop,1000,_pmin,1,_pmin,_pmax)     
  881.      tloop=tloop+thetaint
  882.    Endwhile
  883. Endif

  885. *------------------------
  886. * Plot mixing ratio lines
  887. *------------------------

  889. If (DrawMix = 1)
  890.    If (MixStop < _pmin)
  891.       MixStop = _pmin
  892.    Endif
  893.    "set string 1 l"
  894.    "set z "_zmin " " _zmax
  895.    "set cint 1"
  896.    "set line "MixCol" " MixLine " "MixThk
  897.    cont = 1
  898.    mloop=subwrd(wsclevs,1)
  899.    count = 1
  900.    While (cont = 1)
  901.        BotCoef=log(mloop*_pmax/3801.66)
  902.        BotTval=-245.5*BotCoef/(BotCoef-17.67)
  903.        Botxtemp=GetXLoc(BotTval,_pmax)
  904.        "q w2xy "Botxtemp " " _pmax
  905.        Botxloc=subwrd(result,3)
  906.        Botyloc=byloc            
  907.        TopCoef=log(mloop*MixStop/3801.66)
  908.        TopTval=-245.5*TopCoef/(TopCoef-17.67)
  909.        Topxtemp=GetXLoc(TopTval,MixStop)
  910.        "q w2xy "Topxtemp " " MixStop
  911.        Topxloc=subwrd(result,3)
  912.        Topyloc=subwrd(result,6)
  913.        "set string "MixCol" l 3"
  914.        "set strsiz 0.09"
  915.        If (Botxtemp <= 100 | Topxtemp <= 100)
  916.           If (Topxtemp > 100)
  917.              Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)
  918.              b=Topyloc-Slope*Topxtemp
  919.              Topyloc=Slope*100+b
  920.              Topxloc=rxloc
  921.              "draw string " Topxloc+0.05 " " Topyloc  " " mloop
  922.           Else
  923.              "draw string " Topxloc " " Topyloc+0.1 " " mloop
  924.           Endif
  925.           If (Botxtemp < 0)
  926.              Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)
  927.              b=Botyloc-Slope*Botxtemp
  928.              Botyloc=b
  929.              Botxloc=lxloc        
  930.           Endif
  931.           "draw line "Botxloc " " Botyloc " " Topxloc " " Topyloc
  932.        Endif
  933.        count=count+1
  934.        mloop=subwrd(wsclevs,count)
  935.        If (mloop = "" | count > 50)
  936.           cont = 0
  937.        Endif
  938.    EndWhile
  939.    If (LblAxes = 1)
  940.       "set strsiz 0.15"
  941.       "set string 1 c 3 90"
  942.       If (PageType = "Portrait")
  943.          "draw string 7.40 4.75 Mixing Ratio (g/kg)"
  944.       Else
  945.          "draw string 6.90 4.25 Mixing Ratio (g/kg)"
  946.       Endif
  947.       "set string 1 c 3 0"
  948.    Endif
  949. Endif

  951. *-----------------------------
  952. * Plot moist (pseudo) adiabats
  953. *-----------------------------

  955. If (DrawThtw = 1)
  956.    "set lon 0 100"
  957.    "set y 1"
  958.    "set z 1"
  959.    "set gxout stat"
  960.    tloop=80
  961.    "set line "ThtwCol" "ThtwLine " "ThtwThk
  962.    While (tloop > -80)
  963.      PTemp=LiftWet(tloop,1000,ThtwStop,1,_pmin,_pmax)     
  964.      tloop=tloop-thetwint
  965.    Endwhile
  966. Endif


  969. *-----------------------------------------------------
  970. * Plot transformed user-specified temperature sounding
  971. *-----------------------------------------------------

  973. If (DrawTSnd = 1)
  974.    say "Drawing temperature sounding."
  975.    "set gxout line"
  976.    "set x "_xval
  977.    "set y "_yval
  978.    "set z "_zmin" "_zmax     
  979.    "set ccolor "TSndCol
  980.    "set cstyle "TSndLine
  981.    "set cmark "TSndMrk
  982.    "set digsiz "MrkSize
  983.    "set cthick "TSndThk
  984.    "set missconn on"
  985.    "d tempx"
  986. Endif

  988. *---------------------------------------------------
  989. * Plot transformed user-specified dewpoint sounding
  990. *---------------------------------------------------

  992. If (DrawDSnd = 1)
  993.    say "Drawing dewpoint sounding."
  994.    "set gxout line"
  995.    "set x "_xval
  996.    "set y "_yval
  997.    "set z "_zmin" "_zmax
  998.    "set cmark "DSndMrk
  999.    "set digsiz "MrkSize
  1000.    "set ccolor "DSndCol
  1001.    "set cstyle "DSndLine
  1002.    "set cthick "DSndThk
  1003.    "set missconn on"
  1004.    "d dewpx"
  1005. Endif

  1007. *----------------------------------------
  1008. * Determine lowest level of reported  data
  1009. *----------------------------------------

  1011. If (DrawTSnd = 1)
  1012.    field=sndtemp
  1013. Else
  1014.    field=sndspd
  1015. Endif

  1017. "set gxout stat"
  1018. "set x "_xval
  1019. "set y "_yval
  1020. "set t "_tval
  1021. "set lev " _pmax " " _pmin
  1022. "d maskout(lev,"field"+300)"
  1023. rec=sublin(result,1)
  1024. check=substr(rec,1,6)
  1025. If (check = "Notice")
  1026.     rec=sublin(result,9)
  1027. Else
  1028.     rec=sublin(result,8)
  1029. Endif
  1030. SfcPlev=subwrd(rec,5)

  1032. If (DrawTSnd = 1 & DrawDSnd = 1)
  1033.    "set lev "SfcPlev
  1034.    "d "sndtemp
  1035.    rec=sublin(result,8)
  1036.    Sfctemp=subwrd(rec,4)
  1037.    "d "snddewp
  1038.    rec=sublin(result,8)
  1039.    Sfcdewp=subwrd(rec,4)
  1040.    SfcThee=Thetae(Sfctemp,Sfcdewp,SfcPlev)

  1042. *------------------------------------------
  1043. * Calculate temperature and pressure of LCL
  1044. *------------------------------------------

  1046.    TLcl=Templcl(Sfctemp,Sfcdewp)
  1047.    PLcl=Preslcl(Sfctemp,Sfcdewp,SfcPlev)
  1048. Endif

  1050. *----------------------------------------------------------
  1051. * Plot parcel path from surface to LCL and up moist adiabat
  1052. *----------------------------------------------------------

  1054. If (DrawPrcl = 1)
  1055.    say "Drawing parcel path from surface upward."
  1056.    If (PageType = "Portrait")
  1057.       xloc=7.15
  1058.    Else
  1059.       xloc=6.65
  1060.    Endif
  1061.    "q w2xy 1 "PLcl
  1062.    rec=sublin(result,1)
  1063.    yloc=subwrd(rec,6)
  1064.    "set strsiz 0.1"
  1065.    If (PLcl < _pmax)
  1066.       "set string 1 l"
  1067.       "draw string "xloc" "yloc" LCL"
  1068.       "set line 1 1 1"
  1069.       "draw line "xloc-0.15" "yloc" "xloc-0.05" "yloc
  1070.    Endif
  1071.    "set lon 0 100"
  1072.    "set gxout stat"
  1073.    "set line "PrclCol" "PrclLine " " PrclThk
  1074.    PTemp=LiftDry(Sfctemp,SfcPlev,PLcl,1,_pmin,_pmax)
  1075.    Ptemp=LiftWet(TLcl,PLcl,_pmin,1,_pmin,_pmax)
  1076. Endif

  1078. *-------------------------------------------------------
  1079. * Determine level within lowest 250hPa having highest
  1080. * theta-e value
  1081. *-------------------------------------------------------

  1083. If (DrawTSnd = 1 & DrawDSnd = 1)
  1084.   "set x "_xval
  1085.   "set y "_yval
  1086.   "set t "_tval
  1087.    zz=1
  1088.    MaxThee=-999
  1089.    "set gxout stat"
  1090.    while (zz <= _zmax & pp > _pmax-250)
  1091.        "set z "zz
  1092.        pp=subwrd(result,4)
  1093.        "d "sndtemp
  1094.        rec=sublin(result,8)
  1095.        tt=subwrd(rec,4)
  1096.        "d "snddewp
  1097.        rec=sublin(result,8)
  1098.        dd=subwrd(rec,4)
  1099.        If (abs(tt) < 130 & abs(dd) < 130)
  1100.           Thee=Thetae(tt,dd,pp)
  1101.           If (Thee > MaxThee)
  1102.              MaxThee=Thee
  1103.              TMaxThee=tt
  1104.              DMaxThee=dd
  1105.              PMaxThee=pp
  1106.           Endif
  1107.        endif
  1108.        zz=zz+1
  1109.    Endwhile
  1110.    If (PMaxThee = SfcPlev-250)
  1111.       PMaxThee = SfcPlev
  1112.    Endif
  1113. *------------------------------------------------------
  1114. * Calculate temperature and pressure of LCL from highest   
  1115. * theta-e level
  1116. *------------------------------------------------------
  1117.    If (SfcPlev != PMaxThee)
  1118.       TLclMax=Templcl(TMaxThee,DMaxThee)
  1119.       PLclMax=Preslcl(TMaxThee,DMaxThee,PMaxThee)
  1120.    Endif
  1121. Endif

  1123. *----------------------------------------------------------
  1124. * Plot parcel path from highest theta-e level to LCL and up
  1125. * moist adiabat
  1126. *----------------------------------------------------------

  1128. If (DrawPMax = 1 & SfcPlev != PMaxThee)
  1129.    say "Drawing parcel path from most unstable level upward."
  1130.    If (PageType = "Portrait")
  1131.       xloc=7.15
  1132.    Else
  1133.       xloc=6.65
  1134.    Endif
  1135.    "q w2xy 1 "PLclMax
  1136.    rec=sublin(result,1)
  1137.    yloc=subwrd(rec,6)
  1138.    "set strsiz 0.1"
  1139.    If (PLclMax < _pmax)
  1140.       "set string 1 l"
  1141.       "draw string "xloc" "yloc" LCL"
  1142.       "set line 1 1 1"
  1143.       "draw line "xloc-0.15" "yloc" "xloc-0.05" "yloc
  1144.    Endif
  1145.    "set lon 0 100"
  1146.    "set gxout stat"
  1147.    "set line "PrclCol" "PrclLine " " PrclThk
  1148.    PTemp=LiftDry(TMaxThee,PMaxThee,PLclMax,1,_pmin,_pmax)
  1149.    Ptemp=LiftWet(TLclMax,PLclMax,_pmin,1,_pmin,_pmax)
  1150. Endif

  1152. *--------------------------------
  1153. * Draw thermodynamic indices
  1154. *--------------------------------

  1156. If (DrawIndx = 1)
  1157.    "set string 1 l"
  1158.    "set strsiz 0.10"
  1159.    "set x "_xval
  1160.    "set y "_yval
  1161.    "set t "_tval
  1162.    say "Calculating precipitable water."
  1163.    pw=precipw(sndtemp,snddewp,_pmax,_pmin)
  1164.    say "Calculating thermodynamic indices."
  1165.    Temp850=interp(sndtemp,850)
  1166.    Temp700=interp(sndtemp,700)
  1167.    Temp500=interp(sndtemp,500)
  1168.    Dewp850=interp(snddewp,850)
  1169.    Dewp700=interp(snddewp,700)
  1170.    Dewp500=interp(snddewp,500)
  1171.    If (Temp850>-900 & Dewp850>-900 & Dewp700>-900 & Temp700>-900 & Temp500>-900)
  1172.       K=Temp850+Dewp850+Dewp700-Temp700-Temp500
  1173.    Else
  1174.       K=-999
  1175.    Endif
  1176.    If (Temp850 > -900 & Dewp850 > -900 & Temp500 > -900)
  1177.       tt=Temp850+Dewp850-2*Temp500
  1178.    Else
  1179.       tt=-999
  1180.    Endif
  1181.    Temp500V=virtual2(Temp500+273.15,Dewp500+273.15,500)-273.15
  1182.    PclTemp=LiftWet(TLcl,PLcl,500,0)
  1183.    PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,500)-273.15
  1184.    SLI=Temp500V-PclTempV
  1185.    rec=CAPE(TLcl,PLcl,100,sndtemp,snddewp)
  1186.    Pos=subwrd(rec,1)
  1187.    CIN=subwrd(rec,2)

  1189.    If (SfcPlev != PMaxThee)
  1190.       PclTemp=LiftWet(TLclMax,PLclMax,500,0)
  1191.       PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,500)-273.15
  1192.       LIMax=Temp500V-PclTempV
  1193.       rec=CAPE(TLclMax,PLclMax,100,sndtemp,snddewp)
  1194.       PosMax=subwrd(rec,1)
  1195.       CINMax=subwrd(rec,2)
  1196.    Else
  1197.       LIMax=SLI
  1198.       PosMax=Pos
  1199.       CINMax=CIN
  1200.       MaxThee=SfcThee
  1201.    Endif

  1203.    If (PageType = "Portrait")
  1204.       If (Text1XC = -1)
  1205.          Text1XC=rxloc-0.75
  1206.       Endif
  1207.       If (Text1YC = -1)
  1208.          Text1YC=tyloc-2.25
  1209.       Endif
  1210.       If (Text2XC = -1)
  1211.          Text2XC=rxloc-0.75
  1212.       Endif
  1213.       If (Text2YC = -1)
  1214.          Text2YC=tyloc-3.25
  1215.       Endif
  1216.       If (Text3XC = -1)
  1217.           Text3XC=rxloc-0.75
  1218.       Endif
  1219.       If (Text3YC = -1)
  1220.          Text3YC=tyloc-4.40
  1221.       Endif
  1222.    Else
  1223.       If (Text1XC = -1)
  1224.          Text1XC=rxloc+2.50
  1225.       Endif
  1226.       If (Text1YC = -1)
  1227.          Text1YC=tyloc-3.00
  1228.       Endif
  1229.       If (Text2XC = -1)
  1230.          Text2XC=rxloc+2.50
  1231.       Endif
  1232.       If (Text2YC = -1)
  1233.          Text2YC=tyloc-4.00
  1234.       Endif
  1235.       If (Text3XC = -1)
  1236.          Text3XC=rxloc+2.50
  1237.       Endif
  1238.       If (Text3YC = -1)
  1239.          Text3YC=tyloc-5.10
  1240.       Endif
  1241.    Endif
  1242.    "set string 1 l 3"
  1243.    "set line 0 1 3"
  1244.    "draw recf  "Text1XC-0.75 " " Text1YC-0.40 " " Text1XC+0.75 " " Text1YC+0.25
  1245.    "set line 1 1 3"
  1246.    "draw rec  "Text1XC-0.75 " " Text1YC-0.40 " " Text1XC+0.75 " " Text1YC+0.25
  1247.    "draw string "Text1XC-0.70 " " Text1YC+0.10"  K"
  1248.    "draw string "Text1XC+0.25 " " Text1YC+0.10" " int(K)      
  1249.    "draw string "Text1XC-0.70 " " Text1YC-0.10 "  TT"
  1250.    "draw string "Text1XC+0.25 " " Text1YC-0.10 " " int(tt)
  1251.    "draw string "Text1XC-0.70 " " Text1YC-0.25 "  PW(cm)"
  1252.    "draw string "Text1XC+0.25 " " Text1YC-0.25 " " int(pw*100)/100
  1253.    "set line 0 1 3"
  1254.    "draw recf  "Text2XC-0.75 " " Text2YC-0.60 " " Text2XC+0.75 " " Text2YC+0.60
  1255.    "set line 1 1 3"
  1256.    "draw rec  "Text2XC-0.75 " " Text2YC-0.60 " " Text2XC+0.75 " " Text2YC+0.60
  1257.    "draw string "Text2XC-0.35 " " Text2YC+0.50 " Surface"
  1258.    "draw string "Text2XC-0.70 " " Text2YC+0.30 "  Temp(`3.`0C)"
  1259.    "draw string "Text2XC+0.25 " " Text2YC+0.30 " " int(Sfctemp*10)/10
  1260.    "draw string "Text2XC-0.70 " " Text2YC+0.15 "  Dewp(`3.`0C)"
  1261.    "draw string "Text2XC+0.25 " " Text2YC+0.15 " " int(Sfcdewp*10)/10
  1262.    "draw string "Text2XC-0.70 " " Text2YC "   `3z`0`bE`n(K)"
  1263.    "draw string "Text2XC+0.25 " " Text2YC " " int(SfcThee)
  1264.    "draw string "Text2XC-0.70 " " Text2YC-0.15 "  LI"
  1265.    "draw string "Text2XC+0.25 " " Text2YC-0.15 " " round(SLI)
  1266.    "draw string "Text2XC-0.70 " " Text2YC-0.30 "  CAPE(J)"
  1267.    "draw string "Text2XC+0.25 " " Text2YC-0.30 " " int(Pos)   
  1268.    "draw string "Text2XC-0.70 " " Text2YC-0.45 "  CIN(J)"
  1269.    "draw string "Text2XC+0.25 " " Text2YC-0.45 " " int(CIN)      
  1270.    "set line 0 1 3"
  1271.    "draw recf  "Text3XC-0.75 " " Text3YC-0.55 " "  Text3XC+0.75 " " Text3YC+0.55
  1272.    "set line 1 1 3"
  1273.    "draw rec  "Text3XC-0.75 " " Text3YC-0.55 " "  Text3XC+0.75 " " Text3YC+0.55
  1274.    "draw string "Text3XC-0.60 " " Text3YC+0.45 " Most Unstable"
  1275.    "draw string "Text3XC-0.70 " " Text3YC+0.20 "  Press(hPa)"
  1276.    "draw string "Text3XC+0.25 " " Text3YC+0.20 " " int(PMaxThee)
  1277.    "draw string "Text3XC-0.70 " " Text3YC+0.05 " `3z`0`bE`n(K)"
  1278.    "draw string "Text3XC+0.25 " " Text3YC+0.05 " " int(MaxThee)
  1279.    "draw string "Text3XC-0.70 " " Text3YC-0.10 " LI"
  1280.    "draw string "Text3XC+0.25 " " Text3YC-0.10 " "round(LIMax)
  1281.    "draw string "Text3XC-0.70 " " Text3YC-0.25 " CAPE(J)"
  1282.    "draw string "Text3XC+0.25 " " Text3YC-0.25 " "int(PosMax)
  1283.    "draw string "Text3XC-0.70 " " Text3YC-0.40 " CIN(J)"
  1284.    "draw string "Text3XC+0.25 " " Text3YC-0.40 " " int(CINMax)
  1285. Endif

  1287. *-----------------------------
  1288. * Draw wind profile along side
  1289. *-----------------------------

  1291. If (DrawBarb = 1)
  1292.    say "Drawing Wind Profile."
  1293.    If (poleloc = -1)
  1294.       If (PageType = "Portrait")
  1295.          poleloc = 8.0
  1296.       Else
  1297.          poleloc = 7.5
  1298.       Endif
  1299.    Endif
  1300.    If (barbline = 1)
  1301.       "set line 1 1 3"
  1302.       "draw line "poleloc " " byloc " " poleloc " " tyloc
  1303.    Endif
  1304.    If (BarbCol = -1)
  1305.       'set rgb 41 255 0 132'
  1306.       'set rgb 42 255 0 168'
  1307.       'set rgb 43 255 0 204'
  1308.       'set rgb 44 255 0 240'
  1309.       'set rgb 45 255 0 255'
  1310.       'set rgb 46 204 0 255'
  1311.       'set rgb 47 174 0 255'
  1312.       'set rgb 48 138 0 255'
  1313.       'set rgb 49 108 0 255'
  1314.       'set rgb 50 84 0 255'
  1315.       'set rgb 51 40 0 255'
  1316.       'set rgb 52 0 0  255'
  1317.       'set rgb 53 0 42 255'
  1318.       'set rgb 54 0 84 255'
  1319.       'set rgb 55 0 120 255'
  1320.       'set rgb 56 0 150 255'
  1321.       'set rgb 57 0 192 255'
  1322.       'set rgb 58 0 240 255'
  1323.       'set rgb 59 0 255 210'
  1324.       'set rgb 60 0 255 160'
  1325.       'set rgb 61 0 255 126'
  1326.       'set rgb 62 0 255 78'
  1327.       'set rgb 63 84 255 0'
  1328.       'set rgb 64 138 255 0'
  1329.       'set rgb 65 188 255 0'
  1330.       'set rgb 66 236 255 0'
  1331.       'set rgb 67 255 255 0'
  1332.       'set rgb 68 255 222 0'
  1333.       'set rgb 69 255 192 0'
  1334.       'set rgb 70 255 162 0'
  1335.       'set rgb 71 255 138 0'
  1336.       'set rgb 72 255 108 0'
  1337.       'set rgb 73 255 84 0'
  1338.       'set rgb 74 255 54 0'
  1339.       'set rgb 75 255 12 0'
  1340.       'set rgb 76 255 0 34'
  1341.       'set rgb 77 255 0 70'
  1342.       'set rgb 78 255 0 105'
  1343.       'set rgb 79 255 0 140'
  1344.       'set rgb 80 255 0 175'
  1345.       'set rgb 81 255 0 215'
  1346.       'set rgb 82 255 0 255'
  1347.       'set rgb 83 255 255 255'

  1349.       col1='83 83 83 83 83 83 83 83 83 83 82 81 80 79 78'
  1350.       col2='77 76 75 74 73 72 71 70 69 68 67 66 65 64 63'
  1351.       col3='62 61 60 59 58 57 56 55 54 53 52 51 50 49 48'
  1352.       'set rbcols 'col1' 'col2' 'col3
  1353.    Endif
  1354.    "set z "_zmin" "_zmax
  1355.    "set gxout stat"
  1356.    zz=1
  1357.    wspd=-999
  1358.    cont=1
  1359.    While (cont = 1 & zz < _zmax)
  1360.       "set z "zz
  1361.       pres=subwrd(result,4)
  1362.       "d "sndspd
  1363.       rec=sublin(result,8)
  1364.       wspd=subwrd(rec,4)
  1365.       if (wspd < 0 | pres > _pmax)
  1366.           zz=zz+1
  1367.       else
  1368.           cont=0
  1369.       Endif
  1370.    Endwhile
  1371.    While (zz <= _zmax)
  1372.       "d "sndspd"(z="zz")"
  1373.       rec=sublin(result,8)
  1374.       wspd=subwrd(rec,4)
  1375.       If (BarbCol >= 0)
  1376.          "set line "BarbCol " 1 "BarbThk
  1377.       Else
  1378.          tempbcol=55+wspd/5     
  1379.          If (tempbcol > 83)
  1380.             tempbcol = 83
  1381.          Endif
  1382.          "set line "tempbcol " 1 "BarbThk
  1383.       Endif
  1384.       "d "snddir"(z="zz")"
  1385.       rec=sublin(result,8)
  1386.       wdir=subwrd(rec,4)
  1387.       xwind=GetUWnd(wspd,wdir)
  1388.       ywind=GetVWnd(wspd,wdir)
  1389.       "query gr2xy 5 "zz
  1390.       y1=subwrd(result,6)
  1391.       if (wspd > 0)
  1392.          cc=polelen/wspd
  1393.          xendpole=poleloc-xwind*cc
  1394.          yendpole=y1-ywind*cc
  1395.       endif
  1396.       if (xendpole>0 & wspd >= 0.5)
  1397.         if (flagbase = 1)
  1398.            "draw mark 3 "poleloc " " y1 " 0.05"
  1399.         endif
  1400.         "draw line " poleloc " " y1 "  " xendpole " " yendpole
  1401.         flagloop=wspd/10
  1402.         windcount=wspd
  1403.         flagcount=0
  1404.         xflagstart=xendpole
  1405.         yflagstart=yendpole
  1406.         dx=cos((180-wdir)*_dtr)
  1407.         dy=sin((180-wdir)*_dtr)
  1408.         while (windcount > 47.5)
  1409.            flagcount=flagcount+1
  1410.            dxflag=-len50*dx
  1411.            dyflag=-len50*dy
  1412.            xflagend=xflagstart+dxflag
  1413.            yflagend=yflagstart+dyflag
  1414.            windcount=windcount-50
  1415.            x1=xflagstart+0.5*wid50*xwind/wspd
  1416.            y1=yflagstart+0.5*wid50*ywind/wspd
  1417.            x2=xflagstart-0.5*wid50*xwind/wspd
  1418.            y2=yflagstart-0.5*wid50*ywind/wspd
  1419.            If (Fill50 = 1)
  1420.               "draw polyf "x1" "y1" "x2" "y2" "xflagend" "yflagend" "x1" "y1
  1421.            Else
  1422.               "draw line "x1 " "y1 " " xflagend " " yflagend " "  
  1423.               "draw line "x2 " "y2 " " xflagend " " yflagend
  1424.               "draw line "x1 " "y1 " " x2 " " y2
  1425.            Endif
  1426.            xflagstart=xflagstart+spac50*xwind/wspd
  1427.            yflagstart=yflagstart+spac50*ywind/wspd
  1428.         endwhile
  1429.         while (windcount > 7.5 )
  1430.            flagcount=flagcount+1
  1431.            dxflag=-len10*dx
  1432.            dyflag=-len10*dy
  1433.            xflagend=xflagstart+dxflag
  1434.            yflagend=yflagstart+dyflag
  1435.            windcount=windcount-10
  1436.            "draw line " xflagstart " " yflagstart " " xflagend " " yflagend
  1437.            xflagstart=xflagstart+spac10*xwind/wspd
  1438.            yflagstart=yflagstart+spac10*ywind/wspd
  1439.         endwhile
  1440.         if (windcount > 2.5)
  1441.            flagcount=flagcount+1
  1442.            if (flagcount = 1)
  1443.               xflagstart=xflagstart+spac05*xwind/wspd
  1444.               yflagstart=yflagstart+spac05*ywind/wspd
  1445.            endif
  1446.            dxflag=-len05*dx
  1447.            dyflag=-len05*dy
  1448.            xflagend=xflagstart+dxflag
  1449.            yflagend=yflagstart+dyflag
  1450.            windcount=windcount-5
  1451.            "draw line " xflagstart " " yflagstart " " xflagend " " yflagend
  1452.         endif
  1453.       else
  1454.         if (wspd < 0.5 & wspd >= 0)
  1455.            "draw mark 2 " poleloc " " y1 " 0.08"
  1456.         endif
  1457.       endif
  1458.       zz=zz+barbint
  1459.    endwhile
  1460. Endif

  1462. *----------------
  1463. * Draw Hodograph
  1464. *----------------

  1466. If (DrawHodo = 1)
  1467.    say "Drawing Hodograph."

  1469.    If (HodXcent = -1 | HodYcent = -1)
  1470.       If (PageType = "Portrait")
  1471.          HodXcent=6
  1472.          HodYcent=9.5
  1473.       Else
  1474.          HodXcent=9
  1475.          HodYcent=7.0
  1476.       Endif
  1477.    Endif
  1478.    HodL=HodXcent-HodSize/2.0
  1479.    HodR=HodXcent+HodSize/2.0
  1480.    HodT=HodYcent+HodSize/2.0
  1481.    HodB=HodYcent-HodSize/2.0
  1482.    RingSpac=HodSize/(NumRing*2)
  1483.    "set line 0"
  1484.    "draw recf "HodL" "HodB" "HodR" "HodT
  1485.    "set line "HodoCol" 1 6"
  1486.    "draw rec "HodL" "HodB" "HodR" "HodT
  1487.    "set line 1 1 3"
  1488.    "set string 1 c"
  1489.    "draw mark 1 "HodXcent " "HodYcent " " HodSize
  1490.    i=1
  1491.    While (i <= NumRing)
  1492.      "set strsiz 0.10"
  1493.      "set string 1 c 3 45"
  1494.      uwnd=-i*HodRing*cos(45*_dtr)
  1495.      xloc=HodXcent+uwnd*RingSpac/HodRing
  1496.      yloc=HodYcent+uwnd*RingSpac/HodRing
  1497.   
  1498.      "draw mark 2 "HodXcent " " HodYcent " " i*HodSize/NumRing
  1499.      "draw string "xloc " " yloc " " HodRing*i
  1500.      i=i+1
  1501.    Endwhile
  1502.    "set string 1 l 3 0"
  1503.    If (TickInt > 0)
  1504.       i=0
  1505.       while (i < HodRing*NumRing)
  1506.          dist=i*RingSpac/HodRing
  1507.          hrxloc=HodXcent+dist               
  1508.          hlxloc=HodXcent-dist                     
  1509.          htyloc=HodYcent+dist
  1510.          hbyloc=HodYcent-dist
  1511.          "set line 1 1 3"
  1512.          "draw line "hrxloc " " HodYcent-TickSize/2 " " hrxloc " " HodYcent+TickSize/2
  1513.          "draw line "hlxloc " " HodYcent-TickSize/2 " " hlxloc " " HodYcent+TickSize/2
  1514.          "draw line "HodXcent+TickSize/2 " " htyloc " " HodXcent-TickSize/2 " " htyloc
  1515.          "draw line "HodXcent+TickSize/2 " " hbyloc " " HodXcent-TickSize/2 " " hbyloc
  1516.          i=i+TickInt
  1517.       endwhile
  1518.    Endif
  1519.    "set line "HodoCol " " HodoLine " "HodoThk
  1520.    "draw string "HodL+0.05 " " HodT-0.1 " knots"
  1521.    zloop=_zmin
  1522.    xold=-999
  1523.    yold=-999
  1524.    count=0
  1525.    Depth=0
  1526.    While (zloop < _zmax & Depth < HodoDep)
  1527.       "set z "zloop
  1528.       pres=subwrd(result,4)
  1529.       "d "sndspd
  1530.       rec=sublin(result,8)
  1531.       wspd=subwrd(rec,4)
  1532.       "d "snddir
  1533.       rec=sublin(result,8)
  1534.       wdir=subwrd(rec,4)
  1535.       uwnd=GetUWnd(wspd,wdir)
  1536.       vwnd=GetVWnd(wspd,wdir)
  1537.       If (wspd >= 0)
  1538.          xloc=HodXcent+uwnd*RingSpac/HodRing
  1539.          yloc=HodYcent+vwnd*RingSpac/HodRing
  1540.          If (xloc > 0 & yloc > 0 & xold > 0 & yold > 0)
  1541.             Depth=Depth+pold-pres
  1542.             count=count+1
  1543.             If (count = 1)
  1544.                "draw mark 3 "xold " " yold " 0.05"
  1545.             Endif
  1546.             "draw line "xold" "yold" "xloc" "yloc
  1547.          Endif
  1548.          xold=xloc
  1549.          yold=yloc
  1550.       Endif
  1551.       zloop=zloop+1
  1552.       pold=pres
  1553.    EndWhile

  1555.    If (count > 0)
  1556.       "draw mark 3 "xold " " yold " 0.05"
  1557.    Endif
  1558. Endif

  1560. *----------------------------------------------
  1561. * Calculate and Display Absolute & S-R Helicity
  1562. *----------------------------------------------

  1564. If (DrawHeli = 1)
  1565.    say "Calculating Helicity & SR Helicity."
  1566.    delp=10
  1567.    UTotal=0
  1568.    VTotal=0

  1570. * First, calculate mass-weighted mean wind
  1571. * Since delp is a constant, and mass is proportional to
  1572. * delp, this is a simple sum.

  1574.    "set lev "_pmax " " _pmin
  1575.    "define uwndarr="sndspd"*cos((270-"snddir")*"_dtr")"
  1576.    "define vwndarr="sndspd"*sin((270-"snddir")*"_dtr")"

  1578.    pres=MeanVBot
  1579.    While (pres >= MeanVTop)
  1580.       uwnd=interp(uwndarr,pres)*_ktm
  1581.       vwnd=interp(vwndarr,pres)*_ktm
  1582.       If (uwnd > -900 & vwnd > -900)
  1583.          UTotal=UTotal+uwnd
  1584.          VTotal=VTotal+vwnd
  1585.       Endif
  1586.       pres=pres-delp
  1587.    EndWhile
  1588.    vcount=1+(MeanVBot-MeanVTop)/delp
  1589.    Umean=UTotal/vcount
  1590.    Vmean=VTotal/vcount
  1591.    Spdmean=GetWSpd(Umean,Vmean)
  1592.    MeanDir=GetWDir(Umean,Vmean)

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

  1596.    If (StormMot = 1)
  1597.       If (Spdmean < 15)
  1598.          Reduct=0.25
  1599.          Rotate=30
  1600.       Else
  1601.          Reduct=0.20
  1602.          Rotate=20
  1603.       Endif
  1604.    Else
  1605.       Reduct=0.0
  1606.       Rotate=0.0
  1607.    Endif

  1609.    UReduce=(1-Reduct)*Umean
  1610.    VReduce=(1-Reduct)*Vmean
  1611.    StormSpd=GetWSpd(UReduce,VReduce)

  1613.    StormDir=GetWDir(UReduce,VReduce)+Rotate
  1614.    If (StormDir >= 360)
  1615.       StormDir=StormDir-360
  1616.    Endif

  1618.    StormU=GetUWnd(StormSpd,StormDir)
  1619.    StormV=GetVWnd(StormSpd,StormDir)

  1621. * Draw Storm Motion Vector

  1623.    xloc=HodXcent+_mtk*StormU*RingSpac/HodRing
  1624.    yloc=HodYcent+_mtk*StormV*RingSpac/HodRing

  1626.    "set line 1 1 4"
  1627.    "draw line "HodXcent " " HodYcent " " xloc " " yloc
  1628.    Arr1U=GetUWnd(HodRing/10,StormDir+30)
  1629.    Arr1V=GetVWnd(HodRing/10,StormDir+30)
  1630.    Arr2U=GetUWnd(HodRing/10,StormDir-30)
  1631.    Arr2V=GetVWnd(HodRing/10,StormDir-30)

  1633.    xloc2=xloc-Arr1U/HodRing
  1634.    xloc3=xloc-Arr2U/HodRing
  1635.    yloc2=yloc-Arr1V/HodRing
  1636.    yloc3=yloc-Arr2V/HodRing

  1638.    "set line 1 1 3"

  1640.    If (FillArrw = 0)
  1641.       "draw line "xloc" "yloc" "xloc2" "yloc2
  1642.       "draw line "xloc" "yloc" "xloc3" "yloc3
  1643.    Else
  1644.       "draw polyf "xloc" "yloc" "xloc2" "yloc2" "xloc3" "yloc3" "xloc" "yloc
  1645.    Endif


  1648. * Now, calculate SR and Environmental Helicity

  1650.    helic=0
  1651.    SRhelic=0
  1652.    MinP=SfcPlev-HelicDep
  1653.    pres=SfcPlev
  1654.    uwndold=-999
  1655.    vwndold=-999
  1656.    While (pres >= MinP)
  1657.       uwnd=interp(uwndarr,pres)*_ktm
  1658.       vwnd=interp(vwndarr,pres)*_ktm
  1659.       If (uwnd > -900 & uwndold > -900)
  1660.           du=uwnd-uwndold
  1661.           dv=vwnd-vwndold
  1662.           ubar=0.5*(uwnd+uwndold)
  1663.           vbar=0.5*(vwnd+vwndold)
  1664.           uhelic=-dv*ubar                  
  1665.           vhelic=du*vbar                  
  1666.           SRuhelic=-dv*(ubar-StormU)
  1667.           SRvhelic=du*(vbar-StormV)
  1668.           SRhelic=SRhelic+SRuhelic+SRvhelic
  1669.           helic=helic+uhelic+vhelic
  1670.       Endif
  1671.       uwndold=uwnd
  1672.       vwndold=vwnd
  1673.       pres=pres-delp
  1674.    EndWhile

  1676.    "set strsiz 0.1"
  1677.    "set string 1 l 3"
  1678.    If (PageType = "Portrait")
  1679.       If (Text4XC = -1)
  1680.          Text4XC=rxloc-0.75
  1681.       Endif
  1682.       If (Text4YC = -1)
  1683.          Text4YC=tyloc-5.45
  1684.       Endif
  1685.    Else
  1686.       If (Text4XC = -1)
  1687.          Text4XC=rxloc+2.50
  1688.       Endif
  1689.       If (Text4YC = -1)
  1690.          Text4YC=tyloc-6.10
  1691.       Endif
  1692.    Endif
  1693.    "set line 0 1 3"
  1694.    "draw recf  "Text4XC-0.75 " "Text4YC-0.5 " " Text4XC+0.75 " " Text4YC+0.5
  1695.    "set line 1 1 3"
  1696.    "draw rec  "Text4XC-0.75 " "Text4YC-0.5 " " Text4XC+0.75 " " Text4YC+0.5
  1697.    "draw string "Text4XC-0.45 " " Text4YC+0.40 " Hodograph"
  1698.    "draw string "Text4XC-0.70 " " Text4YC+0.20 " EH"
  1699.    "draw string "Text4XC+0.25 " " Text4YC+0.20 " "int(helic)
  1700.    "draw string "Text4XC-0.70 " " Text4YC+0.05 " SREH"
  1701.    "draw string "Text4XC+0.25 " " Text4YC+0.05 " " int(SRhelic)
  1702.    "draw string "Text4XC-0.70 " " Text4YC-0.20 " StmDir"
  1703.    "draw string "Text4XC+0.25 " " Text4YC-0.20 " " int(StormDir)"`3.`0"
  1704.    "draw string "Text4XC-0.70 " " Text4YC-0.35 " StmSpd(kt)"
  1705.    "draw string "Text4XC+0.25 " " Text4YC-0.35 " " int(_mtk*StormSpd)
  1706. Endif

  1708. *---------------------------------------------------
  1709. * Plot RH profile.
  1710. *---------------------------------------------------

  1712. If (DrawRH = 1)
  1713.   "set z "_zmin" "_zmax
  1714.   "set x "_xval
  1715.   "set y "_yval
  1716.   "set t "_tval
  1717.   "set zlog on"
  1718.   "set gxout line"
  1719.   "set ccolor "RHCol
  1720.   "set cstyle "RHLine
  1721.   "set cmark "RHMrk
  1722.   "set digsiz "MrkSize
  1723.   "set missconn on"
  1724.   "set xlab on"
  1725.   "set frame off"
  1726.   "set vrange 0 350"
  1727.   "set xlpos 0 t"
  1728.   "set xlevs 25 50 75 100"
  1729.   "set grid vertical 5"
  1730.   "define rh=100*exp((17.2694*"snddewp")/("snddewp"+237.3)-(17.2694*"sndtemp")/("sndtemp"+237.3))"
  1731.   "d rh"
  1732.    If (LblAxes = 1)
  1733.      "set string 1 c 3 0"
  1734.      "set strsiz 0.125"
  1735.      If (PageType = "Portrait")
  1736.        "draw string 1.5 10.85 RH (%)"
  1737.      Else
  1738.        "draw string 1.75 8.35 RH (%)"
  1739.      Endif
  1740.    Endif
  1741. Endif

  1743. *------------------------------------------
  1744. * Reset environment to original dimensions
  1745. *------------------------------------------

  1747. "set t "_tval
  1748. "set x "_xval
  1749. "set y "_yval
  1750. "set z "_zmin " "_zmax

  1752. say "Done."

  1754. Return(0)

  1756. *************************************************************************
  1757. function Templcl(temp,dewp)

  1759. *------------------------------------------------------
  1760. * Calculate the temp at the LCL given temp & dewp in C
  1761. *------------------------------------------------------

  1763. tempk=temp+273.15
  1764. dewpk=dewp+273.15
  1765. Parta=1/(dewpk-56)
  1766. Partb=log(tempk/dewpk)/800
  1767. Tlcl=1/(Parta+Partb)+56
  1768. return(Tlcl-273.15)

  1770. **************************************************************************

  1772. function Preslcl(temp,dewp,pres)

  1774. *-------------------------------------------------------
  1775. * Calculate press of LCL given temp & dewp in C and pressure
  1776. *-------------------------------------------------------

  1778. Tlclk=Templcl(temp,dewp)+273.15
  1779. tempk=temp+273.15
  1780. theta=tempk*pow(1000/pres,0.286)
  1781. plcl=1000*pow(Tlclk/theta,3.48)
  1782. return(plcl)

  1784. **************************************************************************
  1785. function LiftWet(startt,startp,endp,display,Pmin,Pmax)

  1787. *--------------------------------------------------------------------
  1788. * Lift a parcel moist adiabatically from startp to endp.
  1789. * Init temp is startt in C.  If you wish to see the parcel's
  1790. * path plotted, display should be 1.  Returns temp of parcel at endp.
  1791. *--------------------------------------------------------------------

  1793. temp=startt
  1794. pres=startp
  1795. cont = 1
  1796. delp=10
  1797. While (pres >= endp & cont = 1)
  1798.     If (display = 1)
  1799.        xtemp=GetXLoc(temp,pres)
  1800.        "q w2xy "xtemp" "pres
  1801.        xloc=subwrd(result,3)
  1802.        yloc=subwrd(result,6)
  1803.        If (xtemp < 0 | xtemp > 100)
  1804.           cont=0
  1805.        Else
  1806.           If (pres >= Pmin & pres < Pmax & pres < startp)  
  1807.              "draw line "xold" "yold" "xloc" "yloc
  1808.           Endif
  1809.        Endif
  1810.     Endif
  1811.     xold=xloc
  1812.     yold=yloc
  1813.     temp=temp-100*delp*gammaw(temp,pres-delp/2,100)
  1814.     pres=pres-delp
  1815. EndWhile
  1816. return(temp)


  1819. **************************************************************************
  1820. function LiftDry(startt,startp,endp,display,Pmin,Pmax)

  1822. *--------------------------------------------------------------------
  1823. * Lift a parcel dry adiabatically from startp to endp.
  1824. * Init temp is startt in C.  If you wish to see the parcel's
  1825. * path plotted, display should be 1.  Returns temp of parcel at endp.
  1826. *--------------------------------------------------------------------

  1828. starttk=startt+273.15
  1829. cont = 1
  1830. delp=10
  1831. round=int(startp/10)*10
  1832. subscr=0.1*round
  1833. powstart=pow(startp,-0.286)
  1834. temp=starttk*_powpres.subscr*powstart-273.15
  1835. pres=round-10
  1836. While (pres >= endp & cont = 1)
  1837.     subscr=0.1*pres
  1838.     temp=starttk*_powpres.subscr*powstart-273.15
  1839.     If (display = 1)
  1840.        xtemp=GetXLoc(temp,pres)
  1841.        "q w2xy "xtemp" "pres
  1842.        xloc=subwrd(result,3)
  1843.        yloc=subwrd(result,6)
  1844.        If (xtempold > 0 & xtempold < 100 & xtemp > 0 & xtemp < 100)
  1845.           If (pres >= Pmin & pres < Pmax & pres < startp)  
  1846.              "draw line "xold" "yold" "xloc" "yloc
  1847.           Endif
  1848.        Endif
  1849.     Endif
  1850.     xold=xloc
  1851.     xtempold=xtemp
  1852.     yold=yloc
  1853.     pres=pres-delp
  1854. EndWhile
  1855. return(temp)

  1857. **************************************************************************
  1858. function CAPE(startt,startp,endp,sndtemp,snddewp)

  1860. *---------------------------------------------------------------------
  1861. * Returns all postive area and convective inhibition above LCL.
  1862. * Parcel is lifted from LCL at startt,startp and is halted
  1863. * at endp.   Integration method used is trapezoid method.
  1864. *---------------------------------------------------------------------

  1866. pres=startp
  1867. PclTemp=startt
  1868. PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15
  1869. delp=10
  1870. Pos=0
  1871. Neg=0
  1872. Neg2=0

  1874. count=0
  1875. While (pres >= endp)
  1876.    EnvTemp=interp(sndtemp,pres)
  1877.    EnvDewp=interp(snddewp,pres)
  1878.    EnvTempV=virtual2(EnvTemp+273.15,EnvDewp+273.15,pres)-273.15
  1879.    DelT=PclTempV-EnvTempV
  1880.    If (abs(EnvTempV) < 130 & abs(PclTempV) < 130)
  1881.      count=count+1
  1882.      If (count > 1)
  1883.        Val=DelT/pres+Prev
  1884.        If (Val > 0)
  1885.           Pos=Pos+Val
  1886.           Neg2=0
  1887.        Else
  1888.           Neg=Neg+abs(Val)
  1889.           Neg2=Neg2+abs(Val)
  1890.        Endif
  1891.      Endif
  1892.      Prev=DelT/pres
  1893.    Endif
  1894.    pres=pres-delp
  1895.    PclTemp=PclTemp-100*delp*gammaw(PclTemp,pres,100)
  1896.    PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15
  1897. Endwhile

  1899. Pos=0.5*Pos*287*delp
  1900. CIN=0.5*(Neg-Neg2)*287*delp

  1902. return(Pos" "CIN)

  1904. ***************************************************************************
  1905. function gammaw(tempc,pres,rh)

  1907. *-----------------------------------------------------------------------
  1908. * Function to calculate the moist adiabatic lapse rate (deg C/Pa) based
  1909. * on the temperature, pressure, and rh of the environment.
  1910. *----------------------------------------------------------------------

  1912. tempk=tempc+273.15
  1913. es=satvap2(tempc)
  1914. ws=mixratio(es,pres)
  1915. w=rh*ws/100
  1916. tempv=virtual(tempk,w)
  1917. latent=latentc(tempc)

  1919. A=1.0+latent*ws/(287*tempk)
  1920. B=1.0+0.622*latent*latent*ws/(1005*287*tempk*tempk)
  1921. Density=100*pres/(287*tempv)
  1922. lapse=(A/B)/(1005*Density)
  1923. return(lapse)

  1925. *************************************************************************
  1926. function latentc(tempc)

  1928. *-----------------------------------------------------------------------
  1929. * Function to return the latent heat of condensation in J/kg given
  1930. * temperature in degrees Celsius.
  1931. *-----------------------------------------------------------------------

  1933. val=2502.2-2.43089*tempc

  1935. return(val*1000)

  1937. *************************************************************************
  1938. function precipw(sndtemp,snddewp,startp,endp)

  1940. *-----------------------------------------------------------------------
  1941. * Function to calculate the precipitable water (cm) in a sounding
  1942. * starting at pressure level startp and ending at pressure level endp.
  1943. *-----------------------------------------------------------------------

  1945. ppold=-999
  1946. ttold=-999
  1947. ddold=-999
  1948. delp=10
  1949. Int=0
  1950. mix=0
  1951. pres=startp
  1952. logpp=log(pres)
  1953. logppm=log(pres-delp)
  1954. while (pres >= endp)
  1955.    tt=interp(sndtemp,pres)
  1956.    dd=interp(snddewp,pres)
  1957.    if (tt>-900 & ttold>-900 & dd>-900 & ddold>-900)
  1958.       e=satvap2(dd)
  1959.       mix=mixratio(e,pres)
  1960.       mixavg=(logpp*mix+logppm*mixold)/(logpp+logppm)
  1961.       Int=Int+1.020408*mixavg*delp
  1962.    endif
  1963.    ttold=tt
  1964.    ddold=dd
  1965.    ppold=pp
  1966.    mixold=mix
  1967.    pres=pres-delp
  1968.    logpp=logppm
  1969.    logppm=log(pres-delp)
  1970. endwhile

  1972. return(Int)

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

  1976. function virtual(temp,mix)

  1978. *------------------------------------------------------------
  1979. * Function to return virtual temperature given temperature in
  1980. * kelvin and mixing ratio in g/g.
  1981. *-------------------------------------------------------------

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

  1985. return (tempv)

  1987. ************************************************************************

  1989. function virtual2(temp,dewp,pres)
  1990.   
  1991. *------------------------------------------------------------
  1992. * Function to return virtual temperature in kelvin given temperature in
  1993. * kelvin and dewpoint in kelvin and pressure in hPa
  1994. *-------------------------------------------------------------

  1996. if (temp > 0 & dewp > 0)
  1997.   vap=satvap2(dewp-273.15)
  1998.   mix=mixratio(vap,pres)
  1999.   tempv=virtual(temp,mix)
  2000. else
  2001.   tempv=-9999
  2002. endif

  2004. return (tempv)
  2005.   
  2006. ************************************************************************

  2008. function satvapor(temp)

  2010. *---------------------------------------------------------------
  2011. * Given temp in Celsius, returns saturation vapor pressure in hPa
  2012. *---------------------------------------------------------------

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

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

  2018. ************************************************************************

  2020. function satvap2(temp)

  2022. *---------------------------------------------------------------
  2023. * Given temp in Celsius, returns saturation vapor pressure in hPa*---------------------------------------------------------------

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

  2027. return(es)

  2029. *************************************************************************

  2031. function mixratio(e,p)

  2033. *------------------------------------------------------
  2034. * Given vapor pressure and pressure, return mixing ratio
  2035. *-------------------------------------------------------

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

  2039. return(mix)

  2041. ************************************************************************

  2043. function getrh(temp,dewp,pres)

  2045. tempk=temp+273.15
  2046. dewpk=dewp+273.15

  2048. es=satvap2(temp)

  2050. If (temp > 0)
  2051.    A=2.53*pow(10,9)
  2052.    B=5420
  2053. Else
  2054.    A=3.41*pow(10,10)
  2055.    B=6130
  2056. Endif

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

  2061. return(100*w/ws)

  2063. ************************************************************************
  2064. function interp(array,pres)

  2066. *------------------------------------------------------------------------
  2067. * Interpolate inside array for pressure level pres.
  2068. * Returns estimated value of array at pressure pres.
  2069. *------------------------------------------------------------------------

  2071. "set gxout stat"
  2072. "set lev "pres
  2073. altpres=subwrd(result,4)
  2074. "q dims"
  2075. rec=sublin(result,4)
  2076. zlev=subwrd(rec,9)

  2078. If (zlev < 2 | zlev > _zmaxfile)
  2079.   Vest = -9999.0
  2080. Else
  2081.   If (altpres > pres)
  2082.     zlev=zlev+1
  2083.   Endif
  2084.   "set z "zlev
  2085.   PAbove=subwrd(result,4)
  2086.   "d "array"(lev="PAbove")"
  2087.   rec=sublin(result,8)
  2088.   VAbove=subwrd(rec,4)
  2089.   "set z "zlev-1
  2090.   PBelow=subwrd(result,4)
  2091.   "d "array"(lev="PBelow")"
  2092.   rec=sublin(result,8)
  2093.   VBelow=subwrd(rec,4)

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

  2097.   If (abs(VAbove) > 130 & zlev > 1 & zlev < _zmaxfile)
  2098.      zz=zlev
  2099.      While (abs(VAbove) > 130 & zz < _zmaxfile)
  2100.        zz=zz+1
  2101.        "set z "zz
  2102.        PAbove=subwrd(result,4)
  2103.        "d "array"(lev="PAbove")"
  2104.        rec=sublin(result,8)
  2105.        VAbove=subwrd(rec,4)
  2106.      EndWhile
  2107.   Endif

  2109.   If (abs(VBelow) > 130 & zlev > 1 & zlev < _zmaxfile)
  2110.       zz=zlev-1
  2111.       While (abs(VBelow) > 130 & zz > 1)
  2112.         zz=zz-1
  2113.         "set z "zz
  2114.         PBelow=subwrd(result,4)
  2115.         "d "array"(lev="PBelow")"
  2116.         rec=sublin(result,8)
  2117.         VBelow=subwrd(rec,4)
  2118.       EndWhile
  2119.   Endif

  2121.   If (abs(VAbove) < 130 & abs(VBelow) < 130)
  2122.      Vest=VBelow+log(PBelow/pres)*(VAbove-VBelow)/log(PBelow/PAbove)
  2123.   Else
  2124.      Vest=-9999.0
  2125.   Endif

  2127. Endif

  2129. Return(Vest)

  2131. ****************************************************************************

  2133. function GetUWnd(wspd,wdir)

  2135. *------------------------
  2136. * Get x-component of wind.
  2137. *------------------------


  2140. If (wspd >= 0)
  2141.    xwind=wspd*cos((270-wdir)*_dtr)
  2142. Else
  2143.    xwind = -9999.0
  2144. Endif
  2145. return(xwind)

  2147. **************************************************************************

  2149. function GetVWnd(wspd,wdir)

  2151. *-----------------------
  2152. * Get y-component of wind
  2153. *------------------------

  2155. If (wspd >= 0)
  2156.    ywind=wspd*sin((270-wdir)*_dtr)
  2157. Else
  2158.    ywind = -9999.0
  2159. Endif
  2160. return(ywind)


  2163. *************************************************************************

  2165. function GetWSpd(xwind,ywind)


  2168. "set gxout stat"
  2169. "d mag("xwind","ywind")"
  2170. rec=sublin(result,8)
  2171. val=subwrd(rec,4)

  2173. return (val)

  2175. *************************************************************************

  2177. function GetWDir(xwind,ywind)

  2179. * Return wind direction given x and y components

  2181. "set gxout stat"
  2182. "define theta=270-"_rtd"*atan2("ywind","xwind")"
  2183. "d theta"
  2184. rec=sublin(result,8)
  2185. Dir=subwrd(rec,4)

  2187. If (Dir > 360)
  2188.    Dir=Dir-360
  2189. Endif

  2191. If (Dir < 0)
  2192.    Dir=360+Dir
  2193. Endif

  2195. return(Dir)

  2197. *************************************************************************

  2199. function GetXLoc(temp,pres)

  2201. *-------------------------------------------------
  2202. * Get x-location on skew-t based on temp, pressure
  2203. *-------------------------------------------------

  2205. xloc=(temp-_m1*log10(pres)-_m3)/_m2
  2206. return(xloc)

  2208. *************************************************************************

  2210. function GetTemp(xloc,pres)

  2212. *-------------------------------------------------
  2213. * Return temperature at location given by xloc,pres
  2214. *-------------------------------------------------

  2216. tempval=_m1*log10(pres)+_m2*xloc+_m3
  2217. return(tempval)

  2219. **************************************************************************

  2221. function GetTheta(temp,pres)         

  2223. *---------------------------------------------------
  2224. * Calculate theta for a given temperature and pressure
  2225. *---------------------------------------------------

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


  2231. *************************************************************************

  2233. function GetThet2(temp,dewp,pres)         

  2235. *---------------------------------------------------
  2236. * Calculate theta for a given temperature,dewp, and pressure
  2237. *---------------------------------------------------

  2239. tempk=273.15+temp
  2240. dewpk=273.15+dewp

  2242. es=satvap2(temp)
  2243. ws=mixratio(es,pres)

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

  2247. exponent=0.2854*(1.0-0.00028*mix)
  2248. theta=(temp+273.15)*pow(1000/pres,exponent)-273.15
  2249. return(theta)

  2251. *************************************************************************

  2253. function Thetae(temp,dewp,pres)

  2255. *--------------------------------------------------------------
  2256. * Return equiv. pot. temp in Kelvin given temp, dewp in celsius
  2257. * From Bolton (1980) Mon Wea Rev
  2258. *--------------------------------------------------------------

  2260. es=satvap2(temp)
  2261. ws=mixratio(es,pres)
  2262. mix=10*getrh(temp,dewp,pres)*ws
  2263. theta=GetThet2(temp,dewp,pres)+273.15
  2264. TLcl=Templcl(temp,dewp)+273.15
  2265. thetae=theta*exp((3.376/TLcl-0.00254)*mix*(1.0+0.00081*mix))

  2267. return(thetae)

  2269. **************************************************************************


  2272. function int(i0)

  2274. *--------------------------
  2275. * Return integer of i0
  2276. *--------------------------
  2277.   i=0
  2278.   while(i<12)
  2279.     i=i+1
  2280.     if(substr(i0,i,1)='.')
  2281.       i0=substr(i0,1,i-1)
  2282.       break
  2283.     endif
  2284.   endwhile
  2285. return(i0)

  2287. *************************************************************************

  2289. function abs(i)

  2291. *----------------------------
  2292. * return absolute value of i
  2293. *----------------------------

  2295.   if (i < 0)
  2296.      absval=-i
  2297.   else
  2298.      absval=i
  2299.   endif

  2301. return(absval)

  2303. *************************************************************************

  2305. function log(i)

  2307. *---------------------------
  2308. * return natural log of i
  2309. *---------------------------

  2311. "set gxout stat"
  2312. "d log("i")"
  2313. rec=sublin(result,8)
  2314. val=subwrd(rec,4)
  2315. return(val)

  2317. *************************************************************************

  2319. function log10(i)

  2321. *--------------------------
  2322. * return log base 10 of i
  2323. *--------------------------

  2325. "set gxout stat"
  2326. "d log10("i")"
  2327. rec=sublin(result,8)
  2328. val=subwrd(rec,4)
  2329. return(val)

  2331. *************************************************************************

  2333. function pow(i,j)

  2335. *-------------------------------
  2336. * return power of i raised to j
  2337. *-------------------------------

  2339. "set gxout stat"
  2340. "d pow("i","j")"
  2341. rec=sublin(result,8)
  2342. val=subwrd(rec,4)
  2343. return(val)

  2345. ************************************************************************

  2347. function cos(i)

  2349. *-----------------------------------------
  2350. * return cosine of i, where i is in radians
  2351. *------------------------------------------

  2353. "set gxout stat"
  2354. "d cos("i")"
  2355. rec=sublin(result,8)
  2356. val=subwrd(rec,4)
  2357. return(val)

  2359. ************************************************************************

  2361. function sin(i)

  2363. *------------------------------------------
  2364. * return sine of i, where i is in radians
  2365. *------------------------------------------

  2367. "set gxout stat"
  2368. "d sin("i")"
  2369. rec=sublin(result,8)
  2370. val=subwrd(rec,4)
  2371. return(val)

  2373. ************************************************************************

  2375. function exp(i)

  2377. *------------------------------------------
  2378. * return exponential of i
  2379. *------------------------------------------

  2381. "set gxout stat"
  2382. "d exp("i")"
  2383. rec=sublin(result,8)
  2384. val=subwrd(rec,4)
  2385. return(val)

  2387. ***********************************************************************
  2388. function round(i)

  2390. rr=abs(1.0*i)
  2391. rr=int(rr+0.5)
  2392. if (i < 0)
  2393.    rr=-1*rr      
  2394. endif
  2395. return(rr)



复制代码


密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

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

能把您计算cape的脚本发给我吗 我也是想利于基础数据计算cape 565653979@qq 万分感谢
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

yn_newyork
发表于 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
回复 支持 反对

使用道具 举报

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

使用道具 举报

123
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|帮助中心|清除痕迹

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

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

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