pro plot_clouds_baum;,input_flag
input_flag='one'
;
;  Common to hold data
;
common clouds_data, base_time, time_offset, height,$
	clutter,cbase,ctop,reflect,velocity,spectrum,artifact
common vceil25k_data, base_timev, time_offsetv, $
	range,detsta,statflg,cbh1,cbh2,cbh3,highsig,lat,lon,alt
;
;  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
get_vceil25k, 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)
;
;  Convert merged moments height from MSL to AGL
;
height=height-318.

;******************************************
;  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 120 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
ijday_ihrfrac_fm_numsec,base_timev,time_offsetv,hrfracv,jdayv,yyv,mmv,ddv,hhv,miv,ssv
;
;  Calculate IDL Julian day for plotting
;
julian_day=make_array(n_elements(jday),/double,value=-9999)
julian_dayv=make_array(n_elements(jdayv),/double,value=-9999)
julian_day=julday(mm,dd,yy,hh,mi,ss)
julian_dayv=julday(mmv,ddv,yyv,hhv,miv,ssv)
;
;  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=3
closetime=hrfrac-stime
closetimev=hrfracv-stime
closest=min(abs(closetime),indst)
closestv=min(abs(closetimev),indstv)
stime=hrfrac[indst]
stimev=hrfracv[indstv]
;***************************
;  Enter end time below
;**************************
etime=6
closetime=hrfrac-etime
closetimev=hrfracv-etime
closest=min(abs(closetime),indet)
closestv=min(abs(closetimev),indetv)
etime=hrfrac[indet]
etimev=hrfracv[indetv]
;***************************
;  Enter start height below
;***************************
sheight=0
closeheight=height-sheight
closest=min(abs(closeheight),indsh)
sheight=height[indsh]
;***************************
;  Enter end height below
;**************************
eheight=15
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]
cbh1=cbh1[indstv:indetv]
cbh2=cbh2[indstv:indetv]
cbh3=cbh3[indstv:indetv]
julian_dayv=julian_dayv[indstv:indetv]
;
;  Set up the plot
;
imagedir='/data/mace4/arm_data/temp_image/'
;plot_type='postscript'  ;postscript
plot_type='gif'
if plot_type eq 'gif' then begin
  set_plot,'X'
  !p.font=-1
  !p.charsize=1.3
  window, 0, xsize=600, ysize=400, TITLE='clouds',/pixmap
  pos1=[0.1,0.15,0.85,0.8]
endif else begin
  !p.font=0
  set_plot,'ps'
  psfilename=imagedir+'mace_baum.20010422.ps'
  device,/inches,xoffset=1.0,yoffset=10.0,/times,/bold,xsize=8,$
    ysize=6,filename=psfilename,/color,bits_per_pixel=8,$
    /landscape
  !p.charsize=1.3
  pos1=[0.1,0.15,0.85,0.8]
endelse

!p.background=!d.n_colors-1
loadct,5
!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'])

;
;  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[254]=200
g_curr[254]=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
;
image=bytscl(reflect,max=25.0,min=-65.0,top=253)
;
;  Put the gray in for missing data
;
result=where(reflect eq -8888,count)
if count gt 0 then image[result]=254
;
;  Put white in for no cloud
;
result=where(reflect eq -9999,count)
if count gt 0 then image[result]=255
;
;  display the image
;
if plot_type eq 'postscript' then image2=congrid(image,sx/20,sy/20) else $
image2=congrid(image,sx,sy)
tv,image2,px[0],py[0],xsize=sx,ysize=sy,/device
;
;  Set up the axes
;
contour,image,julian_day,height,/noerase,/nodata,$
  xtickformat='label_date',xtickunits='time',$
  color=color1,xminor=6,xticklen=-0.02,xtitle='UTC',$
  yrange=[height[0],height[n_elements(height)-1]],$
  ytitle='Height AGL (km)',position=pos1

;  only plot every 5th dot
num_plot=fix(n_elements(cbh1)/5.0)
new_cbh1=fltarr(num_plot)
new_jday=dblarr(num_plot)
for i=0,num_plot-1 do begin
  new_cbh1[i]=cbh1[i*5]
  new_jday[i]=julian_dayv[i*5]
endfor

;oplot,julian_dayv,cbh1/1000.,color=color1,psym=4,symsize=0.5
oplot,new_jday,new_cbh1/1000.,color=color1,psym=4,symsize=0.5
;oplot,julian_dayv,cbh2/1000.,color=254,psym=4,symsize=0.5
;oplot,julian_dayv,cbh2/1000.,color=color1,psym=2
;
;  Create a colorbar
;
colorbar,image,25.0,-65.0,'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.9, alignment=0.5,a_string, /normal, color=0, charsize=1.5
;
;  Create the gif file
;


if plot_type eq 'gif' then begin
  write_gif,imagedir+$
    datastream+'.'+curdate+'_baum.gif',tvrd()
  print,imagedir+datastream+'.'+curdate+'_baum.gif'

endif else begin
  device,/close
endelse


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=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

plot_type='gif'
if plot_type eq 'gif' then begin
  tv, congrid(cb_img,(sx/(35.)),sy),px[1]+sx/7,py[0],xsize=sx/35.,ysize=sy
  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
endif else begin
  tv, congrid(cb_img,(sx/(35.*20.)),sy/20.),px[1]+sx/6,py[0],xsize=sx/35.,ysize=sy
  axis,px[1]+sx/6,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
endelse


return 
end