pro plot_clouds_zyuying,input_flag
;
;  Common to hold data
;
common clouds_data, base_time, time_offset, height,$
	clutter,cbase,ctop,reflect,velocity,spectrum,$
        artifact
;
;  find date and site depending on input flag
;
input_flag=strcompress(input_flag,/remove_all)
if input_flag eq 'one' then begin
  site=' '
  read,site,prompt='Enter ARM site (nsa,sgp,twpc1,twpc2): '
  site=strcompress(site,/remove_all)
  curdate=' '
  read,curdate,prompt='Enter date (yyyymmdd): '
  curdate=strcompress(curdate,/remove_all)
endif
;
;  Format the date to process
;
print,'date  ',curdate,'site',site
curdate=strcompress(string(curdate),/remove_all)
date_string_fm_adtg, curdate,date_string
;
;  Get the data
;
get_clouds,curdate,site,datastream,input_flag
;
;  Transpose the data
;
reflect=transpose(reflect)
;
;  Get the length of the arrays
;
num_heights=n_elements(height)
num_times=n_elements(time_offset)

;******************************************
;  Look for data gaps and fill them in
;******************************************
;
;  Make arrays to hold the new expanded data
;
newreflect=make_array(86400,num_heights,/int,value=-9999)  ;new reflectivity array
newcbase=make_array(86400,/int,value=-9999)
newctop=make_array(86400,/int,value=-9999)
spacing=fltarr(num_times)  ;delta t between two successive time measurements
newtime=fltarr(86400)  ;new time array
;
;  find the missing time intervals and replace with missing data values
;
j=0  ;index for new time array
for i=0,num_times-2 do begin
  spacing[i+1]=time_offset[i+1]-time_offset[i] 
;  print,'********',time_offset[i],time_offset[i+1],spacing[i+1]
;  if spacing[i+1] ge 75 then begin
  if spacing[i+1] ge 150 then begin
    print, 'found a hole',time_offset[i],time_offset[i+1],spacing[i+1]
    x=fix(spacing[i+1]/36.0)
    if x*36.0+time_offset[i] eq time_offset[i+1] then x=x-1.0
      for k=1.0,x do begin
        j=j+1
        newtime[j]=time_offset[i]+k*36.0
        newreflect[j,*]=-8888
        newcbase[j]=-8888
        newctop[j]=-8888
        print,j,newtime[j],newreflect[j,0]
      endfor
	j=j+1
	newtime[j]=time_offset[i+1]
	newreflect[j,*]=reflect[i+1,*]
        newcbase[j]=cbase[i+1]
        newctop[j]=ctop[i+1]
        print,j,newtime[j],newreflect[j,0]
	print, 'filled in hole'
      endif else begin
	j=j+1
	newtime[j]=time_offset[i+1]
	newreflect[j,*]=reflect[i+1,*]
        newcbase[j]=cbase[i+1]
        newctop[j]=ctop[i+1]
        print,j,newtime[j],newreflect[j,0]
      endelse
endfor
;
;  cut the arrays down to their size
;
newtime[0]=time_offset[0]
time_offset=newtime[0:j]
reflect=newreflect[0:j,*] 
cbase=newcbase[0:j]
ctop=newctop[0:j]
;
;  Calculate fractional hour and Julian day
;
ijday_ihrfrac_fm_numsec,base_time,time_offset,hrfrac,jday,yy,mm,dd,hh,mi,ss
;
;  Calculate IDL Julian day for plotting
;
julian_day=make_array(n_elements(jday),/double,value=-9999)
julian_day=julday(mm,dd,yy,hh,mi,ss)
;
;  Change the format of reflectivity
;
if n_elements(reflect) gt 0 then  begin
  reflect=float(reflect)
  result=where(reflect gt -8888,count)
  if count ne 0 then reflect[result]=reflect[result]/100.
endif
;
;  Convert height to km
;
height=height/1000.
;
;  Determine the subsetting indicies
;
;***************************
;  Enter start time below
;***************************
stime=14
closetime=hrfrac-stime
closest=min(abs(closetime),indst)
stime=hrfrac[indst]
;***************************
;  Enter end time below
;**************************
etime=22
closetime=hrfrac-etime
closest=min(abs(closetime),indet)
etime=hrfrac[indet]
;***************************
;  Enter start height below
;***************************
sheight=5
closeheight=height-sheight
closest=min(abs(closeheight),indsh)
sheight=height[indsh]
;***************************
;  Enter end height below
;**************************
eheight=11
closeheight=height-eheight
closest=min(abs(closeheight),indeh)
eheight=height[indeh]

reflect=reflect[indst:indet,indsh:indeh]
julian_day=julian_day[indst:indet]
height=height[indsh:indeh]
cbase=cbase[indst:indet]
ctop=ctop[indst:indet]
hh=hh[indst:indet]
mi=mi[indst:indet]
;
;  Set up the plot
;
!p.background=!d.n_colors-1
loadct,5
!p.charsize=1.8
window, 0, xsize=600, ysize=400, TITLE='clouds',/pixmap
!p.multi = [0,1,1]
!x.range=[julian_day[0],julian_day[n_elements(julian_day)-1]]
!x.style=1
!y.style=1
dummy=label_date(date_format=['%H:%I'])
pos1=[0.1,0.15,0.85,0.8]

;
;  Set up the colors
;  0=black  130=green  25=purple  50=blue  250=red
;
common colors,r_orig,g_orig,b_orig,r_curr,g_curr,b_curr

;  make gray for missing data in color 253
r_curr[!d.n_colors-2]=200   ;r_curr[254]=200
g_curr[!d.n_colors-2]=200   ;g_curr[254]=200
b_curr[!d.n_colors-2]=200   ;b_curr[254]=200

color1=0  ;black
color2=250  ;red
;
;  Plot 1  -  Reflectivity
;

;
;  Set up the plot coordinates
;
contour,reflect,julian_day,height,/nodata,xstyle=4,ystyle=4,$
  position=pos1
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
;
;  Scale and create the image
;
dmax=10.
dmin=-50.
image=bytscl(reflect,max=dmax,min=dmin,top=!d.n_colors-3)
;
;  Put the gray in for missing data
;
result=where(reflect eq -8888,count)
if count gt 0 then image[result]=!d.n_colors-2  ;254
;
;  Put white in for no cloud
;
result=where(reflect eq -9999,count)
if count gt 0 then image[result]=!d.n_colors-1  ;255

tv,congrid(image,sx,sy),px[0],py[0]
;
;  Set up the axes
;
contour,reflect,julian_day,height,/noerase,/nodata,$
  xtickformat='label_date',xtickunits='time',$
  color=color1,xminor=6,xticklen=-0.02,$
  yrange=[height[0],height[n_elements(height)-1]],$
  ytitle='Height (km)',position=pos1,xtitle='Time (UTC)'
colorbar, image,dmax,dmin,'DbZe',px,py,sx,sy
;
;  Add title to plot page
;
if site eq 'twpc1' then capsite='TWP C1'
if site eq 'twpc2' then capsite='TWP C2'
if site eq 'nsa' then capsite='NSA'
if site eq 'sgp' then capsite='SGP'
a_string=capsite+' Merged Moments, '+date_string
;;xyouts, 0.5, 0.88, alignment=0.5,a_string, /normal, color=0, charsize=1.5
;
;  Create the gif file
;

write_gif,'/data/mace4/arm_data/temp_image/'+$
  datastream+'.'+curdate+'_zoom.gif',tvrd()

stop
end
;****************************************************
;  colorbar.pro creates a colorbar for the DBZ images
;****************************************************

pro colorbar, image,dmax,dmin,cb_title,px,py,sx,sy

min_image=0;min(image)
max_image=!d.n_colors-3  ;253;max(image)

cb_arr=fltarr(2,max_image-min_image+1)
for j=0,max_image-min_image do $
  cb_arr(*,j)=dmin+((dmax-dmin)/(max_image-min_image))*j

cb_img=bindgen(2,max_image-min_image+1)

for j=0,max_image-min_image do cb_img(*,j)=min_image+j

tv, congrid(cb_img,(sx/35.),sy),px[1]+sx/7,py[0]
axis,px[1]+sx/7,py[0],yaxis=0,color=0,$
  yrange=[cb_arr[0,0],cb_arr[0,max_image-min_image]],$
  ystyle=1,/device,charsize=1.2,ytitle=cb_title,ticklen=-0.02

return 
end