pro read_mrr_ave_txt_sub, file_name, file_path, reflectivity_attenuated,julian_day,out_height

;  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=file_path+file_name  ;'/Users/u0029340/Documents/data/Capricorn/MRR_DATA/AveData/201603/0314.ave'
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,2000.+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)

reflectivity_attenuated=transpose(reflectivity_attenuated)
reflectivity=transpose(reflectivity)
height=transpose(height)
rain_rate=transpose(rain_rate)
liquid_water_content=transpose(liquid_water_content)
fall_velocity=transpose(fall_velocity)

out_height=height[0,*]

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

;  Color common
;common colors,r_orig,g_orig,b_orig,r_curr,g_curr,b_curr
;
;
;loadct, 11 ;  This is the stern color scheme
;
;
;r_curr[250:255]=255 & r_orig=r_curr
;g_curr[250:255]=255 & g_orig=g_curr
;b_curr[250:255]=255 & b_orig=b_curr
;
;gamma_ct, 1.0
;set_plot, 'Z'
;device, set_resolution=[3000,2000]
;!p.multi = [0,3,6]
;!p.charsize=3.
;!p.font=-1
;!p.background=!d.n_colors
;
;;  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  
;
;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,'/Users/u0029340/Documents/data/Capricorn/MRR_DATA/mrr.timeseries.gif',tvrd()
;
;plot_spectra_flag=0
;if plot_spectra_flag gt 0 then begin
;;***********************************************************
;;  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,'/Users/u0029340/Documents/data/Capricorn/MRR_DATA/mrr.freqt.gif',tvrd()
;
;endif ;if plot_spectra_flag gt 0 then begin


;;*********************************
;;  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,'/Users/u0029340/Documents/data/Capricorn/MRR_DATA/mrr.ave.cfad.gif',tvrd()  
  
return

end