Note how one routine gets one job, numbers are not hard-wired in and keywords are used to set defaults.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro read_all_radiosondes,indir,data,length_files,heights,nv=nv ;+ ;Will read all radiosondes in indir ; ;eg read_all_radiosondes,'~jmarsham/public_html/Teaching/IDL_course/Radiosondes/040720/',data,length_files,heights ;- if not keyword_Set (nv) then nv=6 ;Findfile(IDL routine) is useful! files=findfile(indir+'*',count=n_files) newh=1 ; require interpolated heights for data nvar=6 ;n variables in sonde file ; NOTE: each sonde file has a different number of lines. We cannot of ; put data arrays of different sizes into one big array, so we just ; make it big many=200 ; a large number of lines headers=dblarr(nv,n_files) ; array of headers data=fltarr(nv,many,n_files); array of data heights=fltarr(many,n_files); array of heights length_files=fltarr(n_files); length of each file for i = 0, n_files-1 do begin read_radiosonde,indir+files(i),temphead1,temphead2,tempdata,newheights=newh,nvar=nv headers(*,i)=temphead1 length_files(i)=(size(tempdata))(2);(size(tempdata))(2)=number of lines data(*,0:length_files(i)-1,i)=tempdata heights(0:length_files(i)-1,i)=reform(newh) endfor stop end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro read_radiosonde,infile,header1,header2,data,newheights=newheights,nvar=nvar ;+ ;pro read_radiosonde,infile,header1,header2,data,newheights=newheights ; Reads a stanadrd radiosonde file (infile) into header1,header2,data. ; Will intrepolate pressuer to newheights if newheights is set ; ;INPUTS: infile = file to be read (standardformat) ; nvar=nvar: number of columns in sonde file - default=6 ; ;OUTPUTS header1 = header line 1 (fltarr(6)) ; header2 = header line 2 (strarr(6)) ; data = fltarr(6,n_lines) ; newheights=newheights: heights of data calculated by interpolating ; pressure (logarithmically). ; ;JOHN MARSHAM 13/4/05 ; ;eg read_radiosonde,'/nfs/env-fs-04_u35/jmarsham/BLASIUS/Radiosondes/040720/03743_20040720_11.dat',hedaer1,header2,data,newheights=1 ;- if not keyword_Set(nvar) then nvar=6 header1=dblarr(nvar) ; need double precision to get date correct! ;0 1 2 3 4 5 ;Station code yyymmdd hh ? ? Surface height(m) header2=strarr(nvar) ; PP HT TT TD DD FF ; 0 1 2 3 4 5 i=0 line=fltarr(nvar) openr,12,infile readf,12,header1 readf,12,header2 repeat begin readf, 12, line if i eq 0 then data=[line] else data=[[data],[line]];6 by i array of data i=i+1 endrep until eof(12) close,12 ; Calculate interpolated heights if keyword_Set(newheights) then begin ht=reform(data(1,*));original height data (m) w_ht=where(data(1,*) gt 0,cw_ht);where height is defined ;interpolate newheights using log(pressure) newheights=interpol(ht(w_ht),(alog(data(0,*)))(w_ht),(alog(data(0,*)))) ;TO CHECK THIS WORKS ;w,0 ;r_Set ;plot,data(1,*),alog(data(0,*)),psym=1,xrange=[0,20000],yrange=[9,12],$ ;xtitle='Height (m)',ytitle='ln(pressure/mb)' ;oplot,newheights,alog(data(0,*)),psym=3 ;w,1 ;plot,data(1,*),data(0,*),psym=1,xrange=[0,20000],yrange=[0,1.e5],$ ;xtitle='Height (m)',ytitle='Pressure (mb)' ;oplot,newheights,data(0,*),psym=3 endif end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro plot_sonde_th_td_z,data,sonde,heights,headers ;+ ;plot_sonde_th_td_z,data,sonde,heights,headers ; ;Plot theta vs height and dewpoint vs height for data(6 ;var,nlines,nfiles) ;on one graph ;sonde defines file number, heights are interpolated heights and ;headers are headers for files(See read_all_radiosondes) ;JOHN MARSHAM 13/4/05 ;- ;convert temp to potential temp theta=temp_to_theta(data(2,*,sonde),data(0,*,sonde),ref_press=100000) w,0 ; create window !p.color=0 ;set colour !p.background=256*256*256l-1 ; set background w_theta=where(theta gt 0,cw_theta) ; where pot temp defined w_td=where(data(3,*,sonde) gt -100,cw_td) ; where dew point defined xmin=280;min theta for plot xmax=420; ymin=0;min height for plot ymax=20000 pos=[0.23,0.15,0.95,0.9];position on page - normalised coordinates ;Plot theta and height force axis range, don't draw axes plot,theta(w_theta),(heights(*,sonde))(w_theta),position=pos,$ xrange=[xmin,xmax],yrange=[ymin,ymax],xstyle=5,ystyle=5 ;draw axes axis,xaxis=0,xrange=[xmin,xmax],xstyle=1,xtitle='Potential temperature (K)' axis,yaxis=0,yrange=[ymin,ymax],ystyle=1,ytitle='Height (m)' ;Plot dewpoint data on top xmin2=200 xmax2=300 plot,(data(3,*,sonde))(w_td),(heights(*,sonde))(w_td),xrange=[xmin2,xmax2],$ yrange=[ymin,ymax],/noerase,position=pos,xstyle=5,ystyle=5,linestyle=2 ;draw axis for dew point axis,xaxis=1,xrange=[xmin2,xmax2],xstyle=1,xtitle='Dew point (K)' ; Write time on top xyouts,210,18000,trmfix(fix(headers(1,sonde),type=3))+', '+$ trmfix(fix(headers(2,sonde),type=3))+' UTC',charsize=2 end;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;