Pro readradar,arg,NODATA=nodata,CUT=cut,NO_NAN=no_nan,AVERAGE=average,MAKEFULLDAY=makefullday
;--------------------------------------------------------------------
;Title:		readradar.pro
;Purpose:	Procedure to load merged radar data netCDF files and put
;			them into a usable format.
;			Compatible with new merge files (not SHEBA)
;
;Inputs:		arg - date for the file in "YYYYMMDD" or "YYYYMMDD.HH"
;
;Outputs:		time - time in decimal hours
;			jtime - julian time axis 
;			height - heights  (km) 
;			dbz - reflectivity (dbz) 
;			dopp - velocity (m/s)
;			bheight - cloud base heights (km)
;			theight - cloud top heights (km)
;			spwd - spectral width (m/s)
;			stn - signal to noise ratio (dbz)
;
;Keywords:	nodata - returns 0 if data present and 1 if no data
;			cut - Allows for cutting the data down by the factor given.
;				Takes every "cut"th element.
;			no_nan - if set then cloud boundaries with value 0 are left alone,
;					otherwise they are set to NaN.
;			average - set to the width (in seconds) over which to average the data
;			makefullday - only  works if average is also set.  If set then averages are done over a full 24 hour day with even blocks of time.
;
;I/O Format:	readradar,'20000118'
;
;Author:		Matthew Shupe
;Date:		2/11/99
;Last modified: 8/27/00
;			3/5/01	Added the AVERAGE keyword and MAKEFULLDAY keyword.
;-------------------------------------------------------------------------
;-------------------------------------------------------------------------

;Declare common block variables
;-------------------------------
common com_direct
common com_const
common com_platforms
common com_radar,time,jtime,height,dbz,dopp,spwd,stn,bheight,theight

base_date=julday(1,1,1970,0,0,0)
base=jul_to_dt(base_date)
!dt_base=base

cd,radardir,current=orig_dir
if site eq '' then site='*'
marg=strlowcase(site)+'*'+arg+'*.cdf'
file=findfile(marg,count=nfil)
radarplat=strmid(file[0],0,strpos(file[0],'.'+arg))

if nfil gt 0 then nodata=0 else begin
	print,'No file named ',marg
	nodata=1
	return
endelse

for i=0,nfil-1 do begin

	print,'Loading '+file[i]
	;Load netCDF data
	;--------------------
	fid = ncdf_open(file[i])
	ncdf_varget,fid,ncdf_varid(fid,'base_time'),base_time
	ncdf_varget,fid,ncdf_varid(fid,'time_offset'),time_offset
	if ncdf_varid(fid,'Heights') ne -1 then ncdf_varget,fid,ncdf_varid(fid,'Heights'),height else ncdf_varget,fid,ncdf_varid(fid,'heights'),height
	ncdf_varget,fid,ncdf_varid(fid,'Reflectivity'),db
	ncdf_varget,fid,ncdf_varid(fid,'MeanDopplerVelocity'),dp
	ncdf_varget,fid,ncdf_varid(fid,'SpectralWidth'),sw
	if ncdf_varid(fid,'SignaltoNoiseRatio') ne -1 then  ncdf_varget,fid,ncdf_varid(fid,'SignaltoNoiseRatio'),sn else $
			ncdf_varget,fid,ncdf_varid(fid,'SignalToNoiseRatio'),sn
	if ncdf_varid(fid, 'CloudLayerBottomHeightMplCloth') ne -1 then  ncdf_varget,fid,ncdf_varid(fid,'CloudLayerBottomHeightMplCloth'),bh else $
			 ncdf_varget,fid,ncdf_varid(fid,'CloudLayerBottomHeight'),bh
	if ncdf_varid(fid, 'CloudLayerTopHeightMplCloth') ne -1 then ncdf_varget,fid,ncdf_varid(fid,'CloudLayerTopHeightMplCloth'),th else $
			ncdf_varget,fid,ncdf_varid(fid,'CloudLayerTopHeight'),th
	ncdf_close,fid

	;Convert the data to useful values
	;-----------------------------------
	db=transpose(db)/100.			;puts into dBz
	height=height/1000.				;puts into km
	dp=transpose(dp)/1000.			;puts into m/s
	sw=transpose(sw)/1000.			;puts into m/s
	sn=transpose(sn)/100.			;puts into dBZ
	bh=transpose(bh)/1000.			;puts into km
	th=transpose(th)/1000.			;puts into km

	;Create a julian time axis
	;----------------------------
	starttime=sec_to_dt(base_time)
	ti=starttime.hour+(starttime.minute*60+starttime.second+time_offset)/3600
	jti=float(floor((starttime.julian-julday(1,1,starttime.year,0,0,0)+1)))+ti/24

	if i eq 0 then begin
		jtime=jti
		time=ti
		dbz=db
		dopp=dp
		spwd=sw
		stn=sn
		bheight=bh
		theight=th
	endif else begin
		jtime=[jtime,jti]
		time=[time,ti]
		dbz=[dbz,db]
		dopp=[dopp,dp]
		spwd=[spwd,sw]
		stn=[stn,sn]
		bheight=[bheight,bh]
		theight=[theight,th]
	endelse
endfor

;Cut the height down to only useful range gates (13.5 km for NSA)
;------------------------------------------------
iwh=where(height le float(topheight))
height=height[iwh]
dbz=dbz[*,iwh]
dopp=dopp[*,iwh]
spwd=spwd[*,iwh]
stn=stn[*,iwh]

;Cut the time frequency down if keyword is set
;---------------------------------
if keyword_set(cut) then begin
	iwh=indgen(floor(n_elements(time)/float(cut)))*fix(cut)
	jtime=jtime[iwh]
	time=time[iwh]
	dbz=dbz[iwh,*]
	dopp=dopp[iwh,*]
	spwd=spwd[iwh,*]
	stn=stn[iwh,*]
	bheight=bheight[iwh,*]
	theight=theight[iwh,*]
endif

;Average the data if the AVERAGE keyword is set
;------------------------------------
numheights=n_elements(height)
if keyword_set(AVERAGE) then begin
	width=float(average)/3600.		;puts average width into hours
	mfd=0
	if keyword_set(makefullday) then begin
		newbeams=ceil(24/width)
		mfd=1
	endif else newbeams=ceil((time[n_elements(time)-1]-time[0])/width)      
	tmst=float(fix((average-1)/10))*5.0/3600.0
	jday=floor(jtime[0])
	db=fltarr(newbeams,numheights)*0-999.0
	dp=db
	sw=db
	sn=db
	bh=fltarr(newbeams,10)*0
	th=bh
	tm=fltarr(newbeams)
	jt=tm
	for i=0,newbeams-1 do begin
		if mfd then iwh=where(time ge i*width and time lt (i+1)*width,num1) else iwh=where(time ge i*width+time[0] and time lt (i+1)*width+time[0],num1)
		if num1 gt 0 then begin
			if mfd then begin
				tm[i]=i*width+tmst
				jt[i]=jday+tm[i]/24.0
			endif else begin
				tm[i]=mean(time[iwh])
				jt[i]=mean(jtime[iwh])
			endelse
			numl=0
			incld=0
			for j=0,numheights-1 do begin
				jwh=where(dbz[iwh,j] gt -300, num2)
				if num2 ge 0.33*num1 then begin						;Must be at least 33% of pixels for average.
					db[i,j] = alog10(mean(10.0^(dbz[iwh[jwh],j]/10.0)))*10.0
					dp[i,j] = mean(dopp[iwh[jwh],j])
					sw[i,j] = mean(spwd[iwh[jwh],j])
					sn[i,j] = alog10(mean(10.0^(stn[iwh[jwh],j]/10.0)))*10.0
					if incld eq 0 and numl le 9 then begin

						bh[i,numl]=height[j]
						incld=1
					endif
				endif else begin
					if incld eq 1 then begin
						if height[j-1]-bh[i,numl] gt 0.050 then begin
							th[i,numl]=height[j-1]
							numl=numl+1
						endif
						incld=0
					endif
				endelse
			endfor
		endif else begin
			tm[i]=-999
			if mfd then begin
				tm[i]=i*width+tmst
				jt[i]=jday+tm[i]/24.0
			endif else tm[i]=-999
		endelse
	endfor
	iwh=where(tm ne -999)
	dbz=db[iwh,*]
	dopp=dp[iwh,*]
	spwd=sw[iwh,*]
	stn=sn[iwh,*]
	time=tm[iwh]
	jtime=jt[iwh]
	bheight=bh[iwh,*]
	theight=th[iwh,*]
	iwh=where(theight le 0,num)
	if num gt 0 then bheight[iwh]=0
endif

if not(keyword_set(NO_NAN)) then begin
	iwh=where(bheight eq 0)
	bheight[iwh]=!values.f_nan
	theight[iwh]=!values.f_nan
endif

cd,orig_dir
end