pro cfad_file_write, cfad_file_path, radar_file_prefix, radar_hrfrac, radar_heights, temp_sounding, heights_sounding, dbz, radar_mask, lidar_hrfrac, lidar_cldbase, mwr_hrfrac, mwr_liq,fluxes_hrfrac, latitude, v_wind, pressure

; does the cfad file exist?  Create a new CFAD file for each day.  Then files can be added as desired.

file_name=radar_file_prefix+'_cfad.cdf'

test=file_search(cfad_file_path+file_name, count=count_files)
if count_files eq 0 then begin
  
; file does not exist.  Create it.
; 
;create the dimensions of the cfad.  dbz_bins, height_bins, time_bins, synoptic_bins, cld_base_temp_bins, lwp_bins, latitude_bins
   profile_count=0l & clear_profile_count=0l

   dbz_bins=-45. & while max(dbz_bins) le 20. do dbz_bins=[dbz_bins, max(dbz_bins)+3.]
   height_bins=0.0 & while max(height_bins) le 12. do height_bins=[height_bins, max(height_bins)+0.100]
   time_bins=0.0 & while max(time_bins) lt 24. do time_bins=[time_bins, max(time_bins)+3.]
   latitude_bins=-55. & while max(latitude_bins) lt 45. do latitude_bins=[latitude_bins, max(latitude_bins)+2.5]
   synoptic_bins=[-1., 0., 1.]
   cld_base_temperature_bins=-15. & while max(cld_base_temperature_bins) lt 20. do cld_base_temperature_bins=[cld_base_temperature_bins, max(cld_base_temperature_bins)+1.]
   
   cfad_array=lonarr(n_elements(dbz_bins), n_elements(height_bins),n_elements(time_bins), n_elements(latitude_bins),n_elements(synoptic_bins), n_elements(cld_base_temperature_bins))
  
endif else begin

; file exists.  Read in the data to add then rewrite
cdfid=ncdf_open(cfad_file_path+radar_file_prefix+'_cfad.cdf')
 
cfad_array_id=ncdf_varid(cdfid, 'BASTA_CFAD') & ncdf_varget, cdfid, cfad_array_id, cfad_array
dbz_bins_vector_id=ncdf_varid(cdfid, 'dbz_bins_vector') & ncdf_varget, cdfid, dbz_bins_vector_id, dbz_bins
height_bins_vector_id=ncdf_varid(cdfid, 'height_bins_vector') & ncdf_varget, cdfid, height_bins_vector_id, height_bins
time_bins_vector_id=ncdf_varid(cdfid, 'time_bins_vector') & ncdf_varget, cdfid, time_bins_vector_id, time_bins
latitude_bins_vector_id=ncdf_varid(cdfid, 'latitude_bins_vector') & ncdf_varget, cdfid, latitude_bins_vector_id, latitude_bins
synoptic_bins_vector_id=ncdf_varid(cdfid, 'synoptic_bins_vector') & ncdf_varget, cdfid, synoptic_bins_vector_id, synoptic_bins
cld_base_temp_vector_id=ncdf_varid(cdfid, 'cld_base_temp_vector') & ncdf_varget, cdfid, cld_base_temp_vector_id, cld_base_temperature_bins
profile_count_id=ncdf_varid(cdfid, 'profile_number') &  ncdf_varget, cdfid, profile_count_id, profile_count
cloudless_profile_count_id=ncdf_varid(cdfid, 'cloudless_profile_count') &  ncdf_varget, cdfid, cloudless_profile_count_id, clear_profile_count


ncdf_close, cdfid


endelse ; file written

; populate the arrays
;
;radar_hrfrac, radar_heights, temp_sounding, heights_sounding, dbz, vel_cor, lidar_hrfrac, lidar_cldbase, mwr_hrfrac, mwr_liq,fluxes_hrfrac, latitude, v_wind, pressure

;
;dbz_bins=-45. & while max(dbz_bins) le 20. do dbz_bins=[dbz_bins, max(dbz_bins)+3.]
;height_bins=0.0 & while max(height_bins) le 12. do height_bins=[height_bins, max(height_bins)+0.250
;time_bins=0.0 & while max(time_bins) lt 24. do height_bins=[height_bins, max(height_bins)+3.
;latitude_bins=-55. & while max(latitude_bins) lt 45. do latitude_bins=[latitude_bins, max(latitude_bins)+2.5
;synoptic_bins=[-1., 0., 1.]
;cld_base_temperature_bins=-15. & while max(cld_base_temperature_bins) lt 20. do cld_base_temperature_bins=[cld_base_temperature_bins, max(cld_base_temperature_bins)+1.5
;
;cfad_array=lonarr(n_elements(dbz_bins), n_elements(height_bins),n_elements(time_bins), n_elements(latitude_bins),n_elements(synoptic_bins), n_elements(cld_base_temperature_bins))
;


for j=0,n_elements(radar_hrfrac)-1 do begin
  time_index=0 & latitude_index=0 & synoptic_index=0 & cloud_base_temp_index=0
  profile_count=profile_count+1l
  cloud_flag=0

  lidar_index=-1 & lidar_index=where(abs(lidar_hrfrac-radar_hrfrac[j]) eq min(abs(lidar_hrfrac-radar_hrfrac[j])) and abs(lidar_hrfrac-radar_hrfrac[j]) lt 0.01)
  mwr_index=-1 & mwr_index=where(abs(mwr_hrfrac-radar_hrfrac[j]) eq min(abs(mwr_hrfrac-radar_hrfrac[j])) and abs(mwr_hrfrac-radar_hrfrac[j]) lt 0.02)
  fluxes_index=-1 & fluxes_index=where(abs(fluxes_hrfrac-radar_hrfrac[j]) eq min(abs(fluxes_hrfrac-radar_hrfrac[j])) and abs(fluxes_hrfrac-radar_hrfrac[j]) lt 0.02)

; synoptic_bin

synoptic_bin=1

if fluxes_index[0] gt 5 and fluxes_index[0] lt n_elements(fluxes_hrfrac)-6 then begin
if mean(pressure[fluxes_index[0]:fluxes_index[0]+5])-mean(pressure[fluxes_index[0]-5:fluxes_index[0]]) lt 0.1 or mean(v_wind[fluxes_index[0]-5:fluxes_index[0]+5]) le -2. then begin
; falling pressure and northernly winds - prefrontal
synoptic_index=0
endif
if mean(pressure[fluxes_index[0]:fluxes_index[0]+5])-mean(pressure[fluxes_index[0]-5:fluxes_index[0]]) gt 0.1 or mean(v_wind[fluxes_index[0]-5:fluxes_index[0]+5]) ge 2. then begin
  ; rising pressure and southerly winds - postfrontal
  synoptic_index=2
endif
  
endif ; if (fluxes_index[0] gt 5 and fluxes_index[0] lt n_elements(fluxes_index[0])-6 do begin


;;;;;;;;;;;;;;;;; done with synoptic bin
; time bin
time_index=where(abs(radar_hrfrac[j]-time_bins) eq min(abs(radar_hrfrac[j]-time_bins) ))

; latitude index
latitude_index=where(abs(latitude[fluxes_index[0]]-latitude_bins) eq min(abs(latitude[fluxes_index[0]]-latitude_bins) ))

;cloud base temperature
cloud_base_index=where(abs(lidar_cldbase[lidar_index[0]]-heights_sounding) eq min(abs(lidar_cldbase[lidar_index[0]]-heights_sounding) ))
cloud_base_temp_index=where(abs(temp_sounding[cloud_base_index[0]]-cld_base_temperature_bins) eq min(abs(temp_sounding[cloud_base_index[0]]-cld_base_temperature_bins) ))

cloud_base_index=where(abs(lidar_cldbase[lidar_index[0]]-radar_heights) eq min(abs(lidar_cldbase[lidar_index[0]]-radar_heights) ))


for k=4,n_elements(radar_heights)-1 do begin
dbz_index=0 & height_index=0
height_index=where(abs(radar_heights[k]-height_bins) eq min(abs(radar_heights[k]-height_bins) ))

if dbz[j,k] ge -30. and radar_mask[j,k] eq 1 then begin
  
dbz_index= where(abs(dbz[j,k]-dbz_bins) eq min(abs(dbz[j,k]-dbz_bins) )) 
cloud_flag=1
  
endif else begin  ;if dbz[j,k] gt -30. then begin

dbz_index=0
if cloud_base_index[0] eq k then begin
  dbz_index= where(abs(-35.-dbz_bins) eq min(abs(-35-dbz_bins) ))
  cloud_flag=1
endif


endelse ; if dbz[j,k] gt -30. then begin

;cfad_array=lonarr(n_elements(dbz_bins), n_elements(height_bins),n_elements(time_bins), n_elements(latitude_bins),n_elements(synoptic_bins), n_elements(cld_base_temperature_bins))

cfad_array[dbz_index[0], height_index[0],time_index[0],latitude_index[0],synoptic_index[0],cloud_base_temp_index[0]]=$
  cfad_array[dbz_index[0], height_index[0],time_index[0],latitude_index[0],synoptic_index[0],cloud_base_temp_index[0]]+1l

 
endfor  ;for k=0,n_elements(radar_heights)-1 do begin
  
  if cloud_flag eq 0 then clear_profile_count=clear_profile_count+1l

endfor  ; for j=0,n_elements(radar_hrfrac)-1 do begin
  
;dbz_bins=-45. & while max(dbz_bins) le 20. do dbz_bins=[dbz_bins, max(dbz_bins)+3.]
;height_bins=0.0 & while max(height_bins) le 12. do height_bins=[height_bins, max(height_bins)+0.250
;time_bins=0.0 & while max(time_bins) lt 24. do height_bins=[height_bins, max(height_bins)+3.
;latitude_bins=-55. & while max(latitude_bins) lt 45. do latitude_bins=[latitude_bins, max(latitude_bins)+2.5
;synoptic_bins=[-1., 0., 1.]
;cld_base_temperature_bins=-15. & while max(cld_base_temperature_bins) lt 20. do cld_base_temperature_bins=[cld_base_temperature_bins, max(cld_base_temperature_bins)+1.5
;
;cfad_array=lonarr(n_elements(dbz_bins), n_elements(height_bins),n_elements(time_bins), n_elements(latitude_bins),n_elements(synoptic_bins), n_elements(cld_base_temperature_bins))
 
cdfid=ncdf_create(cfad_file_path+radar_file_prefix+'_cfad.cdf', /clobber)

print, 'created file ',cfad_file_path+radar_file_prefix+'_cfad.cdf'

dbz_bins_did=ncdf_dimdef(cdfid,'dbz_bins', n_elements(dbz_bins))
height_bins_did=ncdf_dimdef(cdfid,'height_bins',  n_elements(height_bins))
time_bins_did=ncdf_dimdef(cdfid,'time_bins',  n_elements(time_bins))
latitude_bins_did=ncdf_dimdef(cdfid,'latitude_bins',  n_elements(latitude_bins))
synoptic_bins_did=ncdf_dimdef(cdfid,'synoptic_bins',  3)
cld_base_temperature_bins_did=ncdf_dimdef(cdfid,'cld_base_temperature_bins',  n_elements(cld_base_temperature_bins))

profile_count_id=ncdf_vardef(cdfid, 'profile_number',  /long)
ncdf_attput, cdfid, profile_count_id, 'long_name','total number of profiles'
ncdf_attput, cdfid, profile_count_id, 'units','none - counts'

cloudless_profile_count_id=ncdf_vardef(cdfid, 'cloudless_profile_count',  /long)
ncdf_attput, cdfid, cloudless_profile_count_id, 'long_name','total number of cloud-free profiles'
ncdf_attput, cdfid, cloudless_profile_count_id, 'units','none - counts'

dbz_bins_vector_id=ncdf_vardef(cdfid, 'dbz_bins_vector', [dbz_bins_did], /float)
ncdf_attput, cdfid, dbz_bins_vector_id, 'long_name','dbz values of the cfad dbz dimension'
ncdf_attput, cdfid, dbz_bins_vector_id, 'units','dbz mm6/m3'

height_bins_vector_id=ncdf_vardef(cdfid, 'height_bins_vector', [height_bins_did], /float)
ncdf_attput, cdfid, height_bins_vector_id, 'long_name','Height values of the height dimension'
ncdf_attput, cdfid, height_bins_vector_id, 'units','height above the ocean in km'

time_bins_vector_id=ncdf_vardef(cdfid, 'time_bins_vector', [time_bins_did], /float)
ncdf_attput, cdfid, time_bins_vector_id, 'long_name','hour values of the time dimension'
ncdf_attput, cdfid, time_bins_vector_id, 'units','hour after UTC midnight'

latitude_bins_vector_id=ncdf_vardef(cdfid, 'latitude_bins_vector', [latitude_bins_did], /float)
ncdf_attput, cdfid, latitude_bins_vector_id, 'long_name','latitude bin values of the latitude dimension'
ncdf_attput, cdfid, latitude_bins_vector_id, 'units','degrees north latitude'

synoptic_bins_vector_id=ncdf_vardef(cdfid, 'synoptic_bins_vector', [synoptic_bins_did], /float)
ncdf_attput, cdfid, synoptic_bins_vector_id, 'long_name','rough estimate of synoptic state of the time of measurment'
ncdf_attput, cdfid, synoptic_bins_vector_id, 'units','index.  -1 means prefrontal - pressure falling and northerly surface wind, 1 post frontal, 0  indeterminant'

cld_base_temp_vector_id=ncdf_vardef(cdfid, 'cld_base_temp_vector', [cld_base_temperature_bins_did], /float)
ncdf_attput, cdfid, cld_base_temp_vector_id, 'long_name','temperatures of the cloud base temperature dimension'
ncdf_attput, cdfid, cld_base_temp_vector_id, 'units','degrees C'

cfad_array_id=ncdf_vardef(cdfid, 'BASTA_CFAD', [dbz_bins_did,height_bins_did, time_bins_did,latitude_bins_did,synoptic_bins_did,cld_base_temperature_bins_did ], /long)
ncdf_attput, cdfid, cfad_array_id, 'long_name','Conditional CFAD of BASTA Observations During Capricorn'
ncdf_attput, cdfid, cfad_array_id, 'units','counts'
ncdf_attput, cdfid, cfad_array_id, 'missing_value', '0'

ncdf_control, cdfid, /endef

ncdf_varput, cdfid, profile_count_id, profile_count
ncdf_varput, cdfid, cloudless_profile_count_id, clear_profile_count
ncdf_varput, cdfid, dbz_bins_vector_id, dbz_bins
ncdf_varput, cdfid, height_bins_vector_id, height_bins
ncdf_varput, cdfid, time_bins_vector_id, time_bins
ncdf_varput, cdfid, latitude_bins_vector_id, latitude_bins
ncdf_varput, cdfid, synoptic_bins_vector_id, synoptic_bins
ncdf_varput, cdfid, cld_base_temp_vector_id, cld_base_temperature_bins

ncdf_varput, cdfid, cfad_array_id, cfad_array

ncdf_close, cdfid

plot_flag=1
if plot_flag eq 1 then begin
  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

  ;;;;;;;;;;;;;;;;;;;;;;;;;;

 
  ;;;;;;;;;;;;;;;;;; plot the track map
  !p.charsize=3.
  !p.font=-1
  !p.background=!d.n_colors-1
  !p.thick=3. 
  

  ;WINDOW, 0, XSIZE=1000, YSIZE=900, TITLE='basta Moments'

  ;plot_color_contour, dbz, min_dbz, max_dbz, height/1000., hrfrac, radar+'Radar Reflectivity '+date, $
  ;'Time (UTC)', 'Height (km)', '!17DbZe',xl, yl, xr, yr
  
  ;cfad_array_id=ncdf_vardef(cdfid, 'BASTA_CFAD', [dbz_bins_did,height_bins_did, time_bins_did,latitude_bins_did,synoptic_bins_did,cld_base_temperature_bins_did ], /long)

plot_cfad_z_dbz=fltarr(n_elements(dbz_bins),n_elements(height_bins))
total_norm=float(total(cfad_array[1:n_elements(dbz_bins)-1,*,*,*,*,*])) 
for k=0,n_elements(height_bins)-1 do begin
for j=1,n_elements(dbz_bins)-1 do begin
plot_cfad_z_dbz[j,k]=float(total(cfad_array[j,k,*,*,*,*]))/total_norm
endfor 
endfor 

  
  ;xl=0.55 & yl=0.85 & xr=0.4 & yr=0.95 & pos_vector=[xl,yl,xr, yr]
  top_top_y=0.93 & dy=0.1 & top_y=0.95 & bot_y=top_y-dy & sep_y=0.05
  plot_num=0
  top_y=top_top_y-(((float(plot_num))*(dy+sep_y)))
  xl=0.1 & yl=top_y-(5.*dy) & xr=0.87 & yr=top_y & pos_vector=[xl,yl,xr, yr]
  color_bar_format_string='(f5.2)' & cb_vertical=1 & cb_charsize=3.5 & cb_pos_vector=[xr+0.035,yl,xr+0.04,yr] & num_contours=20
  
  plot_color_contour_fanning, plot_cfad_z_dbz, 0., max(plot_cfad_z_dbz), height_bins, dbz_bins, ' cfad '+radar_file_prefix,'dBZe', 'Height (km)', 'Frequency',$
    pos_vector,num_contours,color_bar_format_string, cb_vertical,cb_charsize,cb_pos_vector 

;  two_d_array, min_two_d_array, max_two_d_array, y_axis_vector, x_axis_vector, title_string, xtitle_string, ytitle_string, $
;    color_bar_title_string,pos_vector,num_contours,color_bar_format_string, cb_vertical,cb_charsize,cb_pos_vector
    
cld_base_temperature_histo=fltarr(n_elements(cld_base_temperature_bins)) 
for j=0,n_elements(cld_base_temperature_histo)-1 do begin
  cld_base_temperature_histo[j]=float(total(cfad_array[*,*,*,*,*,j]))/float(total(cfad_array))
endfor

time_histo=fltarr(n_elements(time_bins))
for j=0,n_elements(time_histo)-1 do begin
  time_histo[j]=float(total(cfad_array[2:n_elements(dbz_bins)-1,*,j,*,*,*]))/float(total(cfad_array))
endfor

synoptic_histo=fltarr(n_elements(synoptic_bins))
for j=0,n_elements(synoptic_histo)-1 do begin
  synoptic_histo[j]=float(total(cfad_array[2:n_elements(dbz_bins)-1,*,*,*,j,*]))/float(total(cfad_array))
endfor

plot_num=4 & top_y=top_top_y-((float(plot_num))*(dy+sep_y))
xl=0.08 & yl=top_y-dy & xr=0.4 & yr=top_y & pos_vector=[xl,yl,xr, yr]
plot, cld_base_temperature_bins, cld_base_temperature_histo, color=0,xtitle='Cld base temperature',ytitle='frequency',$
  thick=5, charsize=3.0, charthick=3, xthick=5, ythick=5, psym=10, title='Cloud Base Temperature '+radar_file_prefix, $
   pos=pos_vector, xstyle=1 ; yrange=[min([tb23[where(tb23 gt 0.)], tb31[where(tb31 gt 0.)]]), max([tb23[where(tb23 lt 100.)],tb31[where(tb31 lt 100.)]])],
   
 plot_num=4 & top_y=top_top_y-((float(plot_num))*(dy+sep_y))
 xl=0.5 & yl=top_y-dy & xr=0.8 & yr=top_y & pos_vector=[xl,yl,xr, yr]
 plot, time_bins, time_histo, color=0,xtitle='hour bins',ytitle='frequency',$
   thick=5, charsize=3.0, charthick=3, xthick=5, ythick=5, psym=10, title='frequency by hour '+radar_file_prefix, $
   pos=pos_vector, xstyle=1 ; yrange=[min([tb23[where(tb23 gt 0.)], tb31[where(tb31 gt 0.)]]), max([tb23[where(tb23 lt 100.)],tb31[where(tb31 lt 100.)]])], 
   
 plot_num=5 & top_y=top_top_y-((float(plot_num))*(dy+sep_y))
 xl=0.08 & yl=top_y-dy & xr=0.4 & yr=top_y & pos_vector=[xl,yl,xr, yr]
 plot, synoptic_bins, synoptic_histo, color=0,xtitle='synoptic histo',ytitle='frequency',$
   thick=5, charsize=3.0, charthick=3, xthick=5, ythick=5, psym=10, title='-1 pre, 1 post '+radar_file_prefix, $
   pos=pos_vector, xstyle=1     
   
write_gif, cfad_file_path+radar_file_prefix+'_cfad_plot.gif', tvrd()
  
endif ; plot_flag

return
end