;-------------------------------------------------------------------------
;Title:		lai_scs.pro
;
;Purpose:	To load and interpolate all radar, lidar and radiosonde data.
;
;I/O:		lai_scs,jday
;
;Inputs:	jday - the julian day (from IDL>julday(mm,dd,yy) )
;
;Outputs:
;
;Author:	Matthew Shupe
;Date:		2/10/00
;-------------------------------------------------------------------------


;The SCS_FY.PRO program
;----------------------
pro lai_scs,jday

;Set up common block variables
;-----------------------------
common com_rws,rstemp,rsheights,rsjtime
common com_in,tim,jtim,lbase,ltop,ldepol,rbase,rtop
common com_raddn,raddn

;Convert the Julian day into arguments for radar, lidar, and radiosonde
;-----------------------------------------------------------------------
caldat,jday,mo,da,yr
months=['jan','feb','mar','apr','may','jun','jul','aug',$
				'sep','oct','nov','dec']
rmo=months[mo-1]
mo=strtrim(mo,2)
if strlen(mo) eq 1 then mo='0'+mo
da=strtrim(da,2)
if strlen(da) eq 1 then da='0'+da
yr=strtrim(yr,2)

rarg=da+rmo+yr+'*.nc'
larg='*'+yr+'-'+mo+'-'+da+'*.avge'
sarg='*'+strmid(yr,2,2)+mo+da+'.dat'

;Load radar data
; (rtime in decimal hours)
; (rjtime in julian day from 1/1/97)
; (heights in km)
;-----------------------------------
cd,'d:/shebaar/cloudbounds',current=odir
rfile=findfile(rarg,count=rnum)
base_date = julday(1,1,1970,0,0,0)
base=jul_to_dt(base_date)
!dt_base=base
if rnum eq 0 then begin
	rtime=findgen(144)*1/6.
	rbh=fltarr(144)-999.
	rth=rbh
	rjtime=float(floor((jday-julday(1,1,1997,0,0,0)+1)))+rtime/24
endif else begin
	for ir=0,rnum-1 do begin
		print,rfile[ir]
		fid=ncdf_open(rfile[ir])
		ncdf_varget,fid,0,base_time
  		ncdf_varget,fid,1,time_offset
  		ncdf_varget,fid,2,rrb
  		ncdf_varget,fid,3,rrt
  		ncdf_close,fid
  		rrt=transpose(rrt)/1000. 	;puts into km
  		rrb=transpose(rrb)/1000.
  		;Create a julian time axis
  		;----------------------------
  		starttime=sec_to_dt(base_time)
  		rti=starttime.hour+(starttime.minute*60+starttime.second+time_offset)/3600
  		jti=float(floor((starttime.julian-julday(1,1,1997,0,0,0)+1)))+rti/24
		if ir eq 0 then begin
			rbh=rrb
			rth=rrt
			rtime=rti
			rjtime=jti
		endif else begin
			rbh=[rbh,rrb]
			rth=[rth,rrt]
			rtime=[rtime,rti]
			rjtime=[rjtime,jti]
		endelse
	endfor
endelse

;Interpolate radar to 10-min axes
; (tim in decimal hours, 10-min exact)
; (jtim in julian from 1/1/97, 10-min exact)
;--------------------------------
tim=findgen(144)*1/6.
jtim=float(fix(rjtime[0]))+tim/24.
tenmininterp,tim,rtime,rbh,rth,rbase,rtop

;Fill in where radar is down with lookup table
;---------------------------------------------
for r=0,143 do begin
	iwh=where(raddn eq jtim[r])
	if iwh[0] ne -1 then begin
		rbase[r,*]=-999.
		rtop[r,*]=-999.
	endif
endfor

;Load lidar data
; (heights in km)
; (correct for when lidar base=0 and top>0)
;----------------------------------------------
cd,'d:/lidar/cld'
lfile=findfile(larg,count=lnum)
if lnum eq 0 then begin
	ltop=fltarr(144,12)-999.
	lbase=ltop
	ldepol=ltop
endif else begin
	;Read the cld lidar files
	;------------------------
	CLD_read2,lfile[0],lti,ldepol,ltop,lbase,mxl
endelse
iwh=where((ltop le 0.050) and (ltop gt 0))
if iwh[0] ne -1 then begin
	ltop[iwh]=0
	lbase[iwh]=0
	ldepol[iwh]=0
endif
iwh=where((ltop gt 0) and (ltop le 0.090))
if iwh[0] ne -1 then begin
	lbase[iwh]=0.090
	ltop[iwh]=0.095
endif
iwh=where((lbase lt 0.090) and (ltop gt 0.090))
if iwh[0] ne -1 then lbase[iwh]=0.090

for ii=0,143 do begin
	if (lbase[ii,1] ne 0) and (lbase[ii,1] lt lbase[ii,0]) then begin
		iwh=where(lbase[ii,*] gt 0,niwh)
		sot=sort(lbase[ii,iwh])
		lbase[ii,*]=[transpose(lbase[ii,iwh[sot]]),fltarr(12-niwh)*0]
		ltop[ii,*]=[transpose(ltop[ii,iwh[sot]]),fltarr(12-niwh)*0]
		ldepol[ii,*]=[transpose(ldepol[ii,iwh[sot]]),fltarr(12-niwh)*0]
	endif
endfor

;Load radiosonde data and put into same format as radar/lidar data
; i.e. w/ jtim starting 1/1/97 and with heights in km
;------------------------------------------------------------------
cd,'d:/shebaar/rwsonde/raw'
sfile=findfile(sarg)
print,sfile
restore,sfile[0]
rsheights=rsheights/1000.  ;puts into km
if jday ge 2450815L then rsjtime=rsjtime+365.


cd,odir
end ;prog