;
; This program processes a two dimensional variable and writes it to the output
;  netcdf file.  Time averages temperature and then height interpolates it to the 
;  new averaged height grid.
;

Pro two_dimensional_interpolate,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)

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

;
;  Average data in timespacing intervals, keeping height intervals the same
;

avg_data=fltarr(numtimes,num_heights)

for h=0,num_heights-1 do begin  ;loop through old heights
for t=0,numtimes-2 do begin  ;loop through new time intervals
	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
endfor  ;loop through new times
endfor  ;loop through old heights

;
;  Calculate the average for the last time interval
;

for h=0, num_heights-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]+timespacing) 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
;print,newtimes[numtimes-1],newtimes[numtimes-1]+timespacing

;
;  Now interpolate the time averaged temperature data to the new height grid
;

inter_data=fltarr(numtimes,numheights)
for t=0,numtimes-1 do begin  ;loop through new time intervals
	avg_data_hgt=avg_data[t,*]  ;one height column
	result=where(avg_data_hgt ne -9999, count)
	if count ne 0 then begin
		avg_data_hgt=avg_data_hgt[result]
		avg_hgt=height[result]
		inter_data[t,*]=interpol(avg_data_hgt,avg_hgt,newheights)
	endif else begin
		inter_data[t,*]=-9999
	endelse
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
;

inter_data=transpose(inter_data)  ;transpose for writing to netcdf

ncdf_varput, avg_id, avg_data_id, inter_data
inter_data=transpose(inter_data)  ;retrun to original form for plotting

;
;  Close the netcdf file
;

ncdf_close, avg_id

;
;  Surface plots to check the data to make sure the 
;  interpolated 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 interpolated data
;

mav=max(inter_data)  ;max value of the interpolated data
miv=0				;min value of the interpolated data
window,0,title='interpolated  '+dataname
surface,inter_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:
sub_data=data[n1:n2,*]
sub_time_offset=time_offset[n1:n2]
zoomsize=size(data)
if zoomsize[1] ge 2 then surface,sub_data,sub_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(sub_data)  ;max data value
dmin=min(sub_data)  ;null value
result=where(sub_data eq dmin,count)
if count ne 0 then sub_data[result]=1e6
dmin=min(sub_data)  ;min data value
print, 'sub_data min', dmin, 'data max', dmax

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

image=bytarr(num_times,num_heights)  ;initialized to 0
image=bytscl(sub_data,max=dmax,min=dmin,top=255)

;
;  Plot the raw data image
;

window,2,xsize=600,ysize=600
loadct,5
contour, sub_data, sub_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 and interpolated data
;

;
;  Find the interpolated 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(inter_data)  ;null value
;result=where(inter_data eq dmin,count)
;if count ne 0 then inter_data[result]=1e6
;dmin=min(inter_data)  ;min data value
;print, 'inter_data min', dmin, 'newdata max', dmax

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

newimage=bytarr(numtimes,numheights)  ;initialized to 0
newimage=bytscl(inter_data,max=dmax,min=dmin,top=255)

;
;  Plot the new data image
;

window,3,xsize=600,ysize=600
loadct,5
contour, inter_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