;
;  This program processes the two dimensional clutter flag and 
;  writes it to the output netcdf file.  Determines the number
;  of total data points, the number of cloudy data points,
;  mode, minimum, and maximum clutter flag value for each time
;  interval
;

pro two_dimensional_mode4, dataname

common ncdf_file, oldcdfname,twoD,time_dim_name,height_dim_name,$
  height_var_name,avgfilename,curdate

common old_grid, num_times,num_heights,base_time,time_offset,$
  height,site,files,num_files,num_days

common new_grid, numtimes,newtimes, numheights,newheights,$
 timespacing,heightspacing

common bracket, startbst

common data, avg_data

;
;  Loop through all the matching files in the directory
;
m=-1  ;counts number of cdf files without the variable
p=-1  ;counts number of cdf files with variable
j=-1  ;counts number of cdf files
for i = 0, num_files-1 do begin
;
;  Don't try to process any gif files that may be in the directory
;
  if strpos(files[i],'gif') eq -1 and $
  strpos(files[i],'png') eq -1 then begin
    j=j+1
    print,i,j,files[i]
;
;  Open the netcdf file
;
    cdf_id=ncdf_open(files[i])
;
;  Get the dimension id's
;
    time_did=ncdf_dimid(cdf_id,time_dim_name)
    height_did=ncdf_dimid(cdf_id,height_dim_name)
;
;  Get the dimensions
;
    ncdf_diminq,cdf_id,time_did,charstring,ntimes
    ncdf_diminq,cdf_id,height_did,charstring,nheights
;
;  Create array for nodata value
;
    dummy=make_array(ntimes,nheights,/float,value=-8888)
;
;  Get the variable id
;
    data_id=ncdf_varid(cdf_id,dataname)
;
;  Get the variable if it exists
;
    if data_id ne -1 then begin
      print,'Processing 2-D ',dataname
      ncdf_varget, cdf_id, data_id, data
      data=transpose(data)
      p=p+1
    endif else begin
      print, 'Variable not in this file'
      m=m+1
    endelse
;
;  If it is the first real variable get the attributes
;
    if p eq 0 and data_id ne -1 then begin
      print, 'first variable, getting the attributes'
      datainfo=ncdf_varinq(cdf_id, data_id)  ;inquire about the variable
      numatts=datainfo.natts  ;number of attributes assigned to the variable
      parts=strsplit(files[i],'/',/extract)
      part=strsplit(parts[5],'.',/extract)
      origfile=part[0]+'.'+part[1]
;
;  Store the attribute names
;
      dataatr=strarr(numatts)  ;string array to hold attribute names
      for k=0,numatts-1 do begin
;  Returns the name of the attribute
        dataatr[k]=ncdf_attname(cdf_id, data_id, k)  
      endfor
    endif 
;
;  Create the data array
;
    if p eq 0 and m lt 0 and data_id ne -1 then begin
      print,'first file and variable exists' 
      combined_data=data
    endif else if p ge 0 and data_id ne -1 then begin
      print,'adding on existing variable'
      combined_data=[combined_data,data]
    endif else if p lt 0 and m eq 0 and data_id eq -1 then begin
      print,'first file and variable doesnt exist'
      combined_data=dummy
    endif else if m ge 0 and data_id eq -1 then begin
      print,'adding on dummy variable'
      combined_data=[combined_data,dummy]
    endif
;
;  Close the netcdf file
;
    ncdf_close,cdf_id
  endif  ;end of not a gif file
endfor  ;end of looping though the matching files

;
;  Continue processing data because some data exists
;

if j gt -1 and p gt -1 then begin 

;
;  Manipulate the big data array
;

data=combined_data
combined_data=0  ;release memory
data=fix(data)

;
;  This fixes a bug in the Merged Moments file.  An undefined value at each
;  data(*, last height) sets data=-32767.0
;  Reset this value to the nodata value of -9999 for the Merged Moment file
;

result=where(data eq -32767.0, count)
print,'replacing -32767 with -9999',count
if count ne 0 then data [result]=-9999

;
;  Find the mode and the extremes of data in timespacing intervals
;

all_points=fltarr(numtimes,numheights)  ;the number of points in a time interval
cloudy_points=fltarr(numtimes,numheights)  ;cloudy points (has cloudy clutter flag)
mode_data=fltarr(numtimes,numheights)  ;mode value for clutter in the time interval
hi_ext=fltarr(numtimes,numheights)  ;highest data value in the time interval
lo_ext=fltarr(numtimes,numheights)  ;lowest data value in the time interval

for h=0,numheights-1 do begin  ;loop through heights
  for t=0,numtimes-2 do begin  ;loop through time intervals
    data_hgt=data[*,h]  ;data values at a constant height

;   counts existing data, excludes missing data
    result_all=where(time_offset ge newtimes[t] and $
                     time_offset lt newtimes[t+1] and $
                     data_hgt ne -8888, count_all)
    all_points[t,h]=count_all  ;the number of points in a time interval

;   counts truly cloudy data, clutter 1 or 2
    result_cld=where(time_offset ge newtimes[t] and $
                     time_offset lt newtimes[t+1] and $
                     data_hgt ne -9999 and $
                     data_hgt ne -8888 and $
                     (data_hgt eq 1 or data_hgt eq 2), count_cld)
    cloudy_points[t,h]=count_cld  ;cloudy points

;   counts data with a clutter flag value for determining mode,min,max
    result=where(time_offset ge newtimes[t] and $
                 time_offset lt newtimes[t+1] and $
                 data_hgt ne -9999 and $
                 data_hgt ne -8888, count)
    if count ne 0 then begin  ;if there are cloudy points
      sub_data=data_hgt[result]  ;data within the time interval
      ;find the good clutter flags and give them all one number
      gresult=where(sub_data eq 1 or sub_data eq 2, gcount)
      if gcount ne 0 then sub_data[gresult]=1
      ;find the bad clutter flags and give them all one number
      bresult=where(sub_data eq 0 or sub_data eq 3 or sub_data eq 7, bcount)
      if bcount ne 0 then sub_data[bresult]=0
      ;find which has more points, bad or good data
      hist=histogram(sub_data,omax=hi_extreme,omin=lo_extreme)
      bins=findgen(n_elements(hist))+min(sub_data)
      hist_max=max(hist,ind)
      if hist_max eq 1 and count gt 1 then mode_data[t,h]=-9999 else $
      mode_data[t,h]=bins(ind)  ;mode value
        hi_ext[t,h]=hi_extreme
	lo_ext[t,h]=lo_extreme
	;if t ge 10 and t le 50 and h eq 10 then begin
        ;if cloudy_points[t,h] ne 0 then begin
      	  ;print,'*************',count,'************'
	  ;print,'data',data_hgt[result_all]
	  ;print,'sub_data',sub_data
	  ;print,'hist_max',hist_max,bins(ind),'   hi',hi_extreme,$
          ;      '   lo',lo_extreme
	  ;print,'all',all_points[t,h],'cloudy',cloudy_points[t,h]
	  ;window,0,xsize=1200,ysize=400
	  ;mode=make_array(count,/float,value=mode_data[t,h])
	  ;plot_times=time_offset[result]
	  ;plot,time_offset,data_hgt,yrange=[-1,8],ystyle=1, $
          ;xrange=[plot_times[0]-500,plot_times[count-1]+500],xstyle=1
	  ;oplot,time_offset[result],mode,psym=-6
       ;endif
    endif else begin 
      mode_data[t,h]=-9999
      hi_ext[t,h]=-9999
      lo_ext[t,h]=-9999
    endelse
  endfor
endfor

;
;  Find the mode and extremes for the last time interval
;

for h=0, numheights-1 do begin  ;loop through heights
  data_hgt=data[*,h]  ;data values at a constant height

; counts existing data, excludes missing data
  result_all=where(time_offset ge newtimes[numtimes-1] and $
                   time_offset lt (newtimes[numtimes-1]+timespacing) and $
                   data_hgt ne -8888, count_all)
  all_points[numtimes-1,h]=count_all

; counts truly cloudy data, clutter 1 or 2
  result_cld=where(time_offset ge newtimes[numtimes-1] and $
                   time_offset lt (newtimes[numtimes-1]+timespacing) and $
                   data_hgt  ne -9999 and $
                   data_hgt ne -8888 and $
                   (data_hgt eq 1 or data_hgt eq 2), count_cld)
  cloudy_points[numtimes-1,h]=count_cld  ;cloudy points

; counts data with a clutter flag value for determining mode,min,max
  result=where(time_offset ge newtimes[numtimes-1] and $
               time_offset lt (newtimes[numtimes-1]+timespacing) and $
               data_hgt  ne -9999 and $
               data_hgt ne -8888, count)
  if count ne 0 then begin 
    sub_data=data_hgt[result]  ;data within the time interval
    ;find the good clutter flags and give them all one number
    gresult=where(sub_data eq 1 or sub_data eq 2, gcount)
    if gcount ne 0 then sub_data[gresult]=1
    ;find the bad clutter flags and give them all one number
    bresult=where(sub_data eq 0 or sub_data eq 3 or sub_data eq 7, bcount)
    if bcount ne 0 then sub_data[bresult]=0
    ;find which has more points, bad or good data
    hist=histogram(sub_data,omax=hi_extreme,omin=lo_extreme)
    bins=findgen(n_elements(hist))+min(sub_data)
    hist_max=max(hist,ind)
    if hist_max eq 1 then mode_data[numtimes-1,h]=-9999 else $
      mode_data[numtimes-1,h]=bins(ind)  ;mode value
      hi_ext[numtimes-1,h]=hi_extreme
      lo_ext[numtimes-1,h]=lo_extreme
  endif else begin 
    mode_data[numtimes-1,h]=-9999
    hi_ext[t,h]=-9999
    lo_ext[t,h]=-9999
  endelse
endfor


;**************************************
;  Write to the output file
;**************************************

;
;  Open the netcdf file
;

avg_id=ncdf_open(avgfilename, /write)

;
;  Get the variable dimensions
;

time_id=ncdf_dimid(avg_id,'time')
height2_did=ncdf_dimid(avg_id,'height')
datadim=intarr(2)
datadim[1]=time_id
datadim[0]=height2_did

;
;  Put the netcdf file into define mode
;

ncdf_control, avg_id, /redef

;
;  Define variable
;

mode_data_id=ncdf_vardef(avg_id, dataname+'Mode', datadim)  ;mode value
hi_ext_id=ncdf_vardef(avg_id,dataname+'MaxExtreme',datadim)  ;highest value
lo_ext_id=ncdf_vardef(avg_id,dataname+'MinExtreme',datadim)  ;lowest value
all_id=ncdf_vardef(avg_id,dataname+'AllPoints',datadim)  ;total points in time interval
cloudy_id=ncdf_vardef(avg_id,dataname+'CloudyPoints',datadim)  ;truly cloudy data points

;
;  Define attributes
;

ncdf_attput,avg_id,mode_data_id,'long_name','MMCR Reflectivity Clutter Flags - Mode for the time interval'
ncdf_attput,avg_id,mode_data_id,'units','unitless'
ncdf_attput,avg_id,mode_data_id,'comment1','0-not cloud'
ncdf_attput,avg_id,mode_data_id,'comment2','1-cloud'
ncdf_attput,avg_id,mode_data_id,'original data source file ',origfile

ncdf_attput,avg_id,hi_ext_id,'long_name','MMCR Reflectivity Clutter Flags - Max extreme for the time interval'
ncdf_attput,avg_id,hi_ext_id, 'units','unitless'
ncdf_attput,avg_id,hi_ext_id, 'comment1','0-not cloud'
ncdf_attput,avg_id,hi_ext_id, 'comment2','1-cloud'
ncdf_attput,avg_id,hi_ext_id, 'original data source file ',origfile

ncdf_attput,avg_id,lo_ext_id,'long_name','MMCR Reflectivity Clutter Flags - Min extreme for the time interval'
ncdf_attput,avg_id,lo_ext_id,'units','unitless'
ncdf_attput,avg_id,lo_ext_id,'comment1','0-not cloud'
ncdf_attput,avg_id,lo_ext_id,'comment2','1-cloud'
ncdf_attput,avg_id,lo_ext_id,'original data source file ',origfile

ncdf_attput,avg_id,all_id,'long_name','MMCR Reflectivity Clutter Flags - Total data points within the time interval'
ncdf_attput,avg_id,all_id,'units','number of points'
ncdf_attput,avg_id,all_id,'original data source file ',origfile

ncdf_attput,avg_id,cloudy_id,'long_name','MMCR Reflectivity Clutter Flags - Cloudy data points within the time interval'
ncdf_attput,avg_id,cloudy_id,'units','number of points'
ncdf_attput,avg_id,cloudy_id,'original data source file ',origfile

;
;  Change from define mode to data mode
;

ncdf_control, avg_id, /endef

;
;  Write variable
;

mode_data=transpose(mode_data)  ;transpose for writing to netcdf
ncdf_varput, avg_id, mode_data_id, mode_data
mode_data=transpose(mode_data)  ;retrun to original form for plotting

hi_ext=transpose(hi_ext)  ;transpose for writing to netcdf
ncdf_varput, avg_id, hi_ext_id, hi_ext
hi_ext=transpose(hi_ext)  ;retrun to original form for plotting

lo_ext=transpose(lo_ext)  ;transpose for writing to netcdf
ncdf_varput, avg_id, lo_ext_id, lo_ext
lo_ext=transpose(lo_ext)  ;retrun to original form for plotting

all_points=transpose(all_points)  ;transpose for writing to netcdf
ncdf_varput, avg_id, all_id, all_points
all_points=transpose(all_points)  ;retrun to original form for plotting

cloudy_points=transpose(cloudy_points)  ;transpose for writing to netcdf
ncdf_varput, avg_id, cloudy_id, cloudy_points
cloudy_points=transpose(cloudy_points)  ;retrun to original form for plotting

;
;  Close the netcdf file
;

ncdf_close, avg_id

;**********************************
;  Plots to check data
;**********************************

;
;  Surface plots to check the data to make sure the
;  mode data represents the raw data
;

surface_plots='no'  ;Flag this yes if you want surface plots, no if you don't
if surface_plots eq 'yes' then begin

;
;  Surface plot of the mode data
;
mav=max(mode_data)
miv=0
window,0,title='mode '+dataname
surface,mode_data,newtimes,newheights,min_value=miv,max_value=mav,/horizontal

;
;  Surface plot of the raw data
;

window,1,title='raw '+dataname
n1=0
while time_offset[n1] lt newtimes[0] do  n1=n1+1
n2=n1
while time_offset[n2] lt newtimes[numtimes-1] do begin
  n2=n2+1
  if n2 eq num_times then begin
    n2=n2-1
    goto, enough
  endif
endwhile
enough:
old_data=data[n1:n2,*]
old_time_offset=time_offset[n1:n2]
zoomsize=size(data)
if zoomsize[1] ge 2 then surface,old_data,old_time_offset,height,min_value=miv,max_value=mav,/horizontal

endif  ;if you want surface plots

;
;  Image plots as another check of the data
;

image_plots='no'  ;flag this yes if you want image plots, no if you don't
if image_plots eq 'yes' then begin

;
;  Find the data range of the raw data and set the null values to 1e6
;

dmax=max(old_data)  ;max data value
dmin=min(old_data)  ;null value
;result=where(old_data eq dmin,count)
;if count ne 0 then old_data[result]=1e6
;dmin=min(old_data)  ;min data value
;print, 'data min', dmin, 'data max', dmax

;
;  Scale the raw data values and insert them into image
;

image=bytscl(old_data,max=dmax,min=dmin, top=255)

;
;  Plot the old data image
;

window,2,xsize=600,ysize=600
loadct,5
contour, old_data, old_time_offset, height, /nodata,$
	xstyle=1,ystyle=1
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]

;
;  Next plot the regridded data
;

;
;  Find the mode data range and set the null values to 1e6
;  This is commented out so that the raw data dmin and dmax 
;   is used for scaling
;

;dmin=min(mod_data)  ;null value
;result=where(mode_data eq dmin,count)
;if count ne 0 then mode_data[result]=1e6
;dmin=min(mode_data)  ;min data value
;print, 'mode data min', dmin, 'newdata max', dmax

;
;  Scale the data vales and insert then into image
;

newimage=bytscl(mode_data,max=dmax,min=dmin,top=255)

;
;  Plot the new regridded data
;

window,3,xsize=600,ysize=600
loadct,5
contour, mode_data, newtimes, newheights, /nodata,$
 	xstyle=1,ystyle=1
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(newimage,sx,sy),px[0],py[0]

endif  ;if you want image plots

endif  ;variable does exist in one of the files

return
end