function plotskew(sndtemp,snddewp,sndspd,snddir) ************************************************************************* * * GrADS Script to Plot a SkewT/LogP Diagram * * Bob Hart * Penn State University / Dept of Meteorology * Last Update: March 3, 1999 * * Recent Changes: * * 3/3/99 - Made latent heat of vaporization constant a function * of temperature. This improves CAPE calculations * further. Also, bug causing false values of CAPE, CIN * when minimal dewpoint data was present is fixed. * * 2/17/99 - Improved the precision of the precip. water & helicity * calculations by changing them to use the interp function. * Significant accuracy improvement when # of vertical levels * in your dataset is small (< 15). * * 2/16/99 - Fixed a long-lived error in the CAPE calculation that * produced 10% underestimation errors. Also CAPE delp=5 now as * accuracy increase outweighed resulting slight speed slowdown. * * Features: * - All features of standard skewt/logp plot * - 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 * * Please send any further problems, comments, or suggestions to * * * ACKNOWLEDGMENTS: Thanks go to the innumerable users who have helped * fine tune the script from the horrible mess it began as. 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 = -50 * ------------------- 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= 10 thetwint= 5 tempint = 10 wsclevs = "1 1.5 3 4 6 8 10 12 16 20 24 28 32 36 40" * * * ------------------------ 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 * 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 DrawPrcl= 1 DrawPMax= 1 DrawIndx= 1 DrawHeli= 1 DrawHodo= 1 DrawPLev= 1 DrawZLev= 0 DrawZSTD= 0 LblAxes = 1 ThtwStop = 200 MixStop = 600 * * ----------------- 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 barb/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= -1 HodYcent= -1 HodSize = 2 NumRing = 3 HodRing = 15 HodoDep = 300 TickInt = 5 TickSize= 0.05 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= 1 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 * PrclCol = Color of parcel trace * BarbCol = Color of wind barbs (choose -1 for color according to speed) * HodoCol = Color of hodograph trace ThetCol = 2 TempCol = 4 MixCol = 7 ThtwCol = 3 TSndCol = 1 DSndCol = 1 PrclCol = 5 BarbCol = -1 HodoCol = 1 * *-------------------- 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 * PrclLine = Line Style of parcel trace * HodoLine = Line Style of hodograph trace * ThetLine = 1 TempLine = 1 MixLine = 5 ThtwLine = 3 TSndLine = 1 DSndLine = 1 PrclLine = 5 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 * PrclThk = Line thickness of parcel trace * HodoThk = Line thickness of hodograph trace * BarbThk = Line thickness of wind barbs ThetThk = 3 TempThk = 3 MixThk = 3 ThtwThk = 3 TSndThk = 8 DSndThk = 8 PrclThk = 6 HodoThk = 3 BarbThk = 2 * !!!!! 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) *------------------------------------------------------------- * 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) * 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 DrawPrcl = 0 DrawPMax = 0 DrawIndx = 0 Endif If (sndtemp = -1) DrawTSnd = 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.5" 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 3 0.08" 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.06" "set string 1" If (PageType = "Portrait") "draw string 0.5 10.85 mb" Else "draw string 0.5 8.35 mb" 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.08" "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.06" "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.08" 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.06" 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" "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 "draw string " Botxloc-0.2 " " Botyloc-0.2 " " 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 Temperature (`3.`0C)" Else "draw string 3.5 0.15 Temperature (`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 Mixing Ratio (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 0" "set cthick "TSndThk "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 0" "set ccolor "DSndCol "set cstyle "DSndLine "set cthick "DSndThk "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 say PLcl *---------------------------------------------------------- * 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) 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 (tt > -900) 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,100,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=lxloc+0.75 Endif If (Text1YC = -1) Text1YC=tyloc-0.25 Endif If (Text2XC = -1) Text2XC=lxloc+0.75 Endif If (Text2YC = -1) Text2YC=tyloc-1.25 Endif If (Text3XC = -1) Text3XC=lxloc+0.75 Endif If (Text3YC = -1) Text3YC=tyloc-2.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 " PW(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 " Surface" "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 " Most Unstable" "draw string "Text3XC-0.70 " " Text3YC+0.20 " Press(mb)" "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 "sndspd rec=sublin(result,8) wspd=subwrd(rec,4) if (wspd < 0 | pres > pmax) zz=zz+1 else cont=0 Endif Endwhile While (zz <= zmax) "d "sndspd"(z="zz")" rec=sublin(result,8) wspd=subwrd(rec,4) 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 "snddir"(z="zz")" rec=sublin(result,8) wdir=subwrd(rec,4) 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 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 " knots" zloop=zmin xold=-999 yold=-999 count=0 Depth=0 While (zloop < zmax & Depth < HodoDep) "set z "zloop pres=subwrd(result,4) "d "sndspd rec=sublin(result,8) wspd=subwrd(rec,4) "d "snddir rec=sublin(result,8) wdir=subwrd(rec,4) uwnd=GetUWnd(wspd,wdir) vwnd=GetVWnd(wspd,wdir) 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 WgtTotal=0 * First, calculate mass-weighted mean wind "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) weight=log((pres+delp/2)/(pres-delp/2)) WgtTotal=WgtTotal+weight UTotal=UTotal+weight*uwnd VTotal=VTotal+weight*vwnd Endif pres=pres-delp EndWhile Umean=UTotal/WgtTotal Vmean=VTotal/WgtTotal Spdmean=GetWSpd(UTotal/WgtTotal,VTotal/WgtTotal) 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=lxloc+0.75 Endif If (Text4YC = -1) Text4YC=tyloc-3.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 *------------------------------------------ * 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. *--------------------------------------------------------------------- pres=startp PclTemp=startt PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15 delp=10 Pos=0 Neg=0 Neg2=0 logpp=log(pres) logppm=log(pres-delp) While (pres >= endp) EnvTemp=interp(sndtemp,pres) EnvDewp=interp(snddewp,pres) EnvTempV=virtual2(EnvTemp+273.15,EnvDewp+273.15,pres)-273.15 If (abs(EnvTempV) < 130 & abs(PclTempV) < 130) Val=(logpp-logppm)*(PclTempV-EnvTempV) If (Val > 0) Pos=Pos+Val Neg2=0 Else Neg=Neg+abs(Val) Neg2=Neg2+abs(Val) Endif Endif pres=pres-delp logpp=logppm logppm=log(pres-delp) * Need to iterate twice to get decent estimate at lapse rate * Perform first guess at prcl temp at pres PclTempI=PclTemp-100*delp*gammaw(PclTemp,pres+delp/2,100) * Estimate meant temperature of prcl path between pres and pres+delp MeanT=0.5*PclTemp+0.5*PclTempI * Perform second guess of prcl temp at pres using estimated MeanT PclTemp=PclTemp-100*delp*gammaw(MeanT,pres+delp/2,100) PclTempV=virtual2(PclTemp+273.15,PclTemp+273.15,pres)-273.15 Endwhile Pos=Pos*287 CIN=(Neg-Neg2)*287 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 lev "pres altpres=subwrd(result,4) "q dims" rec=sublin(result,4) zlev=subwrd(rec,9) 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) If (abs(VAbove) < 130 & abs(VBelow) < 130) Vest=VBelow+log(PBelow/pres)*(VAbove-VBelow)/log(PBelow/PAbove) Else Vest=-9999.0 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 *-------------------------------------------------------------- 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)