PRO READ_MODIS

COMPILE_OPT IDL2

directory = '../products/'
product = '*.mod04.hdf'
sdsname = 'Optical_Depth_Land_And_Ocean'
scale = 0.001

;- Get list of input files
list = file_search(directory + product)
help, list

;- Set old file date
old_date = 'a1.99001'

;- Get output file unit
get_lun, lun

;- Loop over input files
for i = 0, n_elements(list) - 1 do begin

  ;- Get input file name
  input_file = list[i]
  print, input_file
      
  ;- Read the input file
  hdfid = hdf_sd_start(input_file)
  hdf_sd_varread, hdfid, 'Latitude', lat
  hdf_sd_varread, hdfid, 'Longitude', lon
  hdf_sd_varread, hdfid, sdsname, data
  hdf_sd_end, hdfid

  ;- Apply scale factor
  data = data * scale

  ;- Get file date (e.g., 'a1.11210')
  file_name = basename(input_file)
  date_length = strlen(old_date)
  new_date = strmid(file_name, 0, date_length)
  
  ;- Find valid data values
  loc = where(data gt 0.0, count)
  print, count
  
  ;- Save valid data values
  if (count gt 0) then begin

    ;- Open a new output file if needed
    if (new_date ne old_date) then begin
      free_lun, lun
      output_file = new_date + '.dat'
      openw, lun, output_file, /get_lun
      old_date = new_date
      print, 'Opened new output file ', output_file
    endif

    table = fltarr(3, count)
    table[0, *] = lat[loc]
    table[1, *] = lon[loc]
    table[2, *] = data[loc]
    help, table
    writeu, lun, table

  endif
   
endfor

END