;
;  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
;

pro ctrue_heights

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')
;
;  Get the variable if it exists
;
    if clutter_id ne -1 then begin
      print,'Processing 2-D ','qc_ReflectivityClutterFlag'
      ncdf_varget, cdf_id, clutter_id, clutter
      clutter=transpose(clutter)
      p=p+1
    endif else begin
      print, 'Variable not in this file'
      m=m+1
    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, 'Variable not in this file'
    endelse

;
;  If it is the first real variable get the attributes
;
    if p eq 0 and clutter_id ne -1 then begin
      print, 'first variable, getting the attributes'
      datainfo=ncdf_varinq(cdf_id, clutter_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, clutter_id, k)  
      endfor
    endif 
;
;  Create the data array
;
    if p eq 0 and m lt 0 and clutter_id ne -1 then begin
      print,'first file and variable exists' 
      combined_clutter=clutter
      combined_artifact=artifact
    endif else if p ge 0 and clutter_id ne -1 then begin
      print,'adding on existing variable'
      combined_clutter=[combined_clutter,clutter]
      combined_artifact=[combined_artifact,artifact]
    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_clutter=dummy
      combined_artifact=dummy
    endif else if m ge 0 and clutter_id eq -1 then begin
      print,'adding on dummy variable'
      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
combined_clutter=0  ;release memory
combined_artifact=0  ;release memory
clutter=fix(clutter)
artifact=fix(artifact)

;
;  change the clutter array to reflect the clutter-artifact conditions
;

result=where(clutter eq 2 and artifact eq 1,count)
if count ne 0 then clutter[result]=-9999
result=where(clutter eq 3 and artifact eq 0,count)
if count ne 0 then clutter[result]=-9999
result=where(clutter eq 3 and artifact eq 6,count)
if count ne 0 then clutter[result]=-9999
result=where(clutter eq 0 and artifact eq 0,count)
if count ne 0 then clutter[result]=-9999

;
;  create data arrays to hold the new data
;

cloudy_profile=make_array(1000000,/float,value=0)  
vertcon=make_array(1000000,/float,value=-9999)
deltaz=make_array(1000000,/float,value=-9999)


;
;  Loop through clutter array and put it in the m_height intervals
;
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=-9999)
;ind1=9
;ind2=17
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]
      result=where(one_clutter_hgt gt -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
jstp=4.  ;7lower height

k=-1.0

;
;  Loop through profiles in time index
;

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],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 -8888,count_j)
          result=where(clutter_int[*,m] gt -8888,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 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 -8888,count_layers)
                if count_layers eq 0 then vertcon[k]=0
              endwhile
            endif
           ; 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 -8888,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]
           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),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 -8888,count_j)
        result=where(clutter_int[*,m] gt -8888,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 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 -8888,count_layers)
              if count_layers eq 0 then vertcon[k]=0
            endwhile
          endif
         ; 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 -8888,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]
        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]

;
;  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 file
;**************************************

;
;  Write this to a netcdf file
;
cdfname='robin_hogan_ctrue_'+strcompress(curdate,/remove_all)+'.cdf'
cdf_id=ncdf_create(cdfname,/clobber)

num_id=ncdf_dimdef(cdf_id,'number',/unlimited)

ctrue_id=ncdf_vardef(cdf_id,'ctrue',num_id,/float)
deltaz_id=ncdf_vardef(cdf_id,'deltaz',num_id,/float)
vertcon_id=ncdf_vardef(cdf_id,'vertcon',num_id,/float)

ncdf_control,cdf_id,/endef

ncdf_varput,cdf_id,ctrue_id,cloudy_profile
ncdf_varput,cdf_id,deltaz_id,deltaz
ncdf_varput,cdf_id,vertcon_id,vertcon


ncdf_close,cdf_id

endif  ;the data exits in the file
end