| 1 | 'reinit' |
|---|
| 2 | 'open foo.ctl' |
|---|
| 3 | |
|---|
| 4 | 'set t 1' |
|---|
| 5 | 'q dims' |
|---|
| 6 | rec=sublin(result,5) |
|---|
| 7 | _analysis=subwrd(rec,6) |
|---|
| 8 | say 'Initial Time is ' _analysis |
|---|
| 9 | |
|---|
| 10 | say 'Create gif images as well (1=yes ; 0=no)' |
|---|
| 11 | pull ans |
|---|
| 12 | frame = 1 |
|---|
| 13 | |
|---|
| 14 | 'q file' |
|---|
| 15 | rec=sublin(result,5) |
|---|
| 16 | _endtime=subwrd(rec,12) |
|---|
| 17 | 'q dims' |
|---|
| 18 | rec=sublin(result,2) |
|---|
| 19 | _alon1=subwrd(rec,6) |
|---|
| 20 | _alon2=subwrd(rec,8) |
|---|
| 21 | rec=sublin(result,3) |
|---|
| 22 | _alat1=subwrd(rec,6) |
|---|
| 23 | _alat2=subwrd(rec,8) |
|---|
| 24 | |
|---|
| 25 | BigRun = 1 |
|---|
| 26 | while (BigRun) |
|---|
| 27 | |
|---|
| 28 | irun = 1 |
|---|
| 29 | |
|---|
| 30 | * loc ( 38.58, -77.0 ) |
|---|
| 31 | say 'Please enter coordinated for sounding' |
|---|
| 32 | say ' ' |
|---|
| 33 | |
|---|
| 34 | say 'Latitude values between: ' _alat1 ' and ' _alat2 |
|---|
| 35 | say 'Enter Latitude' |
|---|
| 36 | pull lat1 |
|---|
| 37 | say ' ' |
|---|
| 38 | say 'Longitude values between: ' _alon1 ' and ' _alon2 |
|---|
| 39 | say 'Enter Longitude' |
|---|
| 40 | pull lon1 |
|---|
| 41 | |
|---|
| 42 | say ' ' |
|---|
| 43 | say 'Overlay two time periods 1=yes ; 0=no' |
|---|
| 44 | pull olans |
|---|
| 45 | overlay = 1 |
|---|
| 46 | |
|---|
| 47 | while(irun) |
|---|
| 48 | say ' ' |
|---|
| 49 | say 'Please enter time to plot' |
|---|
| 50 | say 'Available times: 1 to ' _endtime |
|---|
| 51 | pull ans |
|---|
| 52 | 'set t 'ans |
|---|
| 53 | 'set lon 'lon1 |
|---|
| 54 | 'set lat 'lat1 |
|---|
| 55 | 'set lev 1000 100' |
|---|
| 56 | 'define wspd=(sqrt(u*u + v*v))*1.89' |
|---|
| 57 | 'define dir=(atan2(u,v))*57.2957795+180.0' |
|---|
| 58 | if (overlay = 1) |
|---|
| 59 | SCol = 1 |
|---|
| 60 | HCol = 5 |
|---|
| 61 | BCol = -1 |
|---|
| 62 | endif |
|---|
| 63 | if (overlay = 2) |
|---|
| 64 | SCol = 3 |
|---|
| 65 | HCol = 11 |
|---|
| 66 | BCol = 1 |
|---|
| 67 | endif |
|---|
| 68 | |
|---|
| 69 | rc=plotskew(tc,td,wspd,dir,0,SCol,SCol,SCol,HCol,overlay,BCol) |
|---|
| 70 | |
|---|
| 71 | if (olans) |
|---|
| 72 | overlay = 2 |
|---|
| 73 | olans = 0 |
|---|
| 74 | else |
|---|
| 75 | irun = 0 |
|---|
| 76 | endif |
|---|
| 77 | |
|---|
| 78 | endwhile |
|---|
| 79 | |
|---|
| 80 | if(ans) |
|---|
| 81 | 'printim skew'frame'.gif gif' |
|---|
| 82 | frame=frame+1 |
|---|
| 83 | endif |
|---|
| 84 | |
|---|
| 85 | say ' ' |
|---|
| 86 | say 'Continue with another sounding 1=yes ; 0=no' |
|---|
| 87 | pull BigRun |
|---|
| 88 | |
|---|
| 89 | 'reset' |
|---|
| 90 | |
|---|
| 91 | endwhile |
|---|
| 92 | |
|---|
| 93 | |
|---|
| 94 | ************************************************************************* |
|---|
| 95 | function plotskew(sndtemp,snddewp,sndspd,snddir,ClrScrn,TSndCol,DSndCol,PrclCol,RHCol,OL,BarbCol) |
|---|
| 96 | ************************************************************************* |
|---|
| 97 | * |
|---|
| 98 | * GrADS Script to Plot a SkewT/LogP Diagram |
|---|
| 99 | * |
|---|
| 100 | * Bob Hart |
|---|
| 101 | * Penn State University / Dept of Meteorology |
|---|
| 102 | * Last Update: November 10, 1999 |
|---|
| 103 | * |
|---|
| 104 | * Recent Changes: |
|---|
| 105 | * |
|---|
| 106 | * 11/10/99 - Change in calculation method for CAPE/CIN. Trapezoid |
|---|
| 107 | * integration method is now used. Speeds up execution |
|---|
| 108 | * by 25%, and increases accuracy by 5-10%. |
|---|
| 109 | * |
|---|
| 110 | * 10/18/99 - Minor glitch fixed that occasionally caused crash. |
|---|
| 111 | * |
|---|
| 112 | * 8/26/99 - Datasets with missing data can now be used. |
|---|
| 113 | * |
|---|
| 114 | * Features: |
|---|
| 115 | * - All features of standard skewt/logp plot |
|---|
| 116 | * - RH sounding |
|---|
| 117 | * - LCL location |
|---|
| 118 | * - Parcel trajectory for both sfc based convection and elevated from |
|---|
| 119 | * most unstable level (highest theta-e level reported) |
|---|
| 120 | * - Stability indices and precipitable water calculations |
|---|
| 121 | * - CAPE & CIN Calculations |
|---|
| 122 | * - Wind Profile |
|---|
| 123 | * - Hodograph / Hodograph scaling |
|---|
| 124 | * - Helicity and SR Helicity Calculations and Display |
|---|
| 125 | * - Color aspects of output |
|---|
| 126 | * - Line Thickness, style aspects of output |
|---|
| 127 | * - Can be run in either PORTRAIT or LANDSCAPE mode. |
|---|
| 128 | * |
|---|
| 129 | * There are numerous tunable parameters below to change the structure |
|---|
| 130 | * and output for the diagram. |
|---|
| 131 | * |
|---|
| 132 | * Function Arguments: |
|---|
| 133 | * sndtemp - temperature data (Celsius) as a function of pressure |
|---|
| 134 | * snddewp - dewpoint data (Celsius) as a function of pressure |
|---|
| 135 | * sndspd - wind speed data (knots) as a function of pressure |
|---|
| 136 | * snddir - wind direction data as a function of pressure |
|---|
| 137 | * |
|---|
| 138 | * Use '-1' for any of the above 4 arguments to indicate that you |
|---|
| 139 | * are not passing that variable. The appropriate options will |
|---|
| 140 | * be ignored based on your specifying '-1' for that variable. |
|---|
| 141 | * |
|---|
| 142 | * NOTE: Make sure to set the vertical range of the plot before running. |
|---|
| 143 | * I.e., "SET LEV 1050 150", for example. This does not have to |
|---|
| 144 | * be limited to the pressure range of your data. |
|---|
| 145 | * |
|---|
| 146 | * Labelling: Pressure/Height is labelled along left side. Temperature is |
|---|
| 147 | * labelled along bottom. Mixing ratio is labelled along right |
|---|
| 148 | * side/top. |
|---|
| 149 | * |
|---|
| 150 | * |
|---|
| 151 | * PROBLEMS: First check out the web page for the script (which also |
|---|
| 152 | * has a link to a FAQ with answers to many common questions |
|---|
| 153 | * about using the script): |
|---|
| 154 | * http://www.ems.psu.edu/~hart/skew.html |
|---|
| 155 | * |
|---|
| 156 | * Please send any further problems, comments, or suggestions to |
|---|
| 157 | * <hart@ems.psu.edu> |
|---|
| 158 | * |
|---|
| 159 | * ACKNOWLEDGMENTS: Thanks go to the innumerable users who have helped |
|---|
| 160 | * fine tune the script from the horrible mess from which it began. |
|---|
| 161 | * In particular, thanks go out to Steve Lord (NCEP), Mike Fiorino (ECMWF), |
|---|
| 162 | * George Bryan (PSU), Davide Sacchetti (CMIRL), and Enrico Minguzzi (CMIRL). |
|---|
| 163 | * |
|---|
| 164 | ************************************************************************** |
|---|
| 165 | * !!!!! BEGINNING OF USER-SPECIFIED OPTIONS !!!!!! |
|---|
| 166 | ************************************************************************** |
|---|
| 167 | * |
|---|
| 168 | * --------------------- Initialization options ---------------------- |
|---|
| 169 | * |
|---|
| 170 | * ClrScrn = Whether to clear the screen before drawing diagram |
|---|
| 171 | * [1 = yes, 0 = no] |
|---|
| 172 | |
|---|
| 173 | *ClrScrn = 1 |
|---|
| 174 | |
|---|
| 175 | * |
|---|
| 176 | * ------------------- Define Skew-T Diagram Shape/Slope----------------- |
|---|
| 177 | * |
|---|
| 178 | * (P1,T1) = Pres, Temp of some point on left-most side |
|---|
| 179 | * (P2,T2) = Pres, Temp of some point on right-most side |
|---|
| 180 | * (P3,T3) = Pres, Temp of some point in diagram which is mid-point |
|---|
| 181 | * in the horizontal between 1 and 2. |
|---|
| 182 | * |
|---|
| 183 | * P1, P2, P3 are in mb ; T1, T2, T3 are in Celsius |
|---|
| 184 | * |
|---|
| 185 | * These define the SLOPE and WIDTH of the diagram as you see it but DO NOT |
|---|
| 186 | * DEFINE THE HEIGHT of the diagram as you see it. In other words, |
|---|
| 187 | * 1 and 2 do NOT necessarily need to be at the bottom of the diagram and |
|---|
| 188 | * 3 does NOT necessarily need to be at the top. THE VERTICAL PRESSURE |
|---|
| 189 | * RANGE OF THE SKEWT AS YOU SEE IT IS DETERMINED BY YOUR 'SET Z ...' |
|---|
| 190 | * COMMAND OR THE 'SET LEV ...' COMMAND BEFORE RUNNING THIS SCRIPT. |
|---|
| 191 | * |
|---|
| 192 | * _______________________ |
|---|
| 193 | * | | |
|---|
| 194 | * | | |
|---|
| 195 | * | 3 | |
|---|
| 196 | * | | |
|---|
| 197 | * | | |
|---|
| 198 | * | | |
|---|
| 199 | * | | |
|---|
| 200 | * | | |
|---|
| 201 | * | | |
|---|
| 202 | * | | |
|---|
| 203 | * | | |
|---|
| 204 | * |1 2| |
|---|
| 205 | * | | |
|---|
| 206 | * |_______________________| |
|---|
| 207 | * |
|---|
| 208 | * |
|---|
| 209 | * A good set of defining points are given below. Feel free |
|---|
| 210 | * to experiment with variations. |
|---|
| 211 | |
|---|
| 212 | |
|---|
| 213 | P1 = 1000 |
|---|
| 214 | T1 = -40 |
|---|
| 215 | |
|---|
| 216 | P2 = 1000 |
|---|
| 217 | T2 = 40 |
|---|
| 218 | |
|---|
| 219 | P3 = 200 |
|---|
| 220 | T3 = -50 |
|---|
| 221 | |
|---|
| 222 | * Another good set of defining points suggested by George Bryan (PSU) |
|---|
| 223 | * are: |
|---|
| 224 | * |
|---|
| 225 | * P1 = 1000 |
|---|
| 226 | * T1 = -30 |
|---|
| 227 | * |
|---|
| 228 | * P2 = 1000 |
|---|
| 229 | * T2 = 40 |
|---|
| 230 | * |
|---|
| 231 | * P3 = 500 |
|---|
| 232 | * T3 = -18 |
|---|
| 233 | |
|---|
| 234 | * ------------------- Contour Intervals / Levels -------------------------- |
|---|
| 235 | * |
|---|
| 236 | * All variables below are contour intervals/levels for diagram |
|---|
| 237 | * |
|---|
| 238 | * Thetaint = interval for potential temperature lines |
|---|
| 239 | * Thetwint = interval for moist pseudo adiabats |
|---|
| 240 | * tempint = interval for temperature lines |
|---|
| 241 | * wsclevs = contour LEVELS for mixing ratio lines |
|---|
| 242 | * |
|---|
| 243 | * |
|---|
| 244 | thetaint= 10 |
|---|
| 245 | thetwint= 5 |
|---|
| 246 | tempint = 10 |
|---|
| 247 | wsclevs = "1 2 3 4 6 8 10 15 20 25 30 35 40" |
|---|
| 248 | * |
|---|
| 249 | * |
|---|
| 250 | * ------------------------ Output Options -------------------------------- |
|---|
| 251 | * |
|---|
| 252 | * All variables below are logical .. 1=yes, 0=no, unless otherwise |
|---|
| 253 | * specified. |
|---|
| 254 | * |
|---|
| 255 | * DrawBarb = Draw wind barbs along right side of plot |
|---|
| 256 | * DrawThet = Draw dry adiabats |
|---|
| 257 | * DrawThtw = Draw moist pseudo-adiabats |
|---|
| 258 | * DrawTemp = Draw temperature lines |
|---|
| 259 | * DrawMix = Draw mixing ratio lines |
|---|
| 260 | * DrawTSnd = Draw temperature sounding |
|---|
| 261 | * DrawDSnd = Draw dewpoint sounding |
|---|
| 262 | * DrawRH = Draw relative humidity sounding |
|---|
| 263 | * DrawPrcl = Draw parcel path from surface upward |
|---|
| 264 | * DrawPMax = Draw parcel path from most unstable level upward |
|---|
| 265 | * DrawIndx = Display stability indices & CAPE |
|---|
| 266 | * DrawHeli = Calculate and display absolute and storm-relative helicity |
|---|
| 267 | * DrawHodo = Draw hodograph |
|---|
| 268 | * DrawPLev = Draw Pressure Levels |
|---|
| 269 | * DrawZLev = Draw height levels and lines |
|---|
| 270 | * 0 = no lines |
|---|
| 271 | * 1 = above ground level (AGL) |
|---|
| 272 | * 2 = above sea level (ASL) |
|---|
| 273 | * DrawZSTD = Draw Height levels using standard atm lapse rate |
|---|
| 274 | * LblAxes = Label the x,y axes (temperature, pressure,mixing ratio) |
|---|
| 275 | * |
|---|
| 276 | * ThtwStop = Pressure level at which to stop drawing Theta-w lines |
|---|
| 277 | * MixStop = Pressure level at which to stop drawing Mixratio lines |
|---|
| 278 | |
|---|
| 279 | DrawBarb= 1 |
|---|
| 280 | DrawThet= 1 |
|---|
| 281 | DrawThtw= 0 |
|---|
| 282 | DrawTemp= 1 |
|---|
| 283 | DrawMix = 1 |
|---|
| 284 | DrawTSnd= 1 |
|---|
| 285 | DrawDSnd= 1 |
|---|
| 286 | DrawRH = 1 |
|---|
| 287 | DrawPrcl= 1 |
|---|
| 288 | DrawPMax= 0 |
|---|
| 289 | DrawIndx= 1 |
|---|
| 290 | DrawHeli= 1 |
|---|
| 291 | DrawHodo= 1 |
|---|
| 292 | DrawPLev= 1 |
|---|
| 293 | DrawZLev= 0 |
|---|
| 294 | DrawZSTD= 0 |
|---|
| 295 | LblAxes = 1 |
|---|
| 296 | |
|---|
| 297 | ThtwStop = 200 |
|---|
| 298 | MixStop = 600 |
|---|
| 299 | |
|---|
| 300 | |
|---|
| 301 | * |
|---|
| 302 | * ----------------- Sounding Geography options ------------------------ |
|---|
| 303 | * |
|---|
| 304 | * SfcElev = Elevation above sea-level (meters) of lowest level reported |
|---|
| 305 | * in sounding. Used only if DrawZLev = 2 |
|---|
| 306 | |
|---|
| 307 | SfcElev = 0 |
|---|
| 308 | |
|---|
| 309 | |
|---|
| 310 | * |
|---|
| 311 | * ------------------ Thermodynamic Index Options -------------------- |
|---|
| 312 | * |
|---|
| 313 | * All variables here are in inches. Use -1 for the default values. |
|---|
| 314 | * |
|---|
| 315 | * Text1XC = X-location of midpoint of K,TT,PW output box |
|---|
| 316 | * Text1YC = Y-location of midpoint of K,TT,PW output box |
|---|
| 317 | * Text2XC = X-Location of midpoint of surface indices output box |
|---|
| 318 | * Text2YC = Y-location of midpoint of surface indices output box |
|---|
| 319 | * Text3XC = X-Location of midpoint of most unstable level-based indices |
|---|
| 320 | * output box |
|---|
| 321 | * Text3YC = Y-location of midpoint of most unstable level-based indices |
|---|
| 322 | * output box |
|---|
| 323 | |
|---|
| 324 | Text1XC = -1 |
|---|
| 325 | Text1YC = -1 |
|---|
| 326 | Text2XC = -1 |
|---|
| 327 | Text2YC = -1 |
|---|
| 328 | Text3XC = -1 |
|---|
| 329 | Text3YC = -1 |
|---|
| 330 | |
|---|
| 331 | * |
|---|
| 332 | * ----------------- Wind Barb Profile Options ---------------------------- |
|---|
| 333 | * |
|---|
| 334 | * All variables here are in units of inches, unless otherwise specified |
|---|
| 335 | * |
|---|
| 336 | * barbint = Interval for plotting barbs (in units of levels) |
|---|
| 337 | * poleloc = X-Location of profile. Choose -1 for the default. |
|---|
| 338 | * polelen = Length of wind-barb pole |
|---|
| 339 | * Len05 = Length of each 5-knot barb |
|---|
| 340 | * Len10 = Length of each 10-knot barb |
|---|
| 341 | * Len50 = Length of each 50-knot flag |
|---|
| 342 | * Wid50 = Width of base of 50-knot flag |
|---|
| 343 | * Spac50 = Spacing between 50-knot flag and next barb/flag |
|---|
| 344 | * Spac10 = Spacing between 10-knot flag and next flag |
|---|
| 345 | * Spac05 = Spacing between 5-knot flag and next flag |
|---|
| 346 | * Flagbase= Draw flagbase (filled circle) for each windbarb [1=yes, 0 =no] |
|---|
| 347 | * Fill50 = Solid-fill 50-knot flag [1=yes, 0=no] |
|---|
| 348 | * barbline= Draw a vertical line connecting all the wind barbs [1=yes, 0=no] |
|---|
| 349 | * |
|---|
| 350 | barbint = 1 |
|---|
| 351 | poleloc = -1 |
|---|
| 352 | polelen = 0.35 |
|---|
| 353 | len05 = 0.07 |
|---|
| 354 | len10 = 0.15 |
|---|
| 355 | len50 = 0.15 |
|---|
| 356 | wid50 = 0.06 |
|---|
| 357 | spac50 = 0.07 |
|---|
| 358 | spac10 = 0.05 |
|---|
| 359 | spac05 = 0.05 |
|---|
| 360 | Fill50 = 1 |
|---|
| 361 | flagbase= 1 |
|---|
| 362 | barbline= 1 |
|---|
| 363 | |
|---|
| 364 | * |
|---|
| 365 | * |
|---|
| 366 | *---------------- Hodograph Options ------------------------------------- |
|---|
| 367 | * |
|---|
| 368 | * All variables here are in units of inches, unless otherwise specified |
|---|
| 369 | * |
|---|
| 370 | * HodXcent= x-location of hodograph center. Use -1 for default location. |
|---|
| 371 | * HodYcent= y-location of hodograph center. Use -1 for default location. |
|---|
| 372 | * HodSize = Size of hodograph in inches |
|---|
| 373 | * NumRing = Number of rings to place in hodograph (must be at least 1) |
|---|
| 374 | * HodRing = Wind speed increment of each hodograph ring |
|---|
| 375 | * HodoDep = Depth (above lowest level in mb) of end of hodograph trace |
|---|
| 376 | * TickInt = Interval (in kts) at which tick marks are drawn along the axes |
|---|
| 377 | * Use 0 for no tick marks. |
|---|
| 378 | * TickSize= Size of tick mark in inches |
|---|
| 379 | * Text4XC = X-location of midpoint of hodograph text output. Use -1 for default. |
|---|
| 380 | * Text4YC = Y-location of midpoint of hodograph text output. Use -1 for default. |
|---|
| 381 | |
|---|
| 382 | HodXcent= -1 |
|---|
| 383 | HodYcent= -1 |
|---|
| 384 | HodSize = 2 |
|---|
| 385 | NumRing = 4 |
|---|
| 386 | HodRing = 10 |
|---|
| 387 | HodoDep = 300 |
|---|
| 388 | TickInt = 5 |
|---|
| 389 | TickSize= 0.05 |
|---|
| 390 | Text4XC = -1 |
|---|
| 391 | Text4YC = -1 |
|---|
| 392 | |
|---|
| 393 | *--------------- Helicity Options --------------------------------------- |
|---|
| 394 | * |
|---|
| 395 | * MeanVTop = Top pressure level (mb) of mean-wind calculation |
|---|
| 396 | * MeanVBot = Bottom pressure level (mb) of mean-wind calculation |
|---|
| 397 | * HelicDep = Depth in mb (above ground) of helicity integration |
|---|
| 398 | * StormMot = Type of storm motion estimation scheme. Use following: |
|---|
| 399 | * 0 = No departure from mean wind. |
|---|
| 400 | * 1 = Davies-Jones (1990) approach |
|---|
| 401 | * FillArrw = Whether to fill the arrowhead of the storm motion vector |
|---|
| 402 | * [1 = yes, 0 = no] |
|---|
| 403 | |
|---|
| 404 | MeanVTop= 300 |
|---|
| 405 | MeanVBot= 850 |
|---|
| 406 | HelicDep= 300 |
|---|
| 407 | StormMot= 1 |
|---|
| 408 | FillArrw= 1 |
|---|
| 409 | |
|---|
| 410 | * |
|---|
| 411 | *---------------- Color Options ------------------------------------------ |
|---|
| 412 | * |
|---|
| 413 | * ThetCol = Color of dry adiabats |
|---|
| 414 | * TempCol = Color of temperature lines |
|---|
| 415 | * MixCol = Color of mixing ratio lines |
|---|
| 416 | * ThtwCol = Color of moist adiabats |
|---|
| 417 | * TSndCol = Color of Temperature Sounding |
|---|
| 418 | * DSndCol = Color of Dewpoint Sounding |
|---|
| 419 | * RHCol = Color of RH Sounding |
|---|
| 420 | * PrclCol = Color of parcel trace |
|---|
| 421 | * BarbCol = Color of wind barbs (choose -1 for color according to speed) |
|---|
| 422 | * HodoCol = Color of hodograph trace |
|---|
| 423 | |
|---|
| 424 | ThetCol = 2 |
|---|
| 425 | TempCol = 4 |
|---|
| 426 | MixCol = 7 |
|---|
| 427 | ThtwCol = 3 |
|---|
| 428 | *TSndCol = 1 |
|---|
| 429 | *DSndCol = 1 |
|---|
| 430 | *RHCol = 10 |
|---|
| 431 | *PrclCol = 5 |
|---|
| 432 | *BarbCol = -1 |
|---|
| 433 | HodoCol = 1 |
|---|
| 434 | |
|---|
| 435 | * |
|---|
| 436 | *-------------------- Line Style Options ------------------------------------ |
|---|
| 437 | * |
|---|
| 438 | * GrADS Styles: 1=solid;2=long dash;3=short dash;4=long,short dashed; |
|---|
| 439 | * 5=dotted;6=dot dash;7=dot dot dash |
|---|
| 440 | * |
|---|
| 441 | * ThetLine = Line Style of dry adiabats |
|---|
| 442 | * TempLine = Line Style of temperature lines |
|---|
| 443 | * MixLine = Line Style of mixing ratio lines |
|---|
| 444 | * ThtwLine = Line Style of moist adiabats |
|---|
| 445 | * TSndLine = Line Style of Temperature Sounding |
|---|
| 446 | * DSndLine = Line Style of Dewpoint Sounding |
|---|
| 447 | * RHLine = Line Style of RH sounding |
|---|
| 448 | * PrclLine = Line Style of parcel trace |
|---|
| 449 | * HodoLine = Line Style of hodograph trace |
|---|
| 450 | * |
|---|
| 451 | |
|---|
| 452 | ThetLine = 1 |
|---|
| 453 | TempLine = 1 |
|---|
| 454 | MixLine = 5 |
|---|
| 455 | ThtwLine = 3 |
|---|
| 456 | TSndLine = 1 |
|---|
| 457 | DSndLine = 1 |
|---|
| 458 | RHLine = 1 |
|---|
| 459 | PrclLine = 3 |
|---|
| 460 | HodoLine = 1 |
|---|
| 461 | |
|---|
| 462 | * |
|---|
| 463 | *------------------- Line Thickness Options--------------------------------- |
|---|
| 464 | * GrADS Line Thickness: increases with increasing number. Influences |
|---|
| 465 | * hardcopy output more strongly than screen output. |
|---|
| 466 | * |
|---|
| 467 | * |
|---|
| 468 | * ThetThk = Line Thickness of dry adiabats |
|---|
| 469 | * TempThk = Line Thickness of temperature lines |
|---|
| 470 | * MixThk = Line Thickness of mixing ratio lines |
|---|
| 471 | * ThtwThk = Line Thickness of moist adiabats |
|---|
| 472 | * TSndThk = Line Thickness of temperature sounding |
|---|
| 473 | * DSndThk = Line thickness of dewpoint sounding |
|---|
| 474 | * RHThk = Line thickness of RH sounding |
|---|
| 475 | * PrclThk = Line thickness of parcel trace |
|---|
| 476 | * HodoThk = Line thickness of hodograph trace |
|---|
| 477 | * BarbThk = Line thickness of wind barbs |
|---|
| 478 | |
|---|
| 479 | ThetThk = 3 |
|---|
| 480 | TempThk = 3 |
|---|
| 481 | MixThk = 3 |
|---|
| 482 | ThtwThk = 3 |
|---|
| 483 | TSndThk = 8 |
|---|
| 484 | DSndThk = 8 |
|---|
| 485 | RHThk = 8 |
|---|
| 486 | PrclThk = 6 |
|---|
| 487 | HodoThk = 3 |
|---|
| 488 | BarbThk = 2 |
|---|
| 489 | |
|---|
| 490 | * |
|---|
| 491 | *------------------- Data Point Marker Options ----------------------------- |
|---|
| 492 | * GrADS Marker Types: 0 = none ; 1 = cross ; 2 = open circle ; |
|---|
| 493 | * 3 = closed circle ; 4 = open square ; 5 = closed square |
|---|
| 494 | * 6 = X ; 7 = diamond ; 8 = triangle ; 9 = none |
|---|
| 495 | * 10 = open circle with vertical line ; 11 = open oval |
|---|
| 496 | * |
|---|
| 497 | * TSndMrk = Mark type of data point marker for temperature sounding |
|---|
| 498 | * DSndMrk = Mark type of data point marker for dewpoint sounding |
|---|
| 499 | * RHMrk = Mark type of data point marker for relative humidity sounding |
|---|
| 500 | * MrkSize = Mark size (inches) of each data marker |
|---|
| 501 | |
|---|
| 502 | TSndMrk = 0 |
|---|
| 503 | DSndMrk = 0 |
|---|
| 504 | RHMrk = 0 |
|---|
| 505 | MrkSize = 0.1 |
|---|
| 506 | |
|---|
| 507 | |
|---|
| 508 | * !!!!! YOU SHOULD NOT NEED TO CHANGE ANYTHING BELOW HERE !!!!! |
|---|
| 509 | **************************************************************************** |
|---|
| 510 | |
|---|
| 511 | *------------------------------------------- |
|---|
| 512 | * grab user-specified environment dimensions |
|---|
| 513 | *------------------------------------------- |
|---|
| 514 | |
|---|
| 515 | "q dims" |
|---|
| 516 | rec=sublin(result,2) |
|---|
| 517 | _xtype=subwrd(rec,3) |
|---|
| 518 | _xlon=subwrd(rec,6) |
|---|
| 519 | _xval=subwrd(rec,9) |
|---|
| 520 | rec=sublin(result,3) |
|---|
| 521 | _yval=subwrd(rec,9) |
|---|
| 522 | _ylat=subwrd(rec,6) |
|---|
| 523 | _ytype=subwrd(rec,3) |
|---|
| 524 | rec=sublin(result,4) |
|---|
| 525 | _ptype=subwrd(rec,3) |
|---|
| 526 | _pmax=subwrd(rec,6) |
|---|
| 527 | _pmin=subwrd(rec,8) |
|---|
| 528 | _zmin=subwrd(rec,11) |
|---|
| 529 | _zmax=subwrd(rec,13) |
|---|
| 530 | rec=sublin(result,5) |
|---|
| 531 | _ttype=subwrd(rec,3) |
|---|
| 532 | _ttime=subwrd(rec,6) |
|---|
| 533 | _tval=subwrd(rec,9) |
|---|
| 534 | |
|---|
| 535 | "q file" |
|---|
| 536 | rec=sublin(result,5) |
|---|
| 537 | _zmaxfile=subwrd(rec,9) |
|---|
| 538 | |
|---|
| 539 | *------------------------------------------------------------- |
|---|
| 540 | * Check to ensure that dimensions are valid. Warn & exit if not. |
|---|
| 541 | *-------------------------------------------------------------- |
|---|
| 542 | |
|---|
| 543 | dimrc=0 |
|---|
| 544 | If (_xtype != "fixed") |
|---|
| 545 | say "X-Dims Error: Not fixed. Use 'set lon' or 'set x' to specify a value." |
|---|
| 546 | dimrc=-1 |
|---|
| 547 | Endif |
|---|
| 548 | |
|---|
| 549 | If (_ytype != "fixed") |
|---|
| 550 | say "Y-Dims Error: Not fixed. Use 'set lat' or 'set y' to specify a value" |
|---|
| 551 | dimrc=-1 |
|---|
| 552 | Endif |
|---|
| 553 | |
|---|
| 554 | If (_ptype != "varying") |
|---|
| 555 | say "Z-Dims Error: Not varying. Use 'set lev' or 'set z' to specify a range." |
|---|
| 556 | dimrc=-1 |
|---|
| 557 | Endif |
|---|
| 558 | |
|---|
| 559 | If (_ttype != "fixed") |
|---|
| 560 | say "Time Error: Not fixed. Use 'set time' or 'set t' to specify a value" |
|---|
| 561 | dimrc=-1 |
|---|
| 562 | Endif |
|---|
| 563 | |
|---|
| 564 | |
|---|
| 565 | If (dimrc < 0) |
|---|
| 566 | Return(-1) |
|---|
| 567 | Endif |
|---|
| 568 | |
|---|
| 569 | |
|---|
| 570 | * |
|---|
| 571 | * A few global variables used in units conversion |
|---|
| 572 | * |
|---|
| 573 | |
|---|
| 574 | _pi=3.14159265 |
|---|
| 575 | _dtr=_pi/180 |
|---|
| 576 | _rtd=1/_dtr |
|---|
| 577 | _ktm=0.514444 |
|---|
| 578 | _mtk=1/_ktm |
|---|
| 579 | |
|---|
| 580 | * A few global constants used in thermo calcs |
|---|
| 581 | |
|---|
| 582 | _C0=0.99999683 |
|---|
| 583 | _C1=-0.90826951/100 |
|---|
| 584 | _C2= 0.78736169/10000 |
|---|
| 585 | _C3=-0.61117958/1000000 |
|---|
| 586 | _C4= 0.43884187/pow(10,8) |
|---|
| 587 | _C5=-0.29883885/pow(10,10) |
|---|
| 588 | _C6= 0.21874425/pow(10,12) |
|---|
| 589 | _C7=-0.17892321/pow(10,14) |
|---|
| 590 | _C8= 0.11112018/pow(10,16) |
|---|
| 591 | _C9=-0.30994571/pow(10,19) |
|---|
| 592 | |
|---|
| 593 | * A pressure array of power calculations which should be performed |
|---|
| 594 | * only once to reduce execution time. |
|---|
| 595 | |
|---|
| 596 | zz=1100 |
|---|
| 597 | while (zz > 10) |
|---|
| 598 | subscr=zz/10 |
|---|
| 599 | _powpres.subscr=pow(zz,0.286) |
|---|
| 600 | zz=zz-10 |
|---|
| 601 | endwhile |
|---|
| 602 | |
|---|
| 603 | * |
|---|
| 604 | * Turn off options not available due to user data limitations |
|---|
| 605 | * |
|---|
| 606 | |
|---|
| 607 | If (ClrScrn = 1) |
|---|
| 608 | "clear" |
|---|
| 609 | Endif |
|---|
| 610 | |
|---|
| 611 | If (sndspd = -1 | snddir = -1) |
|---|
| 612 | DrawBarb = 0 |
|---|
| 613 | DrawHodo = 0 |
|---|
| 614 | DrawHeli = 0 |
|---|
| 615 | Endif |
|---|
| 616 | |
|---|
| 617 | If (snddewp = -1) |
|---|
| 618 | DrawDSnd = 0 |
|---|
| 619 | DrawRH = 0 |
|---|
| 620 | DrawPrcl = 0 |
|---|
| 621 | DrawPMax = 0 |
|---|
| 622 | DrawIndx = 0 |
|---|
| 623 | Endif |
|---|
| 624 | |
|---|
| 625 | If (sndtemp = -1) |
|---|
| 626 | DrawTSnd = 0 |
|---|
| 627 | DrawRH = 0 |
|---|
| 628 | DrawPrcl = 0 |
|---|
| 629 | DrawPMax = 0 |
|---|
| 630 | DrawIndx = 0 |
|---|
| 631 | DrawZLev = 0 |
|---|
| 632 | Endif |
|---|
| 633 | |
|---|
| 634 | If (NumRing < 1) |
|---|
| 635 | DrawHodo = 0 |
|---|
| 636 | Endif |
|---|
| 637 | |
|---|
| 638 | "q gxinfo" |
|---|
| 639 | rec=sublin(result,2) |
|---|
| 640 | xsize=subwrd(rec,4) |
|---|
| 641 | |
|---|
| 642 | If (xsize = 11) |
|---|
| 643 | PageType = "Landscape" |
|---|
| 644 | Else |
|---|
| 645 | PageType = "Portrait" |
|---|
| 646 | Endif |
|---|
| 647 | |
|---|
| 648 | *------------------------------------------------------ |
|---|
| 649 | * calculate constants determining slope/shape of diagram |
|---|
| 650 | * based on temp/pressure values given by user |
|---|
| 651 | *------------------------------------------------------- |
|---|
| 652 | |
|---|
| 653 | "set x 1" |
|---|
| 654 | "set y 1" |
|---|
| 655 | "set z 1" |
|---|
| 656 | "set t 1" |
|---|
| 657 | _m1=(T1+T2-2*T3)/(2*log10(P2/P3)) |
|---|
| 658 | _m2=(T2-T3-_m1*log10(P2/P3))/50 |
|---|
| 659 | _m3=(T1-_m1*log10(P1)) |
|---|
| 660 | |
|---|
| 661 | "set z "_zmin" "_zmax |
|---|
| 662 | "set zlog on" |
|---|
| 663 | "set xlab off" |
|---|
| 664 | |
|---|
| 665 | *------------------------------------------------- |
|---|
| 666 | * perform coordinate transformation to Skew-T/LogP |
|---|
| 667 | *------------------------------------------------- |
|---|
| 668 | |
|---|
| 669 | "set gxout stat" |
|---|
| 670 | "set x "_xval |
|---|
| 671 | "set y "_yval |
|---|
| 672 | "set t "_tval |
|---|
| 673 | "define tempx=("sndtemp"-"_m1"*log10(lev)-"_m3")/"_m2 |
|---|
| 674 | "define dewpx=("snddewp"-"_m1"*log10(lev)-"_m3")/"_m2 |
|---|
| 675 | |
|---|
| 676 | If (PageType = "Portrait") |
|---|
| 677 | "set parea 0.7 7 0.75 10.5" |
|---|
| 678 | Else |
|---|
| 679 | "set parea 0.7 6.5 0.5 8" |
|---|
| 680 | Endif |
|---|
| 681 | |
|---|
| 682 | "set axlim 0 100" |
|---|
| 683 | "set lon 0 100" |
|---|
| 684 | "set grid on 1 1" |
|---|
| 685 | |
|---|
| 686 | "set z "_zmin " " _zmax |
|---|
| 687 | "set lon 0 100" |
|---|
| 688 | "set clevs -900" |
|---|
| 689 | "set gxout contour" |
|---|
| 690 | |
|---|
| 691 | *------------------------------------- |
|---|
| 692 | * Draw pressure lines |
|---|
| 693 | *------------------------------------- |
|---|
| 694 | |
|---|
| 695 | If (DrawPLev = 0) |
|---|
| 696 | "set ylab off" |
|---|
| 697 | Else |
|---|
| 698 | "set ylab on" |
|---|
| 699 | "set ylopts 1 3 0.10" |
|---|
| 700 | "set xlopts 1 3 0.125" |
|---|
| 701 | Endif |
|---|
| 702 | |
|---|
| 703 | "d lon" |
|---|
| 704 | |
|---|
| 705 | *-------------------------------------- |
|---|
| 706 | * Determine corners of skewt/logp frame |
|---|
| 707 | *-------------------------------------- |
|---|
| 708 | |
|---|
| 709 | "q w2xy 100 "_pmin |
|---|
| 710 | rxloc=subwrd(result,3) |
|---|
| 711 | tyloc=subwrd(result,6) |
|---|
| 712 | "q w2xy 0 "_pmax |
|---|
| 713 | lxloc=subwrd(result,3) |
|---|
| 714 | byloc=subwrd(result,6) |
|---|
| 715 | |
|---|
| 716 | If (DrawPLev = 1 & LblAxes = 1) |
|---|
| 717 | "set strsiz 0.10" |
|---|
| 718 | "set string 1 c 3 0" |
|---|
| 719 | If (PageType = "Portrait") |
|---|
| 720 | "draw string 0.5 10.85 mb" |
|---|
| 721 | Else |
|---|
| 722 | "draw string 0.5 8.35 mb" |
|---|
| 723 | Endif |
|---|
| 724 | Endif |
|---|
| 725 | |
|---|
| 726 | *--------------------------------------------------- |
|---|
| 727 | * Calculate & draw actual height lines using temp data |
|---|
| 728 | *--------------------------------------------------- |
|---|
| 729 | |
|---|
| 730 | If (DrawZLev > 0) |
|---|
| 731 | say "Calculating observed height levels from temp/pressure data." |
|---|
| 732 | zz=1 |
|---|
| 733 | "set gxout stat" |
|---|
| 734 | "set x "_xval |
|---|
| 735 | "set y "_yval |
|---|
| 736 | "set t "_tval |
|---|
| 737 | count=0 |
|---|
| 738 | while (zz < _zmax) |
|---|
| 739 | "set z "zz |
|---|
| 740 | pp.zz=subwrd(result,4) |
|---|
| 741 | lpp.zz=log(pp.zz) |
|---|
| 742 | "d "sndtemp |
|---|
| 743 | rec=sublin(result,1) |
|---|
| 744 | check=substr(rec,1,6) |
|---|
| 745 | If (check = "Notice") |
|---|
| 746 | rec=sublin(result,9) |
|---|
| 747 | Else |
|---|
| 748 | rec=sublin(result,8) |
|---|
| 749 | Endif |
|---|
| 750 | tt=subwrd(rec,4) |
|---|
| 751 | if (tt > -900) |
|---|
| 752 | tk=tt+273.15 |
|---|
| 753 | count=count+1 |
|---|
| 754 | zzm=zz-1 |
|---|
| 755 | If (count = 1) |
|---|
| 756 | If (DrawZLev = 2) |
|---|
| 757 | htlb="ASL" |
|---|
| 758 | height.zz=SfcElev |
|---|
| 759 | Else |
|---|
| 760 | htlb="AGL" |
|---|
| 761 | height.zz=0 |
|---|
| 762 | Endif |
|---|
| 763 | sfcz=height.zz |
|---|
| 764 | Else |
|---|
| 765 | DZ=29.2857*(lpp.zzm-lpp.zz)*(lpp.zz*tk+lpp.zzm*tkold)/(lpp.zz+lpp.zzm) |
|---|
| 766 | height.zz=height.zzm+DZ |
|---|
| 767 | highz=height.zz |
|---|
| 768 | Endif |
|---|
| 769 | else |
|---|
| 770 | height.zz = -9999 |
|---|
| 771 | endif |
|---|
| 772 | tkold=tk |
|---|
| 773 | zz=zz+1 |
|---|
| 774 | endwhile |
|---|
| 775 | |
|---|
| 776 | maxht=int(highz/1000) |
|---|
| 777 | if (int(sfcz/1000) = sfcz/1000) |
|---|
| 778 | minht=int(sfcz/1000) |
|---|
| 779 | else |
|---|
| 780 | minht=1+int(sfcz/1000) |
|---|
| 781 | endif |
|---|
| 782 | |
|---|
| 783 | ht=minht |
|---|
| 784 | "set line 1 3 1" |
|---|
| 785 | "set strsiz 0.10" |
|---|
| 786 | "set string 1 l 3 0" |
|---|
| 787 | while (ht <= maxht) |
|---|
| 788 | zz=1 |
|---|
| 789 | while (height.zz/1000 <= ht) |
|---|
| 790 | zz=zz+1 |
|---|
| 791 | endwhile |
|---|
| 792 | zzm=zz-1 |
|---|
| 793 | PBelow=pp.zzm |
|---|
| 794 | PAbove=pp.zz |
|---|
| 795 | HBelow=height.zzm |
|---|
| 796 | HAbove=height.zz |
|---|
| 797 | DZ=HAbove-HBelow |
|---|
| 798 | DP=PAbove-PBelow |
|---|
| 799 | Del=ht*1000-HBelow |
|---|
| 800 | Est=PBelow+Del*DP/DZ |
|---|
| 801 | If (Est >= _pmin & Est <= _pmax) |
|---|
| 802 | "q w2xy 1 " Est |
|---|
| 803 | yloc=subwrd(result,6) |
|---|
| 804 | "draw line " lxloc " " yloc " " rxloc " " yloc |
|---|
| 805 | "draw string 0.22 "yloc-0.05" "ht |
|---|
| 806 | Endif |
|---|
| 807 | ht=ht+1 |
|---|
| 808 | endwhile |
|---|
| 809 | "set strsiz 0.10" |
|---|
| 810 | "set string 1" |
|---|
| 811 | If (LblAxes = 1) |
|---|
| 812 | If (PageType = "Portrait") |
|---|
| 813 | "draw string 0.25 10.85 km" |
|---|
| 814 | "draw string 0.25 10.75 "htlb |
|---|
| 815 | "draw string 0.25 10.65 OBS" |
|---|
| 816 | Else |
|---|
| 817 | "draw string 0.25 8.35 km" |
|---|
| 818 | "draw string 0.25 8.25 "htlb |
|---|
| 819 | "draw string 0.25 8.15 OBS" |
|---|
| 820 | Endif |
|---|
| 821 | Endif |
|---|
| 822 | Endif |
|---|
| 823 | |
|---|
| 824 | |
|---|
| 825 | *--------------------------------------------------- |
|---|
| 826 | * Draw height levels (height above MSL using Std Atm) |
|---|
| 827 | *--------------------------------------------------- |
|---|
| 828 | |
|---|
| 829 | If (DrawZSTD = 1) |
|---|
| 830 | "set strsiz 0.10" |
|---|
| 831 | minht=30.735*(1-pow(_pmax/1013.26,0.287)) |
|---|
| 832 | minht=int(minht+0.5) |
|---|
| 833 | maxht=30.735*(1-pow(_pmin/1013.26,0.287)) |
|---|
| 834 | maxht=int(maxht) |
|---|
| 835 | "set gxout stat" |
|---|
| 836 | zcount=minht |
|---|
| 837 | while (zcount <= maxht) |
|---|
| 838 | plev=1013.26*pow((1-zcount/30.735),3.4843) |
|---|
| 839 | "q w2xy 0 "plev |
|---|
| 840 | yloc=subwrd(result,6) |
|---|
| 841 | "draw string 0 "yloc-0.05" "zcount |
|---|
| 842 | zcount=zcount+1 |
|---|
| 843 | endwhile |
|---|
| 844 | "set strsiz 0.10" |
|---|
| 845 | If (LblAxes = 1) |
|---|
| 846 | If (PageType = "Portrait") |
|---|
| 847 | "draw string 0 10.85 km" |
|---|
| 848 | "draw string 0 10.75 ASL" |
|---|
| 849 | "draw string 0 10.65 STD" |
|---|
| 850 | Else |
|---|
| 851 | "draw string 0 8.35 km" |
|---|
| 852 | "draw string 0 8.25 ASL" |
|---|
| 853 | "draw string 0 8.15 STD" |
|---|
| 854 | Endif |
|---|
| 855 | Endif |
|---|
| 856 | Endif |
|---|
| 857 | |
|---|
| 858 | |
|---|
| 859 | *----------------------- |
|---|
| 860 | * Plot temperature lines |
|---|
| 861 | *----------------------- |
|---|
| 862 | |
|---|
| 863 | If (DrawTemp = 1) |
|---|
| 864 | "set strsiz 0.1" |
|---|
| 865 | "set z "_zmin " " _zmax |
|---|
| 866 | "set line "TempCol " " TempLine " "TempThk |
|---|
| 867 | "set string 1 c 3 0" |
|---|
| 868 | "set gxout stat" |
|---|
| 869 | maxtline=GetTemp(100,_pmax) |
|---|
| 870 | mintline=GetTemp(0,_pmin) |
|---|
| 871 | |
|---|
| 872 | maxtline=tempint*int(maxtline/tempint) |
|---|
| 873 | mintline=tempint*int(mintline/tempint) |
|---|
| 874 | |
|---|
| 875 | tloop=mintline |
|---|
| 876 | While (tloop <= maxtline) |
|---|
| 877 | Botxtemp=GetXLoc(tloop,_pmax) |
|---|
| 878 | "q w2xy "Botxtemp " " _pmax |
|---|
| 879 | Botxloc=subwrd(result,3) |
|---|
| 880 | Botyloc=byloc |
|---|
| 881 | Topxtemp=GetXLoc(tloop,_pmin) |
|---|
| 882 | "q w2xy "Topxtemp " " _pmin |
|---|
| 883 | Topxloc=subwrd(result,3) |
|---|
| 884 | Topyloc=tyloc |
|---|
| 885 | If (Botxtemp <= 100 | Topxtemp <= 100) |
|---|
| 886 | If (Topxtemp > 100) |
|---|
| 887 | Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp) |
|---|
| 888 | b=Topyloc-Slope*Topxtemp |
|---|
| 889 | Topyloc=Slope*100+b |
|---|
| 890 | Topxloc=rxloc |
|---|
| 891 | Endif |
|---|
| 892 | If (Botxtemp < 0) |
|---|
| 893 | Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp) |
|---|
| 894 | b=Botyloc-Slope*Botxtemp |
|---|
| 895 | Botyloc=b |
|---|
| 896 | Botxloc=lxloc |
|---|
| 897 | Else |
|---|
| 898 | "draw string " Botxloc-0.05 " " Botyloc-0.15 " " tloop |
|---|
| 899 | Endif |
|---|
| 900 | "draw line "Botxloc " " Botyloc " " Topxloc " " Topyloc |
|---|
| 901 | Endif |
|---|
| 902 | tloop=tloop+tempint |
|---|
| 903 | EndWhile |
|---|
| 904 | If (LblAxes = 1) |
|---|
| 905 | "set strsiz 0.15" |
|---|
| 906 | "set string 1 c" |
|---|
| 907 | If (PageType = "Portrait") |
|---|
| 908 | "draw string 4.0 0.35 Temperature (`3.`0C)" |
|---|
| 909 | Else |
|---|
| 910 | "draw string 3.5 0.15 Temperature (`3.`0C)" |
|---|
| 911 | Endif |
|---|
| 912 | Endif |
|---|
| 913 | Endif |
|---|
| 914 | |
|---|
| 915 | |
|---|
| 916 | *------------------ |
|---|
| 917 | * Plot dry adiabats |
|---|
| 918 | *------------------ |
|---|
| 919 | |
|---|
| 920 | If (DrawThet = 1) |
|---|
| 921 | temp=GetTemp(100,_pmin) |
|---|
| 922 | maxtheta=GetThet2(temp,-100,_pmin) |
|---|
| 923 | maxtheta=thetaint*int(maxtheta/thetaint) |
|---|
| 924 | temp=GetTemp(0,_pmax) |
|---|
| 925 | mintheta=GetThet2(temp,-100,_pmax) |
|---|
| 926 | mintheta=thetaint*int(mintheta/thetaint) |
|---|
| 927 | |
|---|
| 928 | "set lon 0 100" |
|---|
| 929 | "set y 1" |
|---|
| 930 | "set z 1" |
|---|
| 931 | tloop=mintheta |
|---|
| 932 | "set line "ThetCol" "ThetLine " "ThetThk |
|---|
| 933 | While (tloop <= maxtheta) |
|---|
| 934 | PTemp=LiftDry(tloop,1000,_pmin,1,_pmin,_pmax) |
|---|
| 935 | tloop=tloop+thetaint |
|---|
| 936 | Endwhile |
|---|
| 937 | Endif |
|---|
| 938 | |
|---|
| 939 | *------------------------ |
|---|
| 940 | * Plot mixing ratio lines |
|---|
| 941 | *------------------------ |
|---|
| 942 | |
|---|
| 943 | If (DrawMix = 1) |
|---|
| 944 | If (MixStop < _pmin) |
|---|
| 945 | MixStop = _pmin |
|---|
| 946 | Endif |
|---|
| 947 | "set string 1 l" |
|---|
| 948 | "set z "_zmin " " _zmax |
|---|
| 949 | "set cint 1" |
|---|
| 950 | "set line "MixCol" " MixLine " "MixThk |
|---|
| 951 | cont = 1 |
|---|
| 952 | mloop=subwrd(wsclevs,1) |
|---|
| 953 | count = 1 |
|---|
| 954 | While (cont = 1) |
|---|
| 955 | BotCoef=log(mloop*_pmax/3801.66) |
|---|
| 956 | BotTval=-245.5*BotCoef/(BotCoef-17.67) |
|---|
| 957 | Botxtemp=GetXLoc(BotTval,_pmax) |
|---|
| 958 | "q w2xy "Botxtemp " " _pmax |
|---|
| 959 | Botxloc=subwrd(result,3) |
|---|
| 960 | Botyloc=byloc |
|---|
| 961 | TopCoef=log(mloop*MixStop/3801.66) |
|---|
| 962 | TopTval=-245.5*TopCoef/(TopCoef-17.67) |
|---|
| 963 | Topxtemp=GetXLoc(TopTval,MixStop) |
|---|
| 964 | "q w2xy "Topxtemp " " MixStop |
|---|
| 965 | Topxloc=subwrd(result,3) |
|---|
| 966 | Topyloc=subwrd(result,6) |
|---|
| 967 | "set string "MixCol" l 3" |
|---|
| 968 | "set strsiz 0.09" |
|---|
| 969 | If (Botxtemp <= 100 | Topxtemp <= 100) |
|---|
| 970 | If (Topxtemp > 100) |
|---|
| 971 | Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp) |
|---|
| 972 | b=Topyloc-Slope*Topxtemp |
|---|
| 973 | Topyloc=Slope*100+b |
|---|
| 974 | Topxloc=rxloc |
|---|
| 975 | "draw string " Topxloc+0.05 " " Topyloc " " mloop |
|---|
| 976 | Else |
|---|
| 977 | "draw string " Topxloc " " Topyloc+0.1 " " mloop |
|---|
| 978 | Endif |
|---|
| 979 | If (Botxtemp < 0) |
|---|
| 980 | Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp) |
|---|
| 981 | b=Botyloc-Slope*Botxtemp |
|---|
| 982 | Botyloc=b |
|---|
| 983 | Botxloc=lxloc |
|---|
| 984 | Endif |
|---|
| 985 | "draw line "Botxloc " " Botyloc " " Topxloc " " Topyloc |
|---|
| 986 | Endif |
|---|
| 987 | count=count+1 |
|---|
| 988 | mloop=subwrd(wsclevs,count) |
|---|
| 989 | If (mloop = "" | count > 50) |
|---|
| 990 | cont = 0 |
|---|
| 991 | Endif |
|---|
| 992 | EndWhile |
|---|
| 993 | If (LblAxes = 1) |
|---|
| 994 | "set strsiz 0.15" |
|---|
| 995 | "set string 1 c 3 90" |
|---|
| 996 | If (PageType = "Portrait") |
|---|
| 997 | "draw string 7.40 4.75 Mixing Ratio (g/kg)" |
|---|
| 998 | Else |
|---|
| 999 | "draw string 6.90 4.25 Mixing Ratio (g/kg)" |
|---|
| 1000 | Endif |
|---|
| 1001 | "set string 1 c 3 0" |
|---|
| 1002 | Endif |
|---|
| 1003 | Endif |
|---|
| 1004 | |
|---|
| 1005 | *----------------------------- |
|---|
| 1006 | * Plot moist (pseudo) adiabats |
|---|
| 1007 | *----------------------------- |
|---|
| 1008 | |
|---|
| 1009 | If (DrawThtw = 1) |
|---|
| 1010 | "set lon 0 100" |
|---|
| 1011 | "set y 1" |
|---|
| 1012 | "set z 1" |
|---|
| 1013 | "set gxout stat" |
|---|
| 1014 | tloop=80 |
|---|
| 1015 | "set line "ThtwCol" "ThtwLine " "ThtwThk |
|---|
| 1016 | While (tloop > -80) |
|---|
| 1017 | PTemp=LiftWet(tloop,1000,ThtwStop,1,_pmin,_pmax) |
|---|
| 1018 | tloop=tloop-thetwint |
|---|
| 1019 | Endwhile |
|---|
| 1020 | Endif |
|---|
| 1021 | |
|---|
| 1022 | |
|---|
| 1023 | *----------------------------------------------------- |
|---|
| 1024 | * Plot transformed user-specified temperature sounding |
|---|
| 1025 | *----------------------------------------------------- |
|---|
| 1026 | |
|---|
| 1027 | If (DrawTSnd = 1) |
|---|
| 1028 | say "Drawing temperature sounding." |
|---|
| 1029 | "set gxout line" |
|---|
| 1030 | "set x "_xval |
|---|
| 1031 | "set y "_yval |
|---|
| 1032 | "set z "_zmin" "_zmax |
|---|
| 1033 | "set ccolor "TSndCol |
|---|
| 1034 | "set cstyle "TSndLine |
|---|
| 1035 | "set cmark "TSndMrk |
|---|
| 1036 | "set digsiz "MrkSize |
|---|
| 1037 | "set cthick "TSndThk |
|---|
| 1038 | "set missconn on" |
|---|
| 1039 | "d tempx" |
|---|
| 1040 | Endif |
|---|
| 1041 | |
|---|
| 1042 | *--------------------------------------------------- |
|---|
| 1043 | * Plot transformed user-specified dewpoint sounding |
|---|
| 1044 | *--------------------------------------------------- |
|---|
| 1045 | |
|---|
| 1046 | If (DrawDSnd = 1) |
|---|
| 1047 | say "Drawing dewpoint sounding." |
|---|
| 1048 | "set gxout line" |
|---|
| 1049 | "set x "_xval |
|---|
| 1050 | "set y "_yval |
|---|
| 1051 | "set z "_zmin" "_zmax |
|---|
| 1052 | "set cmark "DSndMrk |
|---|
| 1053 | "set digsiz "MrkSize |
|---|
| 1054 | "set ccolor "DSndCol |
|---|
| 1055 | "set cstyle "DSndLine |
|---|
| 1056 | "set cthick "DSndThk |
|---|
| 1057 | "set missconn on" |
|---|
| 1058 | "d dewpx" |
|---|
| 1059 | Endif |
|---|
| 1060 | |
|---|
| 1061 | *---------------------------------------- |
|---|
| 1062 | * Determine lowest level of reported data |
|---|
| 1063 | *---------------------------------------- |
|---|
| 1064 | |
|---|
| 1065 | If (DrawTSnd = 1) |
|---|
| 1066 | field=sndtemp |
|---|
| 1067 | Else |
|---|
| 1068 | field=sndspd |
|---|
| 1069 | Endif |
|---|
| 1070 | |
|---|
| 1071 | "set gxout stat" |
|---|
| 1072 | "set x "_xval |
|---|
| 1073 | "set y "_yval |
|---|
| 1074 | "set t "_tval |
|---|
| 1075 | "set lev " _pmax " " _pmin |
|---|
| 1076 | "d maskout(lev,"field"+300)" |
|---|
| 1077 | rec=sublin(result,1) |
|---|
| 1078 | check=substr(rec,1,6) |
|---|
| 1079 | If (check = "Notice") |
|---|
| 1080 | rec=sublin(result,9) |
|---|
| 1081 | Else |
|---|
| 1082 | rec=sublin(result,8) |
|---|
| 1083 | Endif |
|---|
| 1084 | SfcPlev=subwrd(rec,5) |
|---|
| 1085 | |
|---|
| 1086 | If (DrawTSnd = 1 & DrawDSnd = 1) |
|---|
| 1087 | "set lev "SfcPlev |
|---|
| 1088 | "d "sndtemp |
|---|
| 1089 | rec=sublin(result,1) |
|---|
| 1090 | check=substr(rec,1,6) |
|---|
| 1091 | If (check = "Notice") |
|---|
| 1092 | rec=sublin(result,9) |
|---|
| 1093 | Else |
|---|
| 1094 | rec=sublin(result,8) |
|---|
| 1095 | Endif |
|---|
| 1096 | Sfctemp=subwrd(rec,4) |
|---|
| 1097 | "d "snddewp |
|---|
| 1098 | rec=sublin(result,1) |
|---|
| 1099 | check=substr(rec,1,6) |
|---|
| 1100 | If (check = "Notice") |
|---|
| 1101 | rec=sublin(result,9) |
|---|
| 1102 | Else |
|---|
| 1103 | rec=sublin(result,8) |
|---|
| 1104 | Endif |
|---|
| 1105 | Sfcdewp=subwrd(rec,4) |
|---|
| 1106 | SfcThee=Thetae(Sfctemp,Sfcdewp,SfcPlev) |
|---|
| 1107 | |
|---|
| 1108 | *------------------------------------------ |
|---|
| 1109 | * Calculate temperature and pressure of LCL |
|---|
| 1110 | *------------------------------------------ |
|---|
| 1111 | |
|---|
| 1112 | TLcl=Templcl(Sfctemp,Sfcdewp) |
|---|
| 1113 | PLcl=Preslcl(Sfctemp,Sfcdewp,SfcPlev) |
|---|
| 1114 | Endif |
|---|
| 1115 | |
|---|
| 1116 | *---------------------------------------------------------- |
|---|
| 1117 | * Plot parcel path from surface to LCL and up moist adiabat |
|---|
| 1118 | *---------------------------------------------------------- |
|---|
| 1119 | |
|---|
| 1120 | If (DrawPrcl = 1) |
|---|
| 1121 | say "Drawing parcel path from surface upward." |
|---|
| 1122 | If (PageType = "Portrait") |
|---|
| 1123 | xloc=7.15 |
|---|
| 1124 | Else |
|---|
| 1125 | xloc=6.65 |
|---|
| 1126 | Endif |
|---|
| 1127 | "q w2xy 1 "PLcl |
|---|
| 1128 | rec=sublin(result,1) |
|---|
| 1129 | yloc=subwrd(rec,6) |
|---|
| 1130 | "set strsiz 0.1" |
|---|
| 1131 | If (PLcl < _pmax) |
|---|
| 1132 | "set string 1 l" |
|---|
| 1133 | "draw string "xloc" "yloc" LCL" |
|---|
| 1134 | "set line 1 1 1" |
|---|
| 1135 | "draw line "xloc-0.15" "yloc" "xloc-0.05" "yloc |
|---|
| 1136 | Endif |
|---|
| 1137 | "set lon 0 100" |
|---|
| 1138 | "set gxout stat" |
|---|
| 1139 | "set line "PrclCol" "PrclLine " " PrclThk |
|---|
| 1140 | PTemp=LiftDry(Sfctemp,SfcPlev,PLcl,1,_pmin,_pmax) |
|---|
| 1141 | Ptemp=LiftWet(TLcl,PLcl,_pmin,1,_pmin,_pmax) |
|---|
| 1142 | Endif |
|---|
| 1143 | |
|---|
| 1144 | *------------------------------------------------------- |
|---|
| 1145 | * Determine level within lowest 250mb having highest |
|---|
| 1146 | * theta-e value |
|---|
| 1147 | *------------------------------------------------------- |
|---|
| 1148 | |
|---|
| 1149 | If (DrawTSnd = 1 & DrawDSnd = 1) |
|---|
| 1150 | "set x "_xval |
|---|
| 1151 | "set y "_yval |
|---|
| 1152 | "set t "_tval |
|---|
| 1153 | zz=1 |
|---|
| 1154 | MaxThee=-999 |
|---|
| 1155 | "set gxout stat" |
|---|
| 1156 | while (zz <= _zmax & pp > _pmax-250) |
|---|
| 1157 | "set z "zz |
|---|
| 1158 | pp=subwrd(result,4) |
|---|
| 1159 | "d "sndtemp |
|---|
| 1160 | rec=sublin(result,1) |
|---|
| 1161 | check=substr(rec,1,6) |
|---|
| 1162 | If (check = "Notice") |
|---|
| 1163 | rec=sublin(result,9) |
|---|
| 1164 | Else |
|---|
| 1165 | rec=sublin(result,8) |
|---|
| 1166 | Endif |
|---|
| 1167 | tt=subwrd(rec,4) |
|---|
| 1168 | "d "snddewp |
|---|
| 1169 | rec=sublin(result,1) |
|---|
| 1170 | check=substr(rec,1,6) |
|---|
| 1171 | If (check = "Notice") |
|---|
| 1172 | rec=sublin(result,9) |
|---|
| 1173 | Else |
|---|
| 1174 | rec=sublin(result,8) |
|---|
| 1175 | Endif |
|---|
| 1176 | dd=subwrd(rec,4) |
|---|
| 1177 | If (abs(tt) < 130 & abs(dd) < 130) |
|---|
| 1178 | Thee=Thetae(tt,dd,pp) |
|---|
| 1179 | If (Thee > MaxThee) |
|---|
| 1180 | MaxThee=Thee |
|---|
| 1181 | TMaxThee=tt |
|---|
| 1182 | DMaxThee=dd |
|---|
| 1183 | PMaxThee=pp |
|---|
| 1184 | Endif |
|---|
| 1185 | endif |
|---|
| 1186 | zz=zz+1 |
|---|
| 1187 | Endwhile |
|---|
| 1188 | If (PMaxThee = SfcPlev-250) |
|---|
| 1189 | PMaxThee = SfcPlev |
|---|
| 1190 | Endif |
|---|
| 1191 | *------------------------------------------------------ |
|---|
| 1192 | * Calculate temperature and pressure of LCL from highest |
|---|
| 1193 | * theta-e level |
|---|
| 1194 | *------------------------------------------------------ |
|---|
| 1195 | If (SfcPlev != PMaxThee) |
|---|
| 1196 | TLclMax=Templcl(TMaxThee,DMaxThee) |
|---|
| 1197 | PLclMax=Preslcl(TMaxThee,DMaxThee,PMaxThee) |
|---|
| 1198 | Endif |
|---|
| 1199 | Endif |
|---|
| 1200 | |
|---|
| 1201 | *---------------------------------------------------------- |
|---|
| 1202 | * Plot parcel path from highest theta-e level to LCL and up |
|---|
| 1203 | * moist adiabat |
|---|
| 1204 | *---------------------------------------------------------- |
|---|
| 1205 | |
|---|
| 1206 | If (DrawPMax = 1 & SfcPlev != PMaxThee) |
|---|
| 1207 | say "Drawing parcel path from most unstable level upward." |
|---|
| 1208 | If (PageType = "Portrait") |
|---|
| 1209 | xloc=7.15 |
|---|
| 1210 | Else |
|---|
| 1211 | xloc=6.65 |
|---|
| 1212 | Endif |
|---|
| 1213 | "q w2xy 1 "PLclMax |
|---|
| 1214 | rec=sublin(result,1) |
|---|
| 1215 | yloc=subwrd(rec,6) |
|---|
| 1216 | "set strsiz 0.1" |
|---|
| 1217 | If (PLclMax < _pmax) |
|---|
| 1218 | "set string 1 l" |
|---|
| 1219 | "draw string "xloc" "yloc" LCL" |
|---|
| 1220 | "set line 1 1 1" |
|---|
| 1221 | "draw line "xloc-0.15" "yloc" "xloc-0.05" "yloc |
|---|
| 1222 | Endif |
|---|
| 1223 | "set lon 0 100" |
|---|
| 1224 | "set gxout stat" |
|---|
| 1225 | "set line "PrclCol" "PrclLine " " PrclThk |
|---|
| 1226 | PTemp=LiftDry(TMaxThee,PMaxThee,PLclMax,1,_pmin,_pmax) |
|---|
| 1227 | Ptemp=LiftWet(TLclMax,PLclMax,_pmin,1,_pmin,_pmax) |
|---|
| 1228 | Endif |
|---|
| 1229 | |
|---|
| 1230 | *-------------------------------- |
|---|
| 1231 | * Draw thermodynamic indices |
|---|
| 1232 | *-------------------------------- |
|---|
| 1233 | |
|---|
| 1234 | If (DrawIndx = 1) |
|---|
| 1235 | "set string 1 l" |
|---|
| 1236 | "set strsiz 0.10" |
|---|
| 1237 | "set x "_xval |
|---|
| 1238 | "set y "_yval |
|---|
| 1239 | "set t "_tval |
|---|
| 1240 | say "Calculating precipitable water." |
|---|
| 1241 | pw=precipw(sndtemp,snddewp,_pmax,_pmin) |
|---|
| 1242 | say "Calculating thermodynamic indices." |
|---|
| 1243 | Temp850=interp(sndtemp,850) |
|---|
| 1244 | Temp700=interp(sndtemp,700) |
|---|
| 1245 | Temp500=interp(sndtemp,500) |
|---|
| 1246 | Dewp850=interp(snddewp,850) |
|---|
| 1247 | Dewp700=interp(snddewp,700) |
|---|
| 1248 | Dewp500=interp(snddewp,500) |
|---|
| 1249 | If (Temp850>-900 & Dewp850>-900 & Dewp700>-900 & Temp700>-900 & Temp500>-900) |
|---|
| 1250 | K=Temp850+Dewp850+Dewp700-Temp700-Temp500 |
|---|
| 1251 | Else |
|---|
| 1252 | K=-999 |
|---|
| 1253 | Endif |
|---|
| 1254 | If (Temp850 > -900 & Dewp850 > -900 & Temp500 > -900) |
|---|
| 1255 | tt=Temp850+Dewp850-2*Temp500 |
|---|
| 1256 | Else |
|---|
| 1257 | tt=-999 |
|---|
| 1258 | Endif |
|---|
| 1259 | Temp500V=virtual2(Temp500+273.15,Dewp500+273.15,500)-273.15 |
|---|
| 1260 | PclTemp=LiftWet(TLcl,PLcl,500,0) |
|---|
| 1261 | PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,500)-273.15 |
|---|
| 1262 | SLI=Temp500V-PclTempV |
|---|
| 1263 | rec=CAPE(TLcl,PLcl,100,sndtemp,snddewp) |
|---|
| 1264 | Pos=subwrd(rec,1) |
|---|
| 1265 | CIN=subwrd(rec,2) |
|---|
| 1266 | |
|---|
| 1267 | If (SfcPlev != PMaxThee) |
|---|
| 1268 | PclTemp=LiftWet(TLclMax,PLclMax,500,0) |
|---|
| 1269 | PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,500)-273.15 |
|---|
| 1270 | LIMax=Temp500V-PclTempV |
|---|
| 1271 | rec=CAPE(TLclMax,PLclMax,100,sndtemp,snddewp) |
|---|
| 1272 | PosMax=subwrd(rec,1) |
|---|
| 1273 | CINMax=subwrd(rec,2) |
|---|
| 1274 | Else |
|---|
| 1275 | LIMax=SLI |
|---|
| 1276 | PosMax=Pos |
|---|
| 1277 | CINMax=CIN |
|---|
| 1278 | MaxThee=SfcThee |
|---|
| 1279 | Endif |
|---|
| 1280 | |
|---|
| 1281 | If (PageType = "Portrait") |
|---|
| 1282 | If (Text1XC = -1) |
|---|
| 1283 | Text1XC=rxloc-0.75 |
|---|
| 1284 | Endif |
|---|
| 1285 | If (Text1YC = -1) |
|---|
| 1286 | Text1YC=tyloc-2.25 |
|---|
| 1287 | Endif |
|---|
| 1288 | If (Text2XC = -1) |
|---|
| 1289 | Text2XC=rxloc-0.75 |
|---|
| 1290 | Endif |
|---|
| 1291 | If (Text2YC = -1) |
|---|
| 1292 | Text2YC=tyloc-3.25 |
|---|
| 1293 | Endif |
|---|
| 1294 | If (Text3XC = -1) |
|---|
| 1295 | Text3XC=rxloc-0.75 |
|---|
| 1296 | Endif |
|---|
| 1297 | If (Text3YC = -1) |
|---|
| 1298 | Text3YC=tyloc-4.40 |
|---|
| 1299 | Endif |
|---|
| 1300 | Else |
|---|
| 1301 | If (Text1XC = -1) |
|---|
| 1302 | Text1XC=rxloc+2.50 |
|---|
| 1303 | Endif |
|---|
| 1304 | If (Text1YC = -1) |
|---|
| 1305 | Text1YC=tyloc-3.00 |
|---|
| 1306 | Endif |
|---|
| 1307 | If (Text2XC = -1) |
|---|
| 1308 | Text2XC=rxloc+2.50 |
|---|
| 1309 | Endif |
|---|
| 1310 | If (Text2YC = -1) |
|---|
| 1311 | Text2YC=tyloc-4.00 |
|---|
| 1312 | Endif |
|---|
| 1313 | If (Text3XC = -1) |
|---|
| 1314 | Text3XC=rxloc+2.50 |
|---|
| 1315 | Endif |
|---|
| 1316 | If (Text3YC = -1) |
|---|
| 1317 | Text3YC=tyloc-5.10 |
|---|
| 1318 | Endif |
|---|
| 1319 | Endif |
|---|
| 1320 | "set string 1 l 3" |
|---|
| 1321 | "set line 0 1 3" |
|---|
| 1322 | "draw recf "Text1XC-0.75 " " Text1YC-0.40 " " Text1XC+0.75 " " Text1YC+0.25 |
|---|
| 1323 | "set line 1 1 3" |
|---|
| 1324 | "draw rec "Text1XC-0.75 " " Text1YC-0.40 " " Text1XC+0.75 " " Text1YC+0.25 |
|---|
| 1325 | "draw string "Text1XC-0.70 " " Text1YC+0.10" K" |
|---|
| 1326 | "draw string "Text1XC+0.25 " " Text1YC+0.10" " int(K) |
|---|
| 1327 | "draw string "Text1XC-0.70 " " Text1YC-0.10 " TT" |
|---|
| 1328 | "draw string "Text1XC+0.25 " " Text1YC-0.10 " " int(tt) |
|---|
| 1329 | "draw string "Text1XC-0.70 " " Text1YC-0.25 " PW(cm)" |
|---|
| 1330 | "draw string "Text1XC+0.25 " " Text1YC-0.25 " " int(pw*100)/100 |
|---|
| 1331 | "set line 0 1 3" |
|---|
| 1332 | "draw recf "Text2XC-0.75 " " Text2YC-0.60 " " Text2XC+0.75 " " Text2YC+0.60 |
|---|
| 1333 | "set line 1 1 3" |
|---|
| 1334 | "draw rec "Text2XC-0.75 " " Text2YC-0.60 " " Text2XC+0.75 " " Text2YC+0.60 |
|---|
| 1335 | "draw string "Text2XC-0.35 " " Text2YC+0.50 " Surface" |
|---|
| 1336 | "draw string "Text2XC-0.70 " " Text2YC+0.30 " Temp(`3.`0C)" |
|---|
| 1337 | "draw string "Text2XC+0.25 " " Text2YC+0.30 " " int(Sfctemp*10)/10 |
|---|
| 1338 | "draw string "Text2XC-0.70 " " Text2YC+0.15 " Dewp(`3.`0C)" |
|---|
| 1339 | "draw string "Text2XC+0.25 " " Text2YC+0.15 " " int(Sfcdewp*10)/10 |
|---|
| 1340 | "draw string "Text2XC-0.70 " " Text2YC " `3z`0`bE`n(K)" |
|---|
| 1341 | "draw string "Text2XC+0.25 " " Text2YC " " int(SfcThee) |
|---|
| 1342 | "draw string "Text2XC-0.70 " " Text2YC-0.15 " LI" |
|---|
| 1343 | "draw string "Text2XC+0.25 " " Text2YC-0.15 " " round(SLI) |
|---|
| 1344 | "draw string "Text2XC-0.70 " " Text2YC-0.30 " CAPE(J)" |
|---|
| 1345 | "draw string "Text2XC+0.25 " " Text2YC-0.30 " " int(Pos) |
|---|
| 1346 | "draw string "Text2XC-0.70 " " Text2YC-0.45 " CIN(J)" |
|---|
| 1347 | "draw string "Text2XC+0.25 " " Text2YC-0.45 " " int(CIN) |
|---|
| 1348 | "set line 0 1 3" |
|---|
| 1349 | "draw recf "Text3XC-0.75 " " Text3YC-0.55 " " Text3XC+0.75 " " Text3YC+0.55 |
|---|
| 1350 | "set line 1 1 3" |
|---|
| 1351 | "draw rec "Text3XC-0.75 " " Text3YC-0.55 " " Text3XC+0.75 " " Text3YC+0.55 |
|---|
| 1352 | "draw string "Text3XC-0.60 " " Text3YC+0.45 " Most Unstable" |
|---|
| 1353 | "draw string "Text3XC-0.70 " " Text3YC+0.20 " Press(mb)" |
|---|
| 1354 | "draw string "Text3XC+0.25 " " Text3YC+0.20 " " int(PMaxThee) |
|---|
| 1355 | "draw string "Text3XC-0.70 " " Text3YC+0.05 " `3z`0`bE`n(K)" |
|---|
| 1356 | "draw string "Text3XC+0.25 " " Text3YC+0.05 " " int(MaxThee) |
|---|
| 1357 | "draw string "Text3XC-0.70 " " Text3YC-0.10 " LI" |
|---|
| 1358 | "draw string "Text3XC+0.25 " " Text3YC-0.10 " "round(LIMax) |
|---|
| 1359 | "draw string "Text3XC-0.70 " " Text3YC-0.25 " CAPE(J)" |
|---|
| 1360 | "draw string "Text3XC+0.25 " " Text3YC-0.25 " "int(PosMax) |
|---|
| 1361 | "draw string "Text3XC-0.70 " " Text3YC-0.40 " CIN(J)" |
|---|
| 1362 | "draw string "Text3XC+0.25 " " Text3YC-0.40 " " int(CINMax) |
|---|
| 1363 | Endif |
|---|
| 1364 | |
|---|
| 1365 | *----------------------------- |
|---|
| 1366 | * Draw wind profile along side |
|---|
| 1367 | *----------------------------- |
|---|
| 1368 | |
|---|
| 1369 | If (DrawBarb = 1) |
|---|
| 1370 | say "Drawing Wind Profile." |
|---|
| 1371 | If (poleloc = -1) |
|---|
| 1372 | If (PageType = "Portrait") |
|---|
| 1373 | poleloc = 8.0 |
|---|
| 1374 | Else |
|---|
| 1375 | poleloc = 7.5 |
|---|
| 1376 | Endif |
|---|
| 1377 | Endif |
|---|
| 1378 | If (barbline = 1) |
|---|
| 1379 | "set line 1 1 3" |
|---|
| 1380 | "draw line "poleloc " " byloc " " poleloc " " tyloc |
|---|
| 1381 | Endif |
|---|
| 1382 | If (BarbCol = -1) |
|---|
| 1383 | 'set rgb 41 255 0 132' |
|---|
| 1384 | 'set rgb 42 255 0 168' |
|---|
| 1385 | 'set rgb 43 255 0 204' |
|---|
| 1386 | 'set rgb 44 255 0 240' |
|---|
| 1387 | 'set rgb 45 255 0 255' |
|---|
| 1388 | 'set rgb 46 204 0 255' |
|---|
| 1389 | 'set rgb 47 174 0 255' |
|---|
| 1390 | 'set rgb 48 138 0 255' |
|---|
| 1391 | 'set rgb 49 108 0 255' |
|---|
| 1392 | 'set rgb 50 84 0 255' |
|---|
| 1393 | 'set rgb 51 40 0 255' |
|---|
| 1394 | 'set rgb 52 0 0 255' |
|---|
| 1395 | 'set rgb 53 0 42 255' |
|---|
| 1396 | 'set rgb 54 0 84 255' |
|---|
| 1397 | 'set rgb 55 0 120 255' |
|---|
| 1398 | 'set rgb 56 0 150 255' |
|---|
| 1399 | 'set rgb 57 0 192 255' |
|---|
| 1400 | 'set rgb 58 0 240 255' |
|---|
| 1401 | 'set rgb 59 0 255 210' |
|---|
| 1402 | 'set rgb 60 0 255 160' |
|---|
| 1403 | 'set rgb 61 0 255 126' |
|---|
| 1404 | 'set rgb 62 0 255 78' |
|---|
| 1405 | 'set rgb 63 84 255 0' |
|---|
| 1406 | 'set rgb 64 138 255 0' |
|---|
| 1407 | 'set rgb 65 188 255 0' |
|---|
| 1408 | 'set rgb 66 236 255 0' |
|---|
| 1409 | 'set rgb 67 255 255 0' |
|---|
| 1410 | 'set rgb 68 255 222 0' |
|---|
| 1411 | 'set rgb 69 255 192 0' |
|---|
| 1412 | 'set rgb 70 255 162 0' |
|---|
| 1413 | 'set rgb 71 255 138 0' |
|---|
| 1414 | 'set rgb 72 255 108 0' |
|---|
| 1415 | 'set rgb 73 255 84 0' |
|---|
| 1416 | 'set rgb 74 255 54 0' |
|---|
| 1417 | 'set rgb 75 255 12 0' |
|---|
| 1418 | 'set rgb 76 255 0 34' |
|---|
| 1419 | 'set rgb 77 255 0 70' |
|---|
| 1420 | 'set rgb 78 255 0 105' |
|---|
| 1421 | 'set rgb 79 255 0 140' |
|---|
| 1422 | 'set rgb 80 255 0 175' |
|---|
| 1423 | 'set rgb 81 255 0 215' |
|---|
| 1424 | 'set rgb 82 255 0 255' |
|---|
| 1425 | 'set rgb 83 255 255 255' |
|---|
| 1426 | |
|---|
| 1427 | col1='83 83 83 83 83 83 83 83 83 83 82 81 80 79 78' |
|---|
| 1428 | col2='77 76 75 74 73 72 71 70 69 68 67 66 65 64 63' |
|---|
| 1429 | col3='62 61 60 59 58 57 56 55 54 53 52 51 50 49 48' |
|---|
| 1430 | 'set rbcols 'col1' 'col2' 'col3 |
|---|
| 1431 | Endif |
|---|
| 1432 | "set z "_zmin" "_zmax |
|---|
| 1433 | "set gxout stat" |
|---|
| 1434 | zz=1 |
|---|
| 1435 | wspd=-999 |
|---|
| 1436 | cont=1 |
|---|
| 1437 | While (cont = 1 & zz < _zmax) |
|---|
| 1438 | "set z "zz |
|---|
| 1439 | pres=subwrd(result,4) |
|---|
| 1440 | "d "sndspd |
|---|
| 1441 | rec=sublin(result,1) |
|---|
| 1442 | check=substr(rec,1,6) |
|---|
| 1443 | If (check = "Notice") |
|---|
| 1444 | rec=sublin(result,9) |
|---|
| 1445 | Else |
|---|
| 1446 | rec=sublin(result,8) |
|---|
| 1447 | Endif |
|---|
| 1448 | wspd=subwrd(rec,4) |
|---|
| 1449 | if (wspd < 0 | pres > _pmax) |
|---|
| 1450 | zz=zz+1 |
|---|
| 1451 | else |
|---|
| 1452 | cont=0 |
|---|
| 1453 | Endif |
|---|
| 1454 | Endwhile |
|---|
| 1455 | While (zz <= _zmax) |
|---|
| 1456 | "d "sndspd"(z="zz")" |
|---|
| 1457 | rec=sublin(result,1) |
|---|
| 1458 | check=substr(rec,1,6) |
|---|
| 1459 | If (check = "Notice") |
|---|
| 1460 | rec=sublin(result,9) |
|---|
| 1461 | Else |
|---|
| 1462 | rec=sublin(result,8) |
|---|
| 1463 | Endif |
|---|
| 1464 | wspd=subwrd(rec,4) |
|---|
| 1465 | If (BarbCol >= 0) |
|---|
| 1466 | "set line "BarbCol " 1 "BarbThk |
|---|
| 1467 | Else |
|---|
| 1468 | tempbcol=55+wspd/5 |
|---|
| 1469 | If (tempbcol > 83) |
|---|
| 1470 | tempbcol = 83 |
|---|
| 1471 | Endif |
|---|
| 1472 | "set line "tempbcol " 1 "BarbThk |
|---|
| 1473 | Endif |
|---|
| 1474 | "d "snddir"(z="zz")" |
|---|
| 1475 | rec=sublin(result,1) |
|---|
| 1476 | check=substr(rec,1,6) |
|---|
| 1477 | If (check = "Notice") |
|---|
| 1478 | rec=sublin(result,9) |
|---|
| 1479 | Else |
|---|
| 1480 | rec=sublin(result,8) |
|---|
| 1481 | Endif |
|---|
| 1482 | wdir=subwrd(rec,4) |
|---|
| 1483 | xwind=GetUWnd(wspd,wdir) |
|---|
| 1484 | ywind=GetVWnd(wspd,wdir) |
|---|
| 1485 | "query gr2xy 5 "zz |
|---|
| 1486 | y1=subwrd(result,6) |
|---|
| 1487 | if (wspd > 0) |
|---|
| 1488 | cc=polelen/wspd |
|---|
| 1489 | xendpole=poleloc-xwind*cc |
|---|
| 1490 | yendpole=y1-ywind*cc |
|---|
| 1491 | endif |
|---|
| 1492 | if (xendpole>0 & wspd >= 0.5) |
|---|
| 1493 | if (flagbase = 1) |
|---|
| 1494 | "draw mark 3 "poleloc " " y1 " 0.05" |
|---|
| 1495 | endif |
|---|
| 1496 | "draw line " poleloc " " y1 " " xendpole " " yendpole |
|---|
| 1497 | flagloop=wspd/10 |
|---|
| 1498 | windcount=wspd |
|---|
| 1499 | flagcount=0 |
|---|
| 1500 | xflagstart=xendpole |
|---|
| 1501 | yflagstart=yendpole |
|---|
| 1502 | dx=cos((180-wdir)*_dtr) |
|---|
| 1503 | dy=sin((180-wdir)*_dtr) |
|---|
| 1504 | while (windcount > 47.5) |
|---|
| 1505 | flagcount=flagcount+1 |
|---|
| 1506 | dxflag=-len50*dx |
|---|
| 1507 | dyflag=-len50*dy |
|---|
| 1508 | xflagend=xflagstart+dxflag |
|---|
| 1509 | yflagend=yflagstart+dyflag |
|---|
| 1510 | windcount=windcount-50 |
|---|
| 1511 | x1=xflagstart+0.5*wid50*xwind/wspd |
|---|
| 1512 | y1=yflagstart+0.5*wid50*ywind/wspd |
|---|
| 1513 | x2=xflagstart-0.5*wid50*xwind/wspd |
|---|
| 1514 | y2=yflagstart-0.5*wid50*ywind/wspd |
|---|
| 1515 | If (Fill50 = 1) |
|---|
| 1516 | "draw polyf "x1" "y1" "x2" "y2" "xflagend" "yflagend" "x1" "y1 |
|---|
| 1517 | Else |
|---|
| 1518 | "draw line "x1 " "y1 " " xflagend " " yflagend " " |
|---|
| 1519 | "draw line "x2 " "y2 " " xflagend " " yflagend |
|---|
| 1520 | "draw line "x1 " "y1 " " x2 " " y2 |
|---|
| 1521 | Endif |
|---|
| 1522 | xflagstart=xflagstart+spac50*xwind/wspd |
|---|
| 1523 | yflagstart=yflagstart+spac50*ywind/wspd |
|---|
| 1524 | endwhile |
|---|
| 1525 | while (windcount > 7.5 ) |
|---|
| 1526 | flagcount=flagcount+1 |
|---|
| 1527 | dxflag=-len10*dx |
|---|
| 1528 | dyflag=-len10*dy |
|---|
| 1529 | xflagend=xflagstart+dxflag |
|---|
| 1530 | yflagend=yflagstart+dyflag |
|---|
| 1531 | windcount=windcount-10 |
|---|
| 1532 | "draw line " xflagstart " " yflagstart " " xflagend " " yflagend |
|---|
| 1533 | xflagstart=xflagstart+spac10*xwind/wspd |
|---|
| 1534 | yflagstart=yflagstart+spac10*ywind/wspd |
|---|
| 1535 | endwhile |
|---|
| 1536 | if (windcount > 2.5) |
|---|
| 1537 | flagcount=flagcount+1 |
|---|
| 1538 | if (flagcount = 1) |
|---|
| 1539 | xflagstart=xflagstart+spac05*xwind/wspd |
|---|
| 1540 | yflagstart=yflagstart+spac05*ywind/wspd |
|---|
| 1541 | endif |
|---|
| 1542 | dxflag=-len05*dx |
|---|
| 1543 | dyflag=-len05*dy |
|---|
| 1544 | xflagend=xflagstart+dxflag |
|---|
| 1545 | yflagend=yflagstart+dyflag |
|---|
| 1546 | windcount=windcount-5 |
|---|
| 1547 | "draw line " xflagstart " " yflagstart " " xflagend " " yflagend |
|---|
| 1548 | endif |
|---|
| 1549 | else |
|---|
| 1550 | if (wspd < 0.5 & wspd >= 0) |
|---|
| 1551 | "draw mark 2 " poleloc " " y1 " 0.08" |
|---|
| 1552 | endif |
|---|
| 1553 | endif |
|---|
| 1554 | zz=zz+barbint |
|---|
| 1555 | endwhile |
|---|
| 1556 | Endif |
|---|
| 1557 | |
|---|
| 1558 | *---------------- |
|---|
| 1559 | * Draw Hodograph |
|---|
| 1560 | *---------------- |
|---|
| 1561 | |
|---|
| 1562 | If (DrawHodo = 1) |
|---|
| 1563 | say "Drawing Hodograph." |
|---|
| 1564 | |
|---|
| 1565 | If (HodXcent = -1 | HodYcent = -1) |
|---|
| 1566 | If (PageType = "Portrait") |
|---|
| 1567 | HodXcent=6 |
|---|
| 1568 | HodYcent=9.5 |
|---|
| 1569 | Else |
|---|
| 1570 | HodXcent=9 |
|---|
| 1571 | HodYcent=7.0 |
|---|
| 1572 | Endif |
|---|
| 1573 | Endif |
|---|
| 1574 | HodL=HodXcent-HodSize/2.0 |
|---|
| 1575 | HodR=HodXcent+HodSize/2.0 |
|---|
| 1576 | HodT=HodYcent+HodSize/2.0 |
|---|
| 1577 | HodB=HodYcent-HodSize/2.0 |
|---|
| 1578 | RingSpac=HodSize/(NumRing*2) |
|---|
| 1579 | "set line 0" |
|---|
| 1580 | "draw recf "HodL" "HodB" "HodR" "HodT |
|---|
| 1581 | "set line "HodoCol" 1 6" |
|---|
| 1582 | "draw rec "HodL" "HodB" "HodR" "HodT |
|---|
| 1583 | "set line 1 1 3" |
|---|
| 1584 | "set string 1 c" |
|---|
| 1585 | "draw mark 1 "HodXcent " "HodYcent " " HodSize |
|---|
| 1586 | i=1 |
|---|
| 1587 | While (i <= NumRing) |
|---|
| 1588 | "set strsiz 0.10" |
|---|
| 1589 | "set string 1 c 3 45" |
|---|
| 1590 | uwnd=-i*HodRing*cos(45*_dtr) |
|---|
| 1591 | xloc=HodXcent+uwnd*RingSpac/HodRing |
|---|
| 1592 | yloc=HodYcent+uwnd*RingSpac/HodRing |
|---|
| 1593 | |
|---|
| 1594 | "draw mark 2 "HodXcent " " HodYcent " " i*HodSize/NumRing |
|---|
| 1595 | "draw string "xloc " " yloc " " HodRing*i |
|---|
| 1596 | i=i+1 |
|---|
| 1597 | Endwhile |
|---|
| 1598 | "set string 1 l 3 0" |
|---|
| 1599 | If (TickInt > 0) |
|---|
| 1600 | i=0 |
|---|
| 1601 | while (i < HodRing*NumRing) |
|---|
| 1602 | dist=i*RingSpac/HodRing |
|---|
| 1603 | hrxloc=HodXcent+dist |
|---|
| 1604 | hlxloc=HodXcent-dist |
|---|
| 1605 | htyloc=HodYcent+dist |
|---|
| 1606 | hbyloc=HodYcent-dist |
|---|
| 1607 | "set line 1 1 3" |
|---|
| 1608 | "draw line "hrxloc " " HodYcent-TickSize/2 " " hrxloc " " HodYcent+TickSize/2 |
|---|
| 1609 | "draw line "hlxloc " " HodYcent-TickSize/2 " " hlxloc " " HodYcent+TickSize/2 |
|---|
| 1610 | "draw line "HodXcent+TickSize/2 " " htyloc " " HodXcent-TickSize/2 " " htyloc |
|---|
| 1611 | "draw line "HodXcent+TickSize/2 " " hbyloc " " HodXcent-TickSize/2 " " hbyloc |
|---|
| 1612 | i=i+TickInt |
|---|
| 1613 | endwhile |
|---|
| 1614 | Endif |
|---|
| 1615 | "set line "HodoCol " " HodoLine " "HodoThk |
|---|
| 1616 | "draw string "HodL+0.05 " " HodT-0.1 " knots" |
|---|
| 1617 | zloop=_zmin |
|---|
| 1618 | xold=-999 |
|---|
| 1619 | yold=-999 |
|---|
| 1620 | count=0 |
|---|
| 1621 | Depth=0 |
|---|
| 1622 | While (zloop < _zmax & Depth < HodoDep) |
|---|
| 1623 | "set z "zloop |
|---|
| 1624 | pres=subwrd(result,4) |
|---|
| 1625 | "d "sndspd |
|---|
| 1626 | rec=sublin(result,1) |
|---|
| 1627 | check=substr(rec,1,6) |
|---|
| 1628 | If (check = "Notice") |
|---|
| 1629 | rec=sublin(result,9) |
|---|
| 1630 | Else |
|---|
| 1631 | rec=sublin(result,8) |
|---|
| 1632 | Endif |
|---|
| 1633 | wspd=subwrd(rec,4) |
|---|
| 1634 | "d "snddir |
|---|
| 1635 | rec=sublin(result,1) |
|---|
| 1636 | check=substr(rec,1,6) |
|---|
| 1637 | If (check = "Notice") |
|---|
| 1638 | rec=sublin(result,9) |
|---|
| 1639 | Else |
|---|
| 1640 | rec=sublin(result,8) |
|---|
| 1641 | Endif |
|---|
| 1642 | wdir=subwrd(rec,4) |
|---|
| 1643 | uwnd=GetUWnd(wspd,wdir) |
|---|
| 1644 | vwnd=GetVWnd(wspd,wdir) |
|---|
| 1645 | If (wspd >= 0) |
|---|
| 1646 | xloc=HodXcent+uwnd*RingSpac/HodRing |
|---|
| 1647 | yloc=HodYcent+vwnd*RingSpac/HodRing |
|---|
| 1648 | If (xloc > 0 & yloc > 0 & xold > 0 & yold > 0) |
|---|
| 1649 | Depth=Depth+pold-pres |
|---|
| 1650 | count=count+1 |
|---|
| 1651 | If (count = 1) |
|---|
| 1652 | "draw mark 3 "xold " " yold " 0.05" |
|---|
| 1653 | Endif |
|---|
| 1654 | "draw line "xold" "yold" "xloc" "yloc |
|---|
| 1655 | Endif |
|---|
| 1656 | xold=xloc |
|---|
| 1657 | yold=yloc |
|---|
| 1658 | Endif |
|---|
| 1659 | zloop=zloop+1 |
|---|
| 1660 | pold=pres |
|---|
| 1661 | EndWhile |
|---|
| 1662 | |
|---|
| 1663 | If (count > 0) |
|---|
| 1664 | "draw mark 3 "xold " " yold " 0.05" |
|---|
| 1665 | Endif |
|---|
| 1666 | Endif |
|---|
| 1667 | |
|---|
| 1668 | *---------------------------------------------- |
|---|
| 1669 | * Calculate and Display Absolute & S-R Helicity |
|---|
| 1670 | *---------------------------------------------- |
|---|
| 1671 | |
|---|
| 1672 | If (DrawHeli = 1) |
|---|
| 1673 | say "Calculating Helicity & SR Helicity." |
|---|
| 1674 | delp=10 |
|---|
| 1675 | UTotal=0 |
|---|
| 1676 | VTotal=0 |
|---|
| 1677 | |
|---|
| 1678 | * First, calculate mass-weighted mean wind |
|---|
| 1679 | * Since delp is a constant, and mass is proportional to |
|---|
| 1680 | * delp, this is a simple sum. |
|---|
| 1681 | |
|---|
| 1682 | "set lev "_pmax " " _pmin |
|---|
| 1683 | "define uwndarr="sndspd"*cos((270-"snddir")*"_dtr")" |
|---|
| 1684 | "define vwndarr="sndspd"*sin((270-"snddir")*"_dtr")" |
|---|
| 1685 | |
|---|
| 1686 | pres=MeanVBot |
|---|
| 1687 | While (pres >= MeanVTop) |
|---|
| 1688 | uwnd=interp(uwndarr,pres)*_ktm |
|---|
| 1689 | vwnd=interp(vwndarr,pres)*_ktm |
|---|
| 1690 | If (uwnd > -900 & vwnd > -900) |
|---|
| 1691 | UTotal=UTotal+uwnd |
|---|
| 1692 | VTotal=VTotal+vwnd |
|---|
| 1693 | Endif |
|---|
| 1694 | pres=pres-delp |
|---|
| 1695 | EndWhile |
|---|
| 1696 | vcount=1+(MeanVBot-MeanVTop)/delp |
|---|
| 1697 | Umean=UTotal/vcount |
|---|
| 1698 | Vmean=VTotal/vcount |
|---|
| 1699 | Spdmean=GetWSpd(Umean,Vmean) |
|---|
| 1700 | MeanDir=GetWDir(Umean,Vmean) |
|---|
| 1701 | |
|---|
| 1702 | * Now, rotate and reduce mean wind to get storm motion |
|---|
| 1703 | |
|---|
| 1704 | If (StormMot = 1) |
|---|
| 1705 | If (Spdmean < 15) |
|---|
| 1706 | Reduct=0.25 |
|---|
| 1707 | Rotate=30 |
|---|
| 1708 | Else |
|---|
| 1709 | Reduct=0.20 |
|---|
| 1710 | Rotate=20 |
|---|
| 1711 | Endif |
|---|
| 1712 | Else |
|---|
| 1713 | Reduct=0.0 |
|---|
| 1714 | Rotate=0.0 |
|---|
| 1715 | Endif |
|---|
| 1716 | |
|---|
| 1717 | UReduce=(1-Reduct)*Umean |
|---|
| 1718 | VReduce=(1-Reduct)*Vmean |
|---|
| 1719 | StormSpd=GetWSpd(UReduce,VReduce) |
|---|
| 1720 | |
|---|
| 1721 | StormDir=GetWDir(UReduce,VReduce)+Rotate |
|---|
| 1722 | If (StormDir >= 360) |
|---|
| 1723 | StormDir=StormDir-360 |
|---|
| 1724 | Endif |
|---|
| 1725 | |
|---|
| 1726 | StormU=GetUWnd(StormSpd,StormDir) |
|---|
| 1727 | StormV=GetVWnd(StormSpd,StormDir) |
|---|
| 1728 | |
|---|
| 1729 | * Draw Storm Motion Vector |
|---|
| 1730 | |
|---|
| 1731 | xloc=HodXcent+_mtk*StormU*RingSpac/HodRing |
|---|
| 1732 | yloc=HodYcent+_mtk*StormV*RingSpac/HodRing |
|---|
| 1733 | |
|---|
| 1734 | "set line 1 1 4" |
|---|
| 1735 | "draw line "HodXcent " " HodYcent " " xloc " " yloc |
|---|
| 1736 | Arr1U=GetUWnd(HodRing/10,StormDir+30) |
|---|
| 1737 | Arr1V=GetVWnd(HodRing/10,StormDir+30) |
|---|
| 1738 | Arr2U=GetUWnd(HodRing/10,StormDir-30) |
|---|
| 1739 | Arr2V=GetVWnd(HodRing/10,StormDir-30) |
|---|
| 1740 | |
|---|
| 1741 | xloc2=xloc-Arr1U/HodRing |
|---|
| 1742 | xloc3=xloc-Arr2U/HodRing |
|---|
| 1743 | yloc2=yloc-Arr1V/HodRing |
|---|
| 1744 | yloc3=yloc-Arr2V/HodRing |
|---|
| 1745 | |
|---|
| 1746 | "set line 1 1 3" |
|---|
| 1747 | |
|---|
| 1748 | If (FillArrw = 0) |
|---|
| 1749 | "draw line "xloc" "yloc" "xloc2" "yloc2 |
|---|
| 1750 | "draw line "xloc" "yloc" "xloc3" "yloc3 |
|---|
| 1751 | Else |
|---|
| 1752 | "draw polyf "xloc" "yloc" "xloc2" "yloc2" "xloc3" "yloc3" "xloc" "yloc |
|---|
| 1753 | Endif |
|---|
| 1754 | |
|---|
| 1755 | |
|---|
| 1756 | * Now, calculate SR and Environmental Helicity |
|---|
| 1757 | |
|---|
| 1758 | helic=0 |
|---|
| 1759 | SRhelic=0 |
|---|
| 1760 | MinP=SfcPlev-HelicDep |
|---|
| 1761 | pres=SfcPlev |
|---|
| 1762 | uwndold=-999 |
|---|
| 1763 | vwndold=-999 |
|---|
| 1764 | While (pres >= MinP) |
|---|
| 1765 | uwnd=interp(uwndarr,pres)*_ktm |
|---|
| 1766 | vwnd=interp(vwndarr,pres)*_ktm |
|---|
| 1767 | If (uwnd > -900 & uwndold > -900) |
|---|
| 1768 | du=uwnd-uwndold |
|---|
| 1769 | dv=vwnd-vwndold |
|---|
| 1770 | ubar=0.5*(uwnd+uwndold) |
|---|
| 1771 | vbar=0.5*(vwnd+vwndold) |
|---|
| 1772 | uhelic=-dv*ubar |
|---|
| 1773 | vhelic=du*vbar |
|---|
| 1774 | SRuhelic=-dv*(ubar-StormU) |
|---|
| 1775 | SRvhelic=du*(vbar-StormV) |
|---|
| 1776 | SRhelic=SRhelic+SRuhelic+SRvhelic |
|---|
| 1777 | helic=helic+uhelic+vhelic |
|---|
| 1778 | Endif |
|---|
| 1779 | uwndold=uwnd |
|---|
| 1780 | vwndold=vwnd |
|---|
| 1781 | pres=pres-delp |
|---|
| 1782 | EndWhile |
|---|
| 1783 | |
|---|
| 1784 | "set strsiz 0.1" |
|---|
| 1785 | "set string 1 l 3" |
|---|
| 1786 | If (PageType = "Portrait") |
|---|
| 1787 | If (Text4XC = -1) |
|---|
| 1788 | Text4XC=rxloc-0.75 |
|---|
| 1789 | Endif |
|---|
| 1790 | If (Text4YC = -1) |
|---|
| 1791 | Text4YC=tyloc-5.45 |
|---|
| 1792 | Endif |
|---|
| 1793 | Else |
|---|
| 1794 | If (Text4XC = -1) |
|---|
| 1795 | Text4XC=rxloc+2.50 |
|---|
| 1796 | Endif |
|---|
| 1797 | If (Text4YC = -1) |
|---|
| 1798 | Text4YC=tyloc-6.10 |
|---|
| 1799 | Endif |
|---|
| 1800 | Endif |
|---|
| 1801 | "set line 0 1 3" |
|---|
| 1802 | "draw recf "Text4XC-0.75 " "Text4YC-0.5 " " Text4XC+0.75 " " Text4YC+0.5 |
|---|
| 1803 | "set line 1 1 3" |
|---|
| 1804 | "draw rec "Text4XC-0.75 " "Text4YC-0.5 " " Text4XC+0.75 " " Text4YC+0.5 |
|---|
| 1805 | "draw string "Text4XC-0.45 " " Text4YC+0.40 " Hodograph" |
|---|
| 1806 | "draw string "Text4XC-0.70 " " Text4YC+0.20 " EH" |
|---|
| 1807 | "draw string "Text4XC+0.25 " " Text4YC+0.20 " "int(helic) |
|---|
| 1808 | "draw string "Text4XC-0.70 " " Text4YC+0.05 " SREH" |
|---|
| 1809 | "draw string "Text4XC+0.25 " " Text4YC+0.05 " " int(SRhelic) |
|---|
| 1810 | "draw string "Text4XC-0.70 " " Text4YC-0.20 " StmDir" |
|---|
| 1811 | "draw string "Text4XC+0.25 " " Text4YC-0.20 " " int(StormDir)"`3.`0" |
|---|
| 1812 | "draw string "Text4XC-0.70 " " Text4YC-0.35 " StmSpd(kt)" |
|---|
| 1813 | "draw string "Text4XC+0.25 " " Text4YC-0.35 " " int(_mtk*StormSpd) |
|---|
| 1814 | Endif |
|---|
| 1815 | |
|---|
| 1816 | "draw string "Text4XC-0.70 " " Text4YC-0.80 " Longitude" |
|---|
| 1817 | "draw string "Text4XC+0.25 " " Text4YC-0.80 " " _xlon |
|---|
| 1818 | "draw string "Text4XC-0.70 " " Text4YC-1.00 " Latitude" |
|---|
| 1819 | "draw string "Text4XC+0.25 " " Text4YC-1.00 " " _ylat |
|---|
| 1820 | if (OL = 1) |
|---|
| 1821 | "draw string "Text4XC-0.70 " " Text4YC-1.20 " Time" |
|---|
| 1822 | "draw string "Text4XC+0.25 " " Text4YC-1.20 " " _ttime |
|---|
| 1823 | endif |
|---|
| 1824 | if (OL = 2) |
|---|
| 1825 | "set string 3" |
|---|
| 1826 | "draw string "Text4XC+0.25 " " Text4YC-1.40 " " _ttime |
|---|
| 1827 | "set string 1" |
|---|
| 1828 | endif |
|---|
| 1829 | "set string 7" |
|---|
| 1830 | "draw string "Text4XC-0.99 " " Text4YC-1.70 " Initial Time" |
|---|
| 1831 | "draw string "Text4XC+0.25 " " Text4YC-1.70 " " _analysis |
|---|
| 1832 | "set string 1" |
|---|
| 1833 | |
|---|
| 1834 | *--------------------------------------------------- |
|---|
| 1835 | * Plot RH profile. |
|---|
| 1836 | *--------------------------------------------------- |
|---|
| 1837 | |
|---|
| 1838 | If (DrawRH = 1) |
|---|
| 1839 | "set z "_zmin" "_zmax |
|---|
| 1840 | "set x "_xval |
|---|
| 1841 | "set y "_yval |
|---|
| 1842 | "set t "_tval |
|---|
| 1843 | "set zlog on" |
|---|
| 1844 | "set gxout line" |
|---|
| 1845 | "set ccolor "RHCol |
|---|
| 1846 | "set cstyle "RHLine |
|---|
| 1847 | "set cmark "RHMrk |
|---|
| 1848 | "set digsiz "MrkSize |
|---|
| 1849 | "set missconn on" |
|---|
| 1850 | "set xlab on" |
|---|
| 1851 | "set frame off" |
|---|
| 1852 | "set vrange 0 350" |
|---|
| 1853 | "set xlpos 0 t" |
|---|
| 1854 | "set xlevs 25 50 75 100" |
|---|
| 1855 | "set grid vertical 5" |
|---|
| 1856 | "define rh=100*exp((17.2694*"snddewp")/("snddewp"+237.3)-(17.2694*"sndtemp")/("sndtemp"+237.3))" |
|---|
| 1857 | "d rh" |
|---|
| 1858 | If (LblAxes = 1) |
|---|
| 1859 | "set string 1 c 3 0" |
|---|
| 1860 | "set strsiz 0.125" |
|---|
| 1861 | If (PageType = "Portrait") |
|---|
| 1862 | "draw string 1.5 10.85 RH (%)" |
|---|
| 1863 | Else |
|---|
| 1864 | "draw string 1.75 8.35 RH (%)" |
|---|
| 1865 | Endif |
|---|
| 1866 | Endif |
|---|
| 1867 | Endif |
|---|
| 1868 | |
|---|
| 1869 | *------------------------------------------ |
|---|
| 1870 | * Reset environment to original dimensions |
|---|
| 1871 | *------------------------------------------ |
|---|
| 1872 | |
|---|
| 1873 | "set t "_tval |
|---|
| 1874 | "set x "_xval |
|---|
| 1875 | "set y "_yval |
|---|
| 1876 | "set z "_zmin " "_zmax |
|---|
| 1877 | |
|---|
| 1878 | say "Done." |
|---|
| 1879 | |
|---|
| 1880 | Return(0) |
|---|
| 1881 | |
|---|
| 1882 | ************************************************************************* |
|---|
| 1883 | function Templcl(temp,dewp) |
|---|
| 1884 | |
|---|
| 1885 | *------------------------------------------------------ |
|---|
| 1886 | * Calculate the temp at the LCL given temp & dewp in C |
|---|
| 1887 | *------------------------------------------------------ |
|---|
| 1888 | |
|---|
| 1889 | tempk=temp+273.15 |
|---|
| 1890 | dewpk=dewp+273.15 |
|---|
| 1891 | Parta=1/(dewpk-56) |
|---|
| 1892 | Partb=log(tempk/dewpk)/800 |
|---|
| 1893 | Tlcl=1/(Parta+Partb)+56 |
|---|
| 1894 | return(Tlcl-273.15) |
|---|
| 1895 | |
|---|
| 1896 | ************************************************************************** |
|---|
| 1897 | |
|---|
| 1898 | function Preslcl(temp,dewp,pres) |
|---|
| 1899 | |
|---|
| 1900 | *------------------------------------------------------- |
|---|
| 1901 | * Calculate press of LCL given temp & dewp in C and pressure |
|---|
| 1902 | *------------------------------------------------------- |
|---|
| 1903 | |
|---|
| 1904 | Tlclk=Templcl(temp,dewp)+273.15 |
|---|
| 1905 | tempk=temp+273.15 |
|---|
| 1906 | theta=tempk*pow(1000/pres,0.286) |
|---|
| 1907 | plcl=1000*pow(Tlclk/theta,3.48) |
|---|
| 1908 | return(plcl) |
|---|
| 1909 | |
|---|
| 1910 | ************************************************************************** |
|---|
| 1911 | function LiftWet(startt,startp,endp,display,Pmin,Pmax) |
|---|
| 1912 | |
|---|
| 1913 | *-------------------------------------------------------------------- |
|---|
| 1914 | * Lift a parcel moist adiabatically from startp to endp. |
|---|
| 1915 | * Init temp is startt in C. If you wish to see the parcel's |
|---|
| 1916 | * path plotted, display should be 1. Returns temp of parcel at endp. |
|---|
| 1917 | *-------------------------------------------------------------------- |
|---|
| 1918 | |
|---|
| 1919 | temp=startt |
|---|
| 1920 | pres=startp |
|---|
| 1921 | cont = 1 |
|---|
| 1922 | delp=10 |
|---|
| 1923 | While (pres >= endp & cont = 1) |
|---|
| 1924 | If (display = 1) |
|---|
| 1925 | xtemp=GetXLoc(temp,pres) |
|---|
| 1926 | "q w2xy "xtemp" "pres |
|---|
| 1927 | xloc=subwrd(result,3) |
|---|
| 1928 | yloc=subwrd(result,6) |
|---|
| 1929 | If (xtemp < 0 | xtemp > 100) |
|---|
| 1930 | cont=0 |
|---|
| 1931 | Else |
|---|
| 1932 | If (pres >= Pmin & pres < Pmax & pres < startp) |
|---|
| 1933 | "draw line "xold" "yold" "xloc" "yloc |
|---|
| 1934 | Endif |
|---|
| 1935 | Endif |
|---|
| 1936 | Endif |
|---|
| 1937 | xold=xloc |
|---|
| 1938 | yold=yloc |
|---|
| 1939 | temp=temp-100*delp*gammaw(temp,pres-delp/2,100) |
|---|
| 1940 | pres=pres-delp |
|---|
| 1941 | EndWhile |
|---|
| 1942 | return(temp) |
|---|
| 1943 | |
|---|
| 1944 | |
|---|
| 1945 | ************************************************************************** |
|---|
| 1946 | function LiftDry(startt,startp,endp,display,Pmin,Pmax) |
|---|
| 1947 | |
|---|
| 1948 | *-------------------------------------------------------------------- |
|---|
| 1949 | * Lift a parcel dry adiabatically from startp to endp. |
|---|
| 1950 | * Init temp is startt in C. If you wish to see the parcel's |
|---|
| 1951 | * path plotted, display should be 1. Returns temp of parcel at endp. |
|---|
| 1952 | *-------------------------------------------------------------------- |
|---|
| 1953 | |
|---|
| 1954 | starttk=startt+273.15 |
|---|
| 1955 | cont = 1 |
|---|
| 1956 | delp=10 |
|---|
| 1957 | round=int(startp/10)*10 |
|---|
| 1958 | subscr=0.1*round |
|---|
| 1959 | powstart=pow(startp,-0.286) |
|---|
| 1960 | temp=starttk*_powpres.subscr*powstart-273.15 |
|---|
| 1961 | pres=round-10 |
|---|
| 1962 | While (pres >= endp & cont = 1) |
|---|
| 1963 | subscr=0.1*pres |
|---|
| 1964 | temp=starttk*_powpres.subscr*powstart-273.15 |
|---|
| 1965 | If (display = 1) |
|---|
| 1966 | xtemp=GetXLoc(temp,pres) |
|---|
| 1967 | "q w2xy "xtemp" "pres |
|---|
| 1968 | xloc=subwrd(result,3) |
|---|
| 1969 | yloc=subwrd(result,6) |
|---|
| 1970 | If (xtempold > 0 & xtempold < 100 & xtemp > 0 & xtemp < 100) |
|---|
| 1971 | If (pres >= Pmin & pres < Pmax & pres < startp) |
|---|
| 1972 | "draw line "xold" "yold" "xloc" "yloc |
|---|
| 1973 | Endif |
|---|
| 1974 | Endif |
|---|
| 1975 | Endif |
|---|
| 1976 | xold=xloc |
|---|
| 1977 | xtempold=xtemp |
|---|
| 1978 | yold=yloc |
|---|
| 1979 | pres=pres-delp |
|---|
| 1980 | EndWhile |
|---|
| 1981 | return(temp) |
|---|
| 1982 | |
|---|
| 1983 | ************************************************************************** |
|---|
| 1984 | function CAPE(startt,startp,endp,sndtemp,snddewp) |
|---|
| 1985 | |
|---|
| 1986 | *--------------------------------------------------------------------- |
|---|
| 1987 | * Returns all postive area and convective inhibition above LCL. |
|---|
| 1988 | * Parcel is lifted from LCL at startt,startp and is halted |
|---|
| 1989 | * at endp. Integration method used is trapezoid method. |
|---|
| 1990 | *--------------------------------------------------------------------- |
|---|
| 1991 | |
|---|
| 1992 | pres=startp |
|---|
| 1993 | PclTemp=startt |
|---|
| 1994 | PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15 |
|---|
| 1995 | delp=10 |
|---|
| 1996 | Pos=0 |
|---|
| 1997 | Neg=0 |
|---|
| 1998 | Neg2=0 |
|---|
| 1999 | |
|---|
| 2000 | count=0 |
|---|
| 2001 | While (pres >= endp) |
|---|
| 2002 | EnvTemp=interp(sndtemp,pres) |
|---|
| 2003 | EnvDewp=interp(snddewp,pres) |
|---|
| 2004 | EnvTempV=virtual2(EnvTemp+273.15,EnvDewp+273.15,pres)-273.15 |
|---|
| 2005 | DelT=PclTempV-EnvTempV |
|---|
| 2006 | If (abs(EnvTempV) < 130 & abs(PclTempV) < 130) |
|---|
| 2007 | count=count+1 |
|---|
| 2008 | If (count > 1) |
|---|
| 2009 | Val=DelT/pres+Prev |
|---|
| 2010 | If (Val > 0) |
|---|
| 2011 | Pos=Pos+Val |
|---|
| 2012 | Neg2=0 |
|---|
| 2013 | Else |
|---|
| 2014 | Neg=Neg+abs(Val) |
|---|
| 2015 | Neg2=Neg2+abs(Val) |
|---|
| 2016 | Endif |
|---|
| 2017 | Endif |
|---|
| 2018 | Prev=DelT/pres |
|---|
| 2019 | Endif |
|---|
| 2020 | pres=pres-delp |
|---|
| 2021 | PclTemp=PclTemp-100*delp*gammaw(PclTemp,pres,100) |
|---|
| 2022 | PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15 |
|---|
| 2023 | Endwhile |
|---|
| 2024 | |
|---|
| 2025 | Pos=0.5*Pos*287*delp |
|---|
| 2026 | CIN=0.5*(Neg-Neg2)*287*delp |
|---|
| 2027 | |
|---|
| 2028 | return(Pos" "CIN) |
|---|
| 2029 | |
|---|
| 2030 | *************************************************************************** |
|---|
| 2031 | function gammaw(tempc,pres,rh) |
|---|
| 2032 | |
|---|
| 2033 | *----------------------------------------------------------------------- |
|---|
| 2034 | * Function to calculate the moist adiabatic lapse rate (deg C/Pa) based |
|---|
| 2035 | * on the temperature, pressure, and rh of the environment. |
|---|
| 2036 | *---------------------------------------------------------------------- |
|---|
| 2037 | |
|---|
| 2038 | tempk=tempc+273.15 |
|---|
| 2039 | es=satvap2(tempc) |
|---|
| 2040 | ws=mixratio(es,pres) |
|---|
| 2041 | w=rh*ws/100 |
|---|
| 2042 | tempv=virtual(tempk,w) |
|---|
| 2043 | latent=latentc(tempc) |
|---|
| 2044 | |
|---|
| 2045 | A=1.0+latent*ws/(287*tempk) |
|---|
| 2046 | B=1.0+0.622*latent*latent*ws/(1005*287*tempk*tempk) |
|---|
| 2047 | Density=100*pres/(287*tempv) |
|---|
| 2048 | lapse=(A/B)/(1005*Density) |
|---|
| 2049 | return(lapse) |
|---|
| 2050 | |
|---|
| 2051 | ************************************************************************* |
|---|
| 2052 | function latentc(tempc) |
|---|
| 2053 | |
|---|
| 2054 | *----------------------------------------------------------------------- |
|---|
| 2055 | * Function to return the latent heat of condensation in J/kg given |
|---|
| 2056 | * temperature in degrees Celsius. |
|---|
| 2057 | *----------------------------------------------------------------------- |
|---|
| 2058 | |
|---|
| 2059 | val=2502.2-2.43089*tempc |
|---|
| 2060 | |
|---|
| 2061 | return(val*1000) |
|---|
| 2062 | |
|---|
| 2063 | ************************************************************************* |
|---|
| 2064 | function precipw(sndtemp,snddewp,startp,endp) |
|---|
| 2065 | |
|---|
| 2066 | *----------------------------------------------------------------------- |
|---|
| 2067 | * Function to calculate the precipitable water (cm) in a sounding |
|---|
| 2068 | * starting at pressure level startp and ending at pressure level endp. |
|---|
| 2069 | *----------------------------------------------------------------------- |
|---|
| 2070 | |
|---|
| 2071 | ppold=-999 |
|---|
| 2072 | ttold=-999 |
|---|
| 2073 | ddold=-999 |
|---|
| 2074 | delp=10 |
|---|
| 2075 | Int=0 |
|---|
| 2076 | mix=0 |
|---|
| 2077 | pres=startp |
|---|
| 2078 | logpp=log(pres) |
|---|
| 2079 | logppm=log(pres-delp) |
|---|
| 2080 | while (pres >= endp) |
|---|
| 2081 | tt=interp(sndtemp,pres) |
|---|
| 2082 | dd=interp(snddewp,pres) |
|---|
| 2083 | if (tt>-900 & ttold>-900 & dd>-900 & ddold>-900) |
|---|
| 2084 | e=satvap2(dd) |
|---|
| 2085 | mix=mixratio(e,pres) |
|---|
| 2086 | mixavg=(logpp*mix+logppm*mixold)/(logpp+logppm) |
|---|
| 2087 | Int=Int+1.020408*mixavg*delp |
|---|
| 2088 | endif |
|---|
| 2089 | ttold=tt |
|---|
| 2090 | ddold=dd |
|---|
| 2091 | ppold=pp |
|---|
| 2092 | mixold=mix |
|---|
| 2093 | pres=pres-delp |
|---|
| 2094 | logpp=logppm |
|---|
| 2095 | logppm=log(pres-delp) |
|---|
| 2096 | endwhile |
|---|
| 2097 | |
|---|
| 2098 | return(Int) |
|---|
| 2099 | |
|---|
| 2100 | ************************************************************************* |
|---|
| 2101 | |
|---|
| 2102 | function virtual(temp,mix) |
|---|
| 2103 | |
|---|
| 2104 | *------------------------------------------------------------ |
|---|
| 2105 | * Function to return virtual temperature given temperature in |
|---|
| 2106 | * kelvin and mixing ratio in g/g. |
|---|
| 2107 | *------------------------------------------------------------- |
|---|
| 2108 | |
|---|
| 2109 | tempv=temp*(1.0+0.6*mix) |
|---|
| 2110 | |
|---|
| 2111 | return (tempv) |
|---|
| 2112 | |
|---|
| 2113 | ************************************************************************ |
|---|
| 2114 | |
|---|
| 2115 | function virtual2(temp,dewp,pres) |
|---|
| 2116 | |
|---|
| 2117 | *------------------------------------------------------------ |
|---|
| 2118 | * Function to return virtual temperature in kelvin given temperature in |
|---|
| 2119 | * kelvin and dewpoint in kelvin and pressure in mb |
|---|
| 2120 | *------------------------------------------------------------- |
|---|
| 2121 | |
|---|
| 2122 | if (temp > 0 & dewp > 0) |
|---|
| 2123 | vap=satvap2(dewp-273.15) |
|---|
| 2124 | mix=mixratio(vap,pres) |
|---|
| 2125 | tempv=virtual(temp,mix) |
|---|
| 2126 | else |
|---|
| 2127 | tempv=-9999 |
|---|
| 2128 | endif |
|---|
| 2129 | |
|---|
| 2130 | return (tempv) |
|---|
| 2131 | |
|---|
| 2132 | ************************************************************************ |
|---|
| 2133 | |
|---|
| 2134 | function satvapor(temp) |
|---|
| 2135 | |
|---|
| 2136 | *--------------------------------------------------------------- |
|---|
| 2137 | * Given temp in Celsius, returns saturation vapor pressure in mb |
|---|
| 2138 | *--------------------------------------------------------------- |
|---|
| 2139 | |
|---|
| 2140 | pol=_C0+temp*(_C1+temp*(_C2+temp*(_C3+temp*(_C4+temp*(_C5+temp*(_C6+temp*(_C7+temp*(_C8+temp*(_C9))))))))) |
|---|
| 2141 | |
|---|
| 2142 | return(6.1078/pow(pol,8)) |
|---|
| 2143 | |
|---|
| 2144 | ************************************************************************ |
|---|
| 2145 | |
|---|
| 2146 | function satvap2(temp) |
|---|
| 2147 | |
|---|
| 2148 | *--------------------------------------------------------------- |
|---|
| 2149 | * Given temp in Celsius, returns saturation vapor pressure in mb |
|---|
| 2150 | *--------------------------------------------------------------- |
|---|
| 2151 | |
|---|
| 2152 | es=6.112*exp(17.67*temp/(temp+243.5)) |
|---|
| 2153 | |
|---|
| 2154 | return(es) |
|---|
| 2155 | |
|---|
| 2156 | ************************************************************************* |
|---|
| 2157 | |
|---|
| 2158 | function mixratio(e,p) |
|---|
| 2159 | |
|---|
| 2160 | *------------------------------------------------------ |
|---|
| 2161 | * Given vapor pressure and pressure, return mixing ratio |
|---|
| 2162 | *------------------------------------------------------- |
|---|
| 2163 | |
|---|
| 2164 | mix=0.622*e/(p-e) |
|---|
| 2165 | |
|---|
| 2166 | return(mix) |
|---|
| 2167 | |
|---|
| 2168 | ************************************************************************ |
|---|
| 2169 | |
|---|
| 2170 | function getrh(temp,dewp,pres) |
|---|
| 2171 | |
|---|
| 2172 | tempk=temp+273.15 |
|---|
| 2173 | dewpk=dewp+273.15 |
|---|
| 2174 | |
|---|
| 2175 | es=satvap2(temp) |
|---|
| 2176 | |
|---|
| 2177 | If (temp > 0) |
|---|
| 2178 | A=2.53*pow(10,9) |
|---|
| 2179 | B=5420 |
|---|
| 2180 | Else |
|---|
| 2181 | A=3.41*pow(10,10) |
|---|
| 2182 | B=6130 |
|---|
| 2183 | Endif |
|---|
| 2184 | |
|---|
| 2185 | w=A*0.622*exp(-B/dewpk)/pres |
|---|
| 2186 | ws=mixratio(es,pres) |
|---|
| 2187 | |
|---|
| 2188 | return(100*w/ws) |
|---|
| 2189 | |
|---|
| 2190 | ************************************************************************ |
|---|
| 2191 | function interp(array,pres) |
|---|
| 2192 | |
|---|
| 2193 | *------------------------------------------------------------------------ |
|---|
| 2194 | * Interpolate inside array for pressure level pres. |
|---|
| 2195 | * Returns estimated value of array at pressure pres. |
|---|
| 2196 | *------------------------------------------------------------------------ |
|---|
| 2197 | |
|---|
| 2198 | "set gxout stat" |
|---|
| 2199 | "set lev "pres |
|---|
| 2200 | altpres=subwrd(result,4) |
|---|
| 2201 | "q dims" |
|---|
| 2202 | rec=sublin(result,4) |
|---|
| 2203 | zlev=subwrd(rec,9) |
|---|
| 2204 | |
|---|
| 2205 | If (zlev < 2 | zlev > _zmaxfile) |
|---|
| 2206 | Vest = -9999.0 |
|---|
| 2207 | Else |
|---|
| 2208 | If (altpres > pres) |
|---|
| 2209 | zlev=zlev+1 |
|---|
| 2210 | Endif |
|---|
| 2211 | "set z "zlev |
|---|
| 2212 | PAbove=subwrd(result,4) |
|---|
| 2213 | "d "array"(lev="PAbove")" |
|---|
| 2214 | rec=sublin(result,1) |
|---|
| 2215 | check=substr(rec,1,6) |
|---|
| 2216 | If (check = "Notice") |
|---|
| 2217 | rec=sublin(result,9) |
|---|
| 2218 | Else |
|---|
| 2219 | rec=sublin(result,8) |
|---|
| 2220 | Endif |
|---|
| 2221 | VAbove=subwrd(rec,4) |
|---|
| 2222 | "set z "zlev-1 |
|---|
| 2223 | PBelow=subwrd(result,4) |
|---|
| 2224 | "d "array"(lev="PBelow")" |
|---|
| 2225 | rec=sublin(result,1) |
|---|
| 2226 | check=substr(rec,1,6) |
|---|
| 2227 | If (check = "Notice") |
|---|
| 2228 | rec=sublin(result,9) |
|---|
| 2229 | Else |
|---|
| 2230 | rec=sublin(result,8) |
|---|
| 2231 | Endif |
|---|
| 2232 | VBelow=subwrd(rec,4) |
|---|
| 2233 | |
|---|
| 2234 | * Now if we are in a region of missing data, find next good level. |
|---|
| 2235 | |
|---|
| 2236 | If (abs(VAbove) > 130 & zlev > 1 & zlev < _zmaxfile) |
|---|
| 2237 | zz=zlev |
|---|
| 2238 | While (abs(VAbove) > 130 & zz < _zmaxfile) |
|---|
| 2239 | zz=zz+1 |
|---|
| 2240 | "set z "zz |
|---|
| 2241 | PAbove=subwrd(result,4) |
|---|
| 2242 | "d "array"(lev="PAbove")" |
|---|
| 2243 | rec=sublin(result,1) |
|---|
| 2244 | check=substr(rec,1,6) |
|---|
| 2245 | If (check = "Notice") |
|---|
| 2246 | rec=sublin(result,9) |
|---|
| 2247 | Else |
|---|
| 2248 | rec=sublin(result,8) |
|---|
| 2249 | Endif |
|---|
| 2250 | VAbove=subwrd(rec,4) |
|---|
| 2251 | EndWhile |
|---|
| 2252 | Endif |
|---|
| 2253 | |
|---|
| 2254 | If (abs(VBelow) > 130 & zlev > 1 & zlev < _zmaxfile) |
|---|
| 2255 | zz=zlev-1 |
|---|
| 2256 | While (abs(VBelow) > 130 & zz > 1) |
|---|
| 2257 | zz=zz-1 |
|---|
| 2258 | "set z "zz |
|---|
| 2259 | PBelow=subwrd(result,4) |
|---|
| 2260 | "d "array"(lev="PBelow")" |
|---|
| 2261 | rec=sublin(result,1) |
|---|
| 2262 | check=substr(rec,1,6) |
|---|
| 2263 | If (check = "Notice") |
|---|
| 2264 | rec=sublin(result,9) |
|---|
| 2265 | Else |
|---|
| 2266 | rec=sublin(result,8) |
|---|
| 2267 | Endif |
|---|
| 2268 | VBelow=subwrd(rec,4) |
|---|
| 2269 | EndWhile |
|---|
| 2270 | Endif |
|---|
| 2271 | |
|---|
| 2272 | If (abs(VAbove) < 130 & abs(VBelow) < 130) |
|---|
| 2273 | Vest=VBelow+log(PBelow/pres)*(VAbove-VBelow)/log(PBelow/PAbove) |
|---|
| 2274 | Else |
|---|
| 2275 | Vest=-9999.0 |
|---|
| 2276 | Endif |
|---|
| 2277 | |
|---|
| 2278 | Endif |
|---|
| 2279 | |
|---|
| 2280 | Return(Vest) |
|---|
| 2281 | |
|---|
| 2282 | **************************************************************************** |
|---|
| 2283 | |
|---|
| 2284 | function GetUWnd(wspd,wdir) |
|---|
| 2285 | |
|---|
| 2286 | *------------------------ |
|---|
| 2287 | * Get x-component of wind. |
|---|
| 2288 | *------------------------ |
|---|
| 2289 | |
|---|
| 2290 | |
|---|
| 2291 | If (wspd >= 0) |
|---|
| 2292 | xwind=wspd*cos((270-wdir)*_dtr) |
|---|
| 2293 | Else |
|---|
| 2294 | xwind = -9999.0 |
|---|
| 2295 | Endif |
|---|
| 2296 | return(xwind) |
|---|
| 2297 | |
|---|
| 2298 | ************************************************************************** |
|---|
| 2299 | |
|---|
| 2300 | function GetVWnd(wspd,wdir) |
|---|
| 2301 | |
|---|
| 2302 | *----------------------- |
|---|
| 2303 | * Get y-component of wind |
|---|
| 2304 | *------------------------ |
|---|
| 2305 | |
|---|
| 2306 | If (wspd >= 0) |
|---|
| 2307 | ywind=wspd*sin((270-wdir)*_dtr) |
|---|
| 2308 | Else |
|---|
| 2309 | ywind = -9999.0 |
|---|
| 2310 | Endif |
|---|
| 2311 | return(ywind) |
|---|
| 2312 | |
|---|
| 2313 | |
|---|
| 2314 | ************************************************************************* |
|---|
| 2315 | |
|---|
| 2316 | function GetWSpd(xwind,ywind) |
|---|
| 2317 | |
|---|
| 2318 | |
|---|
| 2319 | "set gxout stat" |
|---|
| 2320 | "d mag("xwind","ywind")" |
|---|
| 2321 | rec=sublin(result,1) |
|---|
| 2322 | check=substr(rec,1,6) |
|---|
| 2323 | If (check = "Notice") |
|---|
| 2324 | rec=sublin(result,9) |
|---|
| 2325 | Else |
|---|
| 2326 | rec=sublin(result,8) |
|---|
| 2327 | Endif |
|---|
| 2328 | val=subwrd(rec,4) |
|---|
| 2329 | |
|---|
| 2330 | return (val) |
|---|
| 2331 | |
|---|
| 2332 | ************************************************************************* |
|---|
| 2333 | |
|---|
| 2334 | function GetWDir(xwind,ywind) |
|---|
| 2335 | |
|---|
| 2336 | * Return wind direction given x and y components |
|---|
| 2337 | |
|---|
| 2338 | "set gxout stat" |
|---|
| 2339 | "define theta=270-"_rtd"*atan2("ywind","xwind")" |
|---|
| 2340 | "d theta" |
|---|
| 2341 | rec=sublin(result,1) |
|---|
| 2342 | check=substr(rec,1,6) |
|---|
| 2343 | If (check = "Notice") |
|---|
| 2344 | rec=sublin(result,9) |
|---|
| 2345 | Else |
|---|
| 2346 | rec=sublin(result,8) |
|---|
| 2347 | Endif |
|---|
| 2348 | Dir=subwrd(rec,4) |
|---|
| 2349 | |
|---|
| 2350 | If (Dir > 360) |
|---|
| 2351 | Dir=Dir-360 |
|---|
| 2352 | Endif |
|---|
| 2353 | |
|---|
| 2354 | If (Dir < 0) |
|---|
| 2355 | Dir=360+Dir |
|---|
| 2356 | Endif |
|---|
| 2357 | |
|---|
| 2358 | return(Dir) |
|---|
| 2359 | |
|---|
| 2360 | ************************************************************************* |
|---|
| 2361 | |
|---|
| 2362 | function GetXLoc(temp,pres) |
|---|
| 2363 | |
|---|
| 2364 | *------------------------------------------------- |
|---|
| 2365 | * Get x-location on skew-t based on temp, pressure |
|---|
| 2366 | *------------------------------------------------- |
|---|
| 2367 | |
|---|
| 2368 | xloc=(temp-_m1*log10(pres)-_m3)/_m2 |
|---|
| 2369 | return(xloc) |
|---|
| 2370 | |
|---|
| 2371 | ************************************************************************* |
|---|
| 2372 | |
|---|
| 2373 | function GetTemp(xloc,pres) |
|---|
| 2374 | |
|---|
| 2375 | *------------------------------------------------- |
|---|
| 2376 | * Return temperature at location given by xloc,pres |
|---|
| 2377 | *------------------------------------------------- |
|---|
| 2378 | |
|---|
| 2379 | tempval=_m1*log10(pres)+_m2*xloc+_m3 |
|---|
| 2380 | return(tempval) |
|---|
| 2381 | |
|---|
| 2382 | ************************************************************************** |
|---|
| 2383 | |
|---|
| 2384 | function GetTheta(temp,pres) |
|---|
| 2385 | |
|---|
| 2386 | *--------------------------------------------------- |
|---|
| 2387 | * Calculate theta for a given temperature and pressure |
|---|
| 2388 | *--------------------------------------------------- |
|---|
| 2389 | |
|---|
| 2390 | theta=(temp+273.15)*pow(1000/pres,0.286)-273.15 |
|---|
| 2391 | return(theta) |
|---|
| 2392 | |
|---|
| 2393 | |
|---|
| 2394 | ************************************************************************* |
|---|
| 2395 | |
|---|
| 2396 | function GetThet2(temp,dewp,pres) |
|---|
| 2397 | |
|---|
| 2398 | *--------------------------------------------------- |
|---|
| 2399 | * Calculate theta for a given temperature,dewp, and pressure |
|---|
| 2400 | *--------------------------------------------------- |
|---|
| 2401 | |
|---|
| 2402 | tempk=273.15+temp |
|---|
| 2403 | dewpk=273.15+dewp |
|---|
| 2404 | |
|---|
| 2405 | es=satvap2(temp) |
|---|
| 2406 | ws=mixratio(es,pres) |
|---|
| 2407 | |
|---|
| 2408 | mix=10*getrh(temp,dewp,pres)*ws |
|---|
| 2409 | |
|---|
| 2410 | exponent=0.2854*(1.0-0.00028*mix) |
|---|
| 2411 | theta=(temp+273.15)*pow(1000/pres,exponent)-273.15 |
|---|
| 2412 | return(theta) |
|---|
| 2413 | |
|---|
| 2414 | ************************************************************************* |
|---|
| 2415 | |
|---|
| 2416 | function Thetae(temp,dewp,pres) |
|---|
| 2417 | |
|---|
| 2418 | *-------------------------------------------------------------- |
|---|
| 2419 | * Return equiv. pot. temp in Kelvin given temp, dewp in celsius |
|---|
| 2420 | *-------------------------------------------------------------- |
|---|
| 2421 | |
|---|
| 2422 | es=satvap2(temp) |
|---|
| 2423 | ws=mixratio(es,pres) |
|---|
| 2424 | mix=10*getrh(temp,dewp,pres)*ws |
|---|
| 2425 | theta=GetThet2(temp,dewp,pres)+273.15 |
|---|
| 2426 | TLcl=Templcl(temp,dewp)+273.15 |
|---|
| 2427 | thetae=theta*exp((3.376/TLcl-0.00254)*mix*1.0+0.00081*mix) |
|---|
| 2428 | |
|---|
| 2429 | return(thetae) |
|---|
| 2430 | |
|---|
| 2431 | ************************************************************************** |
|---|
| 2432 | |
|---|
| 2433 | |
|---|
| 2434 | function int(i0) |
|---|
| 2435 | |
|---|
| 2436 | *-------------------------- |
|---|
| 2437 | * Return integer of i0 |
|---|
| 2438 | *-------------------------- |
|---|
| 2439 | i=0 |
|---|
| 2440 | while(i<12) |
|---|
| 2441 | i=i+1 |
|---|
| 2442 | if(substr(i0,i,1)='.') |
|---|
| 2443 | i0=substr(i0,1,i-1) |
|---|
| 2444 | break |
|---|
| 2445 | endif |
|---|
| 2446 | endwhile |
|---|
| 2447 | return(i0) |
|---|
| 2448 | |
|---|
| 2449 | ************************************************************************* |
|---|
| 2450 | |
|---|
| 2451 | function abs(i) |
|---|
| 2452 | |
|---|
| 2453 | *---------------------------- |
|---|
| 2454 | * return absolute value of i |
|---|
| 2455 | *---------------------------- |
|---|
| 2456 | |
|---|
| 2457 | if (i < 0) |
|---|
| 2458 | absval=-i |
|---|
| 2459 | else |
|---|
| 2460 | absval=i |
|---|
| 2461 | endif |
|---|
| 2462 | |
|---|
| 2463 | return(absval) |
|---|
| 2464 | |
|---|
| 2465 | ************************************************************************* |
|---|
| 2466 | |
|---|
| 2467 | function log(i) |
|---|
| 2468 | |
|---|
| 2469 | *--------------------------- |
|---|
| 2470 | * return natural log of i |
|---|
| 2471 | *--------------------------- |
|---|
| 2472 | |
|---|
| 2473 | "set gxout stat" |
|---|
| 2474 | "d log("i")" |
|---|
| 2475 | rec=sublin(result,1) |
|---|
| 2476 | check=substr(rec,1,6) |
|---|
| 2477 | If (check = "Notice") |
|---|
| 2478 | rec=sublin(result,9) |
|---|
| 2479 | Else |
|---|
| 2480 | rec=sublin(result,8) |
|---|
| 2481 | Endif |
|---|
| 2482 | val=subwrd(rec,4) |
|---|
| 2483 | return(val) |
|---|
| 2484 | |
|---|
| 2485 | ************************************************************************* |
|---|
| 2486 | |
|---|
| 2487 | function log10(i) |
|---|
| 2488 | |
|---|
| 2489 | *-------------------------- |
|---|
| 2490 | * return log base 10 of i |
|---|
| 2491 | *-------------------------- |
|---|
| 2492 | |
|---|
| 2493 | "set gxout stat" |
|---|
| 2494 | "d log10("i")" |
|---|
| 2495 | rec=sublin(result,1) |
|---|
| 2496 | check=substr(rec,1,6) |
|---|
| 2497 | If (check = "Notice") |
|---|
| 2498 | rec=sublin(result,9) |
|---|
| 2499 | Else |
|---|
| 2500 | rec=sublin(result,8) |
|---|
| 2501 | Endif |
|---|
| 2502 | val=subwrd(rec,4) |
|---|
| 2503 | return(val) |
|---|
| 2504 | |
|---|
| 2505 | ************************************************************************* |
|---|
| 2506 | |
|---|
| 2507 | function pow(i,j) |
|---|
| 2508 | |
|---|
| 2509 | *------------------------------- |
|---|
| 2510 | * return power of i raised to j |
|---|
| 2511 | *------------------------------- |
|---|
| 2512 | |
|---|
| 2513 | "set gxout stat" |
|---|
| 2514 | "d pow("i","j")" |
|---|
| 2515 | rec=sublin(result,1) |
|---|
| 2516 | check=substr(rec,1,6) |
|---|
| 2517 | If (check = "Notice") |
|---|
| 2518 | rec=sublin(result,9) |
|---|
| 2519 | Else |
|---|
| 2520 | rec=sublin(result,8) |
|---|
| 2521 | Endif |
|---|
| 2522 | val=subwrd(rec,4) |
|---|
| 2523 | return(val) |
|---|
| 2524 | |
|---|
| 2525 | ************************************************************************ |
|---|
| 2526 | |
|---|
| 2527 | function cos(i) |
|---|
| 2528 | |
|---|
| 2529 | *----------------------------------------- |
|---|
| 2530 | * return cosine of i, where i is in radians |
|---|
| 2531 | *------------------------------------------ |
|---|
| 2532 | |
|---|
| 2533 | "set gxout stat" |
|---|
| 2534 | "d cos("i")" |
|---|
| 2535 | rec=sublin(result,1) |
|---|
| 2536 | check=substr(rec,1,6) |
|---|
| 2537 | If (check = "Notice") |
|---|
| 2538 | rec=sublin(result,9) |
|---|
| 2539 | Else |
|---|
| 2540 | rec=sublin(result,8) |
|---|
| 2541 | Endif |
|---|
| 2542 | val=subwrd(rec,4) |
|---|
| 2543 | return(val) |
|---|
| 2544 | |
|---|
| 2545 | ************************************************************************ |
|---|
| 2546 | |
|---|
| 2547 | function sin(i) |
|---|
| 2548 | |
|---|
| 2549 | *------------------------------------------ |
|---|
| 2550 | * return sine of i, where i is in radians |
|---|
| 2551 | *------------------------------------------ |
|---|
| 2552 | |
|---|
| 2553 | "set gxout stat" |
|---|
| 2554 | "d sin("i")" |
|---|
| 2555 | rec=sublin(result,1) |
|---|
| 2556 | check=substr(rec,1,6) |
|---|
| 2557 | If (check = "Notice") |
|---|
| 2558 | rec=sublin(result,9) |
|---|
| 2559 | Else |
|---|
| 2560 | rec=sublin(result,8) |
|---|
| 2561 | Endif |
|---|
| 2562 | val=subwrd(rec,4) |
|---|
| 2563 | return(val) |
|---|
| 2564 | |
|---|
| 2565 | ************************************************************************ |
|---|
| 2566 | |
|---|
| 2567 | function exp(i) |
|---|
| 2568 | |
|---|
| 2569 | *------------------------------------------ |
|---|
| 2570 | * return exponential of i |
|---|
| 2571 | *------------------------------------------ |
|---|
| 2572 | |
|---|
| 2573 | "set gxout stat" |
|---|
| 2574 | "d exp("i")" |
|---|
| 2575 | rec=sublin(result,1) |
|---|
| 2576 | check=substr(rec,1,6) |
|---|
| 2577 | If (check = "Notice") |
|---|
| 2578 | rec=sublin(result,9) |
|---|
| 2579 | Else |
|---|
| 2580 | rec=sublin(result,8) |
|---|
| 2581 | Endif |
|---|
| 2582 | val=subwrd(rec,4) |
|---|
| 2583 | return(val) |
|---|
| 2584 | |
|---|
| 2585 | *********************************************************************** |
|---|
| 2586 | function round(i) |
|---|
| 2587 | |
|---|
| 2588 | rr=abs(1.0*i) |
|---|
| 2589 | rr=int(rr+0.5) |
|---|
| 2590 | if (i < 0) |
|---|
| 2591 | rr=-1*rr |
|---|
| 2592 | endif |
|---|
| 2593 | return(rr) |
|---|
