;
;  Read in the merged moments data
;
pro get_mergedmoments_list,varname,data,site,date,do_clutter,num_files

;***************************************************************
;  Inputs - set here or passed in through the command line
;***************************************************************

;  Site to process
;site='sgp'

;  Date to process
;date=['20050101','20050701']
;date='199701'

;  Variables to get
;varname=['base_time','time_offset','Heights',$
;         'ReflectivityBestEstimateMPLMask']
	 
;  Set to yes if you want to mask the data
;do_clutter='yes'

;********************************************************************
;  If the data is to be masked, then these variables have to be read
;********************************************************************

if do_clutter eq 'yes' then begin
  result=where(varname eq 'qc_ReflectivityClutterFlag',count)
  if count eq 0 then varname=[varname,'qc_ReflectivityClutterFlag']	 
  result=where(varname eq 'qc_RadarArtifacts',count)
  if count eq 0 then varname=[varname,'qc_RadarArtifacts']	 
endif

;**************************************************************
;  Establish some constants
;**************************************************************

;  Missing data flag
missing=-8888

;  Mpl mask
mplmask=-9000

;  Print flag
dp='no'

;**************************************************************
;  Set up strings to build the filename
;**************************************************************

;  Main arm site directory
arm_dir='/uufs/chpc.utah.edu/common/home/mace_data3/data/arm/'+site+'/'


;  Data directory and filename depend on the site
;  True height array - Height array we have chosen for standard.
;  This is because some of 1997 has a different height array. (MSL)
if site eq 'sgp' then begin
  ;  Date when the mmcr data stream changes (20030909)
  limit_date=julday(9,9,2003,0,0,0)
  ;  Strings for the first set of files
  datastream_1='sgpmmcrcalC1.a1/'
  arm_file_1='sgpmmcrMergedMomentsC1.c1'
  ;  Strings for the last set of files
  datastream_2='sgpmmcrmomC1.b1/'
  arm_file_2='sgpmmcrMergedMomentsC1.c1'
  ;  Calculate true height array
  true_height=findgen(167)*90.0+423.0
  true_numheights=167
endif else if site eq 'nsa' then begin
  ;  Date when the mmcr data stream changes (20040415)
  limit_date=julday(4,15,2004,0,0,0)
  ;  Strings for the first set of files
  datastream_1='nsammcrcalC1.a1/'
  arm_file_1='nsammcrMergedMomentsC1.c1'
  ;  Strings for the first set of files
  datastream_2='nsammcrmomC1.b1/'
  arm_file_2='nsammcrMergedMomentsC1.c1'
  ;  Calculate true height array
  true_height=findgen(149)*90.0+112.0
  true_numheights=149
endif else if site eq 'twpc1' then begin
  ;  Date when the mmcr data stream changes (20060616)
  limit_date=julday(6,16,2006,0,0,0)
  ;  Strings for the first set of files
  datastream_1='twpmmcrcalC1.a1/'
  arm_file_1='twpmmcrMergedMomentsC1.c1'
  ;  Strings for the second set of files
  datastream_2='twpmmcrmomC1.b1/'
  arm_file_2='twpmmcrMergedMomentsC1.c1'
  ;  Calculate true height array
  true_height=findgen(226)*90.0+115.0
  true_numheights=226
endif else if site eq 'twpc2' then begin
  ;  Date when the mmcr data stream changes (20061127)
  limit_date=julday(11,27,2006,0,0,0)
  ;  Strings for the first set of files
  datastream_1='twpmmcrcalC2.a1/'
  arm_file_1='twpmmcrMergedMomentsC2.c1'
  ;  Strings for the second set of files
  datastream_2='twpmmcrmomC2.b1/'
  arm_file_2='twpmmcrMergedMomentsC2.c1'
  ;  Calculate true height array
  true_height=findgen(226)*90.0+111.0
  true_numheights=226
endif else if site eq 'twpc3' then begin
  ;  Date when the mmcr data stream changes (20051104)
  limit_date=julday(11,4,2005,0,0,0)
  ;  Strings for the first set of files
  datastream_1='twpmmcrcalC3.a1/'
  arm_file_1='twpmmcrMergedMomentsC3.c1'
  ;  Strings for the second set of files
  datastream_2='twpmmcrmomC3.b1/'
  arm_file_2='twpmmcrMergedMomentsC3.c1'
  ;  Calculate true height array
  true_height=findgen(234)*87.4138+138.6936
  true_numheights=234
endif

;  Height flag = 1 when there is data of different heights (1997 at sgp)
height_flag=0

;  List of 2d variables in the file and their dimensions
var2=[$
'ReflectivityBestEstimate',$
'ReflectivityBestEstimateMPLMask',$
'MeanDopplerVelocity',$
'SpectralWidth',$
'SignaltoNoiseRatio',$
'ModeID',$
'qc_ReflectivityClutterFlag',$
'qc_RadarArtifacts']
vardim2=make_array(n_elements(var2),/int,value=2)

;  List of 1d variables in the file and their dimensions
var1=[$
'number_of_layers',$
'radar_layer_status',$
'CloudBaseBestEstimate',$
'RadarFirstTop',$
'base_height_layer_1',$
'base_temp_layer_1',$
'top_height_layer_1',$
'top_temp_layer_1',$
'data_indicator_flag_layer_1',$
'base_height_layer_2',$
'base_temp_layer_2',$
'top_height_layer_2',$
'top_temp_layer_2',$
'data_indicator_flag_layer_2',$
'base_height_layer_3',$
'base_temp_layer_3',$
'top_height_layer_3',$
'top_temp_layer_3',$
'data_indicator_flag_layer_3']
vardim1=make_array(n_elements(var1),/int,value=1)

;  Variables that aren't to be combined across files
var0=['Heights','base_time','time_offset']
vardim0=make_array(n_elements(var0),/int,value=0)

;  Combine these all into one big array
var=[var2,var1,var0]
vardim=[vardim2,vardim1,vardim0]

;***************************************************************
;  Make sure the dates are in an array
;***************************************************************

;  Number of different dates to pick up
numdates=n_elements(date)

;  Put the search dates in an array of dates called curdate
curdate=make_array(/string,numdates)
if numdates gt 1 then begin
  for i=0,numdates-1 do begin
    curdate[i]=strcompress(string(date[i]),/remove_all)
  endfor
endif else begin
  curdate[0]=strcompress(string(date),/remove_all)
endelse

;***************************************************************
;  If any of the dates have just the month and year, expand
;  to have a date for each day in the month
;***************************************************************

for i=0,numdates-1 do begin
  ;  This is a full date including day
  if strlen(date[i]) eq 8 then begin
    if n_elements(alldates) eq 0 then begin
      alldates=curdate[i]
    endif else begin
      alldates=[alldates,curdate[i]]
    endelse
  ;  Expand the month
  endif else if strlen(date) lt 8 then begin
    year=strmid(date,0,4)
    month=strmid(date,4,2)
    find_number_of_days,year,month,numdays    
    newdates=timegen(numdays,units="Days",start=julday(1,month,year,0,0,0))
    caldat,newdates,mm,dd,yy,hh,mi,ss  
    newdatestr=string(yy,format='(I4)')+string(mm,format='(I02)')+$
               string(dd,format='(I02)')
    if n_elements(alldates) eq 0 then begin
      alldates=newdatestr
    endif else begin
      alldates=[alldates,newdatestr]
    endelse 
  endif
endfor  ;end of loop through dates

;  Get the number of dates again
curdate=alldates
numdates=n_elements(curdate)

;***************************************************************
;  Loop through the dates and make an array of file names
;***************************************************************

;  Loop through the dates
for i=0,numdates-1 do begin

  ;  Date to be analyzed
  year=strmid(curdate[i],0,4)
  month=strmid(curdate[i],4,2)
  day=strmid(curdate[i],6,2)

  ;  Calculate check date
  check_date=julday(month,day,year,0,0,0)

  ;  Choose the arm file search string based on the date
  if check_date lt limit_date then begin
    arm_file=arm_dir+arm_file_1+'/'+year+'/'+arm_file_1+'*'+curdate[i]+'*cdf*'
  endif else if check_date ge limit_date then begin
    arm_file=arm_dir+arm_file_1+'/'+year+'/'+arm_file_1+'*'+curdate[i]+'*cdf'
  endif

  ;  Unzip if necessary
  unix_command='gunzip -fq '+arm_file+'*gz'
  spawn,unix_command

  ;  Get a list of the files for the day
  filesd=file_search(arm_file,count=num_files)  

  ;  Put these files in one big list if files exist for this date
  if num_files gt 0 then begin
    if n_elements(files) eq 0 then begin
      files=filesd
    endif else begin
      files=[files,filesd]
    endelse
  endif
endfor  ;end of loop through the dates

;  Number of files to read and cat together
num_files=n_elements(files)

;***************************************************************
;  Only continue processing if files were found
;***************************************************************

if num_files gt 0 then begin

;***************************************************************
;  Loop through all the matching files and read in the data
;***************************************************************

  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,files[i]

;  Open the netcdf file
      cdf_id=ncdf_open(files[i])

;  Get the dimension id's
      time_did=ncdf_dimid(cdf_id,'time')
      height_did=ncdf_dimid(cdf_id,'nheights')

;  Get the dimensions
      ncdf_diminq,cdf_id,time_did,charstring,numtimes
      ncdf_diminq,cdf_id,height_did,charstring,numheights

;  Create the dummy variables
      dummy1d=make_array(numtimes,/float,value=missing)
      dummy2d=make_array(numheights,numtimes,/float,value=missing)

;  Get the variable ids
      for v=0,n_elements(varname)-1 do begin
        xstr=varname[v]+'_id=ncdf_varid(cdf_id,"'+varname[v]+'")'
        result=execute(xstr)
        if dp eq 'yes' then print,varname[v]
      endfor
;
;  Loop through the variable names.  Get the variable if it exists.
;  If it doesn't exist, set the variable equal to the dummy array.
;

      for v=0,n_elements(varname)-1 do begin

;  Determine the dimension of this variable
        xstr='idx=where(var eq "'+varname[v]+'",count)'
        result=execute(xstr)
        if count eq 0 then begin 
          print,varname[v],'variable not in list of acceptable variables'
          stop
        endif else begin
          dim=vardim[idx]
          if dim eq 1 then dumname='dummy1d' $
          else if dim eq 2 then dumname='dummy2d' $
          else if dim eq 3 then dumname='dummy3d'
        endelse
        if dp eq 'yes' then print,'processing ',varname[v],'dimension',dim

;  Get the variable
        if dim ne 0 then begin
          xstr='if '+varname[v]+'_id gt -1 then ncdf_varget,cdf_id,'+$
            varname[v]+'_id,'+varname[v]+' else '+varname[v]+'='+dumname
        endif else begin
          xstr='if '+varname[v]+'_id gt -1 then ncdf_varget,cdf_id,'+$
            varname[v]+'_id,'+varname[v]    
        endelse
        result=execute(xstr)

;  Transpose 2d variables
        if dim gt 1 then begin
          if dp eq 'yes' then print,'transposing'
          xstr=varname[v]+'=transpose('+varname[v]+')'
          result=execute(xstr)
        endif
      endfor

;  Check the height arrays.  If they are different from the 
;  standard then adjust accordingly.
      if (heights[0] ne true_height[0]) or $
         (numheights ne true_numheights) then begin
	 
        ;  This data has different height data in it
        height_flag=1

        for v=0,n_elements(varname)-1 do begin

          ;  Determine the dimension of this variable
          xstr='idx=where(var eq "'+varname[v]+'",count)'
          result=execute(xstr)
          dim=vardim[idx]
	  if dim gt 1 then begin

            xstr='r_arr=make_array(numtimes,true_numheights,/float,value=-9999)'
            result=execute(xstr)
	    
	    for h=0,true_numheights-1 do begin
	      heightspaces=abs(heights-true_height[h])
	      closest=min(heightspaces,closest_sub)
              if closest lt 100 then begin
	        xstr='r_arr[*,h]='+varname[v]+'[*,closest_sub]'
                result=execute(xstr)
              endif 
	    endfor
            xstr=varname[v]+'=r_arr'
	    result=execute(xstr)
	  endif  ;2dvar
        endfor  ;end of loop through varnames
      endif  ;end of heights array different

;
;  Combine the data 
;

      for v=0,n_elements(varname)-1 do begin

        if varname[v] eq 'time_offset' then begin
          if n_elements(com_btto) eq 0 then begin
            com_btto=long(base_time)+long(time_offset)
          endif else begin
            com_btto=[long(com_btto),long(base_time)+long(time_offset)]
          endelse
        endif
    
        if varname[v] ne 'base_time' and varname[v] ne 'time_offset' and $
           varname[v] ne 'height' then begin
          xstr='if n_elements(com_'+varname[v]+') eq 0 then com_'+$
                varname[v]+'='+varname[v]+' else com_'+varname[v]+$
                '=[com_'+varname[v]+','+varname[v]+']'
          result=execute(xstr)
        endif

      endfor

;  Close the netcdf file
      ncdf_close,cdf_id

    endif  ;only process cdf files
  endfor  ;loop through number of matching files

;***************************************************************
;  Remove the 'com_' prefix from the data array names
;***************************************************************

for v=0,n_elements(varname)-1 do begin

  if varname[v] eq 'base_time' then begin
    base_time=long(com_btto[0])  
  endif else if varname[v] eq 'time_offset' then begin
    time_offset=long(com_btto)-long(base_time)        
  endif else if varname[v] ne 'Heights' then begin
    xstr=varname[v]+'=com_'+varname[v]
    result=execute(xstr)
  endif

endfor  ;end of loop through variables

;***************************************************************
;  Fix the bug where the last reflectivity at the 
;  end of the column is -32767
;***************************************************************

result=where(varname eq 'ReflectivityBestEstimate',count)
if count gt 0 then begin

  result=where(ReflectivityBestEstimate eq -32767,count)
  if count gt 0 then ReflectivityBestEstimate[result]=-9999

endif

;***************************************************************
;  Use clutter and artifact flag to restrict ReflectivityBestEstimate
;***************************************************************

result=where(varname eq 'ReflectivityBestEstimate',count)
if count gt 0 and do_clutter eq 'yes' then begin

  result=where(qc_ReflectivityClutterFlag ne missing,count_clutter)
  result=where(qc_RadarArtifacts ne missing,count_artifact)

  if count_clutter gt 0 and count_artifact gt 0 then begin

    result=where(qc_ReflectivityClutterFlag eq 2 and $
                 qc_RadarArtifacts eq 1,count)
    if count ne 0 then ReflectivityBestEstimate[result]=-9999

    result=where(qc_ReflectivityClutterFlag eq 3 and $
                 qc_RadarArtifacts eq 0,count)
    if count ne 0 then ReflectivityBestEstimate[result]=-9999

    ;result=where(qc_ReflectivityClutterFlag eq 7 and $
    ;             qc_RadarArtifacts eq 7,count);; now good
    ;if count ne 0 then ReflectivityBestEstimate[result]=-9999

    result=where(qc_ReflectivityClutterFlag eq 3 and $
                 qc_RadarArtifacts eq 6 and $
		 ReflectivityBestEstimate ne missing,count)
    if count ne 0 then ReflectivityBestEstimate[result]=-9999

    result=where(qc_ReflectivityClutterFlag eq 0 and $
                 qc_RadarArtifacts eq 0 and $
		 ReflectivityBestEstimate ne missing,count)
    if count ne 0 then ReflectivityBestEstimate[result]=-9999

  endif
endif  ;end of do_clutter

;***************************************************************
;  Use clutter and artifact flag to restrict ReflectivityBestEstimateMPLMask
;***************************************************************

result=where(varname eq 'ReflectivityBestEstimateMPLMask',count)
if count gt 0 and do_clutter eq 'yes' then begin

  result=where(qc_ReflectivityClutterFlag ne missing,count_clutter)
  result=where(qc_RadarArtifacts ne missing,count_artifact)

  if count_clutter gt 0 and count_artifact gt 0 then begin

    result=where(qc_ReflectivityClutterFlag eq 2 and $
                 qc_RadarArtifacts eq 1,count)
    if count ne 0 then ReflectivityBestEstimateMPLMask[result]=-9999

    result=where(qc_ReflectivityClutterFlag eq 3 and $
                 qc_RadarArtifacts eq 0,count)
    if count ne 0 then ReflectivityBestEstimateMPLMask[result]=-9999

    ;result=where(qc_ReflectivityClutterFlag eq 7 and $
    ;             qc_RadarArtifacts eq 7,count);; now good
    ;if count ne 0 then ReflectivityBestEstimateMPLMask[result]=-9999

    result=where(qc_ReflectivityClutterFlag eq 3 and $
                 qc_RadarArtifacts eq 6 and $
		 ReflectivityBestEstimateMPLMask ne missing and $
                 ReflectivityBestEstimateMPLMask ne mplmask,count)
    if count ne 0 then ReflectivityBestEstimateMPLMask[result]=-9999

    result=where(qc_ReflectivityClutterFlag eq 0 and $
                 qc_RadarArtifacts eq 0 and $
		 ReflectivityBestEstimateMPLMask ne missing and $
                 ReflectivityBestEstimateMPLMask ne mplmask,count)
    if count ne 0 then ReflectivityBestEstimateMPLMask[result]=-9999

  endif
endif  ;end of do_clutter

;***************************************************************
;  Zip up the data files
;***************************************************************

;for i=0,num_files-1 do begin
;  unix_command='gzip -fq  '+files[i]  ;zip the data files
;  spawn,unix_command
;endfor

;*******************************************
;  If there are different height arrays, give height the 
;  true height values
;*******************************************

if height_flag eq 1 then heights=true_height

;*****************************************
;  Put the data in a structure to pass back out
;*****************************************

data={data,filename:'string'}
for v=0,n_elements(varname)-1 do begin
  xstr='data=create_struct(data,"'+varname[v]+'",'+varname[v]+')'
  result=execute(xstr)
endfor  

;  Add in a variable to hold the files
xstr='data=create_struct(data,"files",files)'
result=execute(xstr)


endif  ;end of found files

return
end