pro read_mrr_ave_txt

;  Undefine the variable
;pres_o=!NULL
;temp_o=!NULL
;dp_o=!NULL
;rh_o=!NULL
;wspd_o=!NULL
;deg_o=!NULL
;alt_o=!NULL
U=''
SVS=''
Frequency=replicate(!values.f_nan,[64,31])
Diameter=replicate(!values.f_nan,[64,31])
Number=replicate(!values.f_nan,[64,31])
nans=replicate(!values.f_nan,31)
dummy=replicate(!values.F_NAN,[64,31])
dummyd=replicate(!values.F_NAN,[64,31])
dummyn=replicate(!values.F_NAN,[64,31])
first=0

;  Get the data from the txt file
file='1211.ave.txt'
openr,lun,file,/get_lun

while not eof(lun) do begin
  line=''
  readf,lun,line
  ;print,line
  indicator=''
  reads,line,indicator,format='(A1)'
  ;*** M *** TIME
  if indicator eq 'M' then begin
    u=''
    reads,line,y,m,d,h,mm,s,U,ave,stp,seal,smp,svs,dvs,dsn,cc,mdq,$
      format='(4x,I2,I2,I2,I2,I2,I2,1x,A3,4x,I6,4x,I6,4x,I6,4x,E6,4x,A8,4x,F5,4x,I11,3x,I8,4x,I4)'    
    if (first eq 0) then begin
      year=y & month=m & day=d & hour=h & minute=mm & seconds=s & UTC=U
      avg=ave & resolution=stp & asl=seal & sample_rate=smp
      ;Noise_0=nf0 & Noise_1=nf1
      Version_service=svs & Version_firmware=dvs & serial=dsn
      calibration = cc & valid_spectra=mdq          
    endif else begin
      year=[year,y] & month=[month,m] & day=[day,d] & hour=[hour,h] & minute=[minute,mm] & seconds=[seconds,s] & UTC=[UTC,U]
      avg=[avg,ave] & resolution=[resolution,stp] & asl=[asl,seal] & sample_rate=[sample_rate,smp]
      ;Noise_0=[Noise_0,nf0] & Noise_1=[Noise_1,nf1]
      Version_service=[Version_service,svs] & Version_firmware=[Version_firmware,dvs] & serial=[serial,dsn]
      calibration=[calibration,cc] & valid_spectra=[valid_spectra,mdq] 
    endelse  
    ;print,y,m,d,h,mm,s,U,ave,stp,seal,smp,svs,dvs,dsn,cc,mdq     
  ;*** H ***  HEIGHT
  endif else if indicator eq 'H' then begin
    he=fltarr(31)
    reads,line,he,format='(3x,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7,I7)'
    if (first eq 0) then begin
      height = he
    endif else begin
      height = [[height],[he]]
    endelse   
  ;*** T ***  TRANSFER FUNCTION     
  endif else if indicator eq 'T' then begin
    trans=fltarr(31)
    reads,line,trans,format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    if (first eq 0) then begin
      Transfer_function = trans
    endif else begin
      Transfer_function = [[Transfer_function],[trans]]       
    endelse 
  ;*** F ***  Spectral frequancy            
  endif else if indicator eq 'F' then begin
    u=0
    reads,line,u,format='(1x,I2)'
    num=strlen(line)
    num=fix((num-3.)/7.)
    ;num=31
    if (num gt 0) then freq=replicate(!values.f_nan,num) else freq=!values.f_nan
    if (num gt 0) then reads,line,freq,format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    nan_replace=where(freq eq 0.0,count)
    if (count gt 0) then freq[nan_replace]=!values.F_NAN
    ;if (num le 30) then freq=[freq,nans[num:30]]
    ;if (num eq 0) then freq=nans
    if (first eq 0) then begin
      Frequency[u,*]=freq
    endif else begin
      dummy[u,*]=freq       
    endelse  
  ;*** D ***  DROP SIZE 
  endif else if indicator eq 'D' then begin
    reads,line,u,format='(1x,I2)'
    num=strlen(line)
    num=fix((num-3.)/7.)
    if (num gt 0) then diam=replicate(!values.f_nan,num) else diam=!values.f_nan
    if (num gt 0) then reads,line,diam,format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    nan_replace=where(diam eq 0.0,count)
    if (count gt 0) then diam[nan_replace]=!values.F_NAN
    ;if (num le 30) then diam=[diam,nans[num:30]]
    ;if (num eq 0) then diam=nans
    if (first eq 0) then begin
      Diameter[u,*] = diam
    endif else begin
      dummyd[u,*] = diam       
    endelse  
  ;*** N *** Number     
  endif else if indicator eq 'N' then begin
    reads,line,u,format='(1x,I2)'
    num=strlen(line)
    num=fix((num-3.)/7.)
    if (num gt 0) then nums=replicate(!values.f_nan,num) else nums=!values.f_nan
    if (num gt 0) then reads, line, nums, format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    nan_replace=where(nums eq 0.0,count)
    if (count gt 0) then nums[nan_replace]=!values.F_NAN
    if (num le 30) then nums=[nums,nans[num:30]]
    if (num eq 0) then nums=nans
    if (first eq 0) then begin
      Number[u,*] = nums
    endif else begin
      dummyn[u,*] = nums       
    endelse   
  ;*** P *** Path Integrated Attenuation 
  endif else if indicator eq 'P' then begin 
    num=strlen(line)
    num=fix((num-3.)/7.)
    if (num gt 0) then path=replicate(!values.f_nan,num) else path=!values.f_nan
    if (num gt 0) then reads, line, path, format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    if (num le 30) then path=[path,nans[num:30]]
    if (num eq 0) then path=nans
    if (first eq 0) then begin
      Path_attenuation = path
    endif else begin
      Path_attenuation = [[Path_attenuation],[path]]       
    endelse  
  ;*** z *** Attenuated Radar Reflectivity      
  endif else if indicator eq 'z' then begin
    num=strlen(line)
    num=fix((num-3.)/7.)
    if (num gt 0) then a_reflect=replicate(!values.f_nan,num) else a_reflect=!values.f_nan
    if (num gt 0) then reads, line, a_reflect, format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    if (num le 30) then a_reflect=[a_reflect,nans[num:30]]
    if (num eq 0) then a_reflect=nans
    if (first eq 0) then begin
      Reflectivity_attenuated = a_reflect
    endif else begin
      Reflectivity_attenuated = [[Reflectivity_attenuated],[a_reflect]]       
    endelse   
  ;*** Z *** Radar reflectivity 
  endif else if indicator eq 'Z' then begin
    num=strlen(line)
    num=fix((num-3.)/7.)
    if (num gt 0) then reflect=replicate(!values.f_nan,num) else reflect=!values.f_nan
    if (num gt 0) then reads, line, reflect, format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    if (num le 30) then reflect=[reflect,nans[num:30]]
    if (num eq 0) then reflect=nans
    if (first eq 0) then begin
      Reflectivity = reflect
    endif else begin
      Reflectivity = [[Reflectivity],[reflect]]       
    endelse    
  ;*** R ***  rain rate             
  endif else if indicator eq 'R' then begin
    num=strlen(line)
    num=fix((num-3.)/7.)
    if (num gt 0) then rate=replicate(!values.f_nan,num) else rate=!values.f_nan
    if (num gt 0) then reads, line, rate, format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    if (num le 30) then rate=[rate,nans[num:30]]
    if (num eq 0) then rate=nans
    if (first eq 0) then begin
      Rain_rate = rate
    endif else begin
      Rain_rate = [[Rain_rate],[rate]]       
    endelse    
  ;*** L *** liquid water content
  endif else if indicator eq 'L' then begin
    num=strlen(line)
    num=fix((num-3.)/7.)
    if (num gt 0) then lwc=replicate(!values.f_nan,num) else lwc=!values.f_nan
    if (num gt 0) then reads, line, lwc, format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    if (num le 30) then lwc=[lwc,nans[num:30]]
    if (num eq 0) then lwc=nans
    if (first eq 0) then begin
      Liquid_water_content = lwc
    endif else begin
      Liquid_water_content = [[Liquid_water_content],[lwc]]       
    endelse  
  ;*** W *** Fall Velocity        
  endif else if indicator eq 'W' then begin   
    num=strlen(line)
    num=fix((num-3.)/7.)
    if (num gt 0) then velo=replicate(!values.f_nan,num) else velo=!values.f_nan
    if (num gt 0) then reads, line, velo, format='(3x,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7,F7)'
    if (num le 30) then velo=[velo,nans[num:30]]
    if (num eq 0) then velo=nans
    if (first eq 0) then begin
      Fall_velocity = velo
    endif else begin
      Fall_velocity = [[Fall_velocity],[velo]]       
    endelse     
    ;readf, lun,scrap   
    if (first ne 0) then Frequency = [[[Frequency]],[[dummy]]]
    if (first ne 0) then Diameter=[[[Diameter]],[[dummyd]]]
    if (first ne 0) then Number=[[[Number]],[[dummyn]]]
    first=first+1 
    dummy=replicate(!values.F_NAN,[64,31])
    dummyd=replicate(!values.F_NAN,[64,31])
    dummyn=replicate(!values.F_NAN,[64,31])         
  endif    
endwhile
free_lun,lun

rot_num=6
julian_day=julday(month,day,year,hour,minute,seconds)
reflectivity_attenuated=rotate(reflectivity_attenuated,rot_num)
reflectivity=rotate(reflectivity,rot_num)
height=rotate(height,rot_num)
rain_rate=rotate(rain_rate,rot_num)
liquid_water_content=rotate(liquid_water_content,rot_num)
fall_velocity=rotate(fall_velocity,rot_num)

;*********************************************************
;  Plot data
;*********************************************************

;  Color common
common colors,r_orig,g_orig,b_orig,r_curr,g_curr,b_curr

;  Set up the device
set_plot,'z'

;  xsize of the plot window
xwin=800  

;  ysize of the plot window
ywin=900  

;  Size the window
device,set_resolution=[xwin,ywin]

;  Load the color table
loadct,39   

;  Make the background white
!p.background=!d.n_colors-1  

top_color=254
;  Date label format
dummy=label_date(date_format=['%H:%I'])

;  Xticks point out from plot
!x.ticklen=-0.02
!y.ticklen=-0.02

xl=0.08  &  xr=0.85
yb=0.05  &  yt=0.95
sx=0.10        
sy=0.05 
numplots_x=1  ;number of plots in x direction
numplots_y=5  ;number of plots in y direction

position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos

;  String positioning
dr=0.09  ;string offset from the right xcoor
dl=0.025  ;string offset from the left xcoor
dt=0.025  ;string offset from the top ycoor

cbpos=pos
cbpos[*,0]=pos[*,2]+0.09
cbpos[*,2]=cbpos[*,0]+0.005

;**************
;  Reflectivity
;**************
rmax=max([reflectivity,reflectivity_attenuated])
rmin=min([reflectivity,reflectivity_attenuated])
pnum=0
image_var=bytscl(reflectivity,min=rmin,max=rmax,top=top_color)

;  Set up the plot space
contour,reflectivity,julian_day,height[0,*],/nodata,$
  xstyle=4,ystyle=4,position=pos[pnum,*]

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image_var,sx,sy),px[0],py[0]

;  Put up the axes
contour,reflectivity,julian_day,height[0,*]/1000.0,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='label_date',/nodata,$
  ytitle='Height (km)',position=pos[pnum,*],/noerase
xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'reflectivity',$
  color=0,/normal
;  Color bar
;cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
;       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=rmax,minrange=rmin,$
  ncolors=top_color,format='(F6.2)',$
  vertical=1,color=0,position=cbpos[pnum,*],charsize=cbsz

;**************
;  Reflectivity-attenuated
;**************
;dmax=max(reflectivity_attenuated)
;dmin=min(reflectivity_attenuated)
pnum=1
image_var=bytscl(reflectivity_attenuated,min=rmin,max=rmax,top=top_color)

;  Set up the plot space
contour,reflectivity_attenuated,julian_day,height[0,*],/nodata,$
  xstyle=4,ystyle=4,position=pos[pnum,*],/noerase

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image_var,sx,sy),px[0],py[0]

;  Put up the axes
contour,reflectivity_attenuated,julian_day,height[0,*]/1000.0,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='label_date',/nodata,$
  ytitle='Height (km)',position=pos[pnum,*],/noerase
xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'reflectivity attenuated',$
  color=0,/normal
;  Color bar
;cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
;       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=rmax,minrange=rmin,$
  ncolors=top_color,format='(F6.2)',$
  vertical=1,color=0,position=cbpos[pnum,*],charsize=cbsz

;**************
;  Rain-rate
;**************
dmax=max(rain_rate)
dmin=min(rain_rate)
pnum=2
image_var=bytscl(rain_rate,top=top_color)

;  Set up the plot space
contour,rain_rate,julian_day,height[0,*],/nodata,$
  xstyle=4,ystyle=4,position=pos[pnum,*],/noerase

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image_var,sx,sy),px[0],py[0]

;  Put up the axes
contour,rain_rate,julian_day,height[0,*]/1000.0,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='label_date',/nodata,$
  ytitle='Height (km)',position=pos[pnum,*],/noerase
xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'rain rate',$
  color=255,/normal
;  Color bar
;cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
;       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=top_color,format='(F6.2)',$
  vertical=1,color=0,position=cbpos[pnum,*],charsize=cbsz

;**************
;  liquid water content
;**************
dmax=max(liquid_water_content)
dmin=min(liquid_water_content)
pnum=3
image_var=bytscl(liquid_water_content,top=top_color)

;  Set up the plot space
contour,liquid_water_content,julian_day,height[0,*],/nodata,$
  xstyle=4,ystyle=4,position=pos[pnum,*],/noerase

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image_var,sx,sy),px[0],py[0]

;  Put up the axes
contour,liquid_water_content,julian_day,height[0,*]/1000.0,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='label_date',/nodata,$
  ytitle='Height (km)',position=pos[pnum,*],/noerase
xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'liquid water content',$
  color=255,/normal
;  Color bar
;cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
;       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=top_color,format='(F6.2)',$
  vertical=1,color=0,position=cbpos[pnum,*],charsize=cbsz

;**************
;  fall_velocity
;**************
dmax=max(fall_velocity)/9.0
dmin=min(fall_velocity)
pnum=4
image_var=bytscl(fall_velocity,top=top_color)

;  Set up the plot space
contour,fall_velocity,julian_day,height[0,*],/nodata,$
  xstyle=4,ystyle=4,position=pos[pnum,*],/noerase

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image_var,sx,sy),px[0],py[0]

;  Put up the axes
contour,fall_velocity,julian_day,height[0,*]/1000.0,xstyle=1,ystyle=1,$
  color=0,xtickunits='time',xtickformat='label_date',/nodata,$
  ytitle='Height (km)',position=pos[pnum,*],/noerase
xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'fall_velocity',$
  color=255,/normal
;  Color bar
;cbpos=[(pos[pnum,2]+0.050),(pos[pnum,1]+0.02),$
;       (pos[pnum,2]+0.055),(pos[pnum,3]-0.02)]
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=top_color,format='(F6.2)',$
  vertical=1,color=0,position=cbpos[pnum,*],charsize=cbsz

write_gif,'mrr.timeseries.gif',tvrd()

;***********************************************************
;  Plot spectra
;***************************************************
erase
xl=0.1  &  xr=0.85
yb=0.1  &  yt=0.8
sx=0.20
sy=0.15
numplots_x=2  ;number of plots in x direction
numplots_y=2  ;number of plots in y direction

position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
cbpos=pos
cbpos[*,0]=pos[*,2]+0.10
cbpos[*,2]=cbpos[*,0]+0.01

dmax=max(frequency)
dmin=min(frequency)

dopvel=indgen(64)*0.1905

for i=0,3 do begin
freqt=frequency[*,*,i]
 
pnum=i
image_var=bytscl(freqt,min=min,max=dmax,top=top_color)

;  Set up the plot space
s=size(freqt,/dimensions)
contour,freqt,indgen(s[0]),indgen(s[1]),/nodata,$
  xstyle=4,ystyle=4,position=pos[pnum,*],/noerase

;  Grab the plot location
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1

;  Put up the image
tv,congrid(image_var,sx,sy),px[0],py[0]

;  Put up the axes
contour,freqt,dopvel,height[0,*]/1000.0,xstyle=1,ystyle=1,/nodata,$
  color=0,$
  ;xtickunits='time',xtickformat='label_date',$
  xtitle='doppler velocity m/s',$;spectral bin',$
  ytitle='height km',$;'range gate',$
  title='20'+string(year[i],format='(I02)')+string(month[i],format='(I02)')+string(day[i],format='(I02)')$
    +' '+string(hour[i],format='(I02)')+':'+string(minute[i],format='(I02)')+':'+string(seconds[i],format='(I02)'),$
  position=pos[pnum,*],/noerase
;xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'fall_velocity',$
;  color=255,/normal
;  Color bar
colorbar_fanning,maxrange=dmax,minrange=dmin,$
  ncolors=top_color,format='(F8.2)',$
  vertical=1,color=0,position=cbpos[pnum,*]
endfor
write_gif,'mrr.freqt.gif',tvrd()

;*********************************
;  CFAD frequency by altitude diagram
;*********************************
erase
xl=0.1  &  xr=0.85
yb=0.1  &  yt=0.8
sx=0.20
sy=0.15
numplots_x=1  ;number of plots in x direction
numplots_y=1  ;number of plots in y direction

position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
cbpos=pos
cbpos[*,0]=pos[*,2]+0.10
cbpos[*,2]=cbpos[*,0]+0.01

pnum=0

x_start_bin=-20
x_end_bin=25
x_dbin=0.1  ;half bin width

y_start_bin=0
y_end_bin=6000
y_dbin=50  ;half bin width

x_data=reflectivity
y_data=height

histogram_2d,x_data,y_data,x_start_bin,x_end_bin,x_dbin,x_bins,$
  y_start_bin,y_end_bin,y_dbin,y_bins,data_freq,data_counts

image_var=bytscl(data_counts,min=0,max=max(data_counts),top=top_color)
  
  ;  Set up the plot space
  s=size(data_counts,/dimensions)
  contour,data_counts,indgen(s[0]),indgen(s[1]),/nodata,$
    xstyle=4,ystyle=4,position=pos[pnum,*],/noerase

  ;  Grab the plot location
  px=!x.window*!d.x_vsize
  py=!y.window*!d.y_vsize
  sx=px[1]-px[0]+1
  sy=py[1]-py[0]+1

  ;  Put up the image
  tv,congrid(image_var,sx,sy),px[0],py[0]

  ;  Put up the axes
  contour,image_var,x_bins,y_bins,xstyle=1,ystyle=1,/nodata,$
    color=0,$
    ;xtickunits='time',xtickformat='label_date',$
    xtitle='Ze (dBz)',$;spectral bin',$
    ytitle='height (m)',$;'range gate',$
    ;title='20'+string(year[i],format='(I02)')+string(month[i],format='(I02)')+string(day[i],format='(I02)')$
    ;+' '+string(hour[i],format='(I02)')+':'+string(minute[i],format='(I02)')+':'+string(seconds[i],format='(I02)'),$
    position=pos[pnum,*],/noerase
  ;xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'fall_velocity',$
  ;  color=255,/normal
  ;  Color bar
  colorbar_fanning,maxrange=max(data_counts),minrange=0,$
    ncolors=top_color,format='(F8.2)',$
    vertical=1,color=0,position=cbpos[pnum,*]  
  
write_gif,'mrr.ave.cfad.gif',tvrd()  
  
stop

end