;Title:		readspectra.pro
;Purpose:	Read ARM MMCR spectra files, including the moments data
;		and necessary support information about the radar
;		operations.
;
;Inputs:	arg:  Date/time argument for the file(s) of interest:  "YYYYMMDD.HH"
;		mode: String name of the desired mode (i.e., "Stratus", "Cirrus", "Robust", "General", etc.)
;  Keywords:	site: The site acronym for file prefixes (i.e. "nsa", "sgp"...)
;Outputs:
;I/O:
;Author:	Matthew Shupe
;Date:		2/9/05
;-------------------------------------------------------------------------------------

;Temporarily declare variables that may eventually be passed
;-------------------------------------------------------------
arg='20040922.01'
mode='Stratus'
average=30		;Average time in seconds
site='nsa' 		;*****Will have to eventually use a call to keyword_set()

;Get a list of files to read
;----------------------------
farg=site+'mmcrspecmom*'+arg+'*.cdf'
files=findfile(farg,count=numf)
if numf eq 0 then begin
	print,'No spectra data for '+arg
	return
endif

;Read in all files, extract desired mode, and concatenate
;---------------------------------------------------------
for f=0,numf-1 do begin
	fid=ncdf_open(files[f])
	ncdf_varget,fid,'base_time',bt
	ncdf_varget,fid,'time_offset',to
	ncdf_varget,fid,'ModeDescription',md
	ncdf_varget,fid,'ModeNum',mn
	ncdf_varget,fid,'NumFFT',nft
	ncdf_varget,fid,'NyquistVelocity',nyq
	ncdf_varget,fid,'Heights',hts
	ncdf_varget,fid,'Spectra',spec
	ncdf_close,fid
	tm=basetime_to_time(bt,to)			;the time array
	modeind=where(string(md) eq mode)		;the mode index
	modeind=modeind[0]
	height=hts[*,modeind]/1000.			;height in km
	numfft=nft[modeind]				;number of FFTs for the desired mode
	nyquist=nyq[modeind]				;the nyquist velocity [m/s] for the mode
	hwh=where(height gt 0 and height lt 20,numh)	;the valid height indices
	twh=where(mn eq modeind)			;the time indices for the desired mode
	;Cut down the necessary variables in height
	;-------------------------------------------
	height=height[hwh]
	spec=spec[*,hwh,*]
	;Concatenate fields if necessary
	;--------------------------------
	if f eq 0 then begin
		time=tm[twh]
		spectra=spec[*,*,twh]
	endif else begin
		time=[time,tm[twh]]
		spectra=[[[spectra]],[[spec[*,*,twh]]]]
	endelse
	;Create the velocity bins
	;-------------------------
	velbins=findgen(numfft)/numfft*2.0*nyquist-nyquist
	numt=n_elements(time)
endfor

;Average data if requested
;--------------------------
if average gt 0 then begin	;May replace with keyword_set() call
	numta=((time[numt-1]-time[0])*3600.)/average	;Number of averaged beams
	spc=fltarr(numfft,numh,numta)*0.0
	numspecavg=intarr(numta)*0
	ttime=fltarr(numta)
	for i=0,numta-1 do begin
		wh=where(time ge time[0]+(i*average)/3600. and time lt time[0]+((i+1)*average)/3600.,num)
		if num gt 0 then begin
			numspecavg[i]=num
			ttime[i]=mean(time[wh])
			spc[*,*,i]=total(spectra[*,*,wh],3)/num
		endif
	endfor
	spectra=spc
	time=ttime
	spc=0
endif
end