如何利用GRADS和1*1 Ncep再分析资料画CAPE对流有效位能

happyouc · 发表于 2018-3-8 09:21:24

schliezer 发表于 2014-3-10 15:53
这个变量我画过，也拿EC和NCEP的相关数据用skewplot.gs算过，一对比三者量级相差很大，直接画出来的是最 ...

我的qq289619047 谢谢

schliezer · 发表于 2018-3-8 11:04:01

happyouc 发表于 2018-3-8 09:21
我的qq289619047 谢谢

我没有专门算CAPE的程序，传个skewplot的gs程序，可能会输出CAPE，前半截只要改一下路径和变量名就好了，然后直接套用到plotskew()函数里面，后半截是skewplot.gs程序本身：

'reinit'
'open D:\gradsdata\ncepfnl_1_6h\fnl_20120719_00_00_c.ctl'
'set grid off'
'set grads off'
'set z 1 21'
'set lat 37.75'
'set lon 103.92'
tt=1
while(tt<=13)
'c'
'set grads off'
'set t 'tt' '
'define tmp=tmpprs-273.16'
'define RH=rhprs'
'td=tmp-((14.55+0.114*tmp)*(1-0.01*RH) + pow((2.5+0.007*tmp)*(1-0.01*RH),3) + (15.9+0.37*tmp)*pow((1-0.01*RH),14))'
'define uv=mag(vgrdprs,ugrdprs)'
'define pi=3.14159'
if (ugrdprs>=0.)
if (vgrdprs>=0.)
'define dir=180+atan2(ugrdprs,vgrdprs)*180/pi'
endif
endif
if (vgrdprs<=0.)
if (ugrdprs>=0.)
'define dir=270+atan2(ugrdprs,vgrdprs)*180/pi'
endif
endif
if (ugrdprs<=0.)
if (vgrdprs>=0.)
'define dir=180-atan2(ugrdprs,vgrdprs)*180/pi'
endif
endif
if (vgrdprs<0.)
if (ugrdprs<=0.)
'define dir=270-atan2(ugrdprs,vgrdprs)*180/pi'
endif
endif
*if(ugrdprs > 0.);'define dir=90-180/pi*atan2(vgrdprs,ugrdprs)';endif
*if(ugrdprs < 0.);'define dir=270-180/pi*atan2(vgrdprs,ugrdprs)';endif
rc=plotskew(tmp,td,uv,dir)
'printim D:\gradsdata\ncepfnl_1_6h\0719\站点垂直剖面\观测点'%tt%'.gif white'
tt=tt+1
endwhile
;
function plotskew(sndtemp,snddewp,sndspd,snddir)
*************************************************************************
*
* 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
*
* Please send any further problems, comments, or suggestions to
* <hart@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 hPa ; 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
* Another good set of defining points suggested by George Bryan (PSU)
* are:
*
* P1 = 1000
* T1 = -30
*
* P2 = 1000
* T2 = 40
*
* P3 = 500
* T3 = -18
* ------------------- 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 2 3 4 6 8 10 15 20 25 30 35 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
* 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= 0
DrawTemp= 1
DrawMix = 1
DrawTSnd= 1
DrawDSnd= 1
DrawRH = 0
DrawPrcl= 1
DrawPMax= 0
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 hPa) 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 (hPa) of mean-wind calculation
* MeanVBot = Bottom pressure level (hPa) of mean-wind calculation
* HelicDep = Depth in hPa (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
* 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 = 2
TempCol = 4
MixCol = 7
ThtwCol = 3
TSndCol = 1
DSndCol = 1
RHCol = 3
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
* 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 = 3
TSndLine = 1
DSndLine = 2
RHLine = 3
PrclLine = 1
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 = 3
MixThk = 3
ThtwThk = 3
TSndThk = 8
DSndThk = 8
RHThk = 8
PrclThk = 6
HodoThk = 3
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 = 0
DSndMrk = 0
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)
* 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 = 1
DrawPMax = 0
DrawIndx = 0
Endif
If (sndtemp = -1)
DrawTSnd = 0
DrawRH = 0
DrawPrcl = 1
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.10"
"set xlopts 1 3 0.125"
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.10"
"set string 1 c 3 0"
If (PageType = "Portrait")
"draw string 0.5 10.85 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
"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 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 "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 250hPa 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,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=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 " 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(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 "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
* 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 hPa
*-------------------------------------------------------------
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 hPa
*---------------------------------------------------------------
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 hPa*---------------------------------------------------------------
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)

复制代码

yn_newyork · 发表于 2018-3-27 09:30:44

schliezer 发表于 2014-3-10 15:53
这个变量我画过，也拿EC和NCEP的相关数据用skewplot.gs算过，一对比三者量级相差很大，直接画出来的是最 ...

能把您计算cape的脚本发给我吗我也是想利于基础数据计算cape 565653979@qq 万分感谢

yn_newyork · 发表于 2018-4-2 10:58:18

schliezer 发表于 2014-3-10 15:53
这个变量我画过，也拿EC和NCEP的相关数据用skewplot.gs算过，一对比三者量级相差很大，直接画出来的是最 ...

能把cape的计算方法发给我一份吗，我希望通过微波辐射计计算cape，微波辐射计也是只有基础数据

一颗醉甜的草莓 · 发表于 2019-10-17 09:23:02

happyouc 发表于 2018-3-8 09:20
请教能不能把画cape的程序分享一下我想拿ec再分析资料画谢谢

你好你用EC的计算过cape么为什么我计算得到的值和EC提供的cape值相差很大

pefuen · 发表于 2020-3-24 17:04:48

求问1°×1°的怎么下载的！我只找到2.5°×2.5°的再分析资料

		自动登录	找回密码
密码			立即注册

[分享资料] 如何利用GRADS和1*1 Ncep再分析资料画CAPE对流有效位能

浏览过的版块