Jeff,<br>Try using a two-digit color number (i.e., 00) in the printim command as per the recommendation on the index page for the printim command: <a href="http://www.iges.org/grads/gadoc/gradcomdprintim.html">http://www.iges.org/grads/gadoc/gradcomdprintim.html</a><br>
<br>Jeff<br><br><div class="gmail_quote">On Mon, Nov 8, 2010 at 1:41 PM, Jeff Chabot <span dir="ltr"><<a href="mailto:jsc219@gmail.com">jsc219@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
I am using two separate color tables. For precip:<br><br>*Categorical Rain<br>set rgb 30 210 180 140<br>set rgb 16 107 159 91<br>set rgb 17 90 142 75<br>set rgb 18 74 125 60<br>set rgb 19 58 109 45<br>set rgb 20 41 92 19<br>
set rgb 21 25 75 14<br>set rgb 22 253 248 2<br>set rgb 23 229 188 0<br>set rgb 24 253 149 0<br>set rgb 25 253 0 0<br>set rgb 26 212 0 0<br>set rgb 27 255 0 153<br><br>set clevs 0.01 0.05 0.10 0.25 0.50 1.00 2.00 3.00 4.00 5.00<br>
set ccols 0 16 18 19 20 22 23 24 25 26 27<br><br>For precip type:<br><br>*Categorical Rain<br>set rgb 16 107 159 91<br>set rgb 17 90 142 75<br>set rgb 18 74 125 60<br>set rgb 19 58 109 45<br>set rgb 20 41 92 19<br>set rgb 21 25 75 14<br>
set rgb 22 253 248 2<br>set rgb 23 229 188 0<br>set rgb 24 253 149 0<br>set rgb 25 253 0 0<br>set rgb 26 212 0 0<br>set rgb 27 255 0 153<br><br>*Categorical Freezing Rain<br>set rgb 28 255 188 188<br>set rgb 29 255 172 172<br>
set rgb 30 255 156 156<br>set rgb 31 253 140 140<br>set rgb 32 255 124 124<br>set rgb 33 255 108 108<br>set rgb 34 255 96 96<br>set rgb 35 255 80 80<br>set rgb 36 255 64 56<br>set rgb 37 240 32 128<br>set rgb 38 240 16 255<br>
<br>*Categorical Ice Pellets / Mix<br>set rgb 39 255 200 0<br>set rgb 40 255 180 0<br>set rgb 41 255 160 0<br>set rgb 42 255 140 0<br>set rgb 43 255 124 0<br>set rgb 44 255 108 0<br>set rgb 45 255 96 0<br>set rgb 46 255 80 0<br>
set rgb 47 255 64 0<br>set rgb 48 255 32 0<br>set rgb 49 255 16 0<br><br>*Categorical Snow<br>set rgb 50 4 233 231<br>set rgb 51 0 173 255<br>set rgb 52 0 148 255<br>set rgb 53 0 123 255<br>set rgb 54 0 104 255<br>set rgb 55 0 85 255<br>
set rgb 56 4 67 245<br>set rgb 57 0 38 255<br>set rgb 58 0 14 255<br>set rgb 59 0 0 255<br>set rgb 60 0 0 223<br>set rgb 61 255 255 255<br><br>set clevs 0.01 0.05 0.1 0.25 0.5 1 2 3 4 5 10 10.01 10.05 10.1 10.25 10.5 11 12<br>
13 14 15 20 20.01 20.05 20.1 20.25 20.5 21 22 23 24 25 30 30.01 30.05 30.1 30.25<br> 30.5 31 32 33 34 35<br>set ccols 0 16 18 19 20 22 23 24 25 26 27 0 39 40 41 42 43 44 45 46 47 0 0 28 29<br> 30 31 32 33 34 35 36 0 0 50 51 52 53 54 55 56 57 58 59<br>
<br>So, Jeff, I believe that I am using #0 correctly. But, when I do, it just uses white or black as the background as in this example:<br><br><a href="http://jeffsweatherservice.com/grads/nam/ptype_NE_2.gif" target="_blank">http://jeffsweatherservice.com/grads/nam/ptype_NE_2.gif</a> <br>
<br>Thanks,<br><br><br>Jeff<br><br><br><br><div class="gmail_quote">On Mon, Nov 8, 2010 at 10:18 AM, <span dir="ltr"><<a href="mailto:gradsusr-request@gradsusr.org" target="_blank">gradsusr-request@gradsusr.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">Send gradsusr mailing list submissions to<br>
<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a><br>
<br>
To subscribe or unsubscribe via the World Wide Web, visit<div class="im"><br>
<a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br></div>
or, via email, send a message with subject or body 'help' to<br>
<a href="mailto:gradsusr-request@gradsusr.org" target="_blank">gradsusr-request@gradsusr.org</a><br>
<br>
You can reach the person managing the list at<br>
<a href="mailto:gradsusr-owner@gradsusr.org" target="_blank">gradsusr-owner@gradsusr.org</a><br>
<br>
When replying, please edit your Subject line so it is more specific<br>
than "Re: Contents of gradsusr digest..."<br>
<br>
<br>
Today's Topics:<br>
<br>
1. GrADS help, plotting shaded contour precip/ptype products<br>
over a basemap/background (Jeff Chabot)<br>
2. Re: GrADS help, plotting shaded contour precip/ptype products<br>
over a basemap/background (Jeffrey Duda)<br>
3. Problem with skewt plot using WRF output (Cristian ...)<br>
<br>
<br>
----------------------------------------------------------------------<br>
<br>
Message: 1<br>
Date: Mon, 8 Nov 2010 13:51:33 -0500<br>
From: Jeff Chabot <<a href="mailto:jsc219@gmail.com" target="_blank">jsc219@gmail.com</a>><br>
Subject: [gradsusr] GrADS help, plotting shaded contour precip/ptype<br>
products over a basemap/background<br>
To: <a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a><br>
Message-ID:<br>
<AANLkTinKn=<a href="mailto:eTr6myE00LdycRWRxdyy5e8hpBtf%2BFgLPz@mail.gmail.com" target="_blank">eTr6myE00LdycRWRxdyy5e8hpBtf+FgLPz@mail.gmail.com</a>><br>
Content-Type: text/plain; charset="iso-8859-1"<div class="im"><br>
<br>
Dear GrADS Users,<br>
I have been trying to apply a basemap to my precip/ptype shaded contour<br>
products on my personal website for some time now. I have brainstormed and<br>
searched the GrADS online documents and have come up with ideas, but they<br>
have never worked. My latest idea was to use Google Maps terrain data as a<br>
background to the map. I thought it was a good idea, but I discovered that<br>
I still cannot resolve the issue in that precip = 0 shows either all white<br>
or all black depending on what I choose for the background, covering up my<br>
basemap. I really thought the following command would have resolved that<br>
issue, but it did not:<br>
<br>
printim precip_NE_1.png png x570 y435 -b NE.png -t 0 Versions:GrADS: 2.0a9<br>
OS: Fedora 13<br>
Sample: <a href="http://jeffsweatherservice.com/grads/nam/ptype_NE_2.png" target="_blank">http://jeffsweatherservice.com/grads/nam/ptype_NE_2.png</a><br>
<br>
I know this can be done because of the example:<br>
<br>
<a href="http://wxmaps.org/pix/ez.east.html" target="_blank">http://wxmaps.org/pix/ez.east.html</a><br>
<br>
Any assistance here would be much appreciated.Sincerely,<br>
<br>
<br>
Jeff Chabot<br>
Email: <a href="mailto:jsc219@gmail.com" target="_blank">jsc219@gmail.com</a><br>
Web: <a href="http://jeffsweatherservice.com" target="_blank">http://jeffsweatherservice.com</a><br></div>
-------------- next part --------------<br>
An HTML attachment was scrubbed...<br>
URL: <a href="http://gradsusr.org/pipermail/gradsusr/attachments/20101108/93f97642/attachment-0001.html" target="_blank">http://gradsusr.org/pipermail/gradsusr/attachments/20101108/93f97642/attachment-0001.html</a><br>
<br>
------------------------------<br>
<br>
Message: 2<br>
Date: Mon, 8 Nov 2010 12:57:30 -0600<br>
From: Jeffrey Duda <<a href="mailto:jdduda@iastate.edu" target="_blank">jdduda@iastate.edu</a>><br>
Subject: Re: [gradsusr] GrADS help, plotting shaded contour<br>
precip/ptype products over a basemap/background<br>
To: GrADS Users Forum <<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a>><br>
Message-ID:<br>
<<a href="mailto:AANLkTikPmH-Og0-6uNAi4tpUMCoZbAZ3YTqm_HXkndpV@mail.gmail.com" target="_blank">AANLkTikPmH-Og0-6uNAi4tpUMCoZbAZ3YTqm_HXkndpV@mail.gmail.com</a>><br>
Content-Type: text/plain; charset="iso-8859-1"<div><div></div><div class="h5"><br>
<br>
Jeff,<br>
What color table are you using? Are you sure you have a #0 color? I know<br>
that's supposed to be the background, but perhaps you are overwriting it<br>
with some other color that has the same color as your background but a<br>
different color number.<br>
<br>
Jeff Duda<br></div></div></blockquote><div><br> </div><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><div><div></div><div class="h5">
<br>
On Mon, Nov 8, 2010 at 12:51 PM, Jeff Chabot <<a href="mailto:jsc219@gmail.com" target="_blank">jsc219@gmail.com</a>> wrote:<br>
<br>
> Dear GrADS Users,<br>
> I have been trying to apply a basemap to my precip/ptype shaded contour<br>
> products on my personal website for some time now. I have brainstormed and<br>
> searched the GrADS online documents and have come up with ideas, but they<br>
> have never worked. My latest idea was to use Google Maps terrain data as a<br>
> background to the map. I thought it was a good idea, but I discovered that<br>
> I still cannot resolve the issue in that precip = 0 shows either all white<br>
> or all black depending on what I choose for the background, covering up my<br>
> basemap. I really thought the following command would have resolved that<br>
> issue, but it did not:<br>
><br>
> printim precip_NE_1.png png x570 y435 -b NE.png -t 0 Versions:GrADS: 2.0a9<br>
> OS: Fedora 13<br>
> Sample: <a href="http://jeffsweatherservice.com/grads/nam/ptype_NE_2.png" target="_blank">http://jeffsweatherservice.com/grads/nam/ptype_NE_2.png</a><br>
><br>
> I know this can be done because of the example:<br>
><br>
> <a href="http://wxmaps.org/pix/ez.east.html" target="_blank">http://wxmaps.org/pix/ez.east.html</a><br>
><br>
> Any assistance here would be much appreciated. Sincerely,<br>
><br>
><br>
> Jeff Chabot<br>
> Email: <a href="mailto:jsc219@gmail.com" target="_blank">jsc219@gmail.com</a><br>
> Web: <a href="http://jeffsweatherservice.com" target="_blank">http://jeffsweatherservice.com</a><br>
><br>
><br>
> _______________________________________________<br>
> gradsusr mailing list<br>
> <a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a><br>
> <a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
><br>
><br>
<br>
<br>
--<br>
Jeff Duda<br>
Iowa State University<br>
Meteorology Graduate Student<br>
3134 Agronomy Hall<br>
<a href="http://www.meteor.iastate.edu/%7Ejdduda" target="_blank">www.meteor.iastate.edu/~jdduda</a><br></div></div>
-------------- next part --------------<br>
An HTML attachment was scrubbed...<br>
URL: <a href="http://gradsusr.org/pipermail/gradsusr/attachments/20101108/159827b2/attachment-0001.html" target="_blank">http://gradsusr.org/pipermail/gradsusr/attachments/20101108/159827b2/attachment-0001.html</a><br>
<br>
------------------------------<br>
<br>
Message: 3<br>
Date: Mon, 8 Nov 2010 16:23:09 -0300<br>
From: "Cristian ..." <<a href="mailto:galenso85@hotmail.com" target="_blank">galenso85@hotmail.com</a>><br>
Subject: [gradsusr] Problem with skewt plot using WRF output<br>
To: <<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a>><br>
Message-ID: <SNT112-W63DF5D09C452E1C6FB44FBF4F0@phx.gbl><br>
Content-Type: text/plain; charset="iso-8859-1"<br>
<br>
<br>
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>
* <a href="http://www.ems.psu.edu/%7Ehart/skew.html" target="_blank">http://www.ems.psu.edu/~hart/skew.html</a><br>
*<br>
* Please send any further problems, comments, or suggestions to<br>
* <<a href="mailto:hart@ems.psu.edu" target="_blank">hart@ems.psu.edu</a>><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>
<br>
-------------- next part --------------<br>
An HTML attachment was scrubbed...<br>
URL: <a href="http://gradsusr.org/pipermail/gradsusr/attachments/20101108/1c1328b3/attachment.html" target="_blank">http://gradsusr.org/pipermail/gradsusr/attachments/20101108/1c1328b3/attachment.html</a><br>
<br>
------------------------------<div class="im"><br>
<br>
_______________________________________________<br>
gradsusr mailing list<br>
<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a><br>
<a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
<br>
<br></div>
End of gradsusr Digest, Vol 9, Issue 13<br>
***************************************<br>
</blockquote></div><br>
<br>_______________________________________________<br>
gradsusr mailing list<br>
<a href="mailto:gradsusr@gradsusr.org">gradsusr@gradsusr.org</a><br>
<a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
<br></blockquote></div><br><br clear="all"><br>-- <br>Jeff Duda<br>Iowa State University<br>Meteorology Graduate Student<br>3134 Agronomy Hall<br><a href="http://www.meteor.iastate.edu/~jdduda">www.meteor.iastate.edu/~jdduda</a><br>