Context Navigation

skew.gs @ 897

Last change on this file since 897 was 207, checked in by aslmd, 14 years ago
MESOSCALE: A GENERAL CLEAN-UP FOLLOWING UPDATING THE USER MANUAL. EVERYTHING ESSENTIAL IS IN MESOSCALE (much lighter than before). EVERYTHING FOR DEVELOPPERS OR EXPERTS IS IN MESOSCALE_DEV.
Property svn:executable set to ``*
File size: 66.1 KB

Rev	Line
[207]	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(uu + vv))*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-2T3)/(2log10(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.zztk+lpp.zzmtkold)/(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.5wid50xwind/wspd
	1512	y1=yflagstart+0.5wid50ywind/wspd
	1513	x2=xflagstart-0.5wid50xwind/wspd
	1514	y2=yflagstart-0.5wid50ywind/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=-iHodRingcos(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+_mtkStormURingSpac/HodRing
	1732	yloc=HodYcent+_mtkStormVRingSpac/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=100exp((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-100delpgammaw(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.subscrpowstart-273.15
	1961	pres=round-10
	1962	While (pres >= endp & cont = 1)
	1963	subscr=0.1*pres
	1964	temp=starttk_powpres.subscrpowstart-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-100delpgammaw(PclTemp,pres,100)
	2022	PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15
	2023	Endwhile
	2024
	2025	Pos=0.5Pos287*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+latentws/(287tempk)
	2046	B=1.0+0.622latentlatentws/(1005287tempktempk)
	2047	Density=100pres/(287tempv)
	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=(logppmix+logppmmixold)/(logpp+logppm)
	2087	Int=Int+1.020408mixavgdelp
	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.6mix)
	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.112exp(17.67temp/(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=A0.622exp(-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=wspdcos((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=wspdsin((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=_m1log10(pres)+_m2xloc+_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=10getrh(temp,dewp,pres)ws
	2409
	2410	exponent=0.2854(1.0-0.00028mix)
	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=10getrh(temp,dewp,pres)ws
	2425	theta=GetThet2(temp,dewp,pres)+273.15
	2426	TLcl=Templcl(temp,dewp)+273.15
	2427	thetae=thetaexp((3.376/TLcl-0.00254)mix1.0+0.00081mix)
	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)

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format