;
; This program processes a two dimensional variable.  Uses the averaged
;  temperature, clutterflag, and reflectivity to determine the cloud layers for
;  the averaged time interval
;

Pro two_dimensional_derived,layers,reflectivity,temperature,clutterflag

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

;
;  Open the averaged netcdf file
;

avg_id=ncdf_open(avgfilename)

;
;  Get the variable id
;

reflect_id=ncdf_varid(avg_id,reflectivity)
temp_id=ncdf_varid(avg_id,temperature)
clutter_id=ncdf_varid(avg_id,clutterflag)

if reflect_id ne -1 and temp_id ne -1 and clutter_id ne -1 then begin
print,'Processing 2-D ',layers,' from ',reflectivity,',',clutterflag,' and ',temperature

;
;  Get the variable
;

ncdf_varget, avg_id, reflect_id, reflect
ncdf_varget, avg_id, temp_id, temp
ncdf_varget, avg_id, clutter_id, clutter

reflect=transpose(reflect)
temp=transpose(temp)
clutter=transpose(clutter)

;
;  Close the netcdf file
;

ncdf_close, avg_id

;
;  Scan through the reflectivity profile to find the clouds
;

base_temp=make_array(7,numtimes,/float,value=-9999)
base_height=make_array(7,numtimes,/float,value=-9999)
base_reflect=make_array(7,numtimes,/float,value=-9999)
top_temp=make_array(7,numtimes,/float,value=-9999)
top_height=make_array(7,numtimes,/float,value=-9999)
top_reflect=make_array(7,numtimes,/float,value=-9999)
num_layers=make_array(numtimes,/integer,value=-9999)

for t=0,numtimes-1 do begin
	reflect_hgt=reflect[t,*]  ;one reflectivity profile
	temp_hgt=temp[t,*]
	clutter_hgt=clutter[t,*]
	signal=0  ;flags wether we are in cloud or clear sky
	number=0  ;counts the number of consecutive cloud points
	basen=0  ;base number
	topn=0	;top number
	for h=0,numheights-1 do begin
	;for h=0,25 do begin
		if reflect_hgt[h] ne -9999 and signal eq 0 and number le 0 and clutter_hgt[h] gt 0 then begin
			number=number+1
;			print,reflect_hgt[h],newheights[h],temp_hgt[h],clutter_hgt[h],' Testing for Base'
		endif else if reflect_hgt[h] ne -9999 and signal eq 0 and number eq 1 then begin
			signal=1
;			print,reflect_hgt[h-1],newheights[h-1],temp_hgt[h-1],clutter_hgt[h-1],' Found a base'
;			print,reflect_hgt[h],newheights[h],temp_hgt[h],clutter_hgt[h]
			base_temp[basen,t]=temp_hgt[h-1]
			base_height[basen,t]=newheights[h-1]
			base_reflect[basen,t]=reflect_hgt[h-1]
			basen=basen+1
		endif else if reflect_hgt[h] eq -9999 and signal eq 1 then begin
			signal=0
			number=0
;			print,reflect_hgt[h-1],newheights[h-1],temp_hgt[h-1],clutter_hgt[h-1],' Found a top'
;			print,reflect_hgt[h],newheights[h],temp_hgt[h],clutter_hgt[h]
			i=h-1
			while clutter_hgt[i] lt 1 do i=i-1
			if base_height[basen-1,t] eq newheights[i] then begin
				basen=basen-1
;				print, 'bottom equals top, not a layer'
			endif else begin
				top_temp[topn,t]=temp_hgt[i]
				top_height[topn,t]=newheights[i]
				top_reflect[topn,t]=reflect_hgt[i]
				topn=topn+1
;				print,reflect_hgt[h-1],newheights[h-1],clutter_hgt[h-1],' Found a top'
;				print,reflect_hgt[i],newheights[i],clutter_hgt[i],' ***Real top***'
;				print,reflect_hgt[h],newheights[h],clutter_hgt[h]
			endelse
		endif else begin
;			print,reflect_hgt[h],newheights[h],temp_hgt[h],clutter_hgt[h]
			number=0
		endelse
	endfor
	num_layers[t]=topn
;	print,'layers',num_layers[t]
;stop
endfor

;
;  The big arrays have to be separated into smaller arrays for putting into the 
;  netcdf file.  It won't work any other way
;

base_temp_1=fltarr(numtimes)
base_temp_2=fltarr(numtimes)
base_temp_3=fltarr(numtimes)
base_temp_4=fltarr(numtimes)
base_temp_5=fltarr(numtimes)
base_height_1=fltarr(numtimes)
base_height_2=fltarr(numtimes)
base_height_3=fltarr(numtimes)
base_height_4=fltarr(numtimes)
base_height_5=fltarr(numtimes)
top_temp_1=fltarr(numtimes)
top_temp_2=fltarr(numtimes)
top_temp_3=fltarr(numtimes)
top_temp_4=fltarr(numtimes)
top_temp_5=fltarr(numtimes)
top_height_1=fltarr(numtimes)
top_height_2=fltarr(numtimes)
top_height_3=fltarr(numtimes)
top_height_4=fltarr(numtimes)
top_height_5=fltarr(numtimes)

for t=0,numtimes-1 do begin
	base_temp_1[t]=base_temp[0,t]
	base_temp_2[t]=base_temp[1,t]
	base_temp_3[t]=base_temp[2,t]
	base_temp_4[t]=base_temp[3,t]
	base_temp_5[t]=base_temp[4,t]

	base_height_1[t]=base_height[0,t]
	base_height_2[t]=base_height[1,t]
	base_height_3[t]=base_height[2,t]
	base_height_4[t]=base_height[3,t]
	base_height_5[t]=base_height[4,t]

	top_temp_1[t]=top_temp[0,t]
	top_temp_2[t]=top_temp[1,t]
	top_temp_3[t]=top_temp[2,t]
	top_temp_4[t]=top_temp[3,t]
	top_temp_5[t]=top_temp[4,t]

	top_height_1[t]=top_height[0,t]
	top_height_2[t]=top_height[1,t]
	top_height_3[t]=top_height[2,t]
	top_height_4[t]=top_height[3,t]
	top_height_5[t]=top_height[4,t]
endfor

;
;  Open the netcdf file
;

avg_id=ncdf_open(avgfilename, /write)

;
;  Get the variable dimensions
;

time_id=ncdf_dimid(avg_id,'time')

;
;  Put the netcdf file into define mode
;

ncdf_control, avg_id, /redef

;
;  Define variable
;

num_layers_id=ncdf_vardef(avg_id,'number_of_layers',time_id,/short)

base_height_layer_1_id=ncdf_vardef(avg_id,'base_height_layer_1', time_id)
base_height_layer_2_id=ncdf_vardef(avg_id,'base_height_layer_2', time_id)
base_height_layer_3_id=ncdf_vardef(avg_id,'base_height_layer_3', time_id)
base_height_layer_4_id=ncdf_vardef(avg_id,'base_height_layer_4', time_id)
base_height_layer_5_id=ncdf_vardef(avg_id,'base_height_layer_5', time_id)

base_temp_layer_1_id=ncdf_vardef(avg_id,'base_temp_layer_1', time_id)
base_temp_layer_2_id=ncdf_vardef(avg_id,'base_temp_layer_2', time_id)
base_temp_layer_3_id=ncdf_vardef(avg_id,'base_temp_layer_3', time_id)
base_temp_layer_4_id=ncdf_vardef(avg_id,'base_temp_layer_4', time_id)
base_temp_layer_5_id=ncdf_vardef(avg_id,'base_temp_layer_5', time_id)

top_height_layer_1_id=ncdf_vardef(avg_id,'top_height_layer_1', time_id)
top_height_layer_2_id=ncdf_vardef(avg_id,'top_height_layer_2', time_id)
top_height_layer_3_id=ncdf_vardef(avg_id,'top_height_layer_3', time_id)
top_height_layer_4_id=ncdf_vardef(avg_id,'top_height_layer_4', time_id)
top_height_layer_5_id=ncdf_vardef(avg_id,'top_height_layer_5', time_id)

top_temp_layer_1_id=ncdf_vardef(avg_id,'top_temp_layer_1', time_id)
top_temp_layer_2_id=ncdf_vardef(avg_id,'top_temp_layer_2', time_id)
top_temp_layer_3_id=ncdf_vardef(avg_id,'top_temp_layer_3', time_id)
top_temp_layer_4_id=ncdf_vardef(avg_id,'top_temp_layer_4', time_id)
top_temp_layer_5_id=ncdf_vardef(avg_id,'top_temp_layer_5', time_id)

;
;  Define attributes
;

ncdf_attput, avg_id, 'number_of_layers', 'long_name', 'the number of cloud layers'

for i=1,5 do begin
	num=string(i)
	dataname='base_height_layer_'+num
	dataname=strcompress(dataname,/remove_all)
	longname='base height of radar layer'+num
	longname=strcompress(longname)
	theunits='m - agl'
	ncdf_attput, avg_id, dataname, 'long_name', longname
	ncdf_attput, avg_id, dataname, 'units', theunits
endfor

for i=1,5 do begin
	num=string(i)
	dataname='base_temp_layer_'+num
	dataname=strcompress(dataname,/remove_all)
	longname='base temp of radar layer'+num
	longname=strcompress(longname)
	theunits='kelvins'
	ncdf_attput, avg_id, dataname, 'long_name', longname
	ncdf_attput, avg_id, dataname, 'units', theunits
endfor

for i=1,5 do begin
	num=string(i)
	dataname='top_height_layer_'+num
	dataname=strcompress(dataname,/remove_all)
	longname='top height of radar layer'+num
	longname=strcompress(longname)
	theunits='m - agl'
	ncdf_attput, avg_id, dataname, 'long_name', longname
	ncdf_attput, avg_id, dataname, 'units', theunits
endfor

for i=1,5 do begin
	num=string(i)
	dataname='top_temp_layer_'+num
	dataname=strcompress(dataname,/remove_all)
	longname='top temp of radar layer'+num
	longname=strcompress(longname)
	theunits='kelvins'
	ncdf_attput, avg_id, dataname, 'long_name', longname
	ncdf_attput, avg_id, dataname, 'units', theunits
endfor

;
;  Change from define mode to data mode
;

ncdf_control, avg_id, /endef

;
;  Write variable
;

ncdf_varput, avg_id, num_layers_id, num_layers

ncdf_varput, avg_id, base_height_layer_1_id, base_height_1
ncdf_varput, avg_id, base_height_layer_2_id, base_height_2
ncdf_varput, avg_id, base_height_layer_3_id, base_height_3
ncdf_varput, avg_id, base_height_layer_4_id, base_height_4
ncdf_varput, avg_id, base_height_layer_5_id, base_height_5

ncdf_varput, avg_id, base_temp_layer_1_id, base_temp_1
ncdf_varput, avg_id, base_temp_layer_2_id, base_temp_2
ncdf_varput, avg_id, base_temp_layer_3_id, base_temp_3
ncdf_varput, avg_id, base_temp_layer_4_id, base_temp_4
ncdf_varput, avg_id, base_temp_layer_5_id, base_temp_5

ncdf_varput, avg_id, top_height_layer_1_id, top_height_1
ncdf_varput, avg_id, top_height_layer_2_id, top_height_2
ncdf_varput, avg_id, top_height_layer_3_id, top_height_3
ncdf_varput, avg_id, top_height_layer_4_id, top_height_4
ncdf_varput, avg_id, top_height_layer_5_id, top_height_5

ncdf_varput, avg_id, top_temp_layer_1_id, top_temp_1
ncdf_varput, avg_id, top_temp_layer_2_id, top_temp_2
ncdf_varput, avg_id, top_temp_layer_3_id, top_temp_3
ncdf_varput, avg_id, top_temp_layer_4_id, top_temp_4
ncdf_varput, avg_id, top_temp_layer_5_id, top_temp_5

;
;  Close the netcdf file
;

ncdf_close, avg_id

endif  ;variable does exist in the file

return
end