;
;  Read in the average data
;

pro get_average_data_list,varname,data,$
  site,date,avg_int,ret_type,m_or_d,fill_flag,num_files

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

;  List of variables to read in
;varname=['base_time','time_offset','height',$
;         'IWC_best_estimate',$
;         're_ice_best_estimate',$
;         'ice_micro_error_estimate',$
;         'ice_tau_solar',$
;         'ice_omega_solar',$
;         'ice_g_solar',$
;         'ice_tau_IR',$
;         'ice_omega_IR',$
;         'ice_g_IR',$
;         'lwc_best_estimate',$
;         're_liquid_best_estimate',$
;         'liquid_micro_error_estimate',$
;         'liquid_tau_solar',$
;         'liquid_omega_solar',$
;         'liquid_g_solar',$
;         'liquid_tau_IR',$
;         'liquid_omega_IR',$
;         'liquid_g_IR',$
;         'ice_fraction_CCM3',$
;           ]

;site='sgp'
;date='200112'
;avg_int='5min'
;ret_type='ciret4' 
;m_or_d='day'  ;'month' or 'day'
;fill_flag='yes' ;'yes,no' ;need to fill for 2d plots - don't always need data filled

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

;  Printing flag.  Set dp='yes' if you want to see print statements
dp='no'

;  Missing data value
missing=-8888

;  List of 2d variables (cld_type) and their dimensions
var4=['Optical_Depth_LinearSite']
vardim4=make_array(n_elements(var4),/int,value=4)

;  List of 2d variables (scn_type) and their dimensions
var3=['Broadband_LW_FluxSite','Broadband_SW_AlbedoSite']
vardim3=make_array(n_elements(var3),/int,value=3)

;  List of 2d variables in the file and their dimensions
var2=[$
'downwelling_solar_diffuse',$
'downwelling_solar_diffuse_clear',$
'downwelling_solar_direct',$
'downwelling_solar_direct_clear',$
'ice_effective_radius_source',$
'ice_fraction_CCM3',$
'ice_g_IR',$
'ice_g_solar',$
'ice_micro_error_estimate',$
'ice_omega_IR',$
'ice_omega_solar',$
'ice_RadiationParameterization_source_ir',$
'ice_RadiationParameterization_source_solar',$
'ice_tau_IR',$
'ice_tau_solar',$
'IR_flux_down',$
'IR_flux_down_clear',$
'IR_flux_up',$
'IR_flux_up_clear',$
'IWC_best_estimate',$
'iwc_source',$
'liquid_effective_radius_source',$
'liquid_g_IR',$
'liquid_g_solar',$
'liquid_micro_error_estimate',$
'liquid_omega_IR',$
'liquid_omega_solar',$
'liquid_RadiationParameterization_source_ir',$
'liquid_RadiationParameterization_source_solar',$
'liquid_tau_IR',$
'liquid_tau_solar',$
'lwc_best_estimate',$
'lwc_source',$
'MeanDopplerVelocity',$
'mxrat',$
'mxrat_adj',$
'pressure',$
'qc_ReflectivityClutterFlagAllPoints',$
'qc_ReflectivityClutterFlagCloudyPoints',$
'qc_ReflectivityClutterFlagMode',$
'ReflectivityBestEstimate',$
'ReflectivityCloudyPoints',$
'ReflectivityCloudyPoints_mpl',$
'ReflectivityPossiblePoints',$
're_ice_best_estimate',$
're_liquid_best_estimate',$
'solar_flux_down',$
'solar_flux_down_clear',$
'solar_flux_up',$
'solar_flux_up_clear',$
'SpectralWidth',$
'temp',$
'VVEL']     

var2=[var2,$
  'downwelling_solar_diffuse_sdv',$
  'downwelling_solar_diffuse_total',$
  'downwelling_solar_diffuse_valid',$
  'downwelling_solar_diffuse_clear_sdv',$
  'downwelling_solar_diffuse_clear_total',$
  'downwelling_solar_diffuse_clear_valid',$
  'downwelling_solar_direct_sdv',$
  'downwelling_solar_direct_total',$
  'downwelling_solar_direct_valid',$
  'downwelling_solar_direct_clear_sdv',$
  'downwelling_solar_direct_clear_total',$
  'downwelling_solar_direct_clear_valid',$
  'ice_fraction_CCM3_sdv','ice_fraction_CCM3_total','ice_fraction_CCM3_valid',$
  'ice_g_IR_sdv','ice_g_IR_total','ice_g_IR_valid',$
  'ice_g_solar_sdv','ice_g_solar_total','ice_g_solar_valid',$
  'ice_micro_error_estimate_sdv',$
  'ice_micro_error_estimate_total',$
  'ice_micro_error_estimate_valid',$
  'ice_omega_IR_sdv','ice_omega_IR_total','ice_omega_IR_valid',$
  'ice_omega_solar_sdv','ice_omega_solar_total','ice_omega_solar_valid',$
  'ice_tau_IR_sdv','ice_tau_IR_total','ice_tau_IR_valid',$
  'ice_tau_solar_sdv','ice_tau_solar_total','ice_tau_solar_valid',$
  'IR_flux_down_sdv','IR_flux_down_total','IR_flux_down_valid',$
  'IR_flux_down_clear_sdv','IR_flux_down_clear_total','IR_flux_down_clear_valid',$
  'IR_flux_up_sdv','IR_flux_up_total','IR_flux_up_valid',$
  'IR_flux_up_clear_sdv','IR_flux_up_clear_total','IR_flux_up_clear_valid',$
  'IWC_best_estimate_sdv','IWC_best_estimate_total','IWC_best_estimate_valid',$
  'liquid_g_IR_sdv','liquid_g_IR_total','liquid_g_IR_valid',$
  'liquid_g_solar_sdv','liquid_g_solar_total','liquid_g_solar_valid',$
  'liquid_micro_error_estimate_sdv',$
  'liquid_micro_error_estimate_total',$
  'liquid_micro_error_estimate_valid',$
  'liquid_omega_IR_sdv','liquid_omega_IR_total','liquid_omega_IR_valid',$
  'liquid_omega_solar_sdv','liquid_omega_solar_total','liquid_omega_solar_valid',$
  'liquid_tau_IR_sdv','liquid_tau_IR_total','liquid_tau_IR_valid',$
  'liquid_tau_solar_sdv','liquid_tau_solar_total','liquid_tau_solar_valid',$
  'lwc_best_estimate_sdv','lwc_best_estimate_total','lwc_best_estimate_valid',$
  're_ice_best_estimate_sdv','re_ice_best_estimate_total','re_ice_best_estimate_valid',$
  're_liquid_best_estimate_sdv',$
  're_liquid_best_estimate_total',$
  're_liquid_best_estimate_valid',$
  'solar_flux_down_sdv','solar_flux_down_total','solar_flux_down_valid',$
  'solar_flux_down_clear_sdv','solar_flux_down_clear_total','solar_flux_down_clear_valid',$
  'solar_flux_up_sdv','solar_flux_up_total','solar_flux_up_valid',$
  'solar_flux_up_clear_sdv','solar_flux_up_clear_total','solar_flux_up_clear_valid']
vardim2=make_array(n_elements(var2),/int,value=2)

;  List of 1d variables in the file and their dimensions
var1=[$
'base_height_layer_1',$
'base_temp_layer_1',$
'base_press_layer_1',$
'base_height_layer_2',$
'base_temp_layer_2',$
'base_press_layer_2',$
'base_height_layer_3',$
'base_temp_layer_3',$
'base_press_layer_3',$
'base_height_layer_4',$
'base_temp_layer_4',$
'base_press_layer_4',$
'base_height_layer_5',$
'base_temp_layer_5',$
'base_press_layer_5',$
'CloudBaseBestEstimate',$
'clwdn',$
'clwup',$
'column_cfrac',$
'column_cfrac_mpl',$
'confidence_flag',$
'cswup',$
'down_long',$
'down_short',$
'first_cbh',$
'gswfluxdn_clearskyfit',$
'isccp_ptop_1',$
'isccp_ptop_2',$
'isccp_ptop_3',$	 
'isccp_tau_1',$
'isccp_tau_2',$
'isccp_tau_3',$
'liq',$
'longwave_atm_cld_forcing',$
'longwave_atm_cld_forcing_obs',$
'longwave_atm_cld_forcing_obs_num',$
'longwave_sfc_cld_forcing',$
'longwave_sfc_cld_forcing_obs',$
'longwave_sfc_cld_forcing_obs_num',$
'longwave_toa_cld_forcing',$
'longwave_toa_cld_forcing_obs',$
'longwave_toa_cld_forcing_obs_num',$
'mxrat_adj_flag',$
'mwr_confidence_flag',$
'net_atm_cld_forcing',$
'net_atm_cld_forcing_obs',$
'net_atm_cld_forcing_obs_num',$
'net_sfc_cld_forcing',$
'net_sfc_cld_forcing_obs',$
'net_sfc_cld_forcing_obs_num',$
'net_toa_cld_forcing',$
'net_toa_cld_forcing_obs',$
'net_toa_cld_forcing_obs_num',$
'number_of_layers',$
'opt_depth_iwc_flag',$
'precip_flag',$
'precip_rate',$
'RadarFirstTop',$
'sfc_ir_down_clear',$
'sfc_ir_up_clear',$
'sfc_mxrat',$
'sfc_pressure',$
'sfc_solar_down_clear',$
'sfc_solar_up_clear',$
'sfc_temp',$
'shortwave_atm_cld_forcing',$
'shortwave_atm_cld_forcing_obs',$
'shortwave_atm_cld_forcing_obs_num',$
'shortwave_sfc_cld_forcing',$
'shortwave_sfc_cld_forcing_obs',$
'shortwave_sfc_cld_forcing_obs_num',$
'shortwave_toa_cld_forcing',$
'shortwave_toa_cld_forcing_obs',$
'shortwave_toa_cld_forcing_obs_num',$
'skin_temp',$
'sky_ir_temp',$
'solar_cos_z',$
'Solar_Zenith_AngleSite',$
'SolarRatio',$
'SurfaceAlbedo',$
'TOA_IR_up','TOA_IR_up_sdv','TOA_IR_up_total','TOA_IR_up_valid',$
'TOA_IR_up_clear','TOA_IR_up_clear_sdv','TOA_IR_up_clear_total','TOA_IR_up_clear_valid',$
'TOA_solar_down','TOA_solar_down_sdv','TOA_solar_down_total','TOA_solar_down_valid',$
'TOA_solar_down_clear',$
'TOA_solar_down_clear_sdv','TOA_solar_down_clear_total','TOA_solar_down_clear_valid',$
'TOA_solar_up','TOA_solar_up_sdv','TOA_solar_up_total','TOA_solar_up_valid',$
'TOA_solar_up_clear',$
'TOA_solar_up_clear_sdv','TOA_solar_up_clear_total','TOA_solar_up_clear_valid',$
'top_height_layer_1',$
'top_temp_layer_1',$
'top_press_layer_1',$
'top_height_layer_2',$
'top_temp_layer_2',$
'top_press_layer_2',$
'top_height_layer_3',$
'top_temp_layer_3',$
'top_press_layer_3',$
'top_height_layer_4',$
'top_temp_layer_4',$
'top_press_layer_4',$
'top_height_layer_5',$
'top_temp_layer_5',$
'top_press_layer_5',$
'up_long',$
'up_short',$
'vap']
vardim1=make_array(n_elements(var1),/int,value=1)

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

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

;******************************************************
;  Create the name of the file to read
;******************************************************

;  Separate out the year because the 5min files are 
;  in yearly subdirectories
syear=strmid(date,0,4)

;  Arm data directory
arm_dir='/uufs/chpc.utah.edu/common/home/mace-group3/arm/'+site+'/'

;  Selects the intercomparison files that are finally transfered to the archive
if avg_int eq '5min' and ret_type eq 'column' then begin
  file_dir='sgp.Integrate_Column_5min-Mace/'
  arm_file='sgp.Integrate_Column_5min-Mace.'+date+'*.000000.cdf*'
  datastream='sgp.Integrate_Column_5min-Mace'

;  Selects the 5min retrieval file in the basic 5min directory
endif else if avg_int eq '5min' and ret_type eq 'temp' then begin
  file_dir=site+'.average.'+avg_int+'/'+syear+'/'
  arm_file=site+'.average.5min.'+date+'*temp.cdf*'
  datastream=site+'.average.5min.'+ret_type    

;  Selects the 5min retrieval files
endif else if avg_int eq '5min' and $
  (ret_type ne 'mine' and ret_type ne 'archive') then begin
  file_dir=site+'.average.5min.'+ret_type+'/'+syear+'/'
  arm_file=site+'.average.5min.'+date+'*temp.cdf*'
  datastream=site+'.average.5min.'+ret_type    

;  Selects the 5min arm only files
endif else if avg_int eq '5min' and ret_type eq 'mine' then begin
  file_dir=site+'.average.'+avg_int+'/'+syear+'/'
  arm_file=site+'.average.'+avg_int+'.'+date+'*.000000.cdf*'
  datastream=site+'.average.5min'   

;  Selects the 5min files for transfer to the arm archive
endif else if avg_int eq '5min' and ret_type eq 'archive' then begin
  file_dir=site+'.average.'+avg_int+'.'+ret_type+'/'
  arm_file=site+'.average.'+avg_int+'.'+date+'*.000000.cdf*'
  datastream=site+'.average.5min'   

;  Selects the 5min to Xhr arm only files
endif else if avg_int ne '5min' and ret_type eq 'mine' then begin
  file_dir=site+'.average.5min.to.'+avg_int+'.'+ret_type+'/'
  arm_file=site+'.average.5min.to.'+avg_int+'.'+ret_type+'.'+$
    date+'*cd*'
  datastream=site+'.average.5min.to.'+avg_int   

;  Selects the 5min to Xhr intercomparison files for the arm archive
endif else if avg_int ne '5min' and ret_type eq 'column' then begin
  file_dir=site+'.Integrate_Column_'+avg_int+'-Mace/'
  arm_file=site+'.Integrate_Column_'+avg_int+'-Mace.'+date+'*.000000.cdf*'
  datastream=site+'.Integrate_Column_'+avg_int+'-Mace'   

;  Selects the monthly files
endif else if avg_int eq '1mth' then begin
  file_dir=site+'.average.'+avg_int+'.'+ret_type+'/'+syear+'/'
  arm_file=site+'.average.'+avg_int+'.'+ret_type+'.'+date+'*cd*'
  datastream=site+'.average.'+avg_int+'.'+ret_type

;  Selects the year file of monthly averages
endif else if avg_int eq '1yr' then begin
  file_dir=site+'.average.1mth.'+ret_type+'/'
  arm_file=site+'.average.monthly.'+ret_type+'.'+date+'*cd*'
  datastream=site+'.average.monthly.'+ret_type

;  Selects the 5min to Xhr retrieval files
endif else begin
  file_dir=site+'.average.'+avg_int+'.'+ret_type+'/'+syear+'/'
  arm_file=site+'.average.'+avg_int+'.'+ret_type+'.'+date+'*cd*'
  datastream=site+'.average.'+avg_int+'.'+ret_type
endelse

;  Name of the file
armdatadir=arm_dir+file_dir 

;  Unzip the data files
arm_file=armdatadir+arm_file
unix_command='gunzip -fq '+arm_file
spawn,unix_command

;  Get a list of the matching file names
print, 'filename ',arm_file
files=file_search(arm_file,count=num_files)
if dp eq 'yes' then print, 'Found ',num_files,' files'

;*****************************************************
;  Only continue if there are files to process
;*****************************************************

if num_files gt 0 then begin

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

;  Loop through all the matching files in the directory
for f=0,num_files-1 do begin

;  Print the filename
  print,f,files[f]
  
;  Open the netcdf file
  cdf_id=ncdf_open(files[f])

;  Get the dimension ids
  time_did=ncdf_dimid(cdf_id,'time')
  height_did=ncdf_dimid(cdf_id,'height')
  scn_type_did=ncdf_dimid(cdf_id,'scn_type')
  cld_type_did=ncdf_dimid(cdf_id,'cld_type')
  
;  Get the dimensions
  if time_did ne -1 then ncdf_diminq,cdf_id,time_did,charstring,numtimes
  if height_did ne -1 then ncdf_diminq,cdf_id,height_did,charstring,numheights
  if scn_type_did ne -1 then $
    ncdf_diminq,cdf_id,scn_type_did,charstring,numscans $
  else numscans=2
  if cld_type_did ne -1 then $
    ncdf_diminq,cdf_id,cld_type_did,charstring,numcldtypes $
  else numcldtypes=4

;  Create array for missing data.  Value selected at the top of the program
  if time_did ne -1 then begin
    dummy1d=make_array(numtimes,/float,value=missing)
    dummy3d=make_array(numscans,numtimes,/float,value=missing)
    dummy4d=make_array(numcldtypes,numtimes,/float,value=missing)
  endif
  if height_did ne -1 and time_did ne -1 then $
    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' $
      else if dim eq 4 then dumname='dummy4d'
    endelse
    if dp eq 'yes' then print,'processing ',varname[v],'dimension',dim

;  Get the variable
    if dim eq 1 or dim eq 2 or dim eq 3 or dim eq 4 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 eq 2 or dim eq 3 or dim eq 4 then begin
      if dp eq 'yes' then print,'transposing',varname[v]
      xstr=varname[v]+'=transpose('+varname[v]+')'
      result=execute(xstr)
    endif
  endfor

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

endfor  ;end of loop through 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 'height' then begin
    xstr=varname[v]+'=com_'+varname[v]
    result=execute(xstr)
  endif

endfor  ;end of loop through variables

;*****************************************
;  Do some conversions of the retrieval variables
;*****************************************

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

  if varname[v] eq 'IWC_best_estimate' or $
     varname[v] eq 'lwc_best_estimate' or $
     varname[v] eq 're_liquid_best_estimate' or $
     varname[v] eq 're_ice_best_estimate' or $
     varname[v] eq 'ice_fraction_CCM3' or $
     varname[v] eq 'ice_tau_solar' or $
     varname[v] eq 'ice_tau_IR' or $     
     varname[v] eq 'liquid_tau_solar' or $
     varname[v] eq 'liquid_tau_IR' or $
     varname[v] eq 'liquid_g_solar' or $
     varname[v] eq 'liquid_g_IR' or $
     varname[v] eq 'ice_g_solar' or $
     varname[v] eq 'ice_g_IR' or $
     varname[v] eq 'ice_omega_solar' or $
     varname[v] eq 'ice_omega_IR' or $
     varname[v] eq 'liquid_omega_solar' or $
     varname[v] eq 'liquid_omega_IR' then begin

    ;  assign the array to data
    xstr='data='+varname[v]
    result=execute(xstr)

    data=float(data)
    result=where(data gt -8888,count)
    if count gt 0 then begin
      data[result]=10.0^(float(data[result])/1000.0)
    endif
    result_no_cloud=where(data gt 1e-9 and $
                          data lt 9.99e-9,count_no_cloud)
    result_missing=where(data gt 1e-8 and $
                         data lt 9.99e-8,count_missing)
    if count_no_cloud gt 0 then data[result_no_cloud]=-9999
    if count_missing gt 0 then data[result_missing]=-8888

    ;  assign the converted data array to the varname
    xstr=varname[v]+'=data'
    result=execute(xstr)

  endif  ;end of one of these varnames       

endfor  
  
;*****************************************
;  Fill in the blanks of the 1 & 2d arrays
;******************************************

if avg_int eq '5min' and fill_flag eq 'yes' then begin

  ;  Find the 2d variables
  for v=0,n_elements(varname)-1 do begin
    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 ge 2 then begin
      if dim eq 1 or dim eq 2 then begin
        if n_elements(var2d) eq 0 then begin
           var2d=varname[v] 
        endif else begin
           var2d=[var2d,varname[v]]
        endelse
      endif  ;end of this is a 2d var
    endelse  ;end of found this varname in the list
  endfor  ;end of loop through varnames

  ;  If there are variables to fill then continue 
  if n_elements(var2d) gt 0 then begin

  ;  Fill in the internal blanks
    for v=0,n_elements(var2d)-1 do begin
      old_time_offset=time_offset
      xstr='insert_blanks,'+var2d[v]+',old_time_offset,height,missing
      result=execute(xstr)
    endfor

    time_offset=old_time_offset

  ;  Fill in the blanks around the arrays if we are grabbing a month of 5min files
    if strlen(date) eq 6 then begin
      curdate=date+'01'
      for v=0,n_elements(var2d)-1 do begin
        old_time_offset=time_offset & old_base_time=base_time
        xstr='insert_blank_beginning,'+var2d[v]+$
             ',curdate,missing,old_base_time,old_time_offset,height'
        result=execute(xstr)
      endfor
      time_offset=old_time_offset & base_time=old_base_time
    endif

  endif  ;fill the 2d variable
      
endif  ;fill data

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

endif  ;end of there were files to process

return
end