;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