pro find_multiple_layers, site, date, qilong_tau, tau_time, layers_flag, tau_stretch

;layers_flag = 0 if the algorithm only detects one layer.
;layers_flag = 1 if multiple layers are detected with at least XXX
;                height increments between them - determined by the
;                preset variable 'num_spaces'
;layers_flag = 2 if a bright band is detected.

;site = 'sgp'
;date = '20040305'

missing = -9999.
bright_band_flag = 0;Do you want to go looking for a bright band?

plot_flag = 0; Do you want diagnostic plots?

dbz_threshold = -60.; Minimum reflectivity to be considered a cloud
num_spaces = 2; Minimum number of layers that must have no cloud between
              ; recognized clouds.
min_levels = 2; Minimum number of levels that must contain data for
              ; the program to assume a cloud is present. NOTE:
	      ; If this value is changed, the code will have to be
	      ; edited!!
	      
;*****************************************************

sep_date, date, year, month, day
;get qilong tau ONLY if if isn't inputted
if n_elements(qilong_tau) eq 0 then begin
  input_path='/uufs/chpc.utah.edu/common/home/mace_grp/data/arm/sgp/qilong_tau/'+year+'/'
  get_qilong_tau, input_path, date, day, tau_time, qilong_tau, qilong_tau_rms, month
endif

arm_dir = '/uufs/chpc.utah.edu/common/home/mace_grp/data/arm/'
data_dir = site+'/'+site+'.average.5min/'+year+'/'
filename = site+'.average.5min.'+date+'.000000.cdf'

cdf_file = arm_dir+data_dir+filename

spawn, 'gunzip '+cdf_file+'.gz'
cdf_id=ncdf_open(cdf_file)

; Get time and convert it to match tau_time's format:
base_time_id = ncdf_varid(cdf_id,'base_time')
time_offset_id = ncdf_varid(cdf_id,'time_offset')
ncdf_varget, cdf_id, base_time_id, base_time
ncdf_varget, cdf_id, time_offset_id, time_offset

epoch_zero = double(julday(1,1,1970,0,0,0))
;double_base_time = double(base_time)
long_day = (time_offset + base_time)/86400. + epoch_zero

;************************

reflectivity_id = ncdf_varid(cdf_id,'ReflectivityBestEstimate'); dimensions: time,height
if reflectivity_id ne -1 then begin
  ncdf_varget,cdf_id, reflectivity_id, reflectivity
  result=where(reflectivity eq missing,count)
  reflectivity = float(reflectivity)/100.
  if count gt 0 then begin
    reflectivity[result]=missing
  endif
endif else begin
  return
endelse

height_id = ncdf_varid(cdf_id,'height')
if height_id ne -1 then begin
  ncdf_varget, cdf_id, height_id, height
  height = height/1000.; meters to km
endif else begin
  return
endelse

numtimes = n_elements(long_day)
numheights = n_elements(height)
layers_flag = make_array(numtimes,value=0)
cloud_count = make_array(numtimes)
julian_day = make_array(numtimes,/double)

; Find the day of the year
julian_day_calc,year,month,day,0,0,0,jday 


ncdf_close,cdf_id
;************************
;Find spaces
;dummy_flag = 0
for i = 0,numtimes-1 do begin
  refl_column = reflectivity[*,i]
  j_index = where(refl_column gt dbz_threshold,dummy_count)
  cloud_count[i] = dummy_count
  if min_levels eq 1 then begin
    if cloud_count[i] gt 1 then begin
      for j=1,n_elements(j_index)-1 do begin
        if j_index[j]-j_index[j-1] gt num_spaces then begin
	  layers_flag[i] = 1
	  break
	endif
      endfor
    endif
  endif else if min_levels eq 2 then begin
    if cloud_count[i] gt 3 then begin
      for j=3,n_elements(j_index)-1 do begin
        if j_index[j]-j_index[j-1] eq 1 and $
	       j_index[j-1]-j_index[j-2] gt num_spaces and $
	       j_index[j-2]-j_index[j-3] eq 1 then begin
	  layers_flag[i] = 1
	  break
        endif
      endfor
    endif
  endif else if min_levels eq 3 then begin
    if cloud_count[i] gt 5 then begin  
      for j=5,n_elements(j_index)-1 do begin
        if j_index[j]-j_index[j-2] eq 2 and $
	       j_index[j-2]-j_index[j-3] gt num_spaces and $
	       j_index[j-3]-j_index[j-5] eq 3 then begin
	  layers_flag[i] = 1
	  break
	endif
      endfor
    endif
  endif; Code will have to be added here if min_levels is changed.
  
  
;************************
;Find bright bands  
  if cloud_count[i] gt 1 then begin
    do_this = 1
    if layers_flag[i] gt 0 then do_this = 0; if this is the case, then a space was already found, so we don't
                                           ; need to bother finding a bright band
  ;Note that we are still in the loop through time and if there is one or more clouds in this column.
  ;I put this in a subroutine just in case I messed it up.
    if do_this eq 1 then begin
      if bright_band_flag eq 1 then begin
        find_bright_band,refl_column,bright_band_index
        if bright_band_index ge 0 then layers_flag[i] = 3
      endif
    endif
  endif; end if there are reflectivity values in the column (if cloud_count[i] gt 0)
  ;if i eq 170 then stop
  
  ; If a bright band is detected in one place, then likely the entire cloud is mixed-phased
;  do_this = 0; Forget the code below. It screws up.
;  if layers_flag[i] eq 2 and do_this eq 1 then begin
;    k = i-1
;    while k ge 0 do begin
;      if cloud_count[k] gt 0 and layers_flag[k] eq 0 then begin
;	layers_flag[k] = 2
;      endif else begin
;	break
;      endelse
;      k = k-1
;    endwhile
;    k = i+1
;    while k lt numtimes do begin
;      dummy_variable = where(refl_column gt dbz_threshold,another_dummy_count)
;      cloud_count[k] = another_dummy_count
;      if cloud_count[k] gt 0 then begin
;	layers_flag[k] = 2
;      endif else begin
;        break
;      endelse
;      k = k+1
;    endwhile
;  endif; Do this code instead
  do_this = 1
  ; If a bright band is detected, flag its close neighbors as mixed phase as well.
  if layers_flag[i] eq 3 and do_this eq 1 then begin
    for k = 1,10 do begin
      if i-k ge 0 then begin
        if layers_flag[i-k] eq 0 then layers_flag[i-k] = 2
      endif
      if i+k lt numtimes then begin
        if layers_flag[i+k] eq 0 then layers_flag[i+k] = 2
      endif
    endfor
  end
 ;************************ 
 ; find the julian day for this time incriment to be later matched with the time format of qilong tau
  caldat, long_day[i], month, day, year, hour, minute, second
  total_seconds = (double(second)) + (double(minute))*(double(60.)) + (double(hour))*(double(3600.))
  day_fraction = total_seconds / (double(86400.))
  julian_day[i] = double(jday) + day_fraction
 ;************************ 
; if i eq 90 then stop
endfor; end loop through time
result=where(layers_flag eq 3, count)
if count gt 0 then layers_flag[result]=2
;reflectivity[*,90]=20
;reflectivity[16,*]=20
;************************
; Now set qilong tau to missing where appropriate.
;close_enough = (2.5/60.)/24.; 2.5 minutes

numtau = n_elements(qilong_tau)
tau_stretch = make_array(numtimes,value=missing)
old_tau_stretch = tau_stretch
old_tau = qilong_tau
for i=0,numtau-1 do begin
  time_difference = abs(julian_day - tau_time[i])
  min_difference = min(time_difference)
  julian_index = where(time_difference eq min_difference,count)
  if count gt 1 then julian_index = julian_index[0] 
  if plot_flag eq 1 then begin
    old_tau_stretch[julian_index] = old_tau[i]
  endif
  if layers_flag[julian_index] gt 0 then qilong_tau[i] = missing
  tau_stretch[julian_index] = qilong_tau[i]
endfor

;************************************************************************
; Create the diagnostic plots (if desired)

if plot_flag eq 1 then begin
;  Size of the plot window
xwin=800  &  ywin=800  

;  Set up the plotting device
set_plot,'z'
device,set_resolution=[xwin,ywin]
;set_plot,'ps'
;device,/inches,/times,/bold,bits_per_pixel=8,$
; /color,xoffset=0,yoffset=0,xsize=11,ysize=11.1,$
; filename='nsa.'+date+'.eps',/portrait,/encapsulated

;dummy=label_date(date_format=['%H:%I'])
;  Load the color table
loadct,39   

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

;  Date label format
dummy=label_date(date_format=[' %D %M !C %H:%I '])

;  Increase the character size
!p.charsize=0.8  

;  xticks point out from the plot
!x.ticklen=-0.02

;  chose the psym
sys=4

;  Choose the psym size
psz=0.2

;*************************************
;  Position the plots
;*************************************

xl=0.081  &  xr=0.88
yb=0.05  &  yt=0.88
sx=0.07
sy=0.07
numplots_x=1.0
numplots_y=4.0

top_color = 254


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

;*************************************
;  Plot the mmcr reflectivity
;*************************************

dmax = max(reflectivity)
result = where(reflectivity ne missing, count)
if count gt 0 then begin
  dmin = min(reflectivity[result])
endif else begin
  dmin = 0
endelse


reflectivity=transpose(reflectivity)

image=bytscl(reflectivity,max=dmax,min=dmin,top=top_color)
result=where(reflectivity eq missing, count)
if count gt 0 then image[result] = 255

contour,reflectivity,long_day,height,/nodata,xstyle=4,ystyle=4,$
 position=pos[0,*]
px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
sx=px[1]-px[0]+1
sy=py[1]-py[0]+1
tv,congrid(image,sx,sy),px[0],py[0],/device
;  Image
; tv,congrid(image,sx/20.0,sy/20.0),px[0],py[0],$
;  xsize=sx,ysize=sy,/device
;**********************
; Labels and Axes:
;**********************
contour,reflectivity,long_day,height,/noerase,/nodata,$
 xtickformat='label_date',xstyle=1,ystyle=1,$
 xtickunits='hour',xticklen=-0.02,yticklen=-0.02,$
 color=0,title='MMCR Reflectivity', xtickinterval=4,$
 ytitle='Height (km) AGL',position=pos[0,*],xminor=4

;oplot,julian_day,preliminary_cbh,color=255,psym=4,symsize=0.3

;**********************
;  Plot Colorbar
;**********************

colorbar_fanning,maxrange=dmax,minrange=dmin,$
 title='dBz',ncolors=top_color,format='(F6.2)',$
 vertical=1,color=0,position=cbpos[0,*], font=0;, right=1
 
 ;**********************
dmax = max(old_tau)
result = where(old_tau ne missing, count)
if count gt 0 then begin
  dmin = min(old_tau[result])
endif else begin
  dmin = missing
endelse

plot,long_day,old_tau_stretch,color=0,xtickformat='label_date',$
  position=pos[3,*],title='Optical depth calculated by Qilong',$
  yrange=[dmin,dmax],/noerase, xtickinterval=4,xminor=4,$
  xstyle=1,ystyle=1,xtickunits='hour',xticklen=-0.02,yticklen=-0.02
  
;**********************
dmax = max(qilong_tau)
result = where(qilong_tau ne missing, count)
if count gt 0 then begin
  dmin = min(qilong_tau[result])
endif else begin
  dmin = missing
endelse

plot,long_day,tau_stretch,color=0,xtickformat='label_date',$
  position=pos[2,*],title='Qilong tau after multi-phase removal',$
  yrange=[dmin,dmax],/noerase, xtickinterval=4,xminor=4,$
  xstyle=1,ystyle=1,xtickunits='hour',xticklen=-0.02,yticklen=-0.02
;**********************
dmax=3.
dmin=0.

plot,long_day,layers_flag,color=250,xtickformat='label_date',$
  position=pos[1,*],title='Layers Flag',psym=4,$
  yrange=[dmin,dmax],/noerase, xtickinterval=4,xminor=4,$
  xstyle=1,ystyle=1,xtickunits='hour',xticklen=-0.02,yticklen=-0.02
  
image_name = 'multi_phase_removal.'+site+'.'+date+'.gif'
write_gif, image_name, tvrd()

endif
;stop
return
end