; ; This program creates a skewt diagram from the sounding data ; pro pltsondeskewt, site=site,date=date, path=path, desire_gif=desire_gif, $ outpath=outpath, upper_left=upper_left, dostop=dostop, $ t_range=t_range, title=title, facility=facility COMMON RANGES, trange, prange ; ; Establish the location of the data based on the site ; if site eq 'nsa' then begin armdatadir='/data/mace2/arm_data/nsa/' datastream='nsaisssonde10sC1.a1/' armfile='nsaisssonde10sC1.a1'+curdate+'*.cdf*' endif else if site eq 'sgp' then begin armdatadir='/data/mace3/arm_data/sgp/' datastream='sgpsondewnpnC1.a1/' armfile='sgpsondewnpnC1.a1'+curdate+'*.cdf*' endif else if site eq 'twpc1' then begin armdatadir='/data/mace/arm_data/twp/Manus_C1' datastream='twpsondewnpnC1.a1/' armfile='twpsondewnpnC1.a1'+curdate+'*.cdf*' date='' ;openr, 1,'/data/mace/envernon/atmos_o/idl_code/today' openr,1,'today.dat' readf, 1, date close, 1 ;date='20000104' datapath = '/home/grips_y/arm_data/nsa' ;outpath='/data/mace/envernon/atmos_l' outpath='/home/sbenson/Arm_data/Graphics' opath = getenv("SA_QUICKLOOK_HOME")+'/nsa' copypath = '' ; Define platform names and suffixes wrpn_pn = 'nsaisssonde10s' wrpn_sf = '.a1' ls_pn = '' ls_sf = '' pcharsize = 1 tcharsize = 2 ; Manage the keywords ; if there is a date and no facility if ((n_elements(date) ne 0) and (n_elements(facility) eq 0)) then facility = 'C1' ; if you don't want a gif if (n_elements(desire_gif) eq 0) then desire_gif = 0 ; if you want a Z buffer if (!d.name eq 'Z') then begin device,set_resolution=[640,512] desire_gif = 1 endif ; If you don't have a specific title in mind if (n_elements(title) eq 0) then use_default_title=1 $ else use_default_title=0 if (n_elements(path) eq 0) then path = datapath if (keyword_set(t_range) eq 0) then t_range = [-30, 35] ; ; Set up the colors ; loadct,12, ncolors=32 pmain=0 d_red=24 d_dgreen=1 d_blue=15 color=13 ; Input netCDF file of type sgpsondewrpn ; ; If you don't have a specific date in mind, the program ; lets you pick the file from a dialog box ; if (n_elements(date) eq 0) then begin if (n_elements(facility) ne 0) then $ pushd, path + '/' + wrpn_pn + facility + los_sf $ else pushd, path filename = dialog_pickfile(/read, filter='*sondewrpn*cdf') if (filename eq '') then $ message, 'Operation cancelled' $ else begin sepstring = str_sep(filename, '.') date = sepstring[3] endelse if (n_elements(facility) eq 0) then facility = strmid(filename, $ strlen(path)+13, 2) endif $ ; ; If you have a secific date in mind, the filename is created from the pieces ; else begin pushd, path filename=wrpn_pn+facility+wrpn_sf+'/'+wrpn_pn+facility+wrpn_sf+'.*'+date+'.*.cdf' endelse ; ; Stores the names of the files that match what we want ; files = findfile(filename, count=count) ;puts all the files matching the filename into a string if (count eq 0) then begin ; ; No matching wrpn file found ; nowrpn = 1 print,'no wrpn ',filename endif else begin ; ; Found matching wrpn files ; print,'yes wrpn ',filename nowrpn = 0 for j=0, count-1 do begin ;do for all the matching files in the directory ; Read in data and cut array to only good values fid = ncdf_open(files[j]) ;open the netcdf file ncdf_varget, fid, 'base_time', bt ;base time ncdf_varget, fid, 'time_offset', to ;time_offset ncdf_varget, fid, 'press', pres ;pressure presmin=min(pres) ;minimum in the pressure data min_index=where(pres eq presmin) ;index where pressure is at a minimum min_index=min_index[0] to=to[0:min_index] ;cut array to be from the surface to the pressure minimum pres=pres[0:min_index] ncdf_varget, fid, 'dp', dp ;dewpoint dp=dp[0:min_index] ncdf_varget, fid, 'tdry', temp ;temperature temp=temp[0:min_index] ncdf_varget,fid,'wspd',wspd ;wind speed wspd=wspd[0:min_index] ncdf_varget,fid,'deg',wdir ;wind direction wdir=wdir[0:min_index] ncdf_close, fid ;close the netcdf file ; ; Restrict the arrays to only good values for temp and dp ; result=where(temp ne 999 and dp ne 999,count) to=to[result] pres=pres[result] dp=dp[result] temp=temp[result] wspd=wspd[result] wdir=wdir[result] ; Read hhmmss sepstring = str_sep(files[j], '.') hms = sepstring[4] if (j eq 0) then begin ;if it is the first file btto = bt+to[0] num_elem_wrpn = n_elements(to) hms_wrpn = hms pres_wrpn = pres dp_wrpn = dp temp_wrpn = temp wspd_wrpn = wspd wdir_wrpn = wdir endif $ ;end of the first file else begin ;if it is not the first file ;put everthing from the different files into one array btto = [btto, bt+to[0]] num_elem_wrpn = [num_elem_wrpn, n_elements(to)] hms_wrpn = [hms_wrpn, hms] pres_wrpn = [pres_wrpn, pres] dp_wrpn = [dp_wrpn, dp] temp_wrpn = [temp_wrpn, temp] wspd_wrpn = [wspd_wrpn, wspd] wdir_wrpn = [wdir_wrpn, wdir] endelse ;end of not the first file endfor ;end of loop through all matching files endelse ;end of found matching files popd ; Input netCDF file(s) of type sgplssonde ; ; If you have a secific date in mind, the filename is created from the pieces ; pushd, path filename=ls_pn+facility+ls_sf+'/'+ls_pn+facility+ls_sf+'.*'+date+'.*.cdf' files = findfile(filename, count=count) if (count eq 0) then begin ; ; No matching iss file found ; no_ls = 1 print, 'no iss ', filename endif else begin ; ; Found matching iss files ; print, 'yes iss ', filename no_ls = 0 for j=0, count-1 do begin ;do for all the matching files in the directory ; Read in data fid = ncdf_open(files[j]) ;open the netcdf file ncdf_varget, fid, 'time_offset', to ;time offset ncdf_varget, fid, 'press', pres ;pressure ncdf_varget, fid, 'tdry', temp ;temperature ncdf_varget, fid, 'rh', rh ;relative humidity ncdf_varget, fid, 'dp', dp ;dewpoint (I added this. It is in the nsa iss file) ;ncdf_attget, fid, /global, 'MWR_precipitable_water_vapor', pwv ;not in the nsa iss file pwv=0.1 ;some value to make it run ;ncdf_attget, fid, /global, 'H2O_Scale_Factor', h2o ;not in the nsa iss file h2o=0.2 ;some value to make it run ; Convert relative humidity to dew point temperature ;dp = rh2dpt(temp, rh/100.0) ;dp is in the nsa iss file ; Read hhmmss sepstring = str_sep(files[j], '.') hms = sepstring[4] if (j eq 0) then begin ;if it is the first file num_elem_ls = n_elements(to) hms_ls = hms pres_ls = pres dp_ls = dp temp_ls = temp pwv_ls = float(string(pwv)) h2o_ls = float(string(h2o)) endif $ ;end of if it is the first file else begin ;if it isn't the first file ;put everthing from the different files into one array num_elem_ls = [num_elem_ls, n_elements(to)] hms_ls = [hms_ls, hms] pres_ls = [pres_ls, pres] dp_ls = [dp_ls, dp] temp_ls = [temp_ls, temp] pwv_ls = [pwv_ls, float(string(pwv))] h2o_ls = [h2o_ls, float(string(h2o))] endelse ;end of if it isn't the first file endfor ;end of loop through all matching files endelse ;end of found a matching iss file popd ; Manage name of facility to be printed on plot if (facility eq 'C1') then print_facility = 'CF' $ else print_facility = facility ; Exit program if WRPN file is missing if (nowrpn eq 1) then $ message, 'No WRPN files for '+date+', facility '+print_facility+ ' were found.' print, ' A WRPN file was found for date '+date+', facility'+print_facility+', time '+hms_wrpn ; Don't exit if the ISS file is missing if (no_ls eq 0) then $ print, ' An LS SONDE file was found for date '+date+', facility'+print_facility+', time ' + hms_ls ; Match up LS files with corresponding WRPN files ls_exists = make_array(n_elements(hms_wrpn), /integer, value=-1) if (no_ls eq 0) then begin ;there is a correstponding LS file for i=0, n_elements(hms_wrpn)-1 do begin ;loop through the time of the LS file foo=where(hms_ls eq hms_wrpn[i], count) if (count eq 1) then ls_exists[i] = foo ;array of indexes of the matching LS array endfor endif ; ; Establish the outpath for the program ; if ((not keyword_set(outpath)) and (!d.name ne 'X')) then begin outpath = opath+'/'+wrpn_pn+facility+wrpn_sf copypath = opath+'/'+ls_pn+facility+ls_sf endif else outpath = '.' ; Plot skew-T diagrams for the wrpn file for i=0, n_elements(hms_wrpn)-1 do begin century = '' if (strlen(date) lt 8) then century = '19' outname = 'pltsondeskewt'+facility+'.'+century+date+'.'+hms_wrpn[i] yyyy=strmid(date,0,4) mm=strmid(date,4,2) dd=strmid(date,6,2) hh=strmid(hms_wrpn[i],0,2) nn=strmid(hms_wrpn[i],2,2) ss =strmid(hms_wrpn[i],4,2) if (!d.name eq 'PS') then $ device, filename=outpath+'/'+outname+'.ps', /landscape, /color $ else if (!d.name eq 'Z') then erase skewt, t_range ;call skewt subroutine ; Print title ; systime2ymdhms, btto[i], yyyy,mm,dd,hh,nn,ss mnth = [' Jan ',' Feb ',' Mar ',' Apr ',' May ',' Jun ', $ ' Jul ',' Aug ',' Sep ',' Oct ',' Nov ',' Dec '] ; ; Stops if the plot tries to cross over a year boundry ; if(yyyy(0) ne yyyy(n_elements(yyyy)-1)) then $ if(keyword_set(dostop)) then $ stop, 'Crossing the year bounds not handled well - stopping' $ else $ message,'Crossing the year bounds not handled well.' ; ; Make a date and time string ; datestring = string(format='(I0)',dd(0)) + mnth(mm(0)-1) + $ string(format='(I0)',yyyy(0)) timestring = string(format='(I0)',hh(0))+':'+string(format='(I0)', nn(0)) $ +':'+string(format='(I0)',ss(0)) ; ; Title the plot ; if (use_default_title eq 1) then $ title='nsaisssonde10sC1.a1 sonde data for '+datestring+', '+timestring $ ;+ ', facility ' + print_facility $ else begin datestring='date: '+datestring+ $ ', time: '+timestring+ $ ', facility: '+ print_facility xyouts, 0.01, 0.04, datestring, color=pmain, /normal, charsize=pcharsize endelse xyouts, 0.5, 0.95, title, color=pmain, align=0.5, charsize=tcharsize, /normal ; ; Print upper left (string) ; if (n_elements(upper_left) ne 0) then $ xyouts, 0.01, 0.975, upper_left, /normal, color=pmain ; ; Plot LS data, if it exists ; j = ls_exists[i] if (j ge 0) then begin if (j eq 0) then first = 0 $ else first = total(num_elem_ls[0:j-1]) last = total(num_elem_ls[0:j]) - 1 plot_skewt, temp_ls[first:last], dp_ls[first:last], pres_ls[first:last],$ ;call plot_skewt subroutine color=d_red scale_factor='LS Scale Factor = ' + $ strcompress(string(h2o_ls[j]), /remove_all) pwv_string='MWR PWV = ' + $ strcompress(string(pwv_ls[j]), /remove_all) xyouts, 0.99, 0.04, scale_factor, color=pmain, $ align=1, /normal, charsize=pcharsize xyouts, 0.99, 0.02, pwv_string, color=pmain, /normal, align=1, $ charsize=pcharsize endif ; ; Plot WRPN data ; if (i eq 0) then first = 0 $ else first = total(num_elem_wrpn[0:i-1]) last = total(num_elem_wrpn[0:i]) - 1 plot_skewt, temp_wrpn[first:last], dp_wrpn[first:last], pres_wrpn[first:last], $ ;call plot_skewt subroutine color=d_green ;xyouts, 0.01, 0.02, 'blue=WRPN, red=LS', color=pmain, /normal, $ ;charsize=pcharsize ; ; PLOT SONDE WINDS ; windbarb, wdir_wrpn[first:last], wspd_wrpn[first:last], pres_wrpn[first:last] ; ; END PLOTTING WINDS ; ; Create the gif ; ;outpath='/home/atmos/l/envernon/temp_image/' outpath='/home/sbenson/Arm_data/Graphics/' gname=outpath+outname+'.gif' WRITE_GIF, gname, tvrd() ; if (!d.name eq 'PS') then $ ; print, ' Writing PS: '+outpath+'/'+outname+'.ps' $ ; else if (desire_gif eq 1) then begin ; print, ' Writing GIF: '+outpath+'/'+outname+'.gif' ; write_gif, outpath+'/'+outname+'.gif', tvrd(), gred, ggreen, gblue ; endif ; if (copypath ne '') then begin ; if (!d.name eq 'Z') then begin ; print, ' Writing GIF: '+copypath+'/'+outname+'.gif' ; strspawn = 'cp '+outpath+'/'+outname+'.gif '+copypath+'/'+outname+'.gif' ; endif else begin ; print, ' Writing PS: '+copypath+'/'+outname+'.ps' ; strspawn = 'cp '+outpath+'/'+outname+'.ps '+copypath+'/'+outname+'.ps' ; endelse ; spawn,strspawn ; endif ; Allow user to view each diagram successively ; if ((!d.name ne 'Z') and (!d.name ne 'PS')) then begin ; print, 'Press any key to continue:' ; result=get_kbrd(1) ; print, 'Proceeding with execution...' ; endif endfor ;end of loop through all the wrpn data if (!d.name ne 'X') then device,/close ; ; Reset these variables to zero ; !p.multi=0 !p.region=0 !p.position=0 if (keyword_set(dostop)) then stop, 'Stopped in procedure as indicated' if ((!d.name ne 'Z') and (!d.name ne 'PS')) then print, 'Done.' return end ;------------------------------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------------------------------ ; HERE IS WHERE THE SKEWT STUFF STARTS ;------------------------------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------------------------------------------------------------------------------ ; $Id: skewt.pro,v 1.1 1999/11/23 21:24:27 halter Release_1_0 $ ;======================================================================== ; SKEWT.PRO (IDL CODE) ; ; Draw a Skew-T, Log(P) diagram given a temperature range for your data. ; ; Originator: Andrew F. Loughe (afl@cdc.noaa.gov) ; CIRES/NOAA ; Boulder, CO USA ; This code carries no warranty or claim ; as to its usefulness or accuracy! ; ; A Number of the functions found in this file were converted from ; FORTRAN code that was received from NCAR in Boulder, CO USA. ; The original source of the equations is thought to be: ; "Algorithms for Generating a Skew-T, Log P Diagram ; and Computing Selected Meteorological Quantities" ; by G.S. Stipanuk, White Sands Missle Range, Report ECOM-5515. ; ;======================================================================== ; FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE GIVEN TEMP IN KELVIN. ; ESAT_sk(MILLIBARS), T(KELVIN) FUNCTION ESAT_sk, T TC = T - 273.16 ESAT_sk = 6.1078 * EXP( (17.2693882 * TC) / (TC + 237.3) ) RETURN, ESAT_sk END ;======================================================================== ; FUNCTION TO COMPUTE MIXING RATIO GIVEN TEMP. AND PRESS. ; W(GRAMS WATER VAPOR/KILOGRAM DRY AIR), P(MILLIBAR) FUNCTION W_sk, T, P IF (T GE 999.) THEN GOTO, JUMP10 X = ESAT_sk(T) W_sk = 621.97 * X / (P - X) RETURN, W_sk JUMP10: W_sk = 0.0 RETURN, W_sk END ;======================================================================== ; FUNCTION TO COMPUTE SATURATION ADIABATIC TEMP AT 1000 MB GIVEN T & P. ; OS AND T (KELVIN), P (MILLIBARS ) FUNCTION OS_sk, T, P IF (T LT 100.) THEN T1 = T + 273.16 IF (T GE 100.) THEN T1 = T OS_sk = T1 * ((1000./P)^.286) / (EXP( -2.6518986*W_sk(T1,P) / T1) ) RETURN, OS_sk END ;======================================================================== ; FUNCTION TO COMPUTE THE TEMPERATURE (KELVIN) OF AIR AT A GIVEN ; PRESSURE AND WITH A GIVEN MIXING RATIO. ; TMR(KELVIN), W(GRAMS WATER VAPOR/KILOGRAM DRY AIR), P(MILLIBAR) FUNCTION TMR_sk, W, P X = ALOG10 ( W * P / (622.+ W) ) TMR_sk = 10. ^ ( .0498646455 * X + 2.4082965 ) - 7.07475 + $ 38.9114 * ( (10.^( .0915 * X ) - 1.2035 )^2 ) RETURN, TMR_sk END ;======================================================================== ; FUNCTION TO COMPUTE TEMPERATUE (KELVIN) OF A MOIST ADIABAT GIVEN ; OS(KELVIN), P(MILLIBARS) ; SIGN(A,B) REPLACES THE ALGEBRAIC SIGN OF A WITH THE SIGN OF B FUNCTION TSA_sk, OS, P A = OS TQ = 253.16 D = 120 FOR I = 1, 12 DO BEGIN D = D/2. ; IF THE TEMPERATURE DIFFERENCE, X, IS SMALL, EXIT THIS LOOP. X = A * EXP (-2.6518986*W_sk(TQ,P)/TQ)-TQ*((1000./P)^.286) IF ( ABS(X) LT 0.01 ) THEN GOTO, JUMP2 ; TQ = TQ + SIGN(D,X) IF (X LT 0) THEN D = -ABS(D) IF (X GT 0) THEN D = ABS(D) TQ = TQ + D ENDFOR JUMP2: TSA_sk=TQ RETURN, TSA_sk END ;======================================================================== ; Function to determine osition (temp, press) when the isotherms ; in the diagram are rotated (skewed) 45 degrees to the right. ; Used for finding the points needed to connect the dots when ; drawing ALL of the lines (except pressure). ; Originator: Andrew F. Loughe Function Tnew, T, P COMMON RANGES, trange, prange P0 = prange(0) xy1 = convert_coord( [T, P0], /data, /to_device) xy2 = convert_coord( [T, P], /data, /to_device) dy = xy2(1) - xy1(1) dx = dy ; dx = dy for this 45-45-90 triangle xy = convert_coord( [xy2(0)+dx, xy2(1)], /device, /to_data) Tnew = xy(0) return, Tnew end ;======================================================================== ; Function to determine osition (temp, press) in the unskewed ; coordinate system (Opposite of Tnew). ; Used only when placing the labels on various lines. ; Originator: Andrew F. Loughe Function Told, T, METHOD COMMON RANGES, trange, prange P0 = prange(0) P1 = prange(1) T0 = trange(0) T1 = trange(1) if (method eq 1) then begin xy1 = convert_coord( [T, P0], /data, /to_device ) xy2 = convert_coord( [T0, P0], /data, /to_device ) dx = xy2(0) - xy1(0) xy = convert_coord( [xy2(0), xy2(1)+dx], /device, /to_data ) xy1 = convert_coord( [xy(0), xy(1)], /data, /to_device ) xy2 = convert_coord( [xy(0), P1], /data, /to_device ) dy = xy2(1) - xy1(1) xy = convert_coord([xy1(0)+(dy/2.), xy1(1)+(dy/2.)],$ /device, /to_data) endif if (method eq 2) then begin xy1 = convert_coord( [T, P0], /data, /to_device ) xy2 = convert_coord( [T1, P0], /data, /to_device ) dx = xy2(0) - xy1(0) xy = convert_coord( [xy1(0)+dx/2. , xy1(1)+dx/2.], $ /device, /to_data) endif return, xy end ;======================================================================== ; Function to determine the necessary trange for plotting the sounding. Function T_RANGE, p_range, t, td, p COMMON RANGES, trange, prange trange = [-40, 40] ; Default which can be changed if (p_range eq 0) then $ prange = [1050, 100] ; Default which can be changed ; Find number of data levels szd = size(t) nlevels = szd(1) ; Ensure that temperatures are in Celsius. if ( (total(t)/nlevels) gt 100. ) then t = t - 273.16 if ( (total(td)/nlevels) gt 100. ) then td = td - 273.16 ; Set up dummy plot space. 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] plot_io, trange, prange, yrange=prange, /nodata, /xs, /ys, $ color=!p.background, /noerase !p.multi(0) = !p.multi(0) + 1 ; Determine necessary temperature range for the diagram. xx0 = fltarr(nlevels) & yy0=xx0 & xx1=xx0 & yy1=yy0 for i = 0, nlevels-1 do begin xx0(i) = tnew( t(i), p(i) ) yy0(i) = p(i) xx1(i) = tnew( td(i), p(i) ) yy1(i) = p(i) endfor xbegin = fix( (min(xx1)-10.) / 10. ) * 10. xend = fix( (max(xx0)+10.) / 10. ) * 10. return, [xbegin, xend] end ;======================================================================== ; Routine to plot the temperature and dew point temperature sounding ; on top of a skew-T, Log(P) diagram. ; ; Keyword(s): ; COLOR: If set, the color of the line to be plotted over the data line ; for purposes of data visualization pro plot_skewt, t, td, p, color=color pmain=0 d_red=124 d_dgreen=1 d_blue=15 color=13 COMMON RANGES, trange, prange ; Find number of data levels szd = size(t) nlevels = szd(1) ; Define clipping space. clip=[trange(0),prange(0),trange(1),prange(1)] ; Ensure that temperatures are in Celsius. if ( (total(t)/nlevels) gt 100. ) then t = t - 273.16 if ( (total(td)/nlevels) gt 100. ) then td = td - 273.16 ; Overplot the data onto the digram. for i = 0, nlevels-2 do begin ; Plot temperature sounding data. x0 = tnew( t(i), p(i) ) y0 = p(i) x1 = tnew( t(i+1), p(i+1) ) y1 = p(i+1) plots, [x0, x1], [y0, y1], clip=clip, noclip=0, thick=2, color=pmain if (keyword_set(color) eq 1) then plots, [x0, x1], [y0, y1], clip=clip, noclip=0, $ thick=1, color=color ; Plot dew point temperature sounding data. x0 = tnew( td(i), p(i) ) y0 = p(i) x1 = tnew( td(i+1), p(i+1) ) y1 = p(i+1) plots, [x0, x1], [y0, y1], clip=clip, noclip=0, thick=2, color=pmain if (keyword_set(color) eq 1) then plots, [x0, x1], [y0, y1], clip=clip, noclip=0, $ thick=1, color=color endfor ; Close Postscript device, rename output file. ;tdh if ( !d.name eq 'PS') then begin ;tdh device, /close ;tdh spawn, 'mv idl.ps skewt.ps' ;tdh set_plot, 'X' ;tdh endif end ;======================================================================== ; ; PROCEDURE TO DRAW A SKEW-T, Log(P) DIAGRAM GIVEN A DESIRED ; TEMPERATURE RANGE FOR THE DATA. ; ; Originator: Andrew F. Loughe ; PRO SKEWT, t_range, everyT=everyT, everyDA=everyDA, $ everySA=everySA, everyW=everyW, title=title, notitle=notitle on_error, 2 COMMON RANGES, trange, prange pmain=0 d_red=24 d_dgreen=1 d_blue=15 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_elements(prange)) eq 0 then prange = [1050., 100] if (N_params() eq 0) then $ message,$ 'EXAMPLE: skewt, [-20, 20], everyT=10, everyDA=10, everySA=2, everyW=2' if (n_elements(t_range)) eq 1 then t_range=[-40., 40.] ; Set some defaults trange = t_range charsize = 1.25 ; Set default character size window, 0, xsize=800, ysize=800, /pixmap !p.multi = [0,1,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=5, /ys, $ xticklen=.01, ytickname=replicate(' ',30), charsize=charsize, $ title=title, color=pmain, background=!d.n_colors-1 ; 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, color=pmain ; 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, color=pmain ; 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, color=pmain plots, [trange(0), trange(1)], [ytick, ytick], color=pmain ; 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(x0, y0) ; Find rotated temperature position newx1 = tnew(x1, y1) ; Find rotated temperature position plots, [newx0, newx1], [y0, y1], color=pmain, 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(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=pmain, 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(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=pmain, 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(x0, y0) ; Find rotated temperature position newx1 = tnew(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=d_red, charsize=charsize, orientation=-45 endif ; Draw the line. if (icount gt 1) then $ plots, [newx0, newx1], [y0, y1], color=d_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_sk(TK, 1000.) ATSA = TSA_sk(AOS, P) - 273.16 FOR J = 0, 85 DO BEGIN P0 = P T0 = ATSA P = P - 10. ATSA = TSA_sk(AOS, P) - 273.16 if (j gt 0) then begin newx0=tnew(T0,P0) ; Find rotated temperature position newx1=tnew(ATSA,P) ; Find rotated temperature position ; Leave a space for the labels and draw them. if (P gt 730 or P lt 700) then $ plots, [newx0, newx1], [P0, P], $ color=d_dgreen, clip=clip, noclip=0 if ( P eq 730 ) then begin if (newx1 gt trange(0) and newx1 lt trange(1)) then $ xyouts,newx1,P,strcompress(string(fix(TS)),/remove),align=.5,$ color=d_dgreen, 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,3,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) - 20 tmr1 = tmr_sk(Ws(i), press1) - 273.16 press2 = 200. tmr2 = tmr_sk(Ws(i), press2) - 273.16 newx0=tnew(tmr1,press1) ; Find rotated temperature position newx1=tnew(tmr2,press2) ; Find rotated temperature position ; Draw the line. plots, [newx0, newx1], [press1, press2], linestyle=0, $ clip=clip, noclip=0 , color=d_blue ; 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 3.0) then $ xyouts, newx0, press1+70, strcompress(string(fix(Ws(i))),/remove),$ align=.5,color=d_blue, charsize=charsize-.3 if (Ws(i) lt 3.0) then $ xyouts, newx0, press1+70, string(Ws(i),format='(f3.1)'), align=.5,$ color=d_blue, charsize=charsize-.3 endif if (newx1 gt trange(0) and newx1 lt trange(1)) then begin if (Ws(i) ge 2.0) then $ xyouts, newx1, press2-2, strcompress(string(fix(Ws(i))),/remove),$ align=.5, color=d_blue, charsize=charsize-.2 if (Ws(i) lt 2.0) then $ xyouts, newx1, press2-2, string(Ws(i),format='(f3.1)'), align=.5,$ color=d_blue, charsize=charsize-.2 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, color=pmain !p.position = [.05, .05, .95, .95] ; Reset position parameter END