pro find_mixed_layers2, site, date, qilong_tau, tau_time, layers_flag


;This program uses the find_mixed_layers algorithm since in general, it
;was better. However, it failes where either iwp or lwp is missing. 
;If either of those are the case, this program uses the Layers/bright band
;algorithm instead.
;
;This new version also contains a minimum required value for IWP and LWP to
;recognize it. This should cut down on the number of false mixed phase
;readings we get.

;site = 'sgp'
;date = '20041027'

plot_flag = 1
;min_value = 

missing = -9999.

year=strmid(date,0,4)
month=strmid(date,4,2)
day=strmid(date,6,2)



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.temp.cdf'

cdf_file = arm_dir+data_dir+filename

; Get the Qilong Tau only if it is not 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


if site eq 'sgp' then begin
  input_path='/uufs/chpc.utah.edu/common/home/mace_grp/data/arm/sgp/sgp.average.5min/'+year+'/'
endif

get_liq_data2, input_path, site, date, long_day, lwp, iwp
;stop
; Find the day of the year
julian_day_calc,year,month,day,0,0,0,jday

numtimes = n_elements(long_day)
if numtimes eq 0 then return
julian_day = make_array(numtimes,/double)
layers_flag = make_array(numtimes,value=0)
for i = 0,numtimes-1 do begin
  ; 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 lwp[i] gt 0. and iwp[i] gt 0. then layers_flag[i] = 1
endfor

;******************************************************
; This is the new code for version 2 of this program.
missing_result = where(lwp eq missing or iwp eq missing, count)

if count gt 0 then begin
  find_multiple_layers, site, date, qilong_tau, tau_time, multiple_layers_flag
  result = where(multiple_layers_flag gt 1, count)
  if count gt 0 then multiple_layers_flag[result] = 1
  layers_flag[missing_result] = multiple_layers_flag[missing_result]
endif
;******************************************************

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

;************************************************************************
; Get the reflectivity data for plotting

if plot_flag eq 1 then begin

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

  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
  
  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
  ncdf_close, cdf_id
endif

;************************************************************************
; 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=5.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[4,*],$
  title='Optical depth calculated by Qilong for '+date,$
  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[3,*],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=2.
dmin=0.

plot,long_day,layers_flag,color=250,xtickformat='label_date',$
  position=pos[2,*],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
  
dmax1=max(lwp)
dmax2=max(iwp)
dmax=max([dmax1,dmax2])
dmin=0.
plot, long_day,lwp,color=0,xtickformat='label_date',$
  position=pos[1,*],title='LWP (black) and IWP (blue)',$
  yrange=[dmin,dmax],/noerase, xtickinterval=4,xminor=4,$
  xstyle=1,ystyle=1,xtickunits='hour',xticklen=-0.02,yticklen=-0.02

oplot,long_day,iwp,color=70

image_name = 'multi_phase_removal.'+site+'.'+date+'.gif'
write_gif, image_name, tvrd()

endif

;stop
return
end