| Skew-T, Log(P) Diagram [message #3874] |
Tue, 04 April 1995 00:00 |
afl
Messages: 51
Registered: December 1994
|
Member |
|
|
Greetings!
A while back someone asked for an IDL procedure that
is capable of producing a Skew-T, Log(P) diagram. If you
have one, then I have been wasting my time... because I have just
finished beating my head against a wall writing this code
which accomplishes the task (thought it would be nice code
to have lying around). It is not well tested, so please
send comments to afl@cdc.noaa.gov so I can work out any bugs.
Now one must learn to plot sounding data onto the chart.
Use the Tnew function to find Temp. in the "skewed"
coordinate system!
It will also be nice to have the program compute such
quantities as the Wet bulb temp, LCL, LFC, CAPE, etc.
That's for the future.
Andy
; PROCEDURE TO DRAW A SKEW-T, Log(P) DIAGRAM GIVEN A DESIRED
; TEMPERATURE RANGE FOR THE DATA.
;
; Originator: Andrew F. Loughe
;
PRO SKEWT, TRANGE, everyT=everyT, everyDA=everyDA, $
everySA=everySA, everyW=everyW, title=title, notitle=notitle
on_error, 2
if (n_elements(everyT) le 0) then everyT = 10 ; T = Temperature
if (n_elements(everyDA) le 0) then everyDA = 10 ; DA = Dry adiabat
if (n_elements(everySA) le 0) then everySA = 1 ; SA = Saturated adiabat
if (n_elements(everyW) le 0) then everyW = 1 ; W = Mixing ratio
if (not keyword_set(title)) then title='Skew-T, Log(P) Diagram'
if (keyword_set(notitle)) then title=' '
if (N_params() eq 0) then $
message,$
'EXAMPLE: skewt, [-20, 20], everyT=10, everyDA=10, everySA=2, everyW=2'
; Set some defaults
prange = [1050, 100] ; Set default pressure range
charsize = .8 ; Set default character size
RED = 44
GREEN = 22
BLUE = 33
BLACK = 0
WHITE = 1
; Make plot square for arbitrarily chosen trange of 80 degrees.
; Code from Ken Bowman
if (!d.name eq 'PS') then device, /inch, xsize=7, ysize=7
daspect = FLOAT(!D.Y_SIZE)/FLOAT(!D.X_SIZE) * (trange(1)-trange(0))/80.
margin = 0.1
aspect = 1.0 ; A square
x0 = 0.50 - (0.5 - margin)*(daspect/aspect)
y0 = margin
x1 = 0.50 + (0.5 - margin)*(daspect/aspect)
y1 = 1.0 - margin
!P.POSITION = [x0, y0, x1, y1] ; Set value of sytem variable.
; Determine character height and width. Apply charsize.
char_ht = convert_coord([0, !d.y_ch_size], /device, /to_norm)
char_ht = char_ht(1) * 1.0
if (!d.name ne 'X' and charsize gt 1.) then $
char_ht = char_ht * charsize
char_wd = convert_coord([0, !d.x_ch_size], /device, /to_norm)
char_wd = char_wd(1)
; Create the plot space.
plot_io, trange, prange, yrange=prange, /nodata, /xs, /ys, $
xticklen=.01, ytickname=replicate(' ',30), charsize=charsize, $
title=title
; Print PRESSURE title along the y-axis.
lnt=alog(prange(1)) & lnb=alog(prange(0)) & avg=exp(.5*(lnt+lnb))
xy = convert_coord([trange(0), avg],/data,/to_norm)
xyouts, xy(0)-(5.*char_wd), xy(1), 'PRESSURE (hPa)', orient=90, $
/norm, align=.5
; Print TEMPERATURE title along the x-axis.
xy = convert_coord([.5*(trange(0)+trange(1)), prange(0)], /data, /to_norm)
xyouts, xy(0), xy(1)-(3.*char_ht), 'TEMPERATURE (!uo!nC)', align=.5, /norm
; Draw Pressure labels next to tick marks along the y-axis.
pressures = [1050,1000,900,800,700,600,500,400,300,200,100]
for i = 0, 10 do begin
ytick = pressures(i)
if (ytick ge prange(1)) then begin
xy = convert_coord( [trange(0), ytick], /data, /to_norm )
xyouts, xy(0)-(.2*char_wd), xy(1)-(.25*char_ht), $
strcompress(string(ytick),/remove_all), align=1, $
charsize=charsize, /norm
plots, [trange(0), trange(1)], [ytick, ytick] ; Horizontal line.
endif
endfor
clip=[trange(0),prange(0),trange(1),prange(1)] ; Define clipping space.
;=========================================================== =============
; Draw skewed isotherms every "everyT (10C)" (Lines are straight).
for temp = trange(0)-100, trange(1)+5, everyT do begin
x0 = temp
y0 = prange(0)
x1 = temp
y1 = prange(1)
; Draw the line.
newx0 = tnew(prange(0), x0, y0) ; Find rotated temperature position
newx1 = tnew(prange(0), x1, y1) ; Find rotated temperature position
plots, [newx0, newx1], [y0, y1], color=BLUE, clip=clip, noclip=0
; Draw line labels
; Use method #1 in xy function to determine a place for the label.
drew_label = 'no'
xy = Told(prange, trange, temp, 1)
if ( xy(0) gt trange(0) and xy(0) lt trange(1) and $
xy(1) gt prange(1) and xy(1) lt prange(0) ) then begin
drew_label = 'yes'
xyouts, xy(0), xy(1), strcompress(string(fix(temp)), /rem),$
color=BLUE, orient=45, align=.5, charsize=charsize
endif
; Use method #2 in xy function to determine a place for the label.
if (drew_label eq 'no') then xy = Told(prange, trange, temp, 2)
if ( xy(0) gt trange(0) and xy(0) lt trange(1) and $
xy(1) gt prange(1) and xy(1) lt prange(0) and $
drew_label eq 'no') then begin
xyouts, xy(0), xy(1), strcompress(string(fix(temp)), /rem),$
color=BLUE, orient=45, align=.5, charsize=charsize
endif
endfor
;=========================================================== =============
; Draw dry adiabats every "everyDA (10C)" (Lines are curved).
for temp = trange(0), trange(0)+220, everyDA do begin
x1 = float(temp)
y1 = 1050.
inc = -2. ; Lines will be curved, so use a small press. increment.
drew_label='no'
icount = 0
; Dry adiabats from 1050mb up to prange(1).
; For a given temperature and pressure, compute theta and plot a line.
for press = y1, prange(1), inc do begin
icount = icount + 1
x0 = float(x1) ; Orig Temp
y0 = float(press + inc) ; Orig Press
y1 = float(y0 + inc) ; New Press
x1 = (temp+273.16) * ( y1 / 1000. ) ^ (287./1004.) ; New Temp
x1 = x1 - 273.16
newx0 = tnew(prange(0), x0, y0) ; Find rotated temperature position
newx1 = tnew(prange(0), x1, y1) ; Find rotated temperature position
; Draw the labels.
if (fix(x1) eq fix(trange(0)) and drew_label eq 'no') then begin
drew_label='yes'
if ( newx1 gt trange(0) and newx1 lt trange(1) and $
y1 gt prange(1) and y1 lt prange(0) ) then $
xyouts,newx1,y1,strcompress(string(fix(temp)),/remove),$
align=.5, color=RED, charsize=charsize, orientation=-45
endif
; Draw the line.
if (icount gt 1) then $
plots, [newx0, newx1], [y0, y1], color=RED, clip=clip, noclip=0
if (newx1 lt trange(0)) then goto, jump2
endfor
jump2: dummy=0
endfor
;=========================================================== =============
; Draw saturated adiabats. Begin at 40C and step backwards by 4C.
; These lines are curved.
TS = 40.
FOR TS = 40, -64, -everySA*4 DO BEGIN
P = 1060.
TK = TS + 273.16
AOS = OS(TK, 1000.)
ATSA = TSA(AOS, P) - 273.16
FOR J = 0, 85 DO BEGIN
P0 = P
T0 = ATSA
P = P - 10.
ATSA = TSA(AOS, P) - 273.16
if (j gt 0) then begin
newx0=tnew(prange(0),T0,P0) ; Find rotated temperature position
newx1=tnew(prange(0),ATSA,P) ; Find rotated temperature position
; Leave a space for the labels and draw them.
if (P gt 520 or P lt 510) then $
plots, [newx0, newx1], [P0, P], $
color=GREEN, clip=clip, noclip=0
if ( P eq 520 ) then begin
if (newx1 gt trange(0) and newx1 lt trange(1)) then $
xyouts,newx1,P,strcompress(string(fix(TS)),/remove),align=.5 ,$
color=GREEN, charsize=charsize
endif
endif
ENDFOR
ENDFOR
;=========================================================== =============
; Draw mixing ratio lines (Lines are straight).
; Find temperature for a given Ws (g/kg) and Press (mb).
Ws=[ .1,.2,.4,.6,.8,1.,1.5,2.,2.5,4,5,6,7,8,9,10,12, $
14,16,18,20,24,28,32,36,40,44,48,52,56,60,68,76,84 ]
for i = 0, N_elements(Ws)-1, everyW do begin
press1 = prange(0)
tmr1 = tmr(Ws(i), press1) - 273.16
press2 = 200.
tmr2 = tmr(Ws(i), press2) - 273.16
newx0=tnew(prange(0),tmr1,press1) ; Find rotated temperature position
newx1=tnew(prange(0),tmr2,press2) ; Find rotated temperature position
; Draw the line.
plots, [newx0, newx1], [press1, press2], color=22, linestyle=2, $
clip=clip, noclip=0
; Draw the line label.
drew_label='no'
if (newx0 gt trange(0) and newx0 lt trange(1)) then begin
drew_label='yes'
if (Ws(i) ge 1.0) then $
xyouts, newx0, press1-2, strcompress(string(fix(Ws(i))),/remove),$
align=.5,color=GREEN, charsize=charsize
if (Ws(i) lt 1.0) then $
xyouts, newx0, press1-2, string(Ws(i),format='(f3.1)'), align=.5,$
color=GREEN, charsize=charsize
endif
if (newx1 gt trange(0) and newx1 lt trange(1)) then begin
if (Ws(i) ge 1.0) then $
xyouts, newx1, press2-2, strcompress(string(fix(Ws(i))),/remove),$
align=.5, color=GREEN, charsize=charsize
if (Ws(i) lt 1.0) then $
xyouts, newx1, press2-2, string(Ws(i),format='(f3.1)'), align=.5,$
color=GREEN, charsize=charsize
endif
endfor
;=========================================================== =============
; Redraw the plot boundary.
plots, [trange(0),trange(1),trange(1),trange(0),trange(0)], $
[prange(0),prange(0),prange(1),prange(1),prange(0)], thick=2
; Close Postscript device, rename output file, return to calling program.
if ( !d.name eq 'PS') then begin
device, /close
spawn, 'mv idl.ps skewt.ps'
set_plot, 'X'
endif
END
--
Andrew F. Loughe email: afl@cdc.noaa.gov
University of Colorado, CIRES voice: (303) 492-0707
Campus Box 449 fax: (303) 497-7013
Boulder, CO 80309-0449 USA
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
