;Title:		stat_scs.pro
;
;Purpose:	To calculate all desired statistics from the combined
;			radar and lidar datasets.
;
;I/O:
;
;Inputs:
;Outputs:
;
;Author:	Matthew Shupe
;Date:		2/14/00
;Modified:	10/30/00
;--------------------------------------------------------------------

pro stat_scs

;Set up common blocks
;---------------------
common com_rws
common com_in
common com_out

;initialize the output variables (abbreviated names)
; and assign default/bad data values.
;---------------------------------------------------
lbr=fltarr(144)-999.;lowest base radar
lbl=lbr				;lowest base lidar (cloud)
lbe=lbr				;lowest base echo
lbi=intarr(144)-1	;lowest base echo instrument
lbet=lbr			;lowest base echo temperature
lblt=lbr			;lowest base lidar (cloud) temperature
lbb=lbr				;lowest base best (lidar when available, if not then radar)
lbbt=lbr			;lowest base best temperature
htr=lbr				;highest top radar
htl=lbr				;highest top lidar
hte=lbr				;highest top echo
hti=lbi				;highest top instrument
htet=lbr			;highest top echo temperature
nlr=lbi				;max. # layers radar
nll=lbi				;max. # layers lidar
nlc=lbi				;max. # layers combined
nli=lbi				;max. # layers instrument
ctr=lbr				;total cloud thickness radar
ctl=lbr				;total cloud thickness lidar
atl=lbi*0			;attenuation lidar
lpl=lbi				;liquid present lidar (0.11)
lhl=lbr				;height of lowest liquid layer (0.11)
lhlt=lbr			;height of lowest liquid layer temperature (0.11)
lpl2=lbi			;liquid present lidar (0.16)
lhl2=lbr			;height of lowest liquid layer (0.16)
lhlt2=lbr			;height of lowest liquid layer temperature (0.16)
pre=lbi*0			;precipitation
iur=intarr(144)+1	;instrument up radar
iul=iur				;instrument up lidar
cpr=lbi*0			;cloud presence


;Calculate the stats that require looping over time
;---------------------------------------------------
for i=0,143 do begin
	;Cloud base statistics (lbe,lbi)
	; Radar lowest --> lbi=0
	; Lidar lowest --> lbi=1
	; Both w/i 45m --> lbi=2
	; Both no data --> lbi=-1
	;-------------------------------
	if ((lbase[i,0] lt rbase[i,0]) and (lbase[i,0] gt 0)) or $
			((rbase[i,0] le 0) and (lbase[i,0] ge 0)) then begin
		lbe[i]=lbase[i,0]
		lbi[i]=1
	endif else begin
		lbe[i]=rbase[i,0]
		lbi[i]=0
	endelse
	if abs(lbase[i,0]-rbase[i,0]) le 0.045 then lbi[i]=2
	if (rbase[i,0] lt 0) and (lbase[i,0] lt 0) then begin
		lbe[i]=-999.
		lbi[i]=-1
	endif

	;Cloud top statistics (htr,htl,hte,hti)
	; Radar highest --> hti=0
	; Lidar highest --> hti=1
	; Both w/i 45m  --> hti=2
	; Both no data  --> hti=-1
	;--------------------------------------
	htr[i]=max(rtop[i,*])
	htl[i]=max(ltop[i,*])
	if htl[i] gt htr[i] then begin
		hte[i]=htl[i]
		hti[i]=1
	endif else begin
		hte[i]=htr[i]
		hti[i]=0
	endelse
	if abs(htl[i]-htr[i]) le 0.045 then hti[i]=2
	if (htr[i] lt 0) and (htl[i] lt 0) then begin
		hte[i]=-999.
		hti[i]=-1
	endif

	;Number of layers statistics (nlr,nll,nlc,nli)
	; Radar most layers --> nli=0
	; Lidar most layers --> nli=1
	; Both same layers --> nli=2
	; Both no data --> nli=-1
	;---------------------------------------------
	lwh=where(lbase[i,*] gt 0,lnum)
	rwh=where(rbase[i,*] gt 0,rnum)
	if rbase[i,0] ne -999. then nlr[i]=rnum else nlr[i]=-1
	if lbase[i,0] ne -999. then nll[i]=lnum else nll[i]=-1
	case 1 of
		(rnum gt lnum) or ((nlr[i] ge 0) and (nll[i] lt 0)): begin
			nlc[i]=rnum
			nli[i]=0
		end
		(rnum lt lnum) or ((nlr[i] lt 0) and (nll[i] ge 0)): begin
			nlc[i]=lnum
			nli[i]=1
		end
		(rnum eq lnum): begin
			nlc[i]=rnum
			nli[i]=2
			if ((nll[i] lt 0) and (nlr[i] lt 0)) then begin
				nli[i]=-1
				nlc[i]=-1
			endif
		end
	endcase

	;Lidar attunuation (atl)
	; [if high radar top > 1km --> lidar attenuated if seeing >200m lower]
	; [if high radar top < 1km --> lidar attenuated if seeing >80m lower]
	; Not attenuated --> atl=0
	; Attenuated --> atl=1
	; No lidar or radar data --> atl=-1
	;---------------------------------------------------------------------
	if htl[i] gt 0 then begin
		case 1 of
			(htr[i] lt 0): atl[i]=-1
			(htr[i] ge 1.): if htl[i] lt (htr[i]-0.2) then atl[i]=1 else atl[i]=0
			(htr[i] lt 1.): if htl[i] lt (htr[i]-0.08) then atl[i]=1 else atl[i]=0
		endcase
	endif else atl[i]=-1
	if htl[i] eq 0 then atl[i]=0

	;Lidar liquid presence (lpl)
	; [liquid present if depol ratio >0 and <=0.11]
	; [For lpl2, liquid pres if depol >0 and <=0.16]
	; No liquid --> lpl=0
	; Liquid --> lpl=1
	; No data --> lpl=-1
	;Height of liquid layer if present in km
	;----------------------------------------------
	qwh=where((ldepol[i,*] le 0.11) and (ldepol[i,*] gt 0),nliq)
	case 1 of
		(ldepol[i,0] eq -999.): lpl[i]=-1
		(nliq gt 0): begin
			lpl[i]=1
			lhl[i]=lbase[i,qwh[0]]
		end
		(nliq eq 0): lpl[i]=0
	endcase
	qwh=where((ldepol[i,*] le 0.16) and (ldepol[i,*] gt 0),nliq)
	case 1 of
		(ldepol[i,0] eq -999.): lpl2[i]=-1
		(nliq gt 0): begin
			lpl2[i]=1
			lhl2[i]=lbase[i,qwh[0]]
		end
		(nliq eq 0): lpl2[i]=0
	endcase

	;Cloud presence
	; No clouds --> cpr=0
	; Clouds --> cpr=1
	; No data --> cpr=-1
	;--------------------
	case 1 of
		(hte[i] gt 0): cpr[i]=1
		(hte[i] eq 0): cpr[i]=0
		(hte[i] lt 0): cpr[i]=-1
	endcase
endfor

;Calculate the stats that can be done with array operations
;----------------------------------------------------------

	;Instruments up? (iur,iul)
	; instrument functioning --> 1
	; instrument down --> -1
	;-----------------------------
	rdn=where(rbase[*,0] lt 0,nrad)
	if nrad gt 0 then iur[rdn]=-1
	ldn=where(lbase[*,0] lt 0,nlid)
	if nlid gt 0 then iul[ldn]=-1

	;Cloud thickness (ctr,ctl)
	; [Total combined thickness of all layers]
	;-----------------------------------------
	ctr=total(rtop-rbase,2)
	ctl=total(ltop-lbase,2)
	if nrad gt 0 then ctr[rdn]=-999.
	if nlid gt 0 then ctl[ldn]=-999.

	;Lowest cloud bases (lbr,lbl)
	;----------------------------
	lbr=rbase[*,0]
	lbl=lbase[*,0]
	lbb=lbl
	iwh=where(lbl le 0)
	if iwh[0] ne -1 then begin
		jwh=where(lbr[iwh] gt 0)
		if jwh[0] ne -1 then lbb[iwh[jwh]]=lbr[iwh[jwh]]
	end

	;Precipitation present (pre)
	; Precip --> pre=1
	; No precip --> pre=0
	; One/both instruments down --> -1
	; [*****check this out in testing.......]
	; [**What if r sees a lower layer that l doesn't?]
	;--------------------------------------------------
	pwh=where((lbase[*,0]-rbase[*,0]) gt 0.09)
	if pwh[0] ne -1 then pre[pwh]=1
	if nrad gt 0 then pre[rdn]=-1
	if nlid gt 0 then pre[ldn]=-1

	;Base/top temperatures (lbet,lbct,htet,lbbt,lhlt,lhlt2)
	;------------------------------------------------------
	cup=where(lbe ge 0)
	if cup[0] ne -1 then begin
		lbet[cup]=rstemp_ret(lbe[cup],jtim[cup])
		lbbt[cup]=rstemp_ret(lbb[cup],jtim[cup])
		htet[cup]=rstemp_ret(hte[cup],jtim[cup])
	endif
	cup=where(lbl ge 0)
	if cup[0] ne -1 then lblt[cup]=rstemp_ret(lbl[cup],jtim[cup])
	cup=where(lhl ge 0)
	if cup[0] ne -1 then lhlt[cup]=rstemp_ret(lhl[cup],jtim[cup])
	cup=where(lhl2 ge 0)
	if cup[0] ne -1 then lhlt2[cup]=rstemp_ret(lhl2[cup],jtim[cup])


end ;prog