[gradsusr] GrADS help, plotting shaded contour precip/ptype products over a basemap/background

Jeff Chabot jsc219 at gmail.com
Mon Nov 8 15:31:57 EST 2010


Jeff,

I liked the idea, but that didn't help.  I tried the following, both
resulting in a black background for the map area:

printim precip_NE.png png x570 y435 -b /home/jeff/public_html/images/NE.png
-t 00

printim precip_NE.png png x570 y435 -t 00 -b
/home/jeff/public_html/images/NE.png

Example: http://jeffsweatherservice.com/grads/nam/precip_NE_2.png

I even tried using "00" in the color table, but that had the same impact.

Now, if I use the following command, the black background becomes white:

printim precip_NE.png png x570 y435 -b /home/jeff/public_html/images/NE.png
white

Thanks for the idea though.


Jeff








On Mon, Nov 8, 2010 at 10:59 AM, <gradsusr-request at gradsusr.org> wrote:

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


More information about the gradsusr mailing list