;pro filter_radar_precip

; idea is to bring the 5min filename, read in the Robust mode radar file, identify periods with precip
; flip the mwr_flag as needed, and write the adjusted flag to the 5min file.

test_flag=0
if test_flag eq 0 then begin

print, 'opening fmin_files.dat'
openr, 1, 'fivemin_files.dat' & cpc_fname=' '
readf,1, cpc_fname, format='(a100)'
close, 1
cpc_fname=strtrim(cpc_fname,2)
print, cpc_fname

endif else begin

cpc_fname='/home/mace/sgp/sgp.average.5min.ciret4/sgp.average.5min.19970404.000000.temp.cdf'

endelse
if test_flag eq 0 then begin
openr, 1, 'Robust_files.dat' & Robust_fname=' '
readf,1, Robust_fname, format='(a100)'
close, 1
Robust_fname=strtrim(Robust_fname,2)

endif else begin

Robust_fname='/home/mace/sgp/sgpmmcrcalC1.a1/sgpmmcrRobustC1.a1.19970404.000058.cdf'

endelse

; get the 5min times and the mwr flag

cdf_id=ncdf_open(cpc_fname)
print, 'opened ',cpc_fname

time_did=ncdf_dimid(cdf_id, 'time')
baset_id=ncdf_varid(cdf_id, 'base_time')
timeo_id=ncdf_varid(cdf_id, 'time_offset')
height_id=ncdf_varid(cdf_id, 'height')
temp_id=ncdf_varid(cdf_id, 'temp')
precip_flag_id=ncdf_varid(cdf_id, 'precip_flag')
nlayers_id=ncdf_varid(cdf_id, 'number_of_layers')
bzl1_id=ncdf_varid(cdf_id, 'base_height_layer_1')
if bzl1_id ne -1 then ncdf_varget, cdf_id, bzl1_id, bzl1
ncdf_varget, cdf_id, baset_id, base_time
ncdf_varget, cdf_id, timeo_id, time_offset
ncdf_varget, cdf_id, height_id, height
ncdf_varget, cdf_id, temp_id, temp
ncdf_varget, cdf_id, precip_flag_id, precip_flag
ncdf_varget, cdf_id, nlayers_id, n_layers
ncdf_close,cdf_id

ijday_ihrfrac_fm_numsec,base_time,time_offset,hrfrac5,jday5,yy5,mm5,dd5,hh5,mi5,ss5

; get the robust file data

cdf_id=ncdf_open(Robust_fname)

time_did=ncdf_dimid(cdf_id, 'time')
baset_id=ncdf_varid(cdf_id, 'base_time')
timeo_id=ncdf_varid(cdf_id, 'time_offset')
height_id=ncdf_varid(cdf_id, 'Heights')
z_id=ncdf_varid(cdf_id, 'ReflectivityBestEstimate')
v_id=ncdf_varid(cdf_id, 'MeanDopplerVelocity')
ncdf_varget, cdf_id, baset_id, base_time
ncdf_varget, cdf_id, timeo_id, time_offset
ncdf_varget, cdf_id, height_id, heightR
ncdf_varget, cdf_id, z_id, dbz
ncdf_varget, cdf_id, v_id, vel
ncdf_close,cdf_id

dbz[where(dbz eq 0)]=-9999
vel[where(vel eq 0)]=-9999

heightR=heightR-(heightR[0]-height[0])

ijday_ihrfrac_fm_numsec,base_time,time_offset,hrfracR,jdayR,yyR,mmR,ddR,hhR,miR,ssR

; step through the 5min data

for j=0,n_elements(hrfrac5)-1 do begin

    frz_ind=11
    while temp[frz_ind] ge 273. do frz_ind=frz_ind+1   ; the heights of the Robust and 5min files should be the same

    count_precip=0 & brightband=0
    if n_layers[j] gt 0 then begin

       time_dif_index=where(abs(jdayR-jday5[j]) lt 5.d/1440.d)

         for jj=0,n_elements(time_dif_index)-1 do begin
       kk=0
       while heightR[kk] le bzl1[j] and kk lt n_elements(heightR)-2 do kk=kk+1
       for k=0,frz_ind do begin

         if float(dbz[k,time_dif_index[jj]])/100. gt 5.0 and float(vel[k,time_dif_index[jj]])/1000. gt 3.0 then begin
             count_precip=count_precip+1
            ; print, 'precip found at ',hrfrac5[j],heightR[k],float(dbz[k,time_dif_index[jj]])/100.,float(vel[k,time_dif_index[jj]])/1000.
         endif
         if float(dbz[k,time_dif_index[jj]])/100. gt 15.0 then begin
             count_precip=count_precip+1
            ; print, 'precip found at ',hrfrac5[j],heightR[k],float(dbz[k,time_dif_index[jj]])/100.
         endif

       endfor    ; for k=0,frz_ind do begin

       count_dbz_below=0 & count_dbz_abv=0 & count_dbz_in=0 & dbz_below_sum=0. & dbz_abv_sum=0. & dbz_in_sum=0.

       
		k=frz_ind-10
	
	while k le frz_ind-5 do begin
         if float(dbz[k,time_dif_index[jj]])/100. gt -30.0 and heightR[k] lt height[frz_ind-10] then begin
          count_dbz_below=count_dbz_below+1
          dbz_below_sum=dbz_below_sum+float(dbz[k,time_dif_index[jj]])/100.
         endif
         k=k+1
       endwhile
       k=frz_ind+5
       while k le frz_ind+15 do begin
         if float(dbz[k,time_dif_index[jj]])/100. gt -30.0 then begin
          count_dbz_abv=count_dbz_abv+1
          dbz_abv_sum=dbz_abv_sum+float(dbz[k,time_dif_index[jj]])/100.
         endif
         k=k+1
       endwhile
       k=frz_ind-5
       while k le frz_ind+5 do begin
         if float(dbz[k,time_dif_index[jj]])/100. gt -30.0 then begin

          count_dbz_in=count_dbz_in+1
          dbz_in_sum=dbz_in_sum+float(dbz[k,time_dif_index[jj]])/100.
         endif
         k=k+1
       endwhile

       if count_dbz_in gt 0 then begin
         if count_dbz_abv gt 0 then begin
          if dbz_in_sum/float(count_dbz_in) gt 5.+(dbz_abv_sum/float(count_dbz_abv)) and dbz_in_sum/float(count_dbz_in) gt -5. then begin
            brightband=1
           if hrfrac5[j] ge 13.25 and hrfrac5[j] lt 14. then   begin
            print, 'brightband abv ', hrfrac5[j],dbz_in_sum/float(count_dbz_in),(dbz_abv_sum/float(count_dbz_abv))
         endif
         endif
         endif
         if count_dbz_below gt 0 then begin
          if dbz_in_sum/float(count_dbz_in) gt 5.+(dbz_below_sum/float(count_dbz_below)) and dbz_in_sum/float(count_dbz_in) gt -5. then begin
            brightband=1
           print, 'brightband below ', hrfrac5[j],dbz_in_sum/float(count_dbz_in),(dbz_below_sum/float(count_dbz_below))
         endif
         endif
       endif ;   if count_dbz_in gt 0 then begin

         endfor   ; for jj=time_dif_index[0],n_elements[time_dif_index]-1 do begin

    if brightband eq 1 or count_precip gt 5 then begin
        precip_flag[j]=1
        print, 'flag time due to precip ',hrfrac5[j],brightband, count_precip
    endif

    endif  ; if n_layers gt 0 then begin

endfor  ; for j=0,n_elements(hrfracR) do begin

; write out the precip flag back to the 5min file.

print, 'writing to ',cpc_fname
cdf_id=ncdf_open(cpc_fname, /write)


precip_flag_id=ncdf_varid(cdf_id, 'precip_flag')

ncdf_varput, cdf_id, precip_flag_id, precip_flag

ncdf_close,cdf_id
print,  'wrote to ',cpc_fname


end