[gradsusr] GrADS help, plotting shaded contour precip/ptype products over a basemap/background
Jeff Chabot
jsc219 at gmail.com
Mon Nov 8 14:41:17 EST 2010
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>
> >
> 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>
> -------------- 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>
> *
> * 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
> ***************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20101108/0b1712f1/attachment-0001.html
More information about the gradsusr
mailing list