;Title: HRstat.pro ; ;Purpose: Program to calculate heating rate statistics ; ; ;Author: Matthew Shupe ;Date: 2/12/02 ;Modified: 2/12/02 ;----------------------------------------------------------- ;pro HRstat ;Get a list of all the files (key off of the BugsRad output) cd,'/data1/sheba/heatrate' bfiles=findfile('rtoutput_bugsrad*.cdf',count=numf) fargs=strmid(bfiles,17,8) ;Set up distribution vectors and output parameters hrtbins=findgen(201)/20.0 - 5.0 ;-5 to +5 in 0.05 wide bins hrfbins=findgen(201)/20.0 - 5.0 ;-5 to +5 in 0.05 wide bins nhrtb=n_elements(hrtbins) nhrfb=n_elements(hrfbins) dist_hrtbug=fltarr(nhrtb,6) ;HRT dists for "in cloud" (tot,SW,LW) and "all-sky" (tot,SW,LW) dist_hrtstr=dist_hrtbug dist_hrfbug=fltarr(nhrfb,6) ;HRF dists for "in cloud" (tot,SW,LW) and "all-sky" (tot,SW,LW) dist_hrfstr=dist_hrfbug dstat_hrtbug=fltarr(numf,8,6) ;daily mean, stddev, median, num, min, max for each of the above types dstat_hrfbug=dstat_hrtbug dstat_hrtstr=dstat_hrtbug dstat_hrfstr=dstat_hrtbug stat_hrtbug=fltarr(6,6) ;total mean, stddev, median, num, min, max for each of the above types stat_hrfbug=stat_hrtbug stat_hrtstr=stat_hrtbug stat_hrfstr=stat_hrtbug for i=0,numf-1 do begin ;Read in the input data farg='rtinput*'+fargs[i]+'*.cdf' file=findfile(farg) fid=ncdf_open(file[0]) ncdf_varget,fid,ncdf_varid(fid,'time'),time ncdf_varget,fid,ncdf_varid(fid,'hl'),height ncdf_varget,fid,ncdf_varid(fid,'wcil'),iwc ncdf_varget,fid,ncdf_varid(fid,'wcwl'),lwc ncdf_varget,fid,ncdf_varid(fid,'rcil'),ire ncdf_varget,fid,ncdf_varid(fid,'rcwl'),lre ncdf_close,fid ;Read in bugsrad data fid=ncdf_open(bfiles[i]) ncdf_varget,fid,ncdf_varid(fid,'asl'),swhra_b ncdf_varget,fid,ncdf_varid(fid,'aslR'),swhrc_b ncdf_varget,fid,ncdf_varid(fid,'atl'),lwhra_b ncdf_varget,fid,ncdf_varid(fid,'atlR'),lwhrc_b ;ncdf_varget,fid,ncdf_varid(fid,'tclice'),tcice_b ;ncdf_varget,fid,ncdf_varid(fid,'tclwat'),tcwat_b ncdf_close,fid hrt_b=(lwhra_b+swhra_b) hrf_b=(lwhra_b+swhra_b)-(lwhrc_b+swhrc_b) lwhrf_b=(lwhra_b-lwhrc_b) swhrf_b=(swhra_b-swhrc_b) ;Read in Streamer data farg='rtoutput_stream.'+fargs[i]+'.cdf file=findfile(farg,count=nums) if nums gt 0 then begin fid=ncdf_open(file[0]) ncdf_varget,fid,ncdf_varid(fid,'asl'),swhra_s ncdf_varget,fid,ncdf_varid(fid,'aslR'),swhrc_s ncdf_varget,fid,ncdf_varid(fid,'atl'),lwhra_s ncdf_varget,fid,ncdf_varid(fid,'atlR'),lwhrc_s ;ncdf_varget,fid,ncdf_varid(fid,''), ;ncdf_varget,fid,ncdf_varid(fid,''), ;tcice_s=tcice_b*0.0-999 ;tcwat_s=tcwat_b*0.0-999 ncdf_close,fid endif else begin swhra_s=swhra_b*0.0-999. swhrc_s=swhra_s lwhra_s=swhra_s lwhrc_s=swhra_s ;tcice_s=tcice_b*0.0-999 ;tcwat_s=tcwat_b*0.0-999 endelse hrt_s=(lwhra_s+swhra_s) hrf_s=(lwhra_s+swhra_s)-(lwhrc_s+swhrc_s) lwhrf_s=(lwhra_s-lwhrc_s) swhrf_s=(swhra_s-swhrc_s) ;Calculate the distributions stats cwh=where(iwc gt 0.0001 or lwc gt 0.0001 and abs(swhra_b) lt 100 and abs(lwhra_b) lt 100 and abs(swhra_s) lt 100 and abs(lwhra_s) lt 100,num) ;Define in cloud iwh=where(abs(swhra_b) lt 100 and abs(lwhra_b) lt 100 and abs(swhra_s) lt 100 and abs(lwhra_s) lt 100,anum) if num gt 0 then begin dist_hrtbug[*,0]=dist_hrtbug[*,0]+histogram(hrt_b[cwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtbug[*,1]=dist_hrtbug[*,1]+histogram(swhra_b[cwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtbug[*,2]=dist_hrtbug[*,2]+histogram(lwhra_b[cwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtstr[*,0]=dist_hrtstr[*,0]+histogram(hrt_s[cwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtstr[*,1]=dist_hrtstr[*,1]+histogram(swhra_s[cwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtstr[*,2]=dist_hrtstr[*,2]+histogram(lwhra_s[cwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dstat_hrtbug[i,*,0]=[mean(hrt_b[cwh]),stddev(hrt_b[cwh]),median(hrt_b[cwh]),num,min(hrt_b[cwh]),max(hrt_b[cwh]),total(hrt_b[cwh]^2),total(hrt_b[cwh])] dstat_hrtbug[i,*,1]=[mean(swhra_b[cwh]),stddev(swhra_b[cwh]),median(swhra_b[cwh]),num,min(swhra_b[cwh]),max(swhra_b[cwh]),total(swhra_b[cwh]^2),total(swhra_b[cwh])] dstat_hrtbug[i,*,2]=[mean(lwhra_b[cwh]),stddev(lwhra_b[cwh]),median(lwhra_b[cwh]),num,min(lwhra_b[cwh]),max(lwhra_b[cwh]),total(lwhra_b[cwh]^2),total(lwhra_b[cwh])] dstat_hrtstr[i,*,0]=[mean(hrt_s[cwh]),stddev(hrt_s[cwh]),median(hrt_s[cwh]),num,min(hrt_s[cwh]),max(hrt_s[cwh]),total(hrt_s[cwh]^2),total(hrt_s[cwh])] dstat_hrtstr[i,*,1]=[mean(swhra_s[cwh]),stddev(swhra_s[cwh]),median(swhra_s[cwh]),num,min(swhra_s[cwh]),max(swhra_s[cwh]),total(swhra_s[cwh]^2),total(swhra_s[cwh])] dstat_hrtstr[i,*,2]=[mean(lwhra_s[cwh]),stddev(lwhra_s[cwh]),median(lwhra_s[cwh]),num,min(lwhra_s[cwh]),max(lwhra_s[cwh]),total(lwhra_s[cwh]^2),total(lwhra_s[cwh])] dist_hrfbug[*,0]=dist_hrfbug[*,0]+histogram(hrf_b[cwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfbug[*,1]=dist_hrfbug[*,1]+histogram(swhrf_b[cwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfbug[*,2]=dist_hrfbug[*,2]+histogram(lwhrf_b[cwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfstr[*,0]=dist_hrfstr[*,0]+histogram(hrf_s[cwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfstr[*,1]=dist_hrfstr[*,1]+histogram(swhrf_s[cwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfstr[*,2]=dist_hrfstr[*,2]+histogram(lwhrf_s[cwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dstat_hrfbug[i,*,0]=[mean(hrf_b[cwh]),stddev(hrf_b[cwh]),median(hrf_b[cwh]),num,min(hrf_b[cwh]),max(hrf_b[cwh]),total(hrf_b[cwh]^2),total(hrf_b[cwh])] dstat_hrfbug[i,*,1]=[mean(swhrf_b[cwh]),stddev(swhrf_b[cwh]),median(swhrf_b[cwh]),num,min(swhrf_b[cwh]),max(swhrf_b[cwh]),total(swhrf_b[cwh]^2),total(swhrf_b[cwh])] dstat_hrfbug[i,*,2]=[mean(lwhrf_b[cwh]),stddev(lwhrf_b[cwh]),median(lwhrf_b[cwh]),num,min(lwhrf_b[cwh]),max(lwhrf_b[cwh]),total(lwhrf_b[cwh]^2),total(lwhrf_b[cwh])] dstat_hrfstr[i,*,0]=[mean(hrf_s[cwh]),stddev(hrf_s[cwh]),median(hrf_s[cwh]),num,min(hrf_s[cwh]),max(hrf_s[cwh]),total(hrf_s[cwh]^2),total(hrf_s[cwh])] dstat_hrfstr[i,*,1]=[mean(swhrf_s[cwh]),stddev(swhrf_s[cwh]),median(swhrf_s[cwh]),num,min(swhrf_s[cwh]),max(swhrf_s[cwh]),total(swhrf_s[cwh]^2),total(swhrf_s[cwh])] dstat_hrfstr[i,*,2]=[mean(lwhrf_s[cwh]),stddev(lwhrf_s[cwh]),median(lwhrf_s[cwh]),num,min(lwhrf_s[cwh]),max(lwhrf_s[cwh]),total(lwhrf_s[cwh]^2),total(lwhrf_s[cwh])] endif if anum gt 2 then begin dist_hrtbug[*,3]=dist_hrtbug[*,3]+histogram(hrt_b[iwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtbug[*,4]=dist_hrtbug[*,4]+histogram(swhra_b[iwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtbug[*,5]=dist_hrtbug[*,5]+histogram(lwhra_b[iwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtstr[*,3]=dist_hrtstr[*,3]+histogram(hrt_s[iwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtstr[*,4]=dist_hrtstr[*,4]+histogram(swhra_s[iwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dist_hrtstr[*,5]=dist_hrtstr[*,5]+histogram(lwhra_s[iwh],min=hrtbins[0],max=hrtbins[nhrtb-1],binsize=(hrtbins[nhrtb-1]-hrtbins[0])/(nhrtb-1)) dstat_hrtbug[i,*,3]=[mean(hrt_b[iwh]),stddev(hrt_b[iwh]),median(hrt_b[iwh]),anum,min(hrt_b[iwh]),max(hrt_b[iwh]),total(hrt_b[iwh]^2),total(hrt_b[iwh])] dstat_hrtbug[i,*,4]=[mean(swhra_b[iwh]),stddev(swhra_b[iwh]),median(swhra_b[iwh]),anum,min(swhra_b[iwh]),max(swhra_b[iwh]),total(swhra_b[iwh]^2),total(swhra_b[iwh])] dstat_hrtbug[i,*,5]=[mean(lwhra_b[iwh]),stddev(lwhra_b[iwh]),median(lwhra_b[iwh]),anum,min(lwhra_b[iwh]),max(lwhra_b[iwh]),total(lwhra_b[iwh]^2),total(lwhra_b[iwh])] dstat_hrtstr[i,*,3]=[mean(hrt_s[iwh]),stddev(hrt_s[iwh]),median(hrt_s[iwh]),anum,min(hrt_s[iwh]),max(hrt_s[iwh]),total(hrt_s[iwh]^2),total(hrt_s[iwh])] dstat_hrtstr[i,*,4]=[mean(swhra_s[iwh]),stddev(swhra_s[iwh]),median(swhra_s[iwh]),anum,min(swhra_s[iwh]),max(swhra_s[iwh]),total(swhra_s[iwh]^2),total(swhra_s[iwh])] dstat_hrtstr[i,*,5]=[mean(lwhra_s[iwh]),stddev(lwhra_s[iwh]),median(lwhra_s[iwh]),anum,min(lwhra_s[iwh]),max(lwhra_s[iwh]),total(lwhra_s[iwh]^2),total(lwhra_s[iwh])] dist_hrfbug[*,3]=dist_hrfbug[*,3]+histogram(hrf_b[iwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfbug[*,4]=dist_hrfbug[*,4]+histogram(swhrf_b[iwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfbug[*,5]=dist_hrfbug[*,5]+histogram(lwhrf_b[iwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfstr[*,3]=dist_hrfstr[*,3]+histogram(hrf_s[iwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfstr[*,4]=dist_hrfstr[*,4]+histogram(swhrf_s[iwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dist_hrfstr[*,5]=dist_hrfstr[*,5]+histogram(lwhrf_s[iwh],min=hrfbins[0],max=hrfbins[nhrfb-1],binsize=(hrfbins[nhrfb-1]-hrfbins[0])/(nhrfb-1)) dstat_hrfbug[i,*,3]=[mean(hrf_b[iwh]),stddev(hrf_b[iwh]),median(hrf_b[iwh]),anum,min(hrf_b[iwh]),max(hrf_b[iwh]),total(hrf_b[iwh]^2),total(hrf_b[iwh])] dstat_hrfbug[i,*,4]=[mean(swhrf_b[iwh]),stddev(swhrf_b[iwh]),median(swhrf_b[iwh]),anum,min(swhrf_b[iwh]),max(swhrf_b[iwh]),total(swhrf_b[iwh]^2),total(swhrf_b[iwh])] dstat_hrfbug[i,*,5]=[mean(lwhrf_b[iwh]),stddev(lwhrf_b[iwh]),median(lwhrf_b[iwh]),anum,min(lwhrf_b[iwh]),max(lwhrf_b[iwh]),total(lwhrf_b[iwh]^2),total(lwhrf_b[iwh])] dstat_hrfstr[i,*,3]=[mean(hrf_s[iwh]),stddev(hrf_s[iwh]),median(hrf_s[iwh]),anum,min(hrf_s[iwh]),max(hrf_s[iwh]),total(hrf_s[iwh]^2),total(hrf_s[iwh])] dstat_hrfstr[i,*,4]=[mean(swhrf_s[iwh]),stddev(swhrf_s[iwh]),median(swhrf_s[iwh]),anum,min(swhrf_s[iwh]),max(swhrf_s[iwh]),total(swhrf_s[iwh]^2),total(swhrf_s[iwh])] dstat_hrfstr[i,*,5]=[mean(lwhrf_s[iwh]),stddev(lwhrf_s[iwh]),median(lwhrf_s[iwh]),anum,min(lwhrf_s[iwh]),max(lwhrf_s[iwh]),total(lwhrf_s[iwh]^2),total(lwhrf_s[iwh])] endif endfor ;Calculate total stats from the daily stats for i=0,5 do begin sd=((total(dstat_hrtbug[*,6,i])-(total(dstat_hrtbug[*,7,i])^2.0)/total(dstat_hrtbug[*,3,i]))/(total(dstat_hrtbug[*,3,i])-1.0))^0.5 stat_hrtbug[*,i]=[(total(dstat_hrtbug[*,0,i]*dstat_hrtbug[*,3,i]))/total(dstat_hrtbug[*,3,i]),sd,median(dstat_hrtbug[*,2,i]),total(dstat_hrtbug[*,3,i]),min(dstat_hrtbug[*,4,i]),max(dstat_hrtbug[*,5,i])] sd=((total(dstat_hrfbug[*,6,i])-(total(dstat_hrfbug[*,7,i])^2.0)/total(dstat_hrfbug[*,3,i]))/(total(dstat_hrfbug[*,3,i])-1.0))^0.5 stat_hrfbug[*,i]=[(total(dstat_hrfbug[*,0,i]*dstat_hrfbug[*,3,i]))/total(dstat_hrfbug[*,3,i]),sd,median(dstat_hrfbug[*,2,i]),total(dstat_hrfbug[*,3,i]),min(dstat_hrfbug[*,4,i]),max(dstat_hrfbug[*,5,i])] sd=((total(dstat_hrtstr[*,6,i])-(total(dstat_hrtstr[*,7,i])^2.0)/total(dstat_hrtstr[*,3,i]))/(total(dstat_hrtstr[*,3,i])-1.0))^0.5 stat_hrtstr[*,i]=[(total(dstat_hrtstr[*,0,i]*dstat_hrtstr[*,3,i]))/total(dstat_hrtstr[*,3,i]),sd,median(dstat_hrtstr[*,2,i]),total(dstat_hrtstr[*,3,i]),min(dstat_hrtstr[*,4,i]),max(dstat_hrtstr[*,5,i])] sd=((total(dstat_hrfstr[*,6,i])-(total(dstat_hrfstr[*,7,i])^2.0)/total(dstat_hrfstr[*,3,i]))/(total(dstat_hrfstr[*,3,i])-1.0))^0.5 stat_hrfstr[*,i]=[(total(dstat_hrfstr[*,0,i]*dstat_hrfstr[*,3,i]))/total(dstat_hrfstr[*,3,i]),sd,median(dstat_hrfstr[*,2,i]),total(dstat_hrfstr[*,3,i]),min(dstat_hrfstr[*,4,i]),max(dstat_hrfstr[*,5,i])] endfor cd,'/data1/sheba/heatrate/stats' save,hrtbins,hrfbins,dist_hrfbug,dist_hrtbug,dist_hrtstr,dist_hrfstr,stat_hrfbug,stat_hrtbug,stat_hrtstr,stat_hrfstr,$ filename='hrstat_sheba.dat' end ;program