PRO PLOT_MODIS, LAT, LON

COMPILE_OPT IDL2

;- Create equal angle grid
;- (lower left corner of grid is at 180W, 90S, i.e. -180.0, -90.0)
resolution = 1.0
ncol = long(360.0 / resolution)
nrow = long(180.0 / resolution)
grid = fltarr(ncol, nrow)

;- Get list of grid files
list = file_search('a1*.grd')
help, list

;- Create output arrays
data = fltarr(n_elements(list))
time = fltarr(n_elements(list))

;- Compute index of the required grid cell
lonlat_to_index, lon, lat, resolution, index
help, index

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

  ;- Read this file
  input_file = list[i]
  openr, lun, input_file, /get_lun
  readu, lun, grid
  free_lun, lun

  ;- Get data for the desired grid cell
  data[i] = grid[index]

  ;- Get year and day of year from filename
  year = long(strmid(input_file, 3, 2)) + 2000
  day_of_year = long(strmid(input_file, 5, 3))

  ;- Get month and day of month
  ydn2md, year, day_of_year, month, day

  ;- Get Julian date
  time[i] = julday(month, day, year)
      
endfor

help, data
help, time

;- Set missing values to NaN (not a number)
loc = where(data lt 0.0, count)
if (count gt 0) then data[loc] = !VALUES.F_NAN

;- Compute time series, smoothed over 16 days
smoothed_data = smooth(data, 16, /nan, /edge_truncate)

;- Load greyscale color table
loadct, 0

;- Plot the data
result = label_date(date_format=['%M %Y'])
plot, time, smoothed_data, $
  xrange=[min(time), max(time)], $
  xstyle=1, $
  psym=1, $
  color=0, background=255, $
  xtickunits='months', $
  xtickformat='label_date', $
  ytitle='Aerosol Optical Depth', $
  title='MODIS Time Series at' + string(lat, lon, format='(2f8.2)')
 
END