; ================================================================================================
pro l3a_daisy
  ; ================================================================================================
  ; Uses the nearest neighbor approach to get the Level 3A albedo for a single day. That day is
  ; currently set to be the latest day for which Prelim Level 2 data is available in NH18.
  ; This version uses the match_2d routine, which is open source code developed by JD Smith of
  ; the University of Arizone.
  ; After the 3a albedo array is generated, a placeholder actual l3a file is read to fill out the
  ; Level 3 data structure, with the albedo array overwritten by the match_2d results..
  ; The draw_daisy routine is then called to make the png file (NOTE - use the version of
  ; draw_daisy.pro provided in the local directory, as it's been modified to add the correct year).
  ; ------------------------------------------------------------------------------------------------

  set_cips_production_vars
  set_debug_level, 5

  set_production_flag, flag=0

  curr_time=fix(systime(/seconds, /julian, /utc)-0.5, type=3) + 0.5
  start_time=jd2gps(curr_time)*1e6
  aim_day=usec2ad(start_time)
  print, get_routine_name() + "for " + strtrim(string(aim_day -1), 2)

  anc_dir=!COMMON_INPUT_DATA + '/grids/'

  restore,anc_dir+'dlon_grid_7.5km.sav'

  ; L3 grid
  restore,anc_dir+'level_3a_lat_lon_north_v5.sav'
  lat_l3=temporary(latitude)
  lon_l3=temporary(longitude)
  s=size(lat_l3,/dim)
  n3=s[0]

  ; Define albedo threshold
  alb_thresh=0.

  hem='north'
  version='05.20'
  revision='02'
  year='2018'

  ; Path to Level 2 files
  season=hem+'_'+year
  l2_path='/aim/data/cips/'+season+'/level_2/ver_'+version+'/rev_'+revision+'/'
  get_cat_files,l2_path,0,100000,nrev,rev,files,yy,doy,/prelim

  ; Pick latest day
  d0=max(doy)
  x=where(doy eq d0,norbits)
  rev=rev[x]
  files=files[x]

  ; Loop over orbits and save all L2 data

  lat_tot = []
  lon_tot = []
  alb_tot = []

  for irev=0,norbits-1 do begin

    catfile=files[irev]
    cldfile=catfile.replace('cat','cld')
    etcfile=catfile.replace('cat','etc')

    cat=read_cips_file(l2_path+catfile)
    cld=read_cips_file(l2_path+cldfile)
    etc=read_cips_file(l2_path+etcfile)

    ; Select all finite pixels to save

    gd_data=where(finite(cat.zenith_angle_ray_peak))
    lat=cat.latitude[gd_data]
    lon=cat.longitude[gd_data]
    alb=etc.alb_i[gd_data]
    cld_presence=cld.cloud_presence_map[gd_data]

    ; Cloud pixels
    cld=where(cld_presence eq 1 and alb ge alb_thresh,complement=nocld)

    ; Make sure albedo is 0 for non-cloud pixels
    alb[nocld]=0.

    asc=where(abs(lat) gt 90.,complement=dsc)
    if(cat.hemisphere eq 'N' and asc[0] ne -1) then lat[asc]=180.-lat[asc]
    if(cat.hemisphere eq 'S' and asc[0] ne -1) then lat[asc]=-180.-lat[asc]

    lat_tot = [lat_tot, lat ]
    lon_tot = [lon_tot, lon ]
    alb_tot = [alb_tot, alb ]
    heap_gc
  endfor

  ntot=n_elements(lon_tot)

  alb_l2 = temporary(alb_tot)
  lon_l2 = temporary(lon_tot)
  lat_l2 = temporary(lat_tot)

  alb_l3 = replicate(!values.f_nan, 1700L, 1700L)

  del_lon=interpol(dlon,lat_grid,lat_l3)
  b = [dlon[UNIQ(dlon, SORT(dlon))], 7]

  for i =0, n_elements(b)-2 do begin
    ook = where(del_lon gt b[i] and del_lon le b[i+1]) ;; not using exact lon bins here.... selecting bands of bins. This can be improved.

    ;;Match-2d is a modified coyote libraries routine
    ; testing....
    match=match_2d(lon_l3(ook), lat_l3(ook), alb_l2, lon_l2, lat_l2, [B[i]/2.0, dlat/2.], MATCH_DISTANCE=temp, ALB=temp2, missing=!values.f_nan)
    ;      match=match_2d(lon_l3(ook), lat_l3(ook), alb_l2, lon_l2, lat_l2, [B[i]/2.0, 0.05], MATCH_DISTANCE=temp, ALB=temp2, missing=!values.f_nan)

    alb_l3[ook] = temp2

  endfor

  ; ------------------------------------------------------
  ; Now read in dummy L3a file and make the daisy plot
  ; ------------------------------------------------------
  l3_path='/aim/data/cips/north_2014/level_3/ver_05.10/rev_01/'
  f=file_search(l3_path,'*.nc.*')
  ndays=n_elements(f)
  files=strarr(ndays)
  for id=0,ndays-1 do files[id]=strmid(f[id],strlen(l3_path),strlen(f[id])-strlen(l3_path))
  doy=fix(strmid(files,17,3))

  x=where(doy eq d0)
  id=x[0]
  file=files[id]
  p=strpos(file,'.nc')

  png_file=strmid(file,0,p)+'_prelim.png'
  png_file=png_file.replace('2014','2018')
  png_file=png_file.replace('5.10','5.20')
  png_file=png_file.replace('r01', 'r02')

  level_3a_daisy = read_daisy_file(l3_path+file,/net_cdf)
  *level_3a_daisy.albedo=alb_l3

  lower_threshold = 3.0

  albedo_index = WHERE( FINITE(*(level_3a_daisy.albedo)) AND *(level_3a_daisy.albedo) GE lower_threshold, num_pixels )
  IF num_pixels GT 1 THEN BEGIN
    mean_albedo   = MEDIAN( (*(level_3a_daisy.albedo))[albedo_index] )
    stddev_albedo = STDDEV( (*(level_3a_daisy.albedo))[albedo_index] )
    upper_threshold = mean_albedo + (2.0 * stddev_albedo) + 20.0
  ENDIF ELSE BEGIN
    upper_threshold =  20.0
  ENDELSE
  IF upper_threshold LT 20.0 THEN upper_threshold = 20.0

  orbit_location_index = WHERE( finite(*level_3a_daisy.albedo) AND *level_3a_daisy.albedo LT lower_threshold )
  IF orbit_location_index[0] NE -1 THEN (*level_3a_daisy.albedo)[orbit_location_index] = lower_threshold
  p_flag = get_production_flag( )
  IF p_flag EQ 1 THEN preliminary = 0 ELSE preliminary = 1

  temp_fix = '/aim/data/cips/north_2018/level_3/ver_05.20/rev_02/'
  draw_daisy_new,level_3a_daisy,temp_fix + png_file,/west_ct,kmperpix=level_3a_daisy.km_per_pixel, bottom=lower_threshold,  $
    top = upper_threshold,preliminary = preliminary, center_lons = center_longs, orbit_numbers = orbit_numbers

end