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

Pro two_dimensional_mode, 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

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

common bracket, startbst

common data, avg_data

;
;  Open the netcdf file
;

cdf_id=ncdf_open(oldcdfname)

;
;  Get the variable id
;

data_id=ncdf_varid(cdf_id,dataname)
if data_id ne -1 then begin
print,'Processing 2-D ',dataname

;
;  Get the variable
;

ncdf_varget, cdf_id, data_id, data
data=transpose(data)
if dataname eq 'qc_ReflectivityClutterFlag' then data=fix(data)

;
;  Inquire about the variable
;

datainfo=ncdf_varinq(cdf_id, data_id)

;
;  Number of attributes assigned to the variable
;

numatts=datainfo.natts

;
;  Store the attribute names
;

dataatr=strarr(numatts)  ;string array to hold attribute names
for i=0,numatts-1 do begin
  dataatr[i]=ncdf_attname(cdf_id, data_id, i)  ;Returns the name of the attribute
endfor

;
;  Remove bad values
;

result=where (data ge 1e18,count)
if count ne 0 then data[result]=-9999

;
;  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,'count replace in blanks',count
if count ne 0 then data [result]=-9999

;
;  Convert mxrat to relhum
;

if dataname eq 'mxrat' then convert_to_relhum,oldcdfname,time_dim_name,height_dim_name,$
  data,dataname

;
;  Convert ReflectivityBestEstimate to the format for averaging
;

if dataname eq 'ReflectivityBestEstimate' then begin
	print,'inside Reflectivity Best Estimate fix routine'
	data=float(data)
	result=where(data ne -9999, count)
	if count ne 0 then begin 
		data[result]=data[result]/100.  ;convert reflectivity to DbZe
		data[result]=data[result]/10
		data[result]=10^data[result]  ;convert to Ze
	endif
endif

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

possible_points=fltarr(numtimes,numheights)  ;points with clutter flag of 0,3,7
good_points=fltarr(numtimes,numheights)  ;points with clutter flag of 1,2
mode_data=fltarr(numtimes,numheights)  ;mode
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
	result_all=where(time_offset ge newtimes[t] and time_offset lt newtimes[t+1], count_all)
	possible_points[t,h]=count_all  ;the number of points in a time interval
	result=where(time_offset ge newtimes[t] and time_offset lt newtimes[t+1] and data_hgt ne -9999, count)
	good_points[t,h]=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 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
;			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,possible_points[t,h],good_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

;
;  Calculate the average 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
	result_all=where(time_offset ge newtimes[numtimes-1] and time_offset lt (newtimes[numtimes-1]+timespacing), count_all)
	possible_points[numtimes-1,h]=count_all
	result=where(time_offset ge newtimes[numtimes-1] and time_offset lt (newtimes[numtimes-1]+timespacing) and data_hgt  ne -9999, count)
	good_points[numtimes-1,h]=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

;
;  Return the ReflectivityBestEstimate data to the writing form
;

if dataname eq 'ReflectivityBestEstimate' then begin
	result=where(avg_data ne -9999,count)
	if count ne 0 then begin
		avg_data[result]=10*alog10(avg_data[result])
		avg_data[result]=avg_data[result]*100
	endif
	avg_data=fix(avg_data)
endif

;
;  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
possible_id=ncdf_vardef(avg_id,dataname+'PossiblePoints',datadim)  ;total points in time interval
useable_id=ncdf_vardef(avg_id,dataname+'UseablePoints',datadim)  ;good data points
;
;  Define attributes
;

;if dataname eq 'relhum' then begin
;  ncdf_attput, avg_id, mode_data_id, 'long_name', 'relative humidity'
;  ncdf_attput, avg_id, mode_data_id, 'units', 'percent'
;endif else begin
;  for i=0,numatts-1 do begin
;    nowatt=dataatr[i]
;    mytry=ncdf_attcopy(cdf_id, data_id, nowatt, avg_id, mode_data_id)
;    mytry=ncdf_attcopy(cdf_id, data_id, nowatt, avg_id, hi_ext_id)
;    mytry=ncdf_attcopy(cdf_id, data_id, nowatt, avg_id, lo_ext_id)
;   endfor
;endelse
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-no cloud'
ncdf_attput, avg_id, mode_data_id, 'comment2','1-cloud'
ncdf_attput, avg_id, mode_data_id, 'original data source file ',oldcdfname

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-no cloud'
ncdf_attput, avg_id, hi_ext_id, 'comment2','1-cloud'
ncdf_attput, avg_id, hi_ext_id, 'original data source file ',oldcdfname

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-no cloud'
ncdf_attput, avg_id, lo_ext_id, 'comment2','1-cloud'
ncdf_attput, avg_id, lo_ext_id, 'original data source file ',oldcdfname

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

ncdf_attput, avg_id, useable_id, 'long_name','MMCR Reflectivity Clutter Flags - Useable data points within the time interval'
ncdf_attput, avg_id, useable_id, 'units','number of points'
ncdf_attput, avg_id, useable_id, 'original data source file ',oldcdfname

;
;  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
possible_points=transpose(possible_points)  ;transpose for writing to netcdf
ncdf_varput, avg_id, possible_id, possible_points
possible_points=transpose(possible_points)  ;retrun to original form for plotting
good_points=transpose(good_points)  ;transpose for writing to netcdf
ncdf_varput, avg_id, possible_id, good_points
good_points=transpose(good_points)  ;retrun to original form for plotting

;
;  Close the netcdf file
;

ncdf_close, avg_id

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

if dataname eq 'ReflectivityBestEstimate' then begin
	miv=-5000
	mav=max(mode_data)
endif else begin
	mav=max(mode_data)
	miv=0
endelse
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'  ;flat 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=bytarr(num_times,num_heights)  ;initialized to 0
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=bytarr(numtimes,numheights)  ;initialized to 0
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 the file

ncdf_close, cdf_id

return
end