<html>
<head>
<style><!--
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Tahoma
}
--></style>
</head>
<body class='hmmessage'>
Hello,<br><br>My name is Cristian. I am having some troubles to perform some skewt graphics using WRF output.<br>I have created a binary and .ctl archives using ARWpost.<br>The file is OK, and no troubles are expected.<br><br>The problem starts when I want to create a skewt graphic.<br>The temperature and dewpoint profiles are OK, but no information about wind is shown. Neither the Hodograpf, nor the wind speed and direction profiles are shown.<br><br>No message of an error apears, and i have no idea what is happening. <br>Even the undef value is th same as the .ctl archive.<br><br>I am using Grads version 1.a8, in UBUNTU 9.04.<br><br>Here is the script. I only changed parameters at the top of the archive.<br><br>I really appreciate your help.<br><br>Thanks<br><br>-----------------------------------------------------------------------------------------------------------------------------------------<br>* Script para construir radiosondeos<br><br>'reinit'<br>'run /home/cristian/Escritorio/Experimentos/lib/jaecol'<br>'set display color white'<br>'reset'<br><br>* Defino el archivo a utilizar<br><br>'open 20081004_YSU.ctl'<br>*path=''<br><br>* Latitud Longitud y tiempo<br><br>lat=-46<br>lon=-68<br>time=12z04Oct2008<br><br><br>'set time 'time<br>'set lon 'lon<br>'set lat 'lat<br><br>*Entre que niveles hace el perfil hacerlo entre 1000 y 600<br><br>levmax=925<br>levmin=500<br><br>'set lev 'levmax ' ' levmin<br><br>*calculo Td a partir de HR, P y T<br>*'define et=(-2937.4/(t))-(4.9283*log10(t))+22.5518'<br>*'define es=pow(10,et)*10'<br>*'define e=es*0.01*hr'<br>*'define q=0.622*e/(lev-e)'<br>*'define aux1=q*1000*(1.+0.81*q)'<br>*'define aux2=0.2854*(1.-0.28*q)'<br>*'define aux3=(1000/lev)'<br>*'define td=243.5/((17.67/(log(e/6.11)))-1)'<br>*'define tc=t-273'<br><br>'set ylopts 1 5 0.18'<br><br>plotskew(tc,td,u,v)<br><br>'q dims'<br>line1=sublin(result,4)<br>line2=sublin(result,5)<br>itime1=subwrd(line1,6)<br>itit=substr(itime1,1,6)<br>itime1=subwrd(line2,6)<br>itime=substr(itime1,1,12)<br>d=substr(itime1,4,2)<br>m=substr(itime1,6,3)<br>y=substr(itime1,9,4)<br>h=substr(itime1,1,2)<br>lati=-1*lat<br>loni=-1*lon<br><br>'draw title Perfil vertical de T y Td en 'lati'S y 'loni'W\'itime<br>'printim 'd'.'h'sondeo.png png x600 y800 white'<br><br>function plotskew(sndtemp,snddewp,u,v)<br><br>*El valor Undef que se coloca acá debe coincidir con el valor undef del CTL, sinó hay problemas<br>*en la construcción de la hodógrafa y del perfil vertical de viento.<br>undef=1.e30<br>*************************************************************************<br>*<br>* GrADS Script to Plot a SkewT/LogP Diagram <br>*<br>* Bob Hart <br>* Penn State University / Dept of Meteorology<br>* Last Update: January 23, 2001<br>*<br>* Recent Changes:<br>*<br>* 01/23/01 - Fixed a small bug in the theta-e calculation. <br>* Errors averaged 0.5-3K. Thank you George Bryan.<br>*<br>* 11/10/99 - Change in calculation method for CAPE/CIN. Trapezoid <br>* integration method is now used. Speeds up execution<br>* by 25%, and increases accuracy by 5-10%.<br>*<br>* 10/18/99 - Minor glitch fixed that occasionally caused crash.<br>*<br>* 8/26/99 - Datasets with missing data can now be used.<br>*<br>* Features:<br>* - All features of standard skewt/logp plot<br>* - RH sounding <br>* - LCL location<br>* - Parcel trajectory for both sfc based convection and elevated from<br>* most unstable level (highest theta-e level reported)<br>* - Stability indices and precipitable water calculations<br>* - CAPE & CIN Calculations<br>* - Wind Profile<br>* - Hodograph / Hodograph scaling<br>* - Helicity and SR Helicity Calculations and Display<br>* - Color aspects of output<br>* - Line Thickness, style aspects of output<br>* - Can be run in either PORTRAIT or LANDSCAPE mode.<br>*<br>* There are numerous tunable parameters below to change the structure<br>* and output for the diagram.<br>*<br>* Function Arguments:<br>* sndtemp - temperature data (Celsius) as a function of pressure<br>* snddewp - dewpoint data (Celsius) as a function of pressure <br>* sndspd - wind speed data (knots) as a function of pressure<br>* snddir - wind direction data as a function of pressure<br>*<br>* Use '-1' for any of the above 4 arguments to indicate that you<br>* are not passing that variable. The appropriate options will <br>* be ignored based on your specifying '-1' for that variable.<br>*<br>* NOTE: Make sure to set the vertical range of the plot before running.<br>* I.e., "SET LEV 1050 150", for example. This does not have to <br>* be limited to the pressure range of your data.<br>*<br>* Labelling: Pressure/Height is labelled along left side. Temperature is <br>* labelled along bottom. Mixing ratio is labelled along right<br>* side/top.<br>*<br>* <br>* PROBLEMS: First check out the web page for the script (which also<br>* has a link to a FAQ with answers to many common questions<br>* about using the script):<br>* http://www.ems.psu.edu/~hart/skew.html<br>*<br>* Please send any further problems, comments, or suggestions to <br>* <hart@ems.psu.edu><br>*<br>* ACKNOWLEDGMENTS: Thanks go to the innumerable users who have helped<br>* fine tune the script from the horrible mess from which it began.<br>* In particular, thanks go out to Steve Lord (NCEP), Mike Fiorino (ECMWF),<br>* George Bryan (PSU), Davide Sacchetti (CMIRL), and Enrico Minguzzi (CMIRL).<br>*<br>**************************************************************************<br>* !!!!! BEGINNING OF USER-SPECIFIED OPTIONS !!!!!!<br>**************************************************************************<br>*<br>* --------------------- Initialization options ----------------------<br>*<br>* ClrScrn = Whether to clear the screen before drawing diagram<br>* [1 = yes, 0 = no]<br><br>ClrScrn = 1<br><br>*<br>* ------------------- Define Skew-T Diagram Shape/Slope-----------------<br>*<br>* (P1,T1) = Pres, Temp of some point on left-most side<br>* (P2,T2) = Pres, Temp of some point on right-most side<br>* (P3,T3) = Pres, Temp of some point in diagram which is mid-point <br>* in the horizontal between 1 and 2.<br>*<br>* P1, P2, P3 are in mb ; T1, T2, T3 are in Celsius<br>*<br>* These define the SLOPE and WIDTH of the diagram as you see it but DO NOT<br>* DEFINE THE HEIGHT of the diagram as you see it. In other words,<br>* 1 and 2 do NOT necessarily need to be at the bottom of the diagram and<br>* 3 does NOT necessarily need to be at the top. THE VERTICAL PRESSURE <br>* RANGE OF THE SKEWT AS YOU SEE IT IS DETERMINED BY YOUR 'SET Z ...' <br>* COMMAND OR THE 'SET LEV ...' COMMAND BEFORE RUNNING THIS SCRIPT.<br>*<br>* _______________________<br>* | |<br>* | |<br>* | 3 |<br>* | |<br>* | |<br>* | |<br>* | |<br>* | |<br>* | |<br>* | |<br>* | |<br>* |1 2|<br>* | |<br>* |_______________________|<br>* <br>*<br>* A good set of defining points are given below. Feel free<br>* to experiment with variations.<br><br><br>*P1 = 1000<br>*T1 = -40<br><br>*P2 = 1000<br>*T2 = 40<br><br>*P3 = 200<br>*T3 = 0<br><br>* Another good set of defining points suggested by Juan Ruiz (Emagrama)<br>* are:<br>*<br> P1 = 1000<br> T1 = -80<br>*<br> P2 = 1000<br> T2 = 40<br>*<br> P3 = 500<br> T3 = -20<br><br>* ------------------- Contour Intervals / Levels --------------------------<br>*<br>* All variables below are contour intervals/levels for diagram<br>*<br>* Thetaint = interval for potential temperature lines<br>* Thetwint = interval for moist pseudo adiabats<br>* tempint = interval for temperature lines<br>* wsclevs = contour LEVELS for mixing ratio lines<br>*<br>*<br>thetaint= 20<br>thetwint= 10<br>tempint = 10<br>wsclevs = ".1 .5 1 3 6 10 15 20 25 30"<br>*<br>*<br>* ------------------------ Output Options --------------------------------<br>*<br>* All variables below are logical .. 1=yes, 0=no, unless otherwise<br>* specified.<br>*<br>* DrawBarb = Draw wind barbs along right side of plot<br>* DrawThet = Draw dry adiabats<br>* DrawThtw = Draw moist pseudo-adiabats<br>* DrawTemp = Draw temperature lines <br>* DrawMix = Draw mixing ratio lines<br>* DrawTSnd = Draw temperature sounding<br>* DrawDSnd = Draw dewpoint sounding<br>* DrawRH = Draw relative humidity sounding<br>* DrawPrcl = Draw parcel path from surface upward<br>* DrawPMax = Draw parcel path from most unstable level upward<br>* DrawIndx = Display stability indices & CAPE<br>* DrawHeli = Calculate and display absolute and storm-relative helicity<br>* DrawHodo = Draw hodograph<br>* DrawPLev = Draw Pressure Levels <br>* DrawZLev = Draw height levels and lines <br>* 0 = no lines<br>* 1 = above ground level (AGL)<br>* 2 = above sea level (ASL)<br>* DrawZSTD = Draw Height levels using standard atm lapse rate <br>* LblAxes = Label the x,y axes (temperature, pressure,mixing ratio)<br>*<br>* ThtwStop = Pressure level at which to stop drawing Theta-w lines<br>* MixStop = Pressure level at which to stop drawing Mixratio lines<br><br>DrawBarb= 1<br>DrawThet= 1<br>DrawThtw= 1<br>DrawTemp= 1<br>DrawMix = 1<br>DrawTSnd= 1<br>DrawDSnd= 1<br>DrawRH = 0<br>DrawPrcl= 1<br>DrawPMax= 1<br>DrawIndx= 1<br>DrawHeli= 0<br>DrawHodo= 1<br>DrawPLev= 1<br>DrawZLev= 0<br>DrawZSTD= 0<br>LblAxes = 1<br><br>ThtwStop = 200<br>MixStop = 400<br><br><br>* <br>* ----------------- Sounding Geography options ------------------------<br>*<br>* SfcElev = Elevation above sea-level (meters) of lowest level reported<br>* in sounding. Used only if DrawZLev = 2<br><br>SfcElev = 0<br><br><br>*<br>* ------------------ Thermodynamic Index Options --------------------<br>* <br>* All variables here are in inches. Use -1 for the default values.<br>*<br>* Text1XC = X-location of midpoint of K,TT,PW output box<br>* Text1YC = Y-location of midpoint of K,TT,PW output box<br>* Text2XC = X-Location of midpoint of surface indices output box<br>* Text2YC = Y-location of midpoint of surface indices output box<br>* Text3XC = X-Location of midpoint of most unstable level-based indices<br>* output box<br>* Text3YC = Y-location of midpoint of most unstable level-based indices<br>* output box<br><br>Text1XC = -1<br>Text1YC = -1<br>Text2XC = -1<br>Text2YC = -1<br>Text3XC = -1<br>Text3YC = -1<br><br>*<br>* ----------------- Wind Barb Profile Options ----------------------------<br>*<br>* All variables here are in units of inches, unless otherwise specified<br>*<br>* barbint = Interval for plotting barbs (in units of levels)<br>* poleloc = X-Location of profile. Choose -1 for the default.<br>* polelen = Length of wind-barb pole<br>* Len05 = Length of each 5-knot barb <br>* Len10 = Length of each 10-knot barb<br>* Len50 = Length of each 50-knot flag<br>* Wid50 = Width of base of 50-knot flag <br>* Spac50 = Spacing between 50-knot flag and next ror occurred on libarb/flag <br>* Spac10 = Spacing between 10-knot flag and next flag<br>* Spac05 = Spacing between 5-knot flag and next flag<br>* Flagbase= Draw flagbase (filled circle) for each windbarb [1=yes, 0 =no] <br>* Fill50 = Solid-fill 50-knot flag [1=yes, 0=no]<br>* barbline= Draw a vertical line connecting all the wind barbs [1=yes, 0=no]<br>*<br>barbint = 1<br>poleloc = -1<br>polelen = 0.35<br>len05 = 0.07<br>len10 = 0.15<br>len50 = 0.15<br>wid50 = 0.06<br>spac50 = 0.07<br>spac10 = 0.05<br>spac05 = 0.05<br>Fill50 = 1<br>flagbase= 1<br>barbline= 1<br><br>*<br>*<br>*---------------- Hodograph Options -------------------------------------<br>*<br>* All variables here are in units of inches, unless otherwise specified<br>*<br>* HodXcent= x-location of hodograph center. Use -1 for default location.<br>* HodYcent= y-location of hodograph center. Use -1 for default location.<br>* HodSize = Size of hodograph in inches <br>* NumRing = Number of rings to place in hodograph (must be at least 1)<br>* HodRing = Wind speed increment of each hodograph ring<br>* HodoDep = Depth (above lowest level in mb) of end of hodograph trace<br>* TickInt = Interval (in kts) at which tick marks are drawn along the axes<br>* Use 0 for no tick marks.<br>* TickSize= Size of tick mark in inches<br>* Text4XC = X-location of midpoint of hodograph text output. Use -1 for default.<br>* Text4YC = Y-location of midpoint of hodograph text output. Use -1 for default.<br><br>HodXcent= 6<br>HodYcent= 9<br>HodSize = 2<br>NumRing = 3<br>HodRing = 10<br>HodoDep = 500<br>TickInt = 5<br>TickSize= 0.1<br>Text4XC = -1<br>Text4YC = -1<br><br>*--------------- Helicity Options ---------------------------------------<br>*<br>* MeanVTop = Top pressure level (mb) of mean-wind calculation<br>* MeanVBot = Bottom pressure level (mb) of mean-wind calculation<br>* HelicDep = Depth in mb (above ground) of helicity integration<br>* StormMot = Type of storm motion estimation scheme. Use following: <br>* 0 = No departure from mean wind.<br>* 1 = Davies-Jones (1990) approach<br>* FillArrw = Whether to fill the arrowhead of the storm motion vector<br>* [1 = yes, 0 = no]<br><br>MeanVTop= 300 <br>MeanVBot= 850<br>HelicDep= 300<br>StormMot= 0<br>FillArrw= 1<br><br>*<br>*---------------- Color Options ------------------------------------------<br>*<br>* ThetCol = Color of dry adiabats<br>* TempCol = Color of temperature lines<br>* MixCol = Color of mixing ratio lines<br>* ThtwCol = Color of moist adiabats<br>* TSndCol = Color of Temperature Sounding <br>* DSndCol = Color of Dewpoint Sounding<br>* RHCol = Color of RH Sounding<br>* PrclCol = Color of parcel trace<br>* BarbCol = Color of wind barbs (choose -1 for color according to speed)<br>* HodoCol = Color of hodograph trace<br><br>ThetCol = 23<br>TempCol = 79<br>MixCol = 38<br>ThtwCol = 39<br>TSndCol = 29 <br>DSndCol = 49<br>RHCol = 3<br>PrclCol = 1<br>BarbCol = 1<br>HodoCol = 2<br><br>*<br>*-------------------- Line Style Options ------------------------------------<br>*<br>* GrADS Styles: 1=solid;2=long dash;3=short dash;4=long,short dashed;<br>* 5=dotted;6=dot dash;7=dot dot dash<br>*<br>* ThetLine = Line Style of dry adiabats<br>* TempLine = Line Style of temperature lines<br>* MixLine = Line Style of mixing ratio lines<br>* ThtwLine = Line Style of moist adiabats<br>* TSndLine = Line Style of Temperature Sounding<br>* DSndLine = Line Style of Dewpoint Sounding<br>* RHLine = Line Style of RH sounding<br>* PrclLine = Line Style of parcel trace<br>* HodoLine = Line Style of hodograph trace<br>*<br><br>ThetLine = 1<br>TempLine = 1<br>MixLine = 5<br>ThtwLine = 2<br>TSndLine = 1<br>DSndLine = 2<br>RHLine = 1<br>PrclLine = 3<br>HodoLine = 1<br><br>*<br>*------------------- Line Thickness Options---------------------------------<br>* GrADS Line Thickness: increases with increasing number. Influences <br>* hardcopy output more strongly than screen output.<br>*<br>*<br>* ThetThk = Line Thickness of dry adiabats<br>* TempThk = Line Thickness of temperature lines<br>* MixThk = Line Thickness of mixing ratio lines<br>* ThtwThk = Line Thickness of moist adiabats<br>* TSndThk = Line Thickness of temperature sounding<br>* DSndThk = Line thickness of dewpoint sounding<br>* RHThk = Line thickness of RH sounding<br>* PrclThk = Line thickness of parcel trace<br>* HodoThk = Line thickness of hodograph trace<br>* BarbThk = Line thickness of wind barbs<br><br>ThetThk = 3<br>TempThk = 1<br>MixThk = 7<br>ThtwThk = 3<br>TSndThk = 12<br>DSndThk = 12<br>RHThk = 8<br>PrclThk = 6 <br>HodoThk = 6 <br>BarbThk = 2<br><br>*<br>*------------------- Data Point Marker Options -----------------------------<br>* GrADS Marker Types: 0 = none ; 1 = cross ; 2 = open circle ; <br>* 3 = closed circle ; 4 = open square ; 5 = closed square<br>* 6 = X ; 7 = diamond ; 8 = triangle ; 9 = none<br>* 10 = open circle with vertical line ; 11 = open oval<br>*<br>* TSndMrk = Mark type of data point marker for temperature sounding<br>* DSndMrk = Mark type of data point marker for dewpoint sounding<br>* RHMrk = Mark type of data point marker for relative humidity sounding<br>* MrkSize = Mark size (inches) of each data marker<br><br>TSndMrk = 3<br>DSndMrk = 3<br>RHMrk = 0<br>MrkSize = 0.1<br><br><br>* !!!!! YOU SHOULD NOT NEED TO CHANGE ANYTHING BELOW HERE !!!!!<br>****************************************************************************<br><br>*-------------------------------------------<br>* grab user-specified environment dimensions<br>*-------------------------------------------<br><br>"q dims"<br>rec=sublin(result,2)<br>_xtype=subwrd(rec,3)<br>_xval=subwrd(rec,9)<br>rec=sublin(result,3)<br>_yval=subwrd(rec,9)<br>_ytype=subwrd(rec,3)<br>rec=sublin(result,4)<br>_ptype=subwrd(rec,3)<br>_pmax=subwrd(rec,6)<br>_pmin=subwrd(rec,8)<br>_zmin=subwrd(rec,11)<br>_zmax=subwrd(rec,13)<br>rec=sublin(result,5)<br>_ttype=subwrd(rec,3)<br>_tval=subwrd(rec,9)<br><br>"q file"<br>rec=sublin(result,5)<br>_zmaxfile=subwrd(rec,9)<br><br>*-------------------------------------------------------------<br>* Check to ensure that dimensions are valid. Warn & exit if not.<br>*--------------------------------------------------------------<br><br>dimrc=0<br>If (_xtype != "fixed")<br> say "X-Dims Error: Not fixed. Use 'set lon' or 'set x' to specify a value."<br> dimrc=-1<br>Endif<br><br>If (_ytype != "fixed")<br> say "Y-Dims Error: Not fixed. Use 'set lat' or 'set y' to specify a value"<br> dimrc=-1<br>Endif<br><br>If (_ptype != "varying")<br> say "Z-Dims Error: Not varying. Use 'set lev' or 'set z' to specify a range."<br> dimrc=-1<br>Endif<br><br>If (_ttype != "fixed")<br> say "Time Error: Not fixed. Use 'set time' or 'set t' to specify a value"<br> dimrc=-1<br>Endif<br><br><br>If (dimrc < 0)<br> Return(-1)<br>Endif<br><br><br>*<br>* A few global variables used in units conversion<br>*<br><br>_pi=3.14159265<br>_dtr=_pi/180<br>_rtd=1/_dtr <br>_ktm=0.514444<br>_mtk=1/_ktm<br><br>* A few global constants used in thermo calcs<br><br>_C0=0.99999683 <br>_C1=-0.90826951/100<br>_C2= 0.78736169/10000<br>_C3=-0.61117958/1000000<br>_C4= 0.43884187/pow(10,8) <br>_C5=-0.29883885/pow(10,10)<br>_C6= 0.21874425/pow(10,12)<br>_C7=-0.17892321/pow(10,14)<br>_C8= 0.11112018/pow(10,16) <br>_C9=-0.30994571/pow(10,19)<br><br>*Calculo la dirección e intensidad del viento usando las funciones GetWdir y GetWsp<br><br> <br><br><br><br><br>* A pressure array of power calculations which should be performed<br>* only once to reduce execution time.<br><br>zz=1100<br>while (zz > 10)<br> subscr=zz/10<br> _powpres.subscr=pow(zz,0.286)<br> zz=zz-10<br>endwhile<br><br>*<br>* Turn off options not available due to user data limitations<br>*<br><br>If (ClrScrn = 1) <br> "clear"<br>Endif<br><br>If (sndspd = -1 | snddir = -1) <br> DrawBarb = 0<br> DrawHodo = 0<br> DrawHeli = 0<br>Endif<br><br>If (snddewp = -1) <br> DrawDSnd = 0<br> DrawRH = 0<br> DrawPrcl = 0<br> DrawPMax = 0<br> DrawIndx = 0<br>Endif<br><br>If (sndtemp = -1)<br> DrawTSnd = 0<br> DrawRH = 0<br> DrawPrcl = 0<br> DrawPMax = 0<br> DrawIndx = 0<br> DrawZLev = 0<br>Endif<br><br>If (NumRing < 1) <br> DrawHodo = 0<br>Endif <br> <br>"q gxinfo"<br>rec=sublin(result,2)<br>xsize=subwrd(rec,4)<br><br>If (xsize = 11) <br> PageType = "Landscape"<br>Else<br> PageType = "Portrait"<br>Endif<br> <br>*------------------------------------------------------<br>* calculate constants determining slope/shape of diagram<br>* based on temp/pressure values given by user<br>*-------------------------------------------------------<br><br>"set x 1"<br>"set y 1"<br>"set z 1"<br>"set t 1"<br>_m1=(T1+T2-2*T3)/(2*log10(P2/P3))<br>_m2=(T2-T3-_m1*log10(P2/P3))/50<br>_m3=(T1-_m1*log10(P1))<br><br>"set z "_zmin" "_zmax <br>"set zlog on"<br>"set xlab off"<br><br>*-------------------------------------------------<br>* perform coordinate transformation to Skew-T/LogP<br>*-------------------------------------------------<br><br>"set gxout stat"<br>"set x "_xval<br>"set y "_yval<br>"set t "_tval<br>"define tempx=("sndtemp"-"_m1"*log10(lev)-"_m3")/"_m2<br>"define dewpx=("snddewp"-"_m1"*log10(lev)-"_m3")/"_m2<br><br>If (PageType = "Portrait") <br> "set parea 0.7 7 0.75 10"<br>Else<br> "set parea 0.7 6.5 0.5 8"<br>Endif<br> <br>"set axlim 0 100"<br>"set lon 0 100"<br>"set grid on 1 1"<br><br>"set z "_zmin " " _zmax<br>"set lon 0 100"<br>"set clevs -900"<br>"set gxout contour"<br><br>*-------------------------------------<br>* Draw pressure lines <br>*-------------------------------------<br><br>If (DrawPLev = 0) <br> "set ylab off"<br>Else<br> "set ylab on"<br> "set ylopts 1 5 0.18"<br> "set xlopts 1 3 0.18"<br>Endif<br><br>"d lon"<br><br>*--------------------------------------<br>* Determine corners of skewt/logp frame<br>*--------------------------------------<br><br>"q w2xy 100 "_pmin<br>rxloc=subwrd(result,3)<br>tyloc=subwrd(result,6)<br>"q w2xy 0 "_pmax<br>lxloc=subwrd(result,3)<br>byloc=subwrd(result,6)<br><br>If (DrawPLev = 1 & LblAxes = 1)<br> "set strsiz 0.15"<br> "set string 1 c 3 0"<br> If (PageType = "Portrait") <br>* "draw string 0.5 10.5 hPa."<br> Else<br> "draw string 0.5 8.35 hPa."<br> Endif<br>Endif<br><br>*---------------------------------------------------<br>* Calculate & draw actual height lines using temp data<br>*---------------------------------------------------<br><br>If (DrawZLev > 0)<br> say "Calculating observed height levels from temp/pressure data."<br> zz=1<br> "set gxout stat"<br> "set x "_xval<br> "set y "_yval<br> "set t "_tval<br> count=0<br> while (zz < _zmax)<br> "set z "zz<br> pp.zz=subwrd(result,4)<br> lpp.zz=log(pp.zz)<br> "d "sndtemp<br> rec=sublin(result,8)<br> tt=subwrd(rec,4)<br> if (tt > -900) <br> tk=tt+273.15<br> count=count+1<br> zzm=zz-1<br> If (count = 1) <br> If (DrawZLev = 2)<br> htlb="ASL"<br> height.zz=SfcElev<br> Else<br> htlb="AGL"<br> height.zz=0<br> Endif<br> sfcz=height.zz<br> Else<br> DZ=29.2857*(lpp.zzm-lpp.zz)*(lpp.zz*tk+lpp.zzm*tkold)/(lpp.zz+lpp.zzm) <br> height.zz=height.zzm+DZ<br> highz=height.zz<br> Endif<br> else<br> height.zz = -9999<br> endif<br> tkold=tk<br> zz=zz+1<br> endwhile<br><br> maxht=int(highz/1000)<br> if (int(sfcz/1000) = sfcz/1000)<br> minht=int(sfcz/1000)<br> else<br> minht=1+int(sfcz/1000)<br> endif<br><br> ht=minht<br> "set line 1 3 1"<br> "set strsiz 0.10"<br> "set string 1 l 3 0"<br> while (ht <= maxht)<br> zz=1<br> while (height.zz/1000 <= ht)<br> zz=zz+1<br> endwhile<br> zzm=zz-1<br> PBelow=pp.zzm<br> PAbove=pp.zz<br> HBelow=height.zzm<br> HAbove=height.zz<br> DZ=HAbove-HBelow<br> DP=PAbove-PBelow<br> Del=ht*1000-HBelow<br> Est=PBelow+Del*DP/DZ<br> If (Est >= _pmin & Est <= _pmax)<br> "q w2xy 1 " Est<br> yloc=subwrd(result,6)<br> "draw line " lxloc " " yloc " " rxloc " " yloc<br> "draw string 0.22 "yloc-0.05" "ht<br> Endif<br> ht=ht+1<br> endwhile<br> "set strsiz 0.10"<br> "set string 1"<br> If (LblAxes = 1)<br> If (PageType = "Portrait") <br> "draw string 0.25 10.85 km"<br> "draw string 0.25 10.75 "htlb<br> "draw string 0.25 10.65 OBS"<br> Else<br> "draw string 0.25 8.35 km"<br> "draw string 0.25 8.25 "htlb<br> "draw string 0.25 8.15 OBS"<br> Endif<br> Endif<br>Endif<br><br><br>*---------------------------------------------------<br>* Draw height levels (height above MSL using Std Atm)<br>*---------------------------------------------------<br><br>If (DrawZSTD = 1)<br> "set strsiz 0.10"<br> minht=30.735*(1-pow(_pmax/1013.26,0.287))<br> minht=int(minht+0.5)<br> maxht=30.735*(1-pow(_pmin/1013.26,0.287))<br> maxht=int(maxht)<br> "set gxout stat"<br> zcount=minht <br> while (zcount <= maxht) <br> plev=1013.26*pow((1-zcount/30.735),3.4843)<br> "q w2xy 0 "plev <br> yloc=subwrd(result,6)<br> "draw string 0 "yloc-0.05" "zcount<br> zcount=zcount+1<br> endwhile<br> "set strsiz 0.10"<br> If (LblAxes = 1)<br> If (PageType = "Portrait") <br> "draw string 0 10.85 km"<br> "draw string 0 10.75 ASL"<br> "draw string 0 10.65 STD"<br> Else<br> "draw string 0 8.35 km"<br> "draw string 0 8.25 ASL"<br> "draw string 0 8.15 STD"<br> Endif<br> Endif<br>Endif<br><br><br>*-----------------------<br>* Plot temperature lines <br>*-----------------------<br><br>If (DrawTemp = 1)<br> "set strsiz 0.1"<br> "set z "_zmin " " _zmax<br> "set line "TempCol " " TempLine " "TempThk<br> "set string 1 c 3 0"<br> "set gxout stat"<br> maxtline=GetTemp(100,_pmax)<br> mintline=GetTemp(0,_pmin)<br><br> maxtline=tempint*int(maxtline/tempint)<br> mintline=tempint*int(mintline/tempint)<br><br> tloop=mintline<br> While (tloop <= maxtline) <br> Botxtemp=GetXLoc(tloop,_pmax)<br> "q w2xy "Botxtemp " " _pmax<br> Botxloc=subwrd(result,3)<br> Botyloc=byloc <br> Topxtemp=GetXLoc(tloop,_pmin)<br> "q w2xy "Topxtemp " " _pmin<br> Topxloc=subwrd(result,3)<br> Topyloc=tyloc <br> If (Botxtemp <= 100 | Topxtemp <= 100) <br> If (Topxtemp > 100)<br> Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)<br> b=Topyloc-Slope*Topxtemp<br> Topyloc=Slope*100+b<br> Topxloc=rxloc <br> Endif<br> If (Botxtemp < 0)<br> Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)<br> b=Botyloc-Slope*Botxtemp<br> Botyloc=b<br> Botxloc=lxloc <br> Else<br> "set strsiz 0.15"<br> "draw string " Botxloc-0.05 " " Botyloc-0.15 " " tloop<br> Endif<br> "draw line "Botxloc " " Botyloc " " Topxloc " " Topyloc<br> Endif<br> tloop=tloop+tempint<br> EndWhile<br> If (LblAxes = 1) <br> "set strsiz 0.15"<br> "set string 1 c"<br> If (PageType = "Portrait")<br> "draw string 4.0 0.35 Temperatura (`3.`0C)"<br> "draw string 7.7 0.35 Viento m/s "<br> Else<br> "draw string 3.5 0.15 Temperatura (`3.`0C)"<br> Endif<br> Endif<br>Endif<br><br><br>*------------------<br>* Plot dry adiabats<br>*------------------<br><br>If (DrawThet = 1)<br> temp=GetTemp(100,_pmin)<br> maxtheta=GetThet2(temp,-100,_pmin)<br> maxtheta=thetaint*int(maxtheta/thetaint)<br> temp=GetTemp(0,_pmax)<br> mintheta=GetThet2(temp,-100,_pmax)<br> mintheta=thetaint*int(mintheta/thetaint)<br> <br> "set lon 0 100"<br> "set y 1"<br> "set z 1"<br> tloop=mintheta<br> "set line "ThetCol" "ThetLine " "ThetThk<br> While (tloop <= maxtheta)<br> PTemp=LiftDry(tloop,1000,_pmin,1,_pmin,_pmax) <br> tloop=tloop+thetaint<br> Endwhile<br>Endif<br><br>*------------------------<br>* Plot mixing ratio lines<br>*------------------------<br><br>If (DrawMix = 1)<br> If (MixStop < _pmin) <br> MixStop = _pmin<br> Endif<br> "set string 1 l"<br> "set z "_zmin " " _zmax<br> "set cint 1"<br> "set line "MixCol" " MixLine " "MixThk<br> cont = 1<br> mloop=subwrd(wsclevs,1)<br> count = 1<br> While (cont = 1) <br> BotCoef=log(mloop*_pmax/3801.66)<br> BotTval=-245.5*BotCoef/(BotCoef-17.67)<br> Botxtemp=GetXLoc(BotTval,_pmax)<br> "q w2xy "Botxtemp " " _pmax<br> Botxloc=subwrd(result,3)<br> Botyloc=byloc <br> TopCoef=log(mloop*MixStop/3801.66)<br> TopTval=-245.5*TopCoef/(TopCoef-17.67)<br> Topxtemp=GetXLoc(TopTval,MixStop)<br> "q w2xy "Topxtemp " " MixStop<br> Topxloc=subwrd(result,3)<br> Topyloc=subwrd(result,6) <br> "set string "MixCol" l 3"<br> "set strsiz 0.09"<br> If (Botxtemp <= 100 | Topxtemp <= 100) <br> If (Topxtemp > 100)<br> Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)<br> b=Topyloc-Slope*Topxtemp<br> Topyloc=Slope*100+b<br> Topxloc=rxloc <br> "draw string " Topxloc+0.05 " " Topyloc " " mloop<br> Else<br> "draw string " Topxloc " " Topyloc+0.1 " " mloop<br> Endif<br> If (Botxtemp < 0)<br> Slope=(Topyloc-Botyloc)/(Topxtemp-Botxtemp)<br> b=Botyloc-Slope*Botxtemp<br> Botyloc=b<br> Botxloc=lxloc <br> Endif<br> "draw line "Botxloc " " Botyloc " " Topxloc " " Topyloc<br> Endif<br> count=count+1<br> mloop=subwrd(wsclevs,count)<br> If (mloop = "" | count > 50) <br> cont = 0<br> Endif<br> EndWhile<br> If (LblAxes = 1)<br> "set strsiz 0.15"<br> "set string 1 c 3 90"<br> If (PageType = "Portrait")<br>* "draw string 7.40 4.75 Relacion de Mezcla (g/kg)"<br> Else<br>* "draw string 6.90 4.25 Mixing Ratio (g/kg)"<br> Endif<br> "set string 1 c 3 0"<br> Endif<br>Endif<br><br>*-----------------------------<br>* Plot moist (pseudo) adiabats<br>*-----------------------------<br><br>If (DrawThtw = 1)<br> "set lon 0 100"<br> "set y 1"<br> "set z 1"<br> "set gxout stat"<br> tloop=80<br> "set line "ThtwCol" "ThtwLine " "ThtwThk<br> While (tloop > -80)<br> PTemp=LiftWet(tloop,1000,ThtwStop,1,_pmin,_pmax) <br> tloop=tloop-thetwint<br> Endwhile<br>Endif<br><br><br>*-----------------------------------------------------<br>* Plot transformed user-specified temperature sounding<br>*-----------------------------------------------------<br><br>If (DrawTSnd = 1)<br> say "Drawing temperature sounding."<br> "set gxout line"<br> "set x "_xval<br> "set y "_yval<br> "set z "_zmin" "_zmax <br> "set ccolor "TSndCol<br> "set cstyle "TSndLine<br> "set cmark "TSndMrk<br> "set digsiz "MrkSize<br> "set cthick "TSndThk <br> "set missconn on"<br> "d tempx"<br>Endif<br><br>*---------------------------------------------------<br>* Plot transformed user-specified dewpoint sounding<br>*---------------------------------------------------<br><br>If (DrawDSnd = 1)<br> say "Drawing dewpoint sounding."<br> "set gxout line"<br> "set x "_xval<br> "set y "_yval<br> "set z "_zmin" "_zmax<br> "set cmark "DSndMrk <br> "set digsiz "MrkSize<br> "set ccolor "DSndCol<br> "set cstyle "DSndLine<br> "set cthick "DSndThk<br> "set missconn on"<br> "d dewpx"<br>Endif<br><br>*----------------------------------------<br>* Determine lowest level of reported data<br>*----------------------------------------<br><br>If (DrawTSnd = 1) <br> field=sndtemp<br>Else<br> field=sndspd<br>Endif<br><br>"set gxout stat"<br>"set x "_xval<br>"set y "_yval<br>"set t "_tval<br>"set lev " _pmax " " _pmin<br>"d maskout(lev,"field"+300)"<br>rec=sublin(result,1)<br>check=substr(rec,1,6)<br>If (check = "Notice") <br> rec=sublin(result,9)<br>Else<br> rec=sublin(result,8)<br>Endif<br>SfcPlev=subwrd(rec,5)<br><br>If (DrawTSnd = 1 & DrawDSnd = 1)<br> "set lev "SfcPlev<br> "d "sndtemp<br> rec=sublin(result,8)<br> Sfctemp=subwrd(rec,4)<br> "d "snddewp<br> rec=sublin(result,8)<br> Sfcdewp=subwrd(rec,4)<br> SfcThee=Thetae(Sfctemp,Sfcdewp,SfcPlev)<br><br>*------------------------------------------<br>* Calculate temperature and pressure of LCL<br>*------------------------------------------<br><br> TLcl=Templcl(Sfctemp,Sfcdewp)<br> PLcl=Preslcl(Sfctemp,Sfcdewp,SfcPlev)<br>Endif<br><br>*----------------------------------------------------------<br>* Plot parcel path from surface to LCL and up moist adiabat<br>*----------------------------------------------------------<br><br>If (DrawPrcl = 1)<br> say "Drawing parcel path from surface upward."<br> If (PageType = "Portrait")<br> xloc=7.15<br> Else<br> xloc=6.65<br> Endif<br> "q w2xy 1 "PLcl<br> rec=sublin(result,1)<br> yloc=subwrd(rec,6)<br> "set strsiz 0.1"<br> If (PLcl < _pmax) <br> "set string 1 l"<br> "draw string "xloc" "yloc" LCL"<br> "set line 1 1 1"<br> "draw line "xloc-0.15" "yloc" "xloc-0.05" "yloc<br> Endif<br> "set lon 0 100"<br> "set gxout stat"<br> "set line "PrclCol" "PrclLine " " PrclThk<br> PTemp=LiftDry(Sfctemp,SfcPlev,PLcl,1,_pmin,_pmax)<br> Ptemp=LiftWet(TLcl,PLcl,_pmin,1,_pmin,_pmax)<br>Endif<br><br>*-------------------------------------------------------<br>* Determine level within lowest 250mb having highest<br>* theta-e value<br>*-------------------------------------------------------<br><br>If (DrawTSnd = 1 & DrawDSnd = 1)<br> "set x "_xval<br> "set y "_yval<br> "set t "_tval<br> zz=1<br> MaxThee=-999<br> "set gxout stat"<br> while (zz <= _zmax & pp > _pmax-250)<br> "set z "zz<br> pp=subwrd(result,4)<br> "d "sndtemp<br> rec=sublin(result,8)<br> tt=subwrd(rec,4)<br> "d "snddewp<br> rec=sublin(result,8)<br> dd=subwrd(rec,4)<br> If (abs(tt) < 130 & abs(dd) < 130) <br> Thee=Thetae(tt,dd,pp)<br> If (Thee > MaxThee) <br> MaxThee=Thee<br> TMaxThee=tt<br> DMaxThee=dd <br> PMaxThee=pp<br> Endif<br> endif<br> zz=zz+1<br> Endwhile<br> If (PMaxThee = SfcPlev-250) <br> PMaxThee = SfcPlev<br> Endif<br>*------------------------------------------------------<br>* Calculate temperature and pressure of LCL from highest <br>* theta-e level<br>*------------------------------------------------------<br> If (SfcPlev != PMaxThee) <br> TLclMax=Templcl(TMaxThee,DMaxThee)<br> PLclMax=Preslcl(TMaxThee,DMaxThee,PMaxThee)<br> Endif<br>Endif<br><br>*----------------------------------------------------------<br>* Plot parcel path from highest theta-e level to LCL and up<br>* moist adiabat<br>*----------------------------------------------------------<br><br>If (DrawPMax = 1 & SfcPlev != PMaxThee)<br> say "Drawing parcel path from most unstable level upward."<br> If (PageType = "Portrait")<br> xloc=7.15<br> Else<br> xloc=6.65<br> Endif<br> "q w2xy 1 "PLclMax<br> rec=sublin(result,1)<br> yloc=subwrd(rec,6)<br> "set strsiz 0.1"<br> If (PLclMax < _pmax) <br> "set string 1 l"<br> "draw string "xloc" "yloc" LCL"<br> "set line 1 1 1"<br> "draw line "xloc-0.15" "yloc" "xloc-0.05" "yloc<br> Endif<br> "set lon 0 100"<br> "set gxout stat"<br> "set line "PrclCol" "PrclLine " " PrclThk<br> PTemp=LiftDry(TMaxThee,PMaxThee,PLclMax,1,_pmin,_pmax)<br> Ptemp=LiftWet(TLclMax,PLclMax,_pmin,1,_pmin,_pmax)<br>Endif<br><br>*--------------------------------<br>* Draw thermodynamic indices<br>*--------------------------------<br><br>If (DrawIndx = 1) <br> "set string 1 l"<br> "set strsiz 0.10"<br> "set x "_xval<br> "set y "_yval<br> "set t "_tval<br> say "Calculating precipitable water."<br> pw=precipw(sndtemp,snddewp,_pmax,_pmin)<br> say "Calculating thermodynamic indices."<br> Temp850=interp(sndtemp,850)<br> Temp700=interp(sndtemp,700)<br> Temp500=interp(sndtemp,500)<br> Dewp850=interp(snddewp,850)<br> Dewp700=interp(snddewp,700) <br> Dewp500=interp(snddewp,500)<br> If (Temp850>-900 & Dewp850>-900 & Dewp700>-900 & Temp700>-900 & Temp500>-900) <br> K=Temp850+Dewp850+Dewp700-Temp700-Temp500<br> Else <br> K=-999<br> Endif<br> If (Temp850 > -900 & Dewp850 > -900 & Temp500 > -900)<br> tt=Temp850+Dewp850-2*Temp500<br> Else<br> tt=-999<br> Endif<br> Temp500V=virtual2(Temp500+273.15,Dewp500+273.15,500)-273.15<br> PclTemp=LiftWet(TLcl,PLcl,500,0)<br> PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,500)-273.15<br> SLI=Temp500V-PclTempV<br> rec=CAPE(TLcl,PLcl,200,sndtemp,snddewp)<br> Pos=subwrd(rec,1)<br> CIN=subwrd(rec,2)<br><br> If (SfcPlev != PMaxThee) <br> PclTemp=LiftWet(TLclMax,PLclMax,500,0)<br> PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,500)-273.15<br> LIMax=Temp500V-PclTempV<br> rec=CAPE(TLclMax,PLclMax,100,sndtemp,snddewp)<br> PosMax=subwrd(rec,1)<br> CINMax=subwrd(rec,2)<br> Else<br> LIMax=SLI<br> PosMax=Pos<br> CINMax=CIN<br> MaxThee=SfcThee<br> Endif<br> <br> If (PageType = "Portrait") <br> If (Text1XC = -1)<br> Text1XC=rxloc-0.75<br> Endif<br> If (Text1YC = -1)<br> Text1YC=tyloc-2.25<br> Endif<br> If (Text2XC = -1)<br> Text2XC=rxloc-0.75<br> Endif<br> If (Text2YC = -1)<br> Text2YC=tyloc-3.25<br> Endif<br> If (Text3XC = -1)<br> Text3XC=rxloc-0.75<br> Endif<br> If (Text3YC = -1)<br> Text3YC=tyloc-4.40<br> Endif<br> Else<br> If (Text1XC = -1)<br> Text1XC=rxloc+2.50<br> Endif<br> If (Text1YC = -1)<br> Text1YC=tyloc-3.00<br> Endif<br> If (Text2XC = -1)<br> Text2XC=rxloc+2.50<br> Endif<br> If (Text2YC = -1)<br> Text2YC=tyloc-4.00<br> Endif<br> If (Text3XC = -1)<br> Text3XC=rxloc+2.50<br> Endif<br> If (Text3YC = -1)<br> Text3YC=tyloc-5.10<br> Endif<br> Endif<br> "set string 1 l 3"<br> "set line 0 1 3"<br> "draw recf "Text1XC-0.75 " " Text1YC-0.40 " " Text1XC+0.75 " " Text1YC+0.25<br> "set line 1 1 3"<br> "draw rec "Text1XC-0.75 " " Text1YC-0.40 " " Text1XC+0.75 " " Text1YC+0.25<br> "draw string "Text1XC-0.70 " " Text1YC+0.10" K" <br> "draw string "Text1XC+0.25 " " Text1YC+0.10" " int(K) <br> "draw string "Text1XC-0.70 " " Text1YC-0.10 " TT" <br> "draw string "Text1XC+0.25 " " Text1YC-0.10 " " int(tt)<br> "draw string "Text1XC-0.70 " " Text1YC-0.25 " AP(cm)" <br> "draw string "Text1XC+0.25 " " Text1YC-0.25 " " int(pw*100)/100<br> "set line 0 1 3"<br> "draw recf "Text2XC-0.75 " " Text2YC-0.60 " " Text2XC+0.75 " " Text2YC+0.60<br> "set line 1 1 3"<br> "draw rec "Text2XC-0.75 " " Text2YC-0.60 " " Text2XC+0.75 " " Text2YC+0.60<br> "draw string "Text2XC-0.35 " " Text2YC+0.50 " Superficie"<br> "draw string "Text2XC-0.70 " " Text2YC+0.30 " Temp(`3.`0C)" <br> "draw string "Text2XC+0.25 " " Text2YC+0.30 " " int(Sfctemp*10)/10<br> "draw string "Text2XC-0.70 " " Text2YC+0.15 " Dewp(`3.`0C)" <br> "draw string "Text2XC+0.25 " " Text2YC+0.15 " " int(Sfcdewp*10)/10<br> "draw string "Text2XC-0.70 " " Text2YC " `3z`0`bE`n(K)"<br> "draw string "Text2XC+0.25 " " Text2YC " " int(SfcThee) <br> "draw string "Text2XC-0.70 " " Text2YC-0.15 " LI"<br> "draw string "Text2XC+0.25 " " Text2YC-0.15 " " round(SLI)<br> "draw string "Text2XC-0.70 " " Text2YC-0.30 " CAPE(J)"<br> "draw string "Text2XC+0.25 " " Text2YC-0.30 " " int(Pos) <br> "draw string "Text2XC-0.70 " " Text2YC-0.45 " CIN(J)"<br> "draw string "Text2XC+0.25 " " Text2YC-0.45 " " int(CIN) <br> "set line 0 1 3"<br> "draw recf "Text3XC-0.75 " " Text3YC-0.55 " " Text3XC+0.75 " " Text3YC+0.55<br> "set line 1 1 3"<br> "draw rec "Text3XC-0.75 " " Text3YC-0.55 " " Text3XC+0.75 " " Text3YC+0.55<br> "draw string "Text3XC-0.60 " " Text3YC+0.45 " Mas inestable"<br> "draw string "Text3XC-0.70 " " Text3YC+0.20 " Pres(hPa)" <br> "draw string "Text3XC+0.25 " " Text3YC+0.20 " " int(PMaxThee)<br> "draw string "Text3XC-0.70 " " Text3YC+0.05 " `3z`0`bE`n(K)"<br> "draw string "Text3XC+0.25 " " Text3YC+0.05 " " int(MaxThee)<br> "draw string "Text3XC-0.70 " " Text3YC-0.10 " LI" <br> "draw string "Text3XC+0.25 " " Text3YC-0.10 " "round(LIMax)<br> "draw string "Text3XC-0.70 " " Text3YC-0.25 " CAPE(J)" <br> "draw string "Text3XC+0.25 " " Text3YC-0.25 " "int(PosMax) <br> "draw string "Text3XC-0.70 " " Text3YC-0.40 " CIN(J)"<br> "draw string "Text3XC+0.25 " " Text3YC-0.40 " " int(CINMax) <br>Endif<br><br>*-----------------------------<br>* Draw wind profile along side<br>*-----------------------------<br><br>If (DrawBarb = 1) <br> say "Drawing Wind Profile."<br> If (poleloc = -1) <br> If (PageType = "Portrait")<br> poleloc = 8.0<br> Else<br> poleloc = 7.5 <br> Endif<br> Endif<br> If (barbline = 1)<br> "set line 1 1 3"<br> "draw line "poleloc " " byloc " " poleloc " " tyloc<br> Endif<br> If (BarbCol = -1) <br> 'set rgb 41 255 0 132'<br> 'set rgb 42 255 0 168'<br> 'set rgb 43 255 0 204'<br> 'set rgb 44 255 0 240'<br> 'set rgb 45 255 0 255'<br> 'set rgb 46 204 0 255'<br> 'set rgb 47 174 0 255'<br> 'set rgb 48 138 0 255'<br> 'set rgb 49 108 0 255'<br> 'set rgb 50 84 0 255'<br> 'set rgb 51 40 0 255'<br> 'set rgb 52 0 0 255'<br> 'set rgb 53 0 42 255'<br> 'set rgb 54 0 84 255'<br> 'set rgb 55 0 120 255'<br> 'set rgb 56 0 150 255'<br> 'set rgb 57 0 192 255'<br> 'set rgb 58 0 240 255'<br> 'set rgb 59 0 255 210'<br> 'set rgb 60 0 255 160' <br> 'set rgb 61 0 255 126'<br> 'set rgb 62 0 255 78'<br> 'set rgb 63 84 255 0'<br> 'set rgb 64 138 255 0'<br> 'set rgb 65 188 255 0'<br> 'set rgb 66 236 255 0'<br> 'set rgb 67 255 255 0'<br> 'set rgb 68 255 222 0'<br> 'set rgb 69 255 192 0'<br> 'set rgb 70 255 162 0'<br> 'set rgb 71 255 138 0'<br> 'set rgb 72 255 108 0'<br> 'set rgb 73 255 84 0'<br> 'set rgb 74 255 54 0'<br> 'set rgb 75 255 12 0'<br> 'set rgb 76 255 0 34'<br> 'set rgb 77 255 0 70'<br> 'set rgb 78 255 0 105'<br> 'set rgb 79 255 0 140'<br> 'set rgb 80 255 0 175'<br> 'set rgb 81 255 0 215'<br> 'set rgb 82 255 0 255'<br> 'set rgb 83 255 255 255'<br><br> col1='83 83 83 83 83 83 83 83 83 83 82 81 80 79 78'<br> col2='77 76 75 74 73 72 71 70 69 68 67 66 65 64 63'<br> col3='62 61 60 59 58 57 56 55 54 53 52 51 50 49 48'<br> 'set rbcols 'col1' 'col2' 'col3<br> Endif<br> "set z "_zmin" "_zmax<br> "set gxout stat"<br> zz=1<br> wspd=-999<br> cont=1<br> While (cont = 1 & zz < _zmax) <br> "set z "zz<br> pres=subwrd(result,4)<br> "d "u<br> rec=sublin(result,8)<br> uwnd=subwrd(rec,4)<br><br> "d "v<br> rec=sublin(result,8)<br> vwnd=subwrd(rec,4)<br><br> wspd=GetWSpd(uwnd,vwnd)<br><br> if (wspd < 0 | pres > _pmax) <br> zz=zz+1<br> else<br> cont=0<br> Endif<br> Endwhile<br> While (zz <= _zmax)<br><br> "d "u"(z="zz")"<br> rec=sublin(result,8)<br> uwnd=subwrd(rec,4)<br><br> "d "v"(z="zz")"<br> rec=sublin(result,8)<br> vwnd=subwrd(rec,4)<br><br> wspd=GetWSpd(uwnd,vwnd)<br><br> If (BarbCol >= 0)<br> "set line "BarbCol " 1 "BarbThk<br> Else<br> tempbcol=55+wspd/5 <br> If (tempbcol > 83) <br> tempbcol = 83<br> Endif<br> "set line "tempbcol " 1 "BarbThk<br> Endif<br> <br> "d "u"(z="zz")"<br> rec=sublin(result,8)<br> uwnd=subwrd(rec,4)<br><br> "d "v"(z="zz")"<br> rec=sublin(result,8)<br> vwnd=subwrd(rec,4)<br> if(vwnd != undef)<br> wspd=GetWSpd(uwnd,vwnd)<br> wdir=GetWDir(uwnd,vwnd)<br> Else<br> wspd=0<br> wdir=0<br> Endif<br><br> xwind=GetUWnd(wspd,wdir)<br> ywind=GetVWnd(wspd,wdir)<br> "query gr2xy 5 "zz<br> y1=subwrd(result,6) <br> if (wspd > 0) <br> cc=polelen/wspd<br> xendpole=poleloc-xwind*cc<br> yendpole=y1-ywind*cc<br> endif<br> if (xendpole>0 & wspd >= 0.5)<br> if (flagbase = 1) <br> "draw mark 3 "poleloc " " y1 " 0.05"<br> endif<br> "draw line " poleloc " " y1 " " xendpole " " yendpole<br> flagloop=wspd/10<br> windcount=wspd<br> flagcount=0<br> xflagstart=xendpole<br> yflagstart=yendpole<br> dx=cos((180-wdir)*_dtr)<br> dy=sin((180-wdir)*_dtr)<br> while (windcount > 47.5)<br> flagcount=flagcount+1<br> dxflag=-len50*dx<br> dyflag=-len50*dy<br> xflagend=xflagstart+dxflag<br> yflagend=yflagstart+dyflag<br> windcount=windcount-50<br> x1=xflagstart+0.5*wid50*xwind/wspd<br> y1=yflagstart+0.5*wid50*ywind/wspd<br> x2=xflagstart-0.5*wid50*xwind/wspd<br> y2=yflagstart-0.5*wid50*ywind/wspd<br> If (Fill50 = 1) <br> "draw polyf "x1" "y1" "x2" "y2" "xflagend" "yflagend" "x1" "y1<br> Else<br> "draw line "x1 " "y1 " " xflagend " " yflagend " " <br> "draw line "x2 " "y2 " " xflagend " " yflagend<br> "draw line "x1 " "y1 " " x2 " " y2<br> Endif<br> xflagstart=xflagstart+spac50*xwind/wspd<br> yflagstart=yflagstart+spac50*ywind/wspd<br> endwhile<br> while (windcount > 7.5 ) <br> flagcount=flagcount+1<br> dxflag=-len10*dx<br> dyflag=-len10*dy<br> xflagend=xflagstart+dxflag<br> yflagend=yflagstart+dyflag<br> windcount=windcount-10<br> "draw line " xflagstart " " yflagstart " " xflagend " " yflagend<br> xflagstart=xflagstart+spac10*xwind/wspd<br> yflagstart=yflagstart+spac10*ywind/wspd<br> endwhile<br> if (windcount > 2.5) <br> flagcount=flagcount+1<br> if (flagcount = 1) <br> xflagstart=xflagstart+spac05*xwind/wspd<br> yflagstart=yflagstart+spac05*ywind/wspd<br> endif<br> dxflag=-len05*dx<br> dyflag=-len05*dy<br> xflagend=xflagstart+dxflag<br> yflagend=yflagstart+dyflag<br> windcount=windcount-5<br> "draw line " xflagstart " " yflagstart " " xflagend " " yflagend<br> endif<br> else<br> if (wspd < 0.5 & wspd >= 0) <br> "draw mark 2 " poleloc " " y1 " 0.08"<br> endif<br> endif<br> zz=zz+barbint<br> endwhile<br>Endif<br><br>*----------------<br>* Draw Hodograph<br>*----------------<br><br>If (DrawHodo = 1)<br> say "Drawing Hodograph."<br><br> If (HodXcent = -1 | HodYcent = -1) <br> If (PageType = "Portrait") <br> HodXcent=6<br> HodYcent=9.5 <br> rec=sublin(result,8)<br> uwnd=subwrd(rec,4)<br> Else<br> HodXcent=9<br> HodYcent=7.0<br> Endif<br> Endif<br> HodL=HodXcent-HodSize/2.0<br> HodR=HodXcent+HodSize/2.0<br> HodT=HodYcent+HodSize/2.0<br> HodB=HodYcent-HodSize/2.0<br> RingSpac=HodSize/(NumRing*2)<br> "set line 0"<br> "draw recf "HodL" "HodB" "HodR" "HodT<br> "set line "HodoCol" 1 6"<br> "draw rec "HodL" "HodB" "HodR" "HodT<br> "set line 1 1 3"<br> "set string 1 c"<br> "draw mark 1 "HodXcent " "HodYcent " " HodSize<br> i=1<br> While (i <= NumRing)<br> "set strsiz 0.10"<br> "set string 1 c 3 45"<br> uwnd=-i*HodRing*cos(45*_dtr)<br> xloc=HodXcent+uwnd*RingSpac/HodRing<br> yloc=HodYcent+uwnd*RingSpac/HodRing<br> <br> "draw mark 2 "HodXcent " " HodYcent " " i*HodSize/NumRing<br> "draw string "xloc " " yloc " " HodRing*i<br> i=i+1<br> Endwhile<br> "set string 1 l 3 0"<br> If (TickInt > 0) <br> i=0<br> while (i < HodRing*NumRing) <br> dist=i*RingSpac/HodRing<br> hrxloc=HodXcent+dist <br> hlxloc=HodXcent-dist <br> htyloc=HodYcent+dist<br> hbyloc=HodYcent-dist<br> "set line 1 1 3"<br> "draw line "hrxloc " " HodYcent-TickSize/2 " " hrxloc " " HodYcent+TickSize/2<br> "draw line "hlxloc " " HodYcent-TickSize/2 " " hlxloc " " HodYcent+TickSize/2<br> "draw line "HodXcent+TickSize/2 " " htyloc " " HodXcent-TickSize/2 " " htyloc<br> "draw line "HodXcent+TickSize/2 " " hbyloc " " HodXcent-TickSize/2 " " hbyloc<br> i=i+TickInt<br> endwhile<br> Endif<br> "set line "HodoCol " " HodoLine " "HodoThk<br> "draw string "HodL+0.05 " " HodT-0.1 " m/s"<br> zloop=_zmin<br> xold=-999<br> yold=-999<br> count=0<br> Depth=0<br> While (zloop < _zmax & Depth < HodoDep)<br> "set z "zloop<br> pres=subwrd(result,4)<br> <br> "d "u<br> rec=sublin(result,8)<br> uwnd=subwrd(rec,4)<br> "d "v<br> rec=sublin(result,8)<br> vwnd=subwrd(rec,4)<br> if (uwnd = undef)<br> uwnd=0<br> vwnd=0<br> endif<br> <br> If (wspd >= 0) <br> xloc=HodXcent+uwnd*RingSpac/HodRing<br> yloc=HodYcent+vwnd*RingSpac/HodRing<br> If (xloc > 0 & yloc > 0 & xold > 0 & yold > 0) <br> Depth=Depth+pold-pres<br> count=count+1<br> If (count = 1) <br> "draw mark 3 "xold " " yold " 0.05"<br> Endif<br> "draw line "xold" "yold" "xloc" "yloc<br> Endif<br> xold=xloc<br> yold=yloc<br> Endif<br> zloop=zloop+1<br> pold=pres<br> EndWhile<br><br> If (count > 0) <br> "draw mark 3 "xold " " yold " 0.05"<br> Endif<br>Endif<br><br>*----------------------------------------------<br>* Calculate and Display Absolute & S-R Helicity<br>*----------------------------------------------<br><br>If (DrawHeli = 1) <br> say "Calculating Helicity & SR Helicity."<br> delp=10<br> UTotal=0<br> VTotal=0<br><br>* First, calculate mass-weighted mean wind <br>* Since delp is a constant, and mass is proportional to <br>* delp, this is a simple sum.<br><br> "set lev "_pmax " " _pmin<br> "define uwndarr="sndspd"*cos((270-"snddir")*"_dtr")"<br> "define vwndarr="sndspd"*sin((270-"snddir")*"_dtr")"<br><br> pres=MeanVBot<br> While (pres >= MeanVTop)<br> uwnd=interp(uwndarr,pres)*_ktm<br> vwnd=interp(vwndarr,pres)*_ktm<br> If (uwnd > -900 & vwnd > -900) <br> UTotal=UTotal+uwnd<br> VTotal=VTotal+vwnd<br> Endif<br> pres=pres-delp<br> EndWhile<br> vcount=1+(MeanVBot-MeanVTop)/delp <br> Umean=UTotal/vcount<br> Vmean=VTotal/vcount<br> Spdmean=GetWSpd(Umean,Vmean)<br> MeanDir=GetWDir(Umean,Vmean)<br><br>* Now, rotate and reduce mean wind to get storm motion<br><br> If (StormMot = 1) <br> If (Spdmean < 15) <br> Reduct=0.25<br> Rotate=30<br> Else<br> Reduct=0.20<br> Rotate=20<br> Endif<br> Else<br> Reduct=0.0<br> Rotate=0.0<br> Endif<br><br> UReduce=(1-Reduct)*Umean<br> VReduce=(1-Reduct)*Vmean<br> StormSpd=GetWSpd(UReduce,VReduce)<br><br> StormDir=GetWDir(UReduce,VReduce)+Rotate<br> If (StormDir >= 360) <br> StormDir=StormDir-360<br> Endif<br><br> StormU=GetUWnd(StormSpd,StormDir)<br> StormV=GetVWnd(StormSpd,StormDir)<br><br>* Draw Storm Motion Vector<br><br> xloc=HodXcent+_mtk*StormU*RingSpac/HodRing<br> yloc=HodYcent+_mtk*StormV*RingSpac/HodRing<br><br> "set line 1 1 4"<br> "draw line "HodXcent " " HodYcent " " xloc " " yloc<br> Arr1U=GetUWnd(HodRing/10,StormDir+30)<br> Arr1V=GetVWnd(HodRing/10,StormDir+30)<br> Arr2U=GetUWnd(HodRing/10,StormDir-30)<br> Arr2V=GetVWnd(HodRing/10,StormDir-30)<br><br> xloc2=xloc-Arr1U/HodRing<br> xloc3=xloc-Arr2U/HodRing<br> yloc2=yloc-Arr1V/HodRing<br> yloc3=yloc-Arr2V/HodRing<br><br> "set line 1 1 3"<br><br> If (FillArrw = 0) <br> "draw line "xloc" "yloc" "xloc2" "yloc2<br> "draw line "xloc" "yloc" "xloc3" "yloc3<br> Else<br> "draw polyf "xloc" "yloc" "xloc2" "yloc2" "xloc3" "yloc3" "xloc" "yloc<br> Endif<br><br> <br>* Now, calculate SR and Environmental Helicity<br> <br> helic=0<br> SRhelic=0<br> MinP=SfcPlev-HelicDep<br> pres=SfcPlev<br> uwndold=-999<br> vwndold=-999<br> While (pres >= MinP)<br> uwnd=interp(uwndarr,pres)*_ktm<br> vwnd=interp(vwndarr,pres)*_ktm<br> If (uwnd > -900 & uwndold > -900)<br> du=uwnd-uwndold<br> dv=vwnd-vwndold<br> ubar=0.5*(uwnd+uwndold)<br> vbar=0.5*(vwnd+vwndold)<br> uhelic=-dv*ubar <br> vhelic=du*vbar <br> SRuhelic=-dv*(ubar-StormU)<br> SRvhelic=du*(vbar-StormV)<br> SRhelic=SRhelic+SRuhelic+SRvhelic<br> helic=helic+uhelic+vhelic<br> Endif<br> uwndold=uwnd<br> vwndold=vwnd<br> pres=pres-delp<br> EndWhile<br><br> "set strsiz 0.1"<br> "set string 1 l 3"<br> If (PageType = "Portrait") <br> If (Text4XC = -1)<br> Text4XC=rxloc-0.75<br> Endif<br> If (Text4YC = -1) <br> Text4YC=tyloc-5.45<br> Endif<br> Else<br> If (Text4XC = -1)<br> Text4XC=rxloc+2.50<br> Endif<br> If (Text4YC = -1)<br> Text4YC=tyloc-6.10<br> Endif<br> Endif<br> "set line 0 1 3"<br> "draw recf "Text4XC-0.75 " "Text4YC-0.5 " " Text4XC+0.75 " " Text4YC+0.5<br> "set line 1 1 3"<br> "draw rec "Text4XC-0.75 " "Text4YC-0.5 " " Text4XC+0.75 " " Text4YC+0.5<br> "draw string "Text4XC-0.45 " " Text4YC+0.40 " Hodograph"<br> "draw string "Text4XC-0.70 " " Text4YC+0.20 " EH"<br> "draw string "Text4XC+0.25 " " Text4YC+0.20 " "int(helic)<br> "draw string "Text4XC-0.70 " " Text4YC+0.05 " SREH"<br> "draw string "Text4XC+0.25 " " Text4YC+0.05 " " int(SRhelic)<br> "draw string "Text4XC-0.70 " " Text4YC-0.20 " StmDir"<br> "draw string "Text4XC+0.25 " " Text4YC-0.20 " " int(StormDir)"`3.`0"<br> "draw string "Text4XC-0.70 " " Text4YC-0.35 " StmSpd(kt)"<br> "draw string "Text4XC+0.25 " " Text4YC-0.35 " " int(_mtk*StormSpd)<br>Endif<br><br>*---------------------------------------------------<br>* Plot RH profile.<br>*---------------------------------------------------<br><br>If (DrawRH = 1)<br> "set z "_zmin" "_zmax<br> "set x "_xval<br> "set y "_yval<br> "set t "_tval<br> "set zlog on"<br> "set gxout line"<br> "set ccolor "RHCol<br> "set cstyle "RHLine<br> "set cmark "RHMrk<br> "set digsiz "MrkSize<br> "set missconn on"<br> "set xlab on"<br> "set frame off"<br> "set vrange 0 350"<br> "set xlpos 0 t"<br> "set xlevs 25 50 75 100"<br> "set grid vertical 5"<br> "define rh=100*exp((17.2694*"snddewp")/("snddewp"+237.3)-(17.2694*"sndtemp")/("sndtemp"+237.3))"<br> "d rh"<br> If (LblAxes = 1)<br> "set string 1 c 3 0"<br> "set strsiz 0.125"<br> If (PageType = "Portrait")<br> "draw string 1.5 10.85 RH (%)"<br> Else<br> "draw string 1.75 8.35 RH (%)"<br> Endif<br> Endif<br>Endif<br><br>*------------------------------------------<br>* Reset environment to original dimensions<br>*------------------------------------------<br><br>"set t "_tval<br>"set x "_xval <br>"set y "_yval <br>"set z "_zmin " "_zmax<br><br>say "Done."<br><br>*Return(0)<br><br>*************************************************************************<br>function Templcl(temp,dewp)<br><br>*------------------------------------------------------<br>* Calculate the temp at the LCL given temp & dewp in C<br>*------------------------------------------------------<br><br>tempk=temp+273.15<br>dewpk=dewp+273.15<br>Parta=1/(dewpk-56)<br>Partb=log(tempk/dewpk)/800<br>Tlcl=1/(Parta+Partb)+56<br>return(Tlcl-273.15)<br><br>**************************************************************************<br><br>function Preslcl(temp,dewp,pres)<br><br>*-------------------------------------------------------<br>* Calculate press of LCL given temp & dewp in C and pressure<br>*-------------------------------------------------------<br><br>Tlclk=Templcl(temp,dewp)+273.15<br>tempk=temp+273.15<br>theta=tempk*pow(1000/pres,0.286)<br>plcl=1000*pow(Tlclk/theta,3.48)<br>return(plcl)<br><br>**************************************************************************<br>function LiftWet(startt,startp,endp,display,Pmin,Pmax)<br><br>*--------------------------------------------------------------------<br>* Lift a parcel moist adiabatically from startp to endp.<br>* Init temp is startt in C. If you wish to see the parcel's<br>* path plotted, display should be 1. Returns temp of parcel at endp.<br>*--------------------------------------------------------------------<br><br>temp=startt<br>pres=startp<br>cont = 1<br>delp=10<br>While (pres >= endp & cont = 1) <br> If (display = 1) <br> xtemp=GetXLoc(temp,pres)<br> "q w2xy "xtemp" "pres<br> xloc=subwrd(result,3)<br> yloc=subwrd(result,6)<br> If (xtemp < 0 | xtemp > 100)<br> cont=0<br> Else<br> If (pres >= Pmin & pres < Pmax & pres < startp) <br> "draw line "xold" "yold" "xloc" "yloc <br> Endif<br> Endif<br> Endif<br> xold=xloc<br> yold=yloc<br> temp=temp-100*delp*gammaw(temp,pres-delp/2,100)<br> pres=pres-delp<br>EndWhile<br>return(temp)<br><br><br>**************************************************************************<br>function LiftDry(startt,startp,endp,display,Pmin,Pmax)<br><br>*--------------------------------------------------------------------<br>* Lift a parcel dry adiabatically from startp to endp.<br>* Init temp is startt in C. If you wish to see the parcel's<br>* path plotted, display should be 1. Returns temp of parcel at endp.<br>*--------------------------------------------------------------------<br><br>starttk=startt+273.15<br>cont = 1<br>delp=10<br>round=int(startp/10)*10<br>subscr=0.1*round<br>powstart=pow(startp,-0.286)<br>temp=starttk*_powpres.subscr*powstart-273.15<br>pres=round-10<br>While (pres >= endp & cont = 1) <br> subscr=0.1*pres<br> temp=starttk*_powpres.subscr*powstart-273.15<br> If (display = 1) <br> xtemp=GetXLoc(temp,pres)<br> "q w2xy "xtemp" "pres<br> xloc=subwrd(result,3)<br> yloc=subwrd(result,6)<br> If (xtempold > 0 & xtempold < 100 & xtemp > 0 & xtemp < 100) <br> If (pres >= Pmin & pres < Pmax & pres < startp) <br> "draw line "xold" "yold" "xloc" "yloc <br> Endif<br> Endif<br> Endif<br> xold=xloc<br> xtempold=xtemp<br> yold=yloc<br> pres=pres-delp<br>EndWhile<br>return(temp)<br><br>**************************************************************************<br>function CAPE(startt,startp,endp,sndtemp,snddewp)<br><br>*---------------------------------------------------------------------<br>* Returns all postive area and convective inhibition above LCL.<br>* Parcel is lifted from LCL at startt,startp and is halted<br>* at endp. Integration method used is trapezoid method.<br>*---------------------------------------------------------------------<br><br>pres=startp<br>PclTemp=startt<br>PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15<br>delp=10<br>Pos=0<br>Neg=0<br>Neg2=0<br><br>count=0<br>While (pres >= endp)<br> EnvTemp=interp(sndtemp,pres)<br> EnvDewp=interp(snddewp,pres)<br> EnvTempV=virtual2(EnvTemp+273.15,EnvDewp+273.15,pres)-273.15<br> DelT=PclTempV-EnvTempV<br> If (abs(EnvTempV) < 130 & abs(PclTempV) < 130)<br> count=count+1<br> If (count > 1) <br> Val=DelT/pres+Prev <br> If (Val > 0)<br> Pos=Pos+Val<br> Neg2=0<br> Else<br> Neg=Neg+abs(Val)<br> Neg2=Neg2+abs(Val)<br> Endif<br> Endif<br> Prev=DelT/pres<br> Endif<br> pres=pres-delp<br> PclTemp=PclTemp-100*delp*gammaw(PclTemp,pres,100)<br> PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15<br>Endwhile<br><br>Pos=0.5*Pos*287*delp<br>CIN=0.5*(Neg-Neg2)*287*delp<br><br>return(Pos" "CIN)<br><br>***************************************************************************<br>function gammaw(tempc,pres,rh)<br><br>*-----------------------------------------------------------------------<br>* Function to calculate the moist adiabatic lapse rate (deg C/Pa) based<br>* on the temperature, pressure, and rh of the environment.<br>*----------------------------------------------------------------------<br><br>tempk=tempc+273.15<br>es=satvap2(tempc)<br>ws=mixratio(es,pres)<br>w=rh*ws/100<br>tempv=virtual(tempk,w)<br>latent=latentc(tempc)<br><br>A=1.0+latent*ws/(287*tempk)<br>B=1.0+0.622*latent*latent*ws/(1005*287*tempk*tempk)<br>Density=100*pres/(287*tempv)<br>lapse=(A/B)/(1005*Density)<br>return(lapse)<br><br>*************************************************************************<br>function latentc(tempc)<br><br>*-----------------------------------------------------------------------<br>* Function to return the latent heat of condensation in J/kg given<br>* temperature in degrees Celsius.<br>*-----------------------------------------------------------------------<br><br>val=2502.2-2.43089*tempc<br><br>return(val*1000)<br><br>*************************************************************************<br>function precipw(sndtemp,snddewp,startp,endp)<br><br>*-----------------------------------------------------------------------<br>* Function to calculate the precipitable water (cm) in a sounding<br>* starting at pressure level startp and ending at pressure level endp.<br>*-----------------------------------------------------------------------<br><br>ppold=-999<br>ttold=-999<br>ddold=-999<br>delp=10<br>Int=0<br>mix=0 <br>pres=startp<br>logpp=log(pres)<br>logppm=log(pres-delp)<br>while (pres >= endp)<br> tt=interp(sndtemp,pres)<br> dd=interp(snddewp,pres)<br> if (tt>-900 & ttold>-900 & dd>-900 & ddold>-900) <br> e=satvap2(dd)<br> mix=mixratio(e,pres)<br> mixavg=(logpp*mix+logppm*mixold)/(logpp+logppm)<br> Int=Int+1.020408*mixavg*delp<br> endif<br> ttold=tt<br> ddold=dd<br> ppold=pp<br> mixold=mix<br> pres=pres-delp<br> logpp=logppm<br> logppm=log(pres-delp)<br>endwhile<br><br>return(Int)<br><br>*************************************************************************<br><br>function virtual(temp,mix)<br><br>*------------------------------------------------------------<br>* Function to return virtual temperature given temperature in <br>* kelvin and mixing ratio in g/g.<br>*-------------------------------------------------------------<br><br>tempv=temp*(1.0+0.6*mix)<br><br>return (tempv)<br><br>************************************************************************<br><br>function virtual2(temp,dewp,pres)<br> <br>*------------------------------------------------------------<br>* Function to return virtual temperature in kelvin given temperature in<br>* kelvin and dewpoint in kelvin and pressure in mb<br>*-------------------------------------------------------------<br><br>if (temp > 0 & dewp > 0) <br> vap=satvap2(dewp-273.15)<br> mix=mixratio(vap,pres)<br> tempv=virtual(temp,mix)<br>else<br> tempv=-9999<br>endif<br><br>return (tempv)<br> <br>************************************************************************<br><br>function satvapor(temp)<br><br>*---------------------------------------------------------------<br>* Given temp in Celsius, returns saturation vapor pressure in mb<br>*---------------------------------------------------------------<br><br>pol=_C0+temp*(_C1+temp*(_C2+temp*(_C3+temp*(_C4+temp*(_C5+temp*(_C6+temp*(_C7+temp*(_C8+temp*(_C9)))))))))<br><br>return(6.1078/pow(pol,8))<br><br>************************************************************************<br><br>function satvap2(temp)<br><br>*---------------------------------------------------------------<br>* Given temp in Celsius, returns saturation vapor pressure in mb<br>*---------------------------------------------------------------<br><br>es=6.112*exp(17.67*temp/(temp+243.5))<br><br>return(es)<br><br>*************************************************************************<br><br>function mixratio(e,p)<br><br>*------------------------------------------------------<br>* Given vapor pressure and pressure, return mixing ratio<br>*-------------------------------------------------------<br><br>mix=0.622*e/(p-e)<br><br>return(mix)<br><br>************************************************************************<br><br>function getrh(temp,dewp,pres)<br><br>tempk=temp+273.15<br>dewpk=dewp+273.15<br><br>es=satvap2(temp)<br><br>If (temp > 0) <br> A=2.53*pow(10,9)<br> B=5420<br>Else<br> A=3.41*pow(10,10)<br> B=6130<br>Endif<br><br>w=A*0.622*exp(-B/dewpk)/pres<br>ws=mixratio(es,pres)<br><br>return(100*w/ws)<br><br>************************************************************************<br>function interp(array,pres)<br><br>*------------------------------------------------------------------------<br>* Interpolate inside array for pressure level pres.<br>* Returns estimated value of array at pressure pres.<br>*------------------------------------------------------------------------<br><br>"set gxout stat"<br>"set lev "pres<br>altpres=subwrd(result,4)<br>"q dims"<br>rec=sublin(result,4)<br>zlev=subwrd(rec,9)<br><br>If (zlev < 2 | zlev > _zmaxfile)<br> Vest = -9999.0<br>Else<br> If (altpres > pres)<br> zlev=zlev+1<br> Endif<br> "set z "zlev<br> PAbove=subwrd(result,4)<br> "d "array"(lev="PAbove")"<br> rec=sublin(result,8)<br> VAbove=subwrd(rec,4)<br> "set z "zlev-1<br> PBelow=subwrd(result,4)<br> "d "array"(lev="PBelow")"<br> rec=sublin(result,8)<br> VBelow=subwrd(rec,4)<br><br>* Now if we are in a region of missing data, find next good level.<br><br> If (abs(VAbove) > 130 & zlev > 1 & zlev < _zmaxfile)<br> zz=zlev<br> While (abs(VAbove) > 130 & zz < _zmaxfile)<br> zz=zz+1<br> "set z "zz<br> PAbove=subwrd(result,4)<br> "d "array"(lev="PAbove")"<br> rec=sublin(result,8)<br> VAbove=subwrd(rec,4)<br> EndWhile<br> Endif<br><br> If (abs(VBelow) > 130 & zlev > 1 & zlev < _zmaxfile)<br> zz=zlev-1<br> While (abs(VBelow) > 130 & zz > 1) <br> zz=zz-1<br> "set z "zz<br> PBelow=subwrd(result,4)<br> "d "array"(lev="PBelow")"<br> rec=sublin(result,8)<br> VBelow=subwrd(rec,4)<br> EndWhile<br> Endif<br><br> If (abs(VAbove) < 130 & abs(VBelow) < 130)<br> Vest=VBelow+log(PBelow/pres)*(VAbove-VBelow)/log(PBelow/PAbove)<br> Else<br> Vest=-9999.0<br> Endif<br><br>Endif<br><br>Return(Vest)<br><br>****************************************************************************<br><br>function GetUWnd(wspd,wdir)<br><br>*------------------------<br>* Get x-component of wind. <br>*------------------------<br><br><br>If (wspd >= 0) <br> xwind=wspd*cos((270-wdir)*_dtr)<br>Else<br> xwind = -9999.0<br>Endif<br>return(xwind)<br><br>**************************************************************************<br><br>function GetVWnd(wspd,wdir)<br><br>*-----------------------<br>* Get y-component of wind<br>*------------------------<br><br>If (wspd >= 0) <br> ywind=wspd*sin((270-wdir)*_dtr)<br>Else<br> ywind = -9999.0<br>Endif<br>return(ywind)<br><br><br>*************************************************************************<br><br>function GetWSpd(xwind,ywind)<br><br><br>"set gxout stat"<br>"d mag("xwind","ywind")"<br>rec=sublin(result,8)<br>val=subwrd(rec,4)<br><br>return (val)<br><br>*************************************************************************<br><br>function GetWDir(xwind,ywind)<br><br>* Return wind direction given x and y components<br><br>"set gxout stat"<br>"define theta=270-"_rtd"*atan2("ywind","xwind")"<br>"d theta"<br>rec=sublin(result,8)<br>Dir=subwrd(rec,4)<br><br>If (Dir > 360)<br> Dir=Dir-360<br>Endif<br><br>If (Dir < 0)<br> Dir=360+Dir<br>Endif<br><br>return(Dir)<br><br>*************************************************************************<br><br>function GetXLoc(temp,pres)<br><br>*-------------------------------------------------<br>* Get x-location on skew-t based on temp, pressure<br>*-------------------------------------------------<br><br>xloc=(temp-_m1*log10(pres)-_m3)/_m2<br>return(xloc)<br><br>*************************************************************************<br> <br>function GetTemp(xloc,pres) <br><br>*------------------------------------------------- <br>* Return temperature at location given by xloc,pres<br>*-------------------------------------------------<br><br>tempval=_m1*log10(pres)+_m2*xloc+_m3<br>return(tempval)<br><br>**************************************************************************<br><br>function GetTheta(temp,pres) <br><br>*---------------------------------------------------<br>* Calculate theta for a given temperature and pressure<br>*---------------------------------------------------<br><br>theta=(temp+273.15)*pow(1000/pres,0.286)-273.15<br>return(theta)<br><br><br>*************************************************************************<br><br>function GetThet2(temp,dewp,pres) <br><br>*---------------------------------------------------<br>* Calculate theta for a given temperature,dewp, and pressure<br>*---------------------------------------------------<br><br>tempk=273.15+temp<br>dewpk=273.15+dewp<br><br>es=satvap2(temp)<br>ws=mixratio(es,pres)<br><br>mix=10*getrh(temp,dewp,pres)*ws<br><br>exponent=0.2854*(1.0-0.00028*mix)<br>theta=(temp+273.15)*pow(1000/pres,exponent)-273.15<br>return(theta)<br><br>*************************************************************************<br><br>function Thetae(temp,dewp,pres)<br><br>*--------------------------------------------------------------<br>* Return equiv. pot. temp in Kelvin given temp, dewp in celsius<br>* From Bolton (1980) Mon Wea Rev<br>*--------------------------------------------------------------<br><br>es=satvap2(temp)<br>ws=mixratio(es,pres)<br>mix=10*getrh(temp,dewp,pres)*ws<br>theta=GetThet2(temp,dewp,pres)+273.15<br>TLcl=Templcl(temp,dewp)+273.15<br>thetae=theta*exp((3.376/TLcl-0.00254)*mix*(1.0+0.00081*mix))<br><br>return(thetae)<br><br>**************************************************************************<br><br><br>function int(i0)<br><br>*--------------------------<br>* Return integer of i0<br>*--------------------------<br> i=0<br> while(i<12)<br> i=i+1<br> if(substr(i0,i,1)='.')<br> i0=substr(i0,1,i-1)<br> break<br> endif<br> endwhile<br>return(i0)<br><br>*************************************************************************<br><br>function abs(i)<br><br>*----------------------------<br>* return absolute value of i<br>*----------------------------<br><br> if (i < 0) <br> absval=-i<br> else <br> absval=i<br> endif<br><br>return(absval)<br><br>*************************************************************************<br><br>function log(i)<br><br>*---------------------------<br>* return natural log of i<br>*---------------------------<br><br>"set gxout stat"<br>"d log("i")"<br>rec=sublin(result,8)<br>val=subwrd(rec,4)<br>return(val)<br><br>*************************************************************************<br><br>function log10(i)<br><br>*--------------------------<br>* return log base 10 of i<br>*--------------------------<br><br>"set gxout stat"<br>"d log10("i")"<br>rec=sublin(result,8)<br>val=subwrd(rec,4)<br>return(val)<br><br>*************************************************************************<br><br>function pow(i,j)<br><br>*-------------------------------<br>* return power of i raised to j<br>*-------------------------------<br><br>"set gxout stat"<br>"d pow("i","j")"<br>rec=sublin(result,8)<br>val=subwrd(rec,4)<br>return(val)<br><br>************************************************************************<br><br>function cos(i)<br><br>*-----------------------------------------<br>* return cosine of i, where i is in radians<br>*------------------------------------------<br><br>"set gxout stat"<br>"d cos("i")"<br>rec=sublin(result,8)<br>val=subwrd(rec,4)<br>return(val)<br><br>************************************************************************<br><br>function sin(i)<br><br>*------------------------------------------<br>* return sine of i, where i is in radians<br>*------------------------------------------<br><br>"set gxout stat"<br>"d sin("i")"<br>rec=sublin(result,8)<br>val=subwrd(rec,4)<br>return(val)<br><br>************************************************************************<br><br>function exp(i)<br><br>*------------------------------------------<br>* return exponential of i<br>*------------------------------------------<br><br>"set gxout stat"<br>"d exp("i")"<br>rec=sublin(result,8)<br>val=subwrd(rec,4)<br>return(val)<br><br>***********************************************************************<br>function round(i)<br><br>rr=abs(1.0*i)<br>rr=int(rr+0.5)<br>if (i < 0)<br> rr=-1*rr <br>endif<br>return(rr)<br><br><br> </body>
</html>