;   NCL 6.4 0 required
;; If you do not have then uncomment the following. 
;; load "./calculate_monthly_values.640.ncl"
;;
;----------------------------------------
;---COMPUTE MONTHLY MEANS FROM DAILY DATA
;----------------------------------------

;---File information (netCDF)
   diri = "./"      ; input directory
   fili = systemfunc("cd "+diri+" ; ls AOT*20[0-1][1-9].nc")
   pthi = diri + fili
   nfili= dimsizes(fili)

   print(fili)
   print("---")
   print("nfili="+nfili)
   print("---")

;---Open and aggregate all files in specified path
   f    = addfiles(pthi,"r")
   a440   = f[:]->AOT_440
   a500   = f[:]->AOT_500
   a550   = f[:]->AOT_550
   a675   = f[:]->AOT_675
   alpha  = f[:]->ALPHA
   astrom = f[:]->ANGSTROM

   printVarSummary(a500)
   print("---")
   printVarSummary(a675)
   print("---")
   printVarSummary(alpha)
   print("---")
   printVarSummary(astrom)
   print("---")

;---Output convenience; make _FillValue nice
   a440@_FillValue   = -999.0
   a500@_FillValue   = -999.0
   a550@_FillValue   = -999.0
   a675@_FillValue   = -999.0
   alpha@_FillValue  = -999.0
   astrom@_FillValue = -999.0
 
;---Calculate monthly values
;---https://www.ncl.ucar.edu/Document/Functions/Contributed/calculate_monthly_values.shtml

   opt = True 
   opt@nval_crit = 10   ; require at least 10 values for the "avg" is calculated.
                        ; arbitrary set by user
   a440MonAvg   = calculate_monthly_values(a440  , "avg", 0, opt)
   a500MonAvg   = calculate_monthly_values(a500  , "avg", 0, opt)
   a550MonAvg   = calculate_monthly_values(a550  , "avg", 0, opt)
   a675MonAvg   = calculate_monthly_values(a675  , "avg", 0, opt)
   alphaMonAvg  = calculate_monthly_values(alpha , "avg", 0, opt)
   astromMonAvg = calculate_monthly_values(astrom, "avg", 0, opt)

;---time for netCDF
   time     = a500MonAvg&time
   delete(time@_FillValue)
   time&time= time
   printVarSummary(time)
   print("---")

;---https://www.ncl.ucar.edu/Document/Functions/Built-in/cd_calendar.shtml
   yyyymm   = cd_calendar(time,-1)
   yyyymm!0    = "time"
   yyyymm&time =  time
   printVarSummary(yyyymm)
   print("---")
  
;---http://www.ncl.ucar.edu/Document/Functions/Contributed/yyyymm_to_yyyyfrac.shtml
   yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0)
   printVarSummary(yrfrac)
   print("---")
 
   printVarSummary(a500MonAvg)
   print("---")
   printVarSummary(a675MonAvg)
   print("-------------------------------------")
   print(yyyymm+"  "+yrfrac+"  "+a440MonAvg +"   "+a500MonAvg  +"  "+a550MonAvg \
               +"  "+a675MonAvg+"  "+alphaMonAvg+"   "+astromMonAvg)
   print("---")

;---Write to ascii ".txt"
   
   ntxt     = dimsizes(yyyymm)
   txt_data = new (ntxt+1, string, "No_FillValue")
   txt_data(0)  = "yyyymm  yyyy_frac      A440        A500        A550        A675        ALPHA       ASTROM"
   txt_data(1:) = yyyymm + sprintf("%12.5f",yrfrac) \ 
                         + sprintf("%12.5f",a440MonAvg) + sprintf("%12.5f",a500MonAvg) \ 
                         + sprintf("%12.5f",a550MonAvg) + sprintf("%12.5f",a675MonAvg) \
                         + sprintf("%12.5f",alphaMonAvg)+ sprintf("%12.5f",astromMonAvg)
   print(txt_data)
   print("---")

   dtxt = "./"
   ftxt = "AOT.2005-2014.month.txt"
   ptxt = dtxt+ftxt
   system("/bin/rm -f "+ptxt)   ; rm pre-existing file (if present)
   asciiwrite(ptxt, txt_data) 

;---Write to netCDF
   dnc  = "./"
   fnc  = "AOT.2005-2014.month.nc"
   pnc  = dnc +fnc
   system("/bin/rm -f "+pnc)   ; rm pre-existing file (if present)
  
   ncdf = addfile(pnc ,"c")    ; open output netCDF file
;===================================================================
; create global attributes of the file (optional)
;===================================================================
   fAtt               = True            ; assign file attributes
   fAtt@title     = "AOT MONTHLY VALUES"
   if (opt .and. isatt(opt,"nval_crit")) then
       fAtt@information = "Monthly values calculated only if "+opt@nval_crit+" daily values available"
   else
       fAtt@information = "Monthly values calculate even if only one daily value is available"
   end if
   fAtt@Conventions   = "None"
   fAtt@creation_date = systemfunc ("date")
   fileattdef( ncdf, fAtt )            ; copy file attributes

;===================================================================
; make time an UNLIMITED dimension
;===================================================================
   filedimdef(ncdf,"time",-1,True) 

;===================================================================
; output variables directly; NCL will call appropriate functions
; to write the meta data associated with each variable
;===================================================================
    ncdf->time    = time
    ncdf->yyyymm  = yyyymm
    ncdf->yrfrac  = yrfrac
    ncdf->AOT_440_MonMean    = a440MonAvg
    ncdf->AOT_500_MonMean    = a500MonAvg
    ncdf->AOT_550_MonMean    = a550MonAvg
    ncdf->AOT_675_MonMean    = a675MonAvg
    ncdf->AOT_alpha_MonMean  = alphaMonAvg
    ncdf->AOT_astrom_MonMean = astromMonAvg