;Title:			acse_readrds
;Purpose:		Read ACSE radiosonde files
;
;Inputs:		sarg:  date/time string, YYYYMMDDHH
;				    earg:  data/time string, YYYYMMDDHH
;Outputs:		time:  time of sounding(s) [hours]
;				    height: height of measurement [km]
;				    pres:  pressure [mb]
;				    temp:  temperature [C]
;				    dewpt: dew point temperature [C]
;				    rh: relative humidity [%]
;				    speed:  wind speed [m/s]
;				    direction:  wind direction [deg]
;Keywords:  directory:  directory containing the data
;				    htinterp:  height array to which to interpolate data [km]
;				    htlim:  upper height limit above which data are removed [km]
;				    jtime:  time in decimal julian days [days]
;I/O:
;Author:		Matthew Shupe
;Date:			7/15/08, adapted 7/7/2014
;---------------------------------------------------

pro acse_readrds,sarg,earg,time,height,pres,temp,dewpt,rh,speed,direction,directory=directory,htinterp=htinterp,htlim=htlim,jtime=jtime

if keyword_set(directory) then cd, directory, current=orig_dir
if not keyword_set(htlim) then htlim=50.

;Determine the files to open
sjd=julday(strmid(sarg,4,2),strmid(sarg,6,2),strmid(sarg,0,4))
ejd=julday(strmid(earg,4,2),strmid(earg,6,2),strmid(earg,0,4))
if strmid(sarg,6,2) ne strmid(earg,6,2) then begin
	numa=ejd-sjd+1
	darg=strarr(numa)
	for i=0,numa-1 do begin
		caldat,sjd+i,mo,da,yr
		if mo lt 10 then mo='0'+strtrim(mo,2) else mo=strtrim(mo,2)
		if da lt 10 then da='0'+strtrim(da,2) else da=strtrim(da,2)
		yr=strtrim(yr,2)
		darg[i]=yr+mo+da
	endfor
endif else begin
	darg=strmid(sarg,0,8)
	numa=n_elements(darg)
endelse
;Add the file from the end of the day
caldat,ejd+1,mo,da,yr
if mo lt 10 then mo='0'+strtrim(mo,2) else mo=strtrim(mo,2)
if da lt 10 then da='0'+strtrim(da,2) else da=strtrim(da,2)
yr=strtrim(yr,2)
darg=[darg,yr+mo+da]
numa=numa+1
sjd=sjd-julday(1,1,2014)+1
ejd=ejd-julday(1,1,2014)+1

;Initialize arrays
st1=strarr(20)
jnk=' '
time=fltarr(6*numa)
temp=fltarr(6*numa,10000)*0-999
dewpt=temp
height=temp
rh=temp
pres=temp
speed=temp
direction=temp

;Loop over all daily arguments
indf=0
mxind=0
for d=0,numa-1 do begin
	;Get the list of files for the daily argument
	arg='FLEDT_*'+darg[d]+'*.tsv'
	file=findfile(arg,count=numf)
	file=file[sort(long(strmid(file,6,8))*100+long(strmid(file,15,2)))]
	;Loop over all files in the day
	for f=0,numf-1 do begin
		ind=0
		time[indf]=fix(strmid(file[f],15,2))+d*24
		print,'Reading ',file[f]
		openr,lun,/get_lun,file[f]
		;print,file[f]
		;read header
		for i=0,39 do readf,lun,jnk
		;read rest of file
		while not eof(lun) do begin
		  readf,lun,jnk
		  st1=strsplit(jnk,/regex,/extract)
		  temp[indf,ind]=st1[2]
		  height[indf,ind]=st1[6]
		  dewpt[indf,ind]=st1[8]
		  rh[indf,ind]=st1[3]
		  pres[indf,ind]=st1[7]
		  speed[indf,ind]=st1[11]
		  direction[indf,ind]=st1[10]
			ind=ind+1
		endwhile
		free_lun,lun
		if ind gt mxind then mxind=ind
		indf=indf+1
	endfor
endfor
if indf gt 0 then begin
	time=time[0:indf-1]
	temp=temp[0:indf-1,*]
	dewpt=dewpt[0:indf-1,*]
	height=height[0:indf-1,*]
	rh=rh[0:indf-1,*]
	pres=pres[0:indf-1,*]
	speed=speed[0:indf-1,*]
	direction=direction[0:indf-1,*]
	temp=temp[*,0:mxind-1]-273.16
	dewpt=dewpt[*,0:mxind-1]-273.16
	height=height[*,0:mxind-1]/1000.
	rh=rh[*,0:mxind-1]
	pres=pres[*,0:mxind-1]
	speed=speed[*,0:mxind-1]
	direction=direction[*,0:mxind-1]
	wh=where(height lt 0 or height gt htlim)
	height[wh]=!values.f_nan
	temp[wh]=!values.f_nan
	rh[wh]=!values.f_nan
	dewpt[wh]=!values.f_nan
	pres[wh]=!values.f_nan
	speed[wh]=!values.f_nan
	direction[wh]=!values.f_nan
endif else begin
	time=[0,6,12,18]
	height=fltarr(4,20)
	temp=height*0-999. & dewpt=height*0-999. & rh=height*0-999. & pres=height*0-999. & speed=height*0-999. & direction=height*0-999.
	for i=0,3 do height[i,*]=findgen(20)/2.0
endelse

;Interpolate the data to height grid if requested
;  (make values constant if they are outside  of the time range
if keyword_set(htinterp) then begin
	indf=indf>1   ;If indf=0 (i.e. no soundings) then there is one dummy)
	nht=n_elements(htinterp)
	temp0=fltarr(indf,nht)*0-999.
	dewpt0=fltarr(indf,nht)*0-999.
	rh0=fltarr(indf,nht)*0-999.
	pres0=fltarr(indf,nht)*0-999.
	speed0=fltarr(indf,nht)*0-999.
	direction0=fltarr(indf,nht)*0-999.
	for i=0,indf-1 do begin
		wh=where(height[i,*] gt 0,nm)
		if nm gt 0 then begin
			temp0[i,*]=interpol(temp[i,wh],height[i,wh],htinterp)
			dewpt0[i,*]=interpol(dewpt[i,wh],height[i,wh],htinterp)
			rh0[i,*]=interpol(rh[i,wh],height[i,wh],htinterp)
			pres0[i,*]=interpol(pres[i,wh],height[i,wh],htinterp)
			speed0[i,*]=interpol(speed[i,wh],height[i,wh],htinterp)
			direction0[i,*]=interpol(direction[i,wh],height[i,wh],htinterp)
		endif
	endfor
	temp=temp0
	dewpt=dewpt0
	rh=rh0
	pres=pres0
	speed=speed0
	direction=direction0
	height=htinterp
endif

;Trim data to the hour range requested
shr=fix(strmid(sarg,8,2))
ehr=fix(strmid(earg,8,2))
sda=fix(strmid(sarg,6,2))
eda=fix(strmid(earg,6,2))
wh=where(time ge shr and time le ehr+24.0*(ejd-sjd),nm)
if nm gt 0 then begin
	time=time[wh]
	if not keyword_set(htinterp) then height=height[wh,*]
	temp=temp[wh,*]
	dewpt=dewpt[wh,*]
	rh=rh[wh,*]
	pres=pres[wh,*]
	speed=speed[wh,*]
	direction=direction[wh,*]
endif else begin
	time=float(shr)
	if not keyword_set(htinterp) then height=transpose(findgen(20)/2.0)
	temp=height*0-999. & dewpt=height*0-999. & rh=height*0-999. & pres=height*0-999. & speed=height*0-999. & direction=height*0-999.
endelse

;compute julian time
jtime=floor(sjd)+time/24.0

if keyword_set(directory) then cd, orig_dir
end