;;++++++++++++++++++++++++++++++++++++++++++++++++++++++++

undef("parse_MODIS")
function parse_MODIS (fnam[1]:string)
; MODIS file names are standardized:
; MOD04_L2.AYYYYDDD.HHMM.VVV.YYYYDDDHHMMSS.hdf 
; MOD04_L2.A2014305.0000.006.2015077213911.pscs_000500937810.hdf
; 0123456789012345678901
begin
   yyyyddd = todouble( str_get_cols(fnam,10,16) )
   hhmm    = todouble( str_get_cols(fnam,18,21) )
   ydhm    = yyyyddd*10000 + hhmm
   return( ydhm )
end

undef("addMeta_MODIS")
procedure addMeta_MODIS (x[*][*][*]:numeric, time[*]:numeric
                        ,LongName[1]:string, Units[1]:string)
; add meta data to existing variable
begin
  x!0 = "time"
  x!1 = "Cell_Along_Swath"
  x!2 = "Cell_ACross_Swath"
  x&time = time

  x@long_name = LongName
  if (Units.ne."") then
      x@units = Units
  end if
end

;;++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;;                        MAIN
;;++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  netCDF = True ; create netCDF
  var   = "Optical_Depth_Land_And_Ocean"
  
  diri  = "./"
  fili  = systemfunc("cd "+diri+"; ls MOD04_L2*.hdf")
  nfili = dimsizes(fili)
  print(fili)


; give one extra dimension cuz not all swaths are the same size.
       
  mAcross  = 135    
  nAlong   = 203
  nAlong1  = nAlong + 1

  time  = new (   nfili                  , "double", "No_FillValue"   )
  data  = new ( (/nfili,nAlong1,mAcross/), "float" , -9999.0)
  lat2d = new ( (/nfili,nAlong1,mAcross/), "float" ,  -999.0)
  lon2d = new ( (/nfili,nAlong1,mAcross/), "float" ,  -999.0)

  do nf=0,nfili-1                  ; loop over files
     a = addfile(diri+fili(nf),"r")
     time(nf)  = (/ parse_MODIS( fili(nf) ) /)

     x        := short2flt(a->$var$) 
     dimx      = dimsizes(x)

     if (dimx(0).eq.203 .and. dimx(1).eq.135) then
         data(nf,0:202,:)  = (/ x /)
         lat2d(nf,0:202,:) = (/ a->Latitude  /)
         lon2d(nf,0:202,:) = (/ a->Longitude /)
     else
         data(nf,:,:)  = (/ x /)
         lat2d(nf,:,:) = (/ a->Latitude  /)
         lon2d(nf,:,:) = (/ a->Longitude /)
     end if

  end do

; manually add meta data

  time!0 = "time"
  time@long_name = "current time as YYYYDDDHHMM"
  time&time      = time
  printVarSummary(time)
 
;;LongName  = "AOT at 0.55 micron for both ocean (Average) (Quality flag=1,2,3) and land (corrected) (Quality flag=3)"
  LongName  = " AOT at 0.55 micron for both ocean and land"
  addMeta_MODIS(data, time, LongName, "")
  addMeta_MODIS(lat2d, time, "Geodetic Latitude" , "degrees_north")
  addMeta_MODIS(lon2d, time, "Geodetic Longitude", "degrees_east" )

  printVarSummary(data)
  print("data: min="+min(data)+"   max="+max(data))
  printVarSummary(lat2d)
  print("lat2d: min="+min(lat2d)+"   max="+max(lat2d))
  printVarSummary(lon2d)
  print("lon2d: min="+min(lon2d)+"   max="+max(lon2d))
  print("-----")
  
;===================================================================
; netCDF
;===================================================================
  
    ;diro  = "/home/bliujuss/Plots/"
     diro  = "./"
     filo  = "MODIS_L2.Swath.nc"
     system("rm -f " + diro+filo)
     ncdf = addfile (diro + filo, "c")

;===================================================================
; create global attributes of the file (optional)
;===================================================================
     fAtt               = True            ; assign file attributes
     fAtt@title         = "MODIS L2 Swath"
   ;;fAtt@source_file   =  .........
     fAtt@Conventions   = "None"
     fAtt@creation_date = systemfunc ("date")
     fileattdef( ncdf, fAtt )            ; copy file attributes

;===================================================================
; Write variables with meta data
;===================================================================
     ncdf->AOT          = data
     ncdf->LAT2D        = lat2d
     ncdf->LON2D        = lon2d