;
; This program processes a two dimensional variable
;

Pro profile,dataname

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

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)

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

;
;  Convert mxrat to relhum
;

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

;
;  Average data in timespacing intervals 
;

avg_data=fltarr(numtimes,numheights)

for h=0,numheights-1 do begin  ;loop through heights
;h=4
for t=0,numtimes-2 do begin  ;loop through time intervals
;for t=0,4 do begin
	data_hgt=data[*,h]  ;data values at a constant height
	result=where(time_offset ge newtimes[t] and time_offset lt newtimes[t+1] and data_hgt ne -9999, count)
	if count ge 3 then avg_data[t,h]=total(data_hgt[result])/count $
	else avg_data[t,h]=-9999

;	if count ge 3 then begin
;		print,'data ',t,count,data_hgt[result]
;		print,'time ',t,count,time_offset[result]
;		print,avg_data[t,h]
;	endif else print, t, count, '******************not enough values'
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=where(time_offset ge newtimes[numtimes-1] and time_offset lt (newtimes[numtimes-1]+300) and data_hgt  ne -9999, count)
	if count ge 3 then avg_data[numtimes-1,h]=total(data_hgt[result])/count $
	else avg_data[numtimes-1,h]=-9999
endfor

;  Get reflectivity best estimate from the alread averaged file

;
;  Open the averaged netcdf file
;

cdf_id=ncdf_open(avgfilename)

;
;  Get the variable id
;

reflectivity_id=ncdf_varid(cdf_id,'ReflectivityBestEstimate')

if reflectivity_id ne -1 then begin
print,'Processing 2-D Reflectivity Best Estimate'

;
;  Get the variable
;

ncdf_varget, cdf_id, reflectivity_id, reflectivity
reflectivity=transpose(reflectivity)

for t=0, numtimes-1 do begin
	for h=0, numheights-1 do begin
		if reflectivity[t,h] ne -9999 then print,reflectivity[t,h], newheights[h],newtimes[t],'yes' else print,reflectivity[t,h],newheights[h]
	endfor
	stop
endfor

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

avg_data_id=ncdf_vardef(avg_id, dataname, datadim)

;
;  Define attributes
;

if dataname eq 'relhum' then begin
  ncdf_attput, avg_id, avg_data_id, 'long_name', 'relative humidity'
  ncdf_attput, avg_id, avg_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, avg_data_id)
  endfor
endelse
ncdf_attput, avg_id, avg_data_id, 'original data source file ',oldcdfname

;
;  Change from define mode to data mode
;

ncdf_control, avg_id, /endef

;
;  Write variable
;

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

;
;  Close the netcdf file
;

ncdf_close, avg_id

;
;  Plot to check data
;

;
;  Scatter plot to check the data
;

;if dataname eq 'ReflectivityBestEstimate' then begin
	;miv=-5000
	;mav=max(avg_data)
;endif else begin
;	mav=max(avg_data)
;	miv=0
;endelse
;window,2,title='avg '+dataname
;surface,avg_data,newtimes,newheights,min_value=miv,max_value=mav,/horizontal

;window,3,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:
;data=data[n1:n2,*]
;sub_time_offset=time_offset[n1:n2]
;zoomsize=size(data)
;if zoomsize[1] ge 2 then surface,data,sub_time_offset,height,min_value=miv,max_value=mav,/horizontal

;
;  First plot old data
;

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

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

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

;image=bytarr(num_times,num_heights)  ;initialized to 0
;result=where(data ne 1e6, count)
;if count ne 0 then image[result]=((255)/(dmax-dmin))*(data[result]-dmin)
;print,max(image),min(image)

;
;  Switch the max and min values of image
;

;image(where(image eq max(image)))=0
;image(where(image eq min(image)))=255

;
;  Plot the old data
;

;window,3,xsize=600,ysize=600
;loadct,5
;contour, newdata, 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
;n1=0
;while time_offset[n1] lt newtimes[0] do  n1=n1+1
;n2=n1
;while time_offset[n2] lt newtimes[numtimes-1] do n2=n2+1
;image=image[n1:n2,*]
;tv,congrid(image,sx,sy),px[0],py[0]

;
;  Next plot the regridded data
;

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

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

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

;newimage=bytarr(numtimes,newnumheights)  ;initialized to 0
;result=where(newdata ne 1e6, count)
;if count ne 0 then newimage[result]=((255)/(dmax-dmin))*(newdata[result]-dmin)
;print,max(newimage),min(newimage)

;
;  Switch the max and min values of image
;

;newimage(where(newimage eq max(newimage)))=0
;newimage(where(newimage eq min(newimage)))=255

;
;  Plot the new regridded data
;

;window,4,xsize=600,ysize=600
;loadct,5
;contour, newdata, 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  ;ReflectivityBestEstimate doesn't exist
ncdf_close, avg_id
endif  ;variable does exist in the file

ncdf_close, cdf_id

return
end