;-------------------------------------------------------
;Title:		ascos_readwbanddbz.pro
;Purpose:	Program to read WBAND reflectivity only.
;Note:    For use with single mode files
;
;Input:		sarg:  date/time string, YYYYMMDDHH
;			    earg:  data/time string, YYYYMMDDHH
;Output:	time: radar time axis, decimal hours
;			    height: radar height axis, km
;			    dbz: radar reflectivity, dBZ
;		
;Keywords:  directory: directory containing the data
;				  default is current working directory
;				  factor:  factor to cut down time array by
;				  jtime: decimal year day
;				  stnfilter:  Set to the filter level for applying to reflectivity
;Author:	Matthew Shupe
;Date:		8/5/05, Adapted 7/7/2014
;-------------------------------------------------------

pro acse_readwbanddbz,sarg,earg,time,height,dbz,directory=directory,factor=factor,jtime=jtime,stnfilter=stnfilter

if keyword_set(directory) then cd,directory,current=orig_dir 
base_date=julday(1,1,1970,0,0,0)
base=jul_to_dt(base_date)
!dt_base=base

;----------------------------------
;Determine a list of files to read
;----------------------------------
sjd=julday(strmid(sarg,4,2),strmid(sarg,6,2),strmid(sarg,0,4))-julday(1,1,strmid(sarg,0,4))+1
ejd=julday(strmid(earg,4,2),strmid(earg,6,2),strmid(earg,0,4))-julday(1,1,strmid(earg,0,4))+1
if strmid(sarg,6,2) ne strmid(earg,6,2) then begin
	if ejd lt sjd then begin
		print,'acse_readwband:  cannot go across a new year.'
		stop
	endif
	numa=ejd-sjd+1
	darg=strarr(numa)
	for i=0,numa-1 do begin
		darg[i]=strmid(sarg,0,4)+strtrim(sjd+i,2)
	endfor
endif else begin
	darg=strmid(sarg,0,4)+strtrim(sjd,2)
	numa=n_elements(darg)
endelse

;----------------------------
;Read in the daily *mom*.nc
;----------------------------
;Get the height array
farg=darg[0]+'*MMCRMom.nc'
momfile=findfile(farg,count=numf)
if numf eq 0 then begin 
  time=-999
  print,'acse_readwband:  No W-band data for requested interval.'
  return
endif
fid=ncdf_open(momfile[0])
ncdf_varget,fid,'Heights',ht
ncdf_close,fid
height=ht[*,0]/1000.
hwh=where(height ge 0 and height lt 13,numh)
height=height[hwh]

len=numa*200000L
time=fltarr(len)*0-999.
dbz=fltarr(len,numh)

fill=0
ind=0L
for d=0,numa-1 do begin
	farg=darg[d]+'*MMCRMom.nc'
	momfile=findfile(farg,count=numf)
	momfile=momfile[sort(long(strmid(momfile,0,11)))]
	for f=0,numf-1 do begin
	  print,'Reading ',momfile[f]
		fid=ncdf_open(momfile[f])
		ncdf_varget,fid,'base_time',bt
		ncdf_varget,fid,'time_offset',to			;sec
		ncdf_varget,fid,'Reflectivity',Refl			;dBZ
		ncdf_varget,fid,'SignalToNoiseRatio',Sig		;dB
		ncdf_close,fid
		;Make a time array in hours
		;---------------------------
		starttime=sec_to_dt(bt)
		tm=starttime.hour+(starttime.minute*60+starttime.second+to)/3600
		numt=n_elements(tm)
		
		if keyword_set(stnfilter) then begin
		  wh=where(Sig le stnfilter,nm)
		  if nm gt 0 then Refl[wh]=-999.
		endif
		
		;Add data to variables
		;----------------------
		time[ind:ind+numt-1]=tm+24.0*d
		dbz[ind:ind+numt-1,*]=transpose(Refl[hwh,*])
		ind=ind+numt
	endfor
endfor

;Trim data to the hour range requested
;----------------------------------------
shr=fix(strmid(sarg,8,2))
ehr=fix(strmid(earg,8,2))
wh=where(time ge shr and time le ehr+24.0*(ejd-sjd),nm)
if nm gt 0 then begin
	time=time[wh]
	dbz=dbz[wh,*]
endif

;Cut down the data rate by a specified factor
;----------------------------------------------
if keyword_set(factor) then begin
  numt=n_elements(time)
  wh=lindgen(numt/factor)*factor
  dbz=dbz[wh,*]
  time=time[wh]
endif

;make yearday time
jtime=floor(sjd)+time/24.0

if keyword_set(directory) then cd,orig_dir
end ;prog