;***************************************************************
; Utility functions to get (1) date information (2) facilitate netCDF 
;********************************************************************
; MCD19A2.AYYYYDDD.hHHvVV.CCC.YYYYDDDHHMMSS.hdf
; YYYYDDD = Year and Day of Year of acquisition
; hHH = Horizontal tile number (0-35)
; vVV = Vertical tile number (0-17)
; CCC = Collection number
; YYYYDDDHHMMSS = Production Date and Time
; Example:
; MCD19A2.A2019161.h32v11.006.2019163201550.hdf
;********************************************************************
 undef("parseFileName")
function parseFileName ( fNam[*]:string )
; parse the file name(s) and extract data information
; https://www.ncl.ucar.edu/Document/Functions/Built-in/str_get_cols.shtml
; 012345678901234567890123456789012345678901234
;           1         2         3         4
; MCD19A2.AYYYYDDD.hHHvVV.CCC.YYYYDDDHHMMSS.hdf
;          YYYYDDD  tile
local n, nNam, time_info, YYYY, DDD, yyyyddd, yyyymmdd, mmdd, mm, dd  
begin  
   nNam = dimsizes( fNam )

   YYYY      = toint(str_get_cols(fNam, 9,12))     ; YYYY 
   DDD       = toint(str_get_cols(fNam,13,15))     ; DDD

   yyyyddd   = YYYY*1000 + DDD
   yyyymmdd  = (/ yyyyddd_to_yyyymmdd(yyyyddd) /)  ; (/ ... /) no meta data
   mmdd      = yyyymmdd-(yyyymmdd/10000)*10000
   mm        = mmdd/100                            
   dd        = mmdd-mm*100

   time_info    = new( (/nNam,6/), integer,"No_FillValue")

   do n=0,nNam-1
      time_info(n,0) = yyyyddd(n) 
      time_info(n,1) = YYYY(n)      
      time_info(n,2) = DDD(n)
      time_info(n,3) = yyyymmdd(n)
      time_info(n,4) = mm(n)           
      time_info(n,5) = dd(n)           
   end do
   
   return (time_info)   
end

;**************************************************************
;                     MAIN
;**************************************************************

  vNam   = "Optical_Depth_047_grid1km"
  latloc =  40.821                           ; can be more than on location
  lonloc = -73.949  
  print("latloc"+latloc+"  lonloc="+lonloc )
  print("=============================")

  diri = "./"
  fili = systemfunc("cd "+diri+" ; ls MCD19A2.A*.hdf")
  print(fili)
  nfil = dimsizes( fili )
  print("nfil="+nfil)  
  print("=============================")
  time_info = parseFileName ( fili )
  print(time_info)
  print("=============================")
  print(time_info(:,0)+" "+time_info(:,1)+" "+time_info(:,2)+" "+ \
        time_info(:,3)+" "+time_info(:,4)+" "+time_info(:,5) )
  print("=============================")
  
;*****************************************************************
; Read over the files: Sum the quantities
; *NO* data filtering is performed: No range checking nor use of AOD_QA_grid1km
;*****************************************************************
  f     = addfiles(diri+fili+".he2", "r")   ;  .he2 causes NCL to add lat/lon arrays 
                                     ; read data
  lat2d = f[0]->GridLat_grid1km 
  lon2d = f[0]->GridLon_grid1km 
  nm    = getind_latlon2d (lat2d,lon2d, latloc, lonloc)
  nl    = nm(0,0)
  ml    = nm(0,1)
  print("nl="+nl+"  ml="+ml+"  lat2d(nl,ml)="+lat2d(nl,ml)+"  lon2d(nl,ml)="+lon2d(nl,ml) )
  print("=============================")
  X     = short2flt( f[:]->$vNam$(:,nl,ml) )   ; ( Orbits_grid1km )
  print(X)
  print("=============================")