;***************************************************************
; binning_2.ncl
;
; Concepts illustrated:
;   - Reading multiple fileslwith data from MODIS "MCD19A2"
;   - Binning and summing the data from each pass [swath]
;   - Averaging binned data
;   - Plotting
;   - Creating a netCDF file
;
;***************************************************************
; 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) 
   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 
      time_info(n,1) = YYYY      
      time_info(n,2) = DDD
      time_info(n,3) = yyyymmdd
      time_info(n,4) = mm           
      time_info(n,5) = dd           
   end do
   
   return (time_info)   
end

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

  vNam  = "Optical_Depth_047_grid1km"
  PLOT  = True

;*****************************************************************
; binned grid definition     
; The grid must be constant spacing in lat and in lon.
; The spacing in the lat and lon directions may be different.
; so: 2x2, 3x5, etc
;*****************************************************************

  nlat  = 3600                     ; grid coordinates
  mlon  = 7200      
  lat   = latGlobeFo(nlat,"lat","latitude","degrees_north")
  lon   = lonGlobeFo(mlon,"lon","longitude","degrees_east")
  lon   = (/ lon - 180. /)  ; subtract 180 from all values 
  lon&lon = lon             ; update coordinates

;*****************************************************************
; Variables to hold binned quantities
;*****************************************************************

  GBIN  = new ( (/nlat,mlon/), float ) 
  GKNT  = new ( (/nlat,mlon/), integer ) 

  GBIN  = 0.0                      ; initialization
  GKNT  = 0

;*****************************************************************
; File info                          
;     HDFEOSVersion : HDFEOS_V2.19
;     StructMetadata_0 : GROUP=SwathStructure
;*****************************************************************

  diri = "./"
  fili = systemfunc("cd "+diri+" ; ls MCD19A2.A2019071*")
  nfil = dimsizes( fili )
  print("nfil="+nfil)  
  f0   = addfile(diri+fili(0)+".he2", "r")   ;  .he2 causes NCL to add lat/lon arrays 
  x    = short2flt( f0->$vNam$   )   ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km )
 
;*****************************************************************
; Loop over the files: Sum the quantities
; *NO* data filtering is performed: No range checking nor use of AOD_QA_grid1km
;*****************************************************************
  tStrt = systemfunc("date")         ; time the loop (wall clock)

  do nf=0,nfil-1                     ; loop over all files
     print(nf+"   "+fili(nf))
     f       = addfile(diri+fili(nf)+".he2", "r")   ;  .he2 causes NCL to add lat/lon arrays 
                                     ; read data
     lat2d   = f->GridLat_grid1km 
     lon2d   = f->GridLon_grid1km 
     x       = short2flt( f->$vNam$   )   ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km )
  printMinMax(x,0)
     dimx    = dimsizes(x)
     norbit  = dimx(0)

    do no=0,norbit-1                 ; loop over all orbits within the current file  
        bin_sum(GBIN,GKNT,lon,lat   \
               ,ndtooned(lon2d), ndtooned(lat2d),ndtooned(x(no,:,:)) )
    end do

     delete(  x  )                   ; size may change
  end do
  
  wallClockElapseTime(tStrt, "Main Sum Loop: nlat="+nlat+", mlon="+mlon, 0)

;*****************************************************************
; Perform averaging
;*****************************************************************
                                    ; avoid division by 0.0
  GKNT     = where(GKNT.eq.0 , GKNT@_FillValue, GKNT)
  GBIN     = GBIN/GKNT
 
  GBIN!0   = "lat"
  GBIN!1   = "lon"
  GBIN&lat =  lat
  GBIN&lon =  lon

  copy_VarCoords(GBIN, GKNT)        ; copy coords

  if (isfilevaratt(f, vNam, "long_name")) then
      GBIN@long_name = "BINNED AVERAGE: "+vNam
      GKNT@long_name = "BINNED COUNT: "+vNam
  end if

  if (isfilevaratt(f, vNam, "units")) then
      GBIN@units     = f->$vNam$@units
  end if

;*****************************************************************
; time/date
;*****************************************************************
  dateinfo= parseFileName(fili(0))

  yyyyddd = dateinfo(0,0)    ; this is the date info for the 1st file
  yyyy    = dateinfo(0,1)
  ddd     = dateinfo(0,2)
  yyyymmdd= dateinfo(0,3)
  month   = dateinfo(0,4)
  day     = dateinfo(0,5)

;*****************************************************************
; PLOTS   
;*****************************************************************

  if (PLOT) then

      setvalues NhlGetWorkspaceObjectId()
                "wsMaximumSize" : 200000000
      end setvalues   

      plot    = new ( 2, "graphic")

     
      wks     = gsn_open_wks("png", "MCD19A2_binning_"+yyyymmdd); send graphics to PNG file
      
      res                     = True          ; Use plot options
      res@gsnDraw             = False
      res@gsnFrame            = False

      res@cnFillOn            = True          ; Turn on color fill
    ;;res@cnFillPalette       = "wh-bl-gr-ye-re"  ; set color map
      res@cnFillPalette       = "BlAqGrYeOrReVi200"  ; set color map
      res@cnFillMode          = "RasterFill"  ; Turn on raster color
      res@cnLinesOn           = False         ; Turn off contourn lines

      res@mpCenterLonF        = 180           ; set map center at 180
      res@lbOrientation       = "vertical"

      plot(0) = gsn_csm_contour_map(wks,GBIN, res)
      plot(1) = gsn_csm_contour_map(wks,GKNT, res)

;*****************************************************************
; Panel   
;*****************************************************************
      dlat    = lat(1)-lat(0)

      resP                     = True
      resP@gsnPanelMainString  = "MCD19A2: "+yyyy+"/"+month+"/"+day+ \
                                 ": "+nlat+"x"+mlon+" dlat="+sprintf("%5.3f", dlat)
      resP@gsnMaximize         = True          ; max size 
     ;resP@gsnPaperOrientation = "portrait"
      gsn_panel(wks,plot,(/2,1/),    resP)

  end if