;
;  This program processes the two dimensional clutter flag and 
;  writes it to the output netcdf file.  Determines the number
;  of total data points, the number of cloudy data points,
;  mode, minimum, and maximum clutter flag value for each time
;  interval
;
;  note: missing  = -8888
;        no cloud = -7777
;

pro ctrue_cmin_cmax_crand3

common ncdf_file, oldcdfname,twoD,time_dim_name,height_dim_name,$
  height_var_name,avgfilename,curdate

common old_grid, num_times,num_heights,base_time,time_offset,$
  height,site,files,num_files,num_days

common new_grid, numtimes,newtimes, numheights,newheights,$
 timespacing,heightspacing

common bracket, startbst

common data, avg_data

;
;  Loop through all the matching files in the directory
;
m=-1  ;counts number of cdf files without the variable
p=-1  ;counts number of cdf files with variable
j=-1  ;counts number of cdf files
for i = 0, num_files-1 do begin
;
;  Don't try to process any gif files that may be in the directory
;
  if strpos(files[i],'gif') eq -1 and $
  strpos(files[i],'png') eq -1 then begin
    j=j+1
    print,i,j,files[i]
;
;  Open the netcdf file
;
    cdf_id=ncdf_open(files[i])
;
;  Get the dimension id's
;
    time_did=ncdf_dimid(cdf_id,time_dim_name)
    height_did=ncdf_dimid(cdf_id,height_dim_name)
;
;  Get the dimensions
;
    ncdf_diminq,cdf_id,time_did,charstring,ntimes
    ncdf_diminq,cdf_id,height_did,charstring,nheights
;
;  Create array for nodata value
;
    dummy=make_array(ntimes,nheights,/float,value=-8888)
;
;  Get the variable id
;
    clutter_id=ncdf_varid(cdf_id,'qc_ReflectivityClutterFlag')
    artifact_id=ncdf_varid(cdf_id,'qc_RadarArtifacts')
    reflect_id=ncdf_varid(cdf_id,'ReflectivityBestEstimate')
;
;  Get the variable if it exists
;
    if reflect_id ne -1 then begin
      print,'Processing 2-D ','qc_ReflectivityBestEstimate'
      ncdf_varget, cdf_id, reflect_id, reflect
      reflect=transpose(reflect)
      p=p+1
    endif else begin
      print, 'reflect not in this file'
      m=m+1
    endelse
    if clutter_id ne -1 then begin
      print,'Processing 2-D ','qc_ReflectivityClutterFlag'
      ncdf_varget, cdf_id, clutter_id, clutter
      clutter=transpose(clutter)
    endif else begin
      print, 'clutter not in this file'
    endelse
    if artifact_id ne -1 then begin
      print,'Processing 2-D ','qc_RadarArtifacts'
      ncdf_varget, cdf_id, artifact_id, artifact
      artifact=transpose(artifact)
    endif else begin
      print, 'artifact not in this file'
    endelse
;
;  If it is the first real variable get the attributes
;
    if p eq 0 and reflect_id ne -1 then begin
      print, 'first variable, getting the attributes'
      datainfo=ncdf_varinq(cdf_id, reflect_id)  ;inquire about the variable
      numatts=datainfo.natts  ;number of attributes assigned to the variable
      parts=strsplit(files[i],'/',/extract)
      part=strsplit(parts[5],'.',/extract)
      origfile=part[0]+'.'+part[1]
;
;  Store the attribute names
;
      dataatr=strarr(numatts)  ;string array to hold attribute names
      for k=0,numatts-1 do begin
;  Returns the name of the attribute
        dataatr[k]=ncdf_attname(cdf_id, reflect_id, k)  
      endfor
    endif 
;
;  Create the data array
;
    if p eq 0 and m lt 0 and reflect_id ne -1 then begin
      print,'first file and variable exists' 
      combined_reflect=reflect

      if clutter_id ne -1 then begin
        combined_clutter=clutter 
        combined_artifact=artifact
      endif else begin
        combined_clutter=dummy
        combined_artifact=dummy
      endelse

    endif else if p ge 0 and reflect_id ne -1 then begin
      print,'adding on existing variable'
      combined_reflect=[combined_reflect,reflect]

      if clutter_id ne -1 then begin      
        combined_clutter=[combined_clutter,clutter]
        combined_artifact=[combined_artifact,artifact]
      endif else begin
        combined_clutter=[combined_clutter,dummy]
        combined_artifact=[combined_artifact,dummy]
      endelse

    endif else if p lt 0 and m eq 0 and clutter_id eq -1 then begin
      print,'first file and variable doesnt exist'
      combined_reflect=dummy
      combined_clutter=dummy
      combined_artifact=dummy
    endif else if m ge 0 and clutter_id eq -1 then begin
      print,'adding on dummy variable'
      combined_reflect=[combined_reflect,dummy]
      combined_clutter=[combined_clutter,dummy]
      combined_artifact=[combined_artifact,dummy]
    endif
;
;  Close the netcdf file
;
    ncdf_close,cdf_id
  endif  ;end of not a gif file
endfor  ;end of looping though the matching files
print,'finished reading in the data'
;
;  Continue processing data because some data exists
;

if j gt -1 and p gt -1 then begin 

;
;  Manipulate the big data array
;

clutter=combined_clutter
artifact=combined_artifact
reflect=combined_reflect
combined_clutter=0  ;release memory
combined_artifact=0  ;release memory
combined_reflect=0
clutter=fix(clutter)
artifact=fix(artifact)

;
;  For the sgp site, change the clutter array to reflect 
;  the clutter-artifact conditions
;
print,site
if site eq 'sgp' then begin
  result=where(clutter eq 2 and artifact eq 1,count)
  if count ne 0 then clutter[result]=-7777
  result=where(clutter eq 3 and artifact eq 0,count)
  if count ne 0 then clutter[result]=-7777
  result=where(clutter eq 3 and artifact eq 6,count)
  if count ne 0 then clutter[result]=-7777
  result=where(clutter eq 0 and artifact eq 0,count)
  if count ne 0 then clutter[result]=-7777
  reflect=0  ;release memory
  artifact=0  ;release memory
endif else begin
  result=where(reflect ne -9999,count)
  if count ne 0 then reflect[result]=1
  result=where(reflect eq -9999,count)
  if count ne 0 then reflect[result]=-7777
  clutter=reflect
  reflect=0
  artifact=0
endelse
;
;  create data arrays to hold the new data
;

cloudy_profile=make_array(1000000,/float,value=0)  
vertcon=make_array(1000000,/float,value=-7777)
deltaz=make_array(1000000,/float,value=-7777)
cmin=make_array(1000000,/float,value=-7777)
cmax=make_array(1000000,/float,value=-7777)
crand=make_array(1000000,/float,value=-7777)

all_points=make_array(1000000,/float,value=0) 
cloudy_points=make_array(1000000,/float,value=0)
upper_hgt=make_array(1000000,/float,value=0)
lower_hgt=make_array(1000000,/float,value=0)
cloudy_upper=make_array(1000000,/float,value=0)
cloudy_lower=make_array(1000000,/float,value=0)

;year=make_array(1000000,/float,value=0)
;month=make_array(1000000,/float,value=0)
;day=make_array(1000000,/float,value=0)
;hour=make_array(1000000,/float,value=0)
;minute=make_array(1000000,/float,value=0)
;second=make_array(1000000,/float,value=0)
numsec=make_array(1000000,/long,value=0)

;*******************************
;  Loop through clutter array and put it in the m_height intervals
;  In other words, 350m intervals which start at 750m and end at 10500m
;*******************************
do_height='yes'
if do_height eq 'yes' then begin
mheight=750+(findgen(29)*360)
mclutter=make_array(n_elements(time_offset),n_elements(mheight),/float,value=-7777)
ind1=21
ind2=28
for h=0,n_elements(mheight)-2 do begin
  result_hgt=where(height ge mheight[h] and height lt mheight[h+1],count_hgt)
  ;print, mheight[h],height[result_hgt],mheight[h+1]
  if count_hgt gt 0 then begin
    clutter_hgt=clutter[*,result_hgt]
    ;print, clutter_hgt[ind1:ind2,*]
    for i=0,count_hgt-1 do begin
      one_clutter_hgt=clutter_hgt[*,i]
;  put in cloudy values
      result=where(one_clutter_hgt gt -7777,count)
      if count gt 0 then mclutter[result,h]=one_clutter_hgt[result]
;  put in missing values
      result=where(one_clutter_hgt eq -8888,count)
      if count gt 0 then mclutter[result,h]=one_clutter_hgt[result]
    endfor
    ;print,'mclutter'
    ;print, fix(mclutter[ind1:ind2,h])
  endif
endfor
clutter=mclutter[*,0:27]
height=mheight[0:27]
endif  ;do_height

;*******************************************
;  Find the column cloud fraction of data in timespacing intervals
;*******************************************


istp=1.  ;5time stop for printing
jstp=4.  ;7lower height stop for printing

k=-1.0   ;counts number of data points

;
;  Loop through profiles in time index
;

clutter_first_level=clutter[*,0]
for t=0,numtimes-2 do begin  ;loop through time intervals
  result_int=where(time_offset ge newtimes[t] and $
                   time_offset lt newtimes[t+1] and $
                   clutter_first_level ne -8888,count_int)
  if count_int gt 0 then begin
    clutter_int=clutter[result_int,*]
    ;openw,1,'file.txt'
    ;for p=0,n_elements(height)-1 do begin
    ;  printf,1,height[p],clutter_int[*,p]
    ;endfor
    ;close,1
;  Loop throught the heights in the profiles
    for j=0,n_elements(height)-2 do begin  ;lower layer
      for m=0,n_elements(height)-1 do begin;upper layer
        if height[m] gt height[j] then begin
          result=where(clutter_int[*,j] gt -7777,count_j)
          result=where(clutter_int[*,m] gt -7777,count_m)
          ;print,t,j,m,k,height[j],height[m]
          ;print,clutter_int[*,j],clutter_int[*,m]
          ;print,count_j,count_m,count_int 
          if (count_j gt 0 and count_j lt count_int) and $
             (count_m gt 0 and count_m lt count_int) then begin
            k=k+1
        ; determine cfrac for the two layers
            cfrac_j=float(count_j)/float(count_int)
            cfrac_m=float(count_m)/float(count_int)
        ; determine min overlap cloud cover   
            cmin[k]=min([1,cfrac_j+cfrac_m])
        ; determine max overlap cloud cover
            cmax[k]=max([cfrac_j,cfrac_m])
        ; determine random cloud cover
            crand[k]=cfrac_j+cfrac_m-(cfrac_j*cfrac_m)
        ; determine the layer separation
            deltaz[k]=height[m]-height[j]
        ; determine if they are vertically continuous
            vertcon[k]=1
            result_vert=where(height ge height[j] and $
                              height le height[m],count_vert)
            if count_vert gt 0 then begin
              layers=clutter_int[*,result_vert]
              q=-1.0
              while vertcon[k] eq 1 and q lt count_vert-1.0 do begin
                q=q+1.0 
                ;print,fix(height[result_vert[q]]),layers[*,q]
                result=where(layers[*,q] gt -7777,count_layers)
                if count_layers eq 0 then vertcon[k]=0
              endwhile
            endif else stop
        ; loop through to find colum cfrac between the layers
            both_layers=make_array(count_int,2,/float,value=-9999)
            both_layers[*,1]=clutter_int[*,m]
            both_layers[*,0]=clutter_int[*,j] 
            ctrue=0
            for r=0,count_int-1 do begin
              result=where(both_layers[r,*] gt -7777,count_bl)
              if count_bl gt 0 then ctrue=ctrue+1
            endfor
            cloudy_profile[k]=float(ctrue)/float(count_int)
            ;print,'vertcon',vertcon[k],'ctrue',cloudy_profile[k]
        ; other variables we want to keep track off
            all_points[k]=count_int  ;num data points 
            cloudy_points[k]=ctrue   ;num cloudy points
            upper_hgt[k]=height[m]   ;upper layer height
            lower_hgt[k]=height[j]   ;lower layer height
            cloudy_upper[k]=count_m  ;upper layer cloudy points
            cloudy_lower[k]=count_j  ;lower layer cloudy points
        ; time information
            ;standard_date_time,(startbst+newtimes[t]),yy,mm,dd,hh,mi,ss
            ;year[k]=yy
            ;month[k]=mm
            ;day[k]=dd
            ;hour[k]=hh
            ;minute[k]=mi
            ;second[k]=ss
            numsec[k]=newtimes[t]

          endif   ;if the layer isn't rejected
        endif   ;if the upper layer is above the lower
      endfor  ;loop through upper layer
    endfor  ;loop through lower layer
  endif   ;data is in the time interval
endfor  ;loop through time intervals

;
;  Find the column cloud fraction for the last time interval
;

result_int=where(time_offset ge newtimes[numtimes-1] and $
                 time_offset lt (newtimes[numtimes-1]+timespacing) and $
                 clutter_first_level ne -8888,count_int)
if count_int gt 0 then begin
  clutter_int=clutter[result_int,*]
  ;openw,1,'file.txt'
  ;for p=0,n_elements(height)-1 do begin
  ;  printf,1,height[p],clutter_int[*,p]
  ;endfor
  ;close,1
;  Loop throught the heights in the profiles
  for j=0,n_elements(height)-2 do begin  ;lower layer
    for m=0,n_elements(height)-1 do begin;upper layer
      if height[m] gt height[j] then begin
        result=where(clutter_int[*,j] gt -7777,count_j)
        result=where(clutter_int[*,m] gt -7777,count_m)
        ;print,t,j,m,k,height[j],height[m]
        ;print,clutter_int[*,j],clutter_int[*,m]
        ;print,count_j,count_m,count_int 
        if (count_j gt 0 and count_j lt count_int) and $
           (count_m gt 0 and count_m lt count_int) then begin
          k=k+1
       ; determine cfrac for the two layers
          cfrac_j=float(count_j)/float(count_int)
          cfrac_m=float(count_m)/float(count_int)
       ; determine min overlap cloud cover   
          cmin[k]=min([1,cfrac_j+cfrac_m])
       ; determine max overlap cloud cover
          cmax[k]=max(cfrac_j,cfrac_m)
       ; determine random cloud cover
          crand[k]=cfrac_j+cfrac_m-(cfrac_j*cfrac_m)
       ; determine the layer separation
          deltaz[k]=height[m]-height[j]
       ; determine if they are vertically continuous
          vertcon[k]=1
          result_vert=where(height ge height[j] and $
                            height le height[m],count_vert)
          if count_vert gt 0 then begin
            layers=clutter_int[*,result_vert]
            q=-1.0
            while vertcon[k] eq 1 and q lt count_vert-1.0 do begin
              q=q+1.0 
              ;print,fix(height[result_vert[q]]),layers[*,q]
              result=where(layers[*,q] gt -7777,count_layers)
              if count_layers eq 0 then vertcon[k]=0
            endwhile
          endif else stop
         ; loop through to find colum cfrac between the layers
          both_layers=make_array(count_int,2,/float,value=-9999)
          both_layers[*,1]=clutter_int[*,m]
          both_layers[*,0]=clutter_int[*,j] 
          ctrue=0
          for r=0,count_int-1 do begin
            result=where(both_layers[r,*] gt -7777,count_bl)
            if count_bl gt 0 then ctrue=ctrue+1
          endfor
          cloudy_profile[k]=float(ctrue)/float(count_int)
          ;print,'vertcon',vertcon[k],'ctrue',cloudy_profile[k]

       ; other variables we want to keep track off
          all_points[k]=count_int  ;num data points 
          cloudy_points[k]=ctrue   ;num cloudy points
          upper_hgt[k]=height[m]   ;upper layer height
          lower_hgt[k]=height[j]   ;lower layer height
          cloudy_upper[k]=count_m  ;upper layer cloudy points
          cloudy_lower[k]=count_j  ;lower layer cloudy points
       ; time information
          ;standard_date_time,(startbst+newtimes[t]),yy,mm,dd,hh,mi,ss
          ;year[k]=yy
          ;month[k]=mm
          ;day[k]=dd
          ;hour[k]=hh
          ;minute[k]=mi
          ;second[k]=ss
          numsec[k]=newtimes[t]

        endif   ;if the layer isn't rejected
      endif   ;if the upper layer is above the lower
    endfor  ;loop through upper layer
  endfor  ;loop through lower layer
endif   ;data is in the time interval

;
;  Cut the arrays down to size
;
cloudy_profile=cloudy_profile[0:k]
deltaz=deltaz[0:k]
vertcon=vertcon[0:k]
cmin=cmin[0:k]
cmax=cmax[0:k]
crand=crand[0:k]

all_points=all_points[0:k]
cloudy_points=cloudy_points[0:k]
upper_hgt=upper_hgt[0:k]
lower_hgt=lower_hgt[0:k]
cloudy_upper=cloudy_upper[0:k]
cloudy_lower=cloudy_lower[0:k]

numsec=numsec[0:k]

;
;  Convert to percent
;

;result=where(cloudy_profile ne -9999,count)
;if count gt 0 then cloudy_profile[result]=cloudy_profile[result]*100.

;**************************************
;  Write to the output netcdf file
;**************************************

;
;  Create the netcdf file
;
cdf_id=ncdf_create(avgfilename,/clobber)
;
;  Define the dimensions
;
num_id=ncdf_dimdef(cdf_id,'number',/unlimited)
height_grid_id=ncdf_dimdef(cdf_id,'height_grid',$
  n_elements(height))
;
;  Define the variables
;
base_time_id=ncdf_vardef(cdf_id,'base_time',/long)
ncdf_attput,cdf_id,base_time_id,'long_name','base time in epoch'
ncdf_attput,cdf_id,base_time_id,'units','seconds since 1970-01-01 00:00:00 00:00'

height_id=ncdf_vardef(cdf_id,'height',height_grid_id,/float)
ncdf_attput,cdf_id,height_id,'long_name','AGL'
ncdf_attput,cdf_id,height_id,'units','meters'

time_offset_id=ncdf_vardef(cdf_id,'time_offset',num_id,/long)
ncdf_attput,cdf_id,time_offset_id,'long_name','time offset from base_time'
ncdf_attput,cdf_id,time_offset_id,'units','seconds'

deltaz_id=ncdf_vardef(cdf_id,'deltaz',num_id,/float)
ncdf_attput,cdf_id,deltaz_id,'long_name','layer separation'
ncdf_attput,cdf_id,deltaz_id,'units','meters'

vertcon_id=ncdf_vardef(cdf_id,'vertcon',num_id,/float)
ncdf_attput,cdf_id,vertcon_id,'long_name','vertically continuous layers'
ncdf_attput,cdf_id,vertcon_id,'units','1=yes  0=no'

ctrue_id=ncdf_vardef(cdf_id,'ctrue',num_id,/float)
ncdf_attput,cdf_id,ctrue_id,'long_name','observed overlap'
ncdf_attput,cdf_id,ctrue_id,'units','%'

cmin_id=ncdf_vardef(cdf_id,'cmin',num_id,/float)
ncdf_attput,cdf_id,cmin_id,'long_name','minimum overlap'
ncdf_attput,cdf_id,cmin_id,'units','%'

cmax_id=ncdf_vardef(cdf_id,'cmax',num_id,/float)
ncdf_attput,cdf_id,cmax_id,'long_name','maximum overlap'
ncdf_attput,cdf_id,cmax_id,'units','%'

crand_id=ncdf_vardef(cdf_id,'crand',num_id,/float)
ncdf_attput,cdf_id,crand_id,'long_name','random overlap'
ncdf_attput,cdf_id,crand_id,'units','%'

all_points_id=ncdf_vardef(cdf_id,'all_points',num_id,/float)
ncdf_attput,cdf_id,all_points_id,'long_name','data points in time interval'
ncdf_attput,cdf_id,all_points_id,'units','number of points'

cloudy_points_id=ncdf_vardef(cdf_id,'cloudy_points',num_id,/float)
ncdf_attput,cdf_id,cloudy_points_id,'long_name','cloudy data points in time interval'
ncdf_attput,cdf_id,cloudy_points_id,'units','number of points'

cloudy_upper_id=ncdf_vardef(cdf_id,'cloudy_upper',num_id,/float)
ncdf_attput,cdf_id,cloudy_upper_id,'long_name','cloudy data points in time interval'
ncdf_attput,cdf_id,cloudy_upper_id,'long_name','for the upper layer'
ncdf_attput,cdf_id,cloudy_upper_id,'units','number of points'

cloudy_lower_id=ncdf_vardef(cdf_id,'cloudy_lower',num_id,/float)
ncdf_attput,cdf_id,cloudy_lower_id,'long_name','cloudy data points in time interval'
ncdf_attput,cdf_id,cloudy_lower_id,'long_name','for the lower layer'
ncdf_attput,cdf_id,cloudy_lower_id,'units','number of points'

upper_hgt_id=ncdf_vardef(cdf_id,'upper_hgt',num_id,/float)
ncdf_attput,cdf_id,upper_hgt_id,'long_name','upper layer height'
ncdf_attput,cdf_id,upper_hgt_id,'units','meters'

lower_hgt_id=ncdf_vardef(cdf_id,'lower_hgt',num_id,/float)
ncdf_attput,cdf_id,lower_hgt_id,'long_name','lower layer height'
ncdf_attput,cdf_id,lower_hgt_id,'units','meters'

;year_id=ncdf_vardef(cdf_id,'year',num_id,/float)
;month_id=ncdf_vardef(cdf_id,'month',num_id,/float)
;day_id=ncdf_vardef(cdf_id,'day',num_id,/float)
;hour_id=ncdf_vardef(cdf_id,'hour',num_id,/float)
;minute_id=ncdf_vardef(cdf_id,'minute',num_id,/float)
;second_id=ncdf_vardef(cdf_id,'second',num_id,/float)

;
;  Change from define to data mode
;
ncdf_control,cdf_id,/endef
; 
;  Put the variables in the netcdf file
;
ncdf_varput,cdf_id,base_time_id,startbst
ncdf_varput,cdf_id,height_id,height
ncdf_varput,cdf_id,ctrue_id,cloudy_profile
ncdf_varput,cdf_id,deltaz_id,deltaz
ncdf_varput,cdf_id,vertcon_id,vertcon
ncdf_varput,cdf_id,cmin_id,cmin
ncdf_varput,cdf_id,cmax_id,cmax
ncdf_varput,cdf_id,crand_id,crand

ncdf_varput,cdf_id,all_points_id,all_points
ncdf_varput,cdf_id,cloudy_points_id,cloudy_points
ncdf_varput,cdf_id,upper_hgt_id,upper_hgt
ncdf_varput,cdf_id,lower_hgt_id,lower_hgt
ncdf_varput,cdf_id,cloudy_upper_id,cloudy_upper
ncdf_varput,cdf_id,cloudy_lower_id,cloudy_lower

;ncdf_varput,cdf_id,year_id,year
;ncdf_varput,cdf_id,month_id,month
;ncdf_varput,cdf_id,day_id,day
;ncdf_varput,cdf_id,hour_id,hour
;ncdf_varput,cdf_id,minute_id,minute
;ncdf_varput,cdf_id,second_id,second
ncdf_varput,cdf_id,time_offset_id,numsec
;
;  Close the netcdf file
;
ncdf_close,cdf_id

endif  ;the data exits in the file
end