[gradsusr] Problem with skewt plot using WRF output
Jeffrey Duda
jdduda at iastate.edu
Mon Nov 8 15:00:16 EST 2010
Cristian,
A couple of things I noticed.
1) If you are running WRF-ARW, then your u and v wind are in m/s, but the
script assumes your winds will be in kts. Just know that unless you adjust
u and v before calling the main function, the results you get will be in m/s
instead of kts. However, it's possible that this is causing the hodograph
not to be plotted, but I would think an error would appear if that was the
case.
2) I would consider setting HodXcent and HodYcent to -1. Since the script
redefines the plotting window, your choice of 6 and 9 for the center of the
hodograph may not be on the virtual page/plotting area and thus you wouldn't
even see it.
3) Try setting BarbCol to -1 instead of 1
See if any of this helps.
Jeff Duda
On Mon, Nov 8, 2010 at 1:23 PM, Cristian ... <galenso85 at hotmail.com> wrote:
> 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)
>
>
>
> _______________________________________________
> 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/69d616b2/attachment-0003.html
More information about the gradsusr
mailing list