load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/calendar_decode2.ncl"

begin

;-- define file name

; read data ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  wkdir = "/home/model-user/"
  fil_dat  = "AOT_550"
; list all AOT files for different years
  files = systemfunc("ls " + wkdir +"AOT_550*.txt") 
  print(files)
; number of files
  nfiles = dimsizes(files)
  nday = 365  ; number of days in a year

; define fields with number of years (files) and number days per year
  AOT_550_all = new((/nfiles,nday/),"float")
  day_all = new((/nfiles,nday/),"float")

; loop through different years to read the files
 do n= 0,nfiles-1
  print(files(n))
  z1 = asciiread(files(n),(/nday+1/),"string")
  printVarSummary(z1)
  new_str = str_split_csv(z1,",",1)
  printVarSummary(new_str)
  variables = new_str(0,0:6)
  print(variables)
  Day = stringtoint(new_str(1:nday,0))
  AOT_675 = stringtofloat(new_str(1:nday,1))
  AOT_500 = stringtofloat(new_str(1:nday,2))
  AOT_440 = stringtofloat(new_str(1:nday,3))
  AOT_440_675 = stringtofloat(new_str(1:nday,4))
  Alpha = stringtoint(new_str(1:nday,5))
  AOT_550 = stringtofloat(new_str(1:nday,6))

; write to the fields
  day_all(n,0:nday-1) = Day
  AOT_550_all(n,0:nday-1) =  AOT_550
 end do
 AOT_550_all!0 = "years"
 AOT_550_all!1 = "days"
 printVarSummary(AOT_550_all)
 
 ;compute monthly mean from daily values julian days

  wkdir = "/home/model-user/"
  fil_dat  = "AOT_550"

  files = systemfunc("ls " + wkdir +"AOT_550*.txt") 
  print(files)
  nfiles = dimsizes(files)
  nday = 365  ; number of days in a year Julian days
   
   nyrs = max(2014)-min(2005)+1
   ntim = nyrs*12
   
    
   monavg = ((/ntim/) typeof(AOT_550), getFillValue(AOT_550)) 

   nt = -1
   do nyr=0,nyrs-1
      do nmo=0,11
           i = ind(years(nyr).eq.year .and. month(nmo).eq.(nmo)+1)
           

    MONAVG(nyr,nmo,:,:) = dim_avg(:, time|i)

      end do
    end do 


;**************************************************
; create plot
;**************************************************
  wks = gsn_open_wks("png","AOT_550_")

  res = True
  res@gsnDraw = False
  res@gsnFrame = False
  res@tiMainString = "Figure 2, Ilorin"
  res@tiYAxisString = "AOT_550 (nm)"
  res@tiXAxisString = "Days"                 ; x axis title

  res@xyMarkerColors      = (/"black","black","black","black","black","black","black","black","black","black"/)  ; line colors
  res@xyLineThicknesses = (/4,2,2,2,2,2,2,2,2,2/)        ; line thicknesses
  res@xyMarkLineModes   = (/"Markers","Markers","Markers","Markers","Markers","Markers","Markers","Markers","Markers","Markers","Markers"/)

  res@trXMinF       = 0.          ; Control limits of X
  res@trXMaxF       = 365.     ; and Y axes.
  res@trYMinF       = 0.          ; Control limits of X
  res@trYMaxF       = 2.          ; Control limits of X
  plot = gsn_xy(wks,day_all,AOT_550,res)          ; Draw an XY plot with 1 curve.
  res@xyMarkLineModes   = (/"Lines","Markers","Markers","Markers","Markers","Markers","Markers","Markers","Markers","Markers","Markers"/)
  res@xyLineColors      = (/"blue"/)  ; line colors
  plot2 = gsn_xy(wks,day,AOT_550,res)      ; Draw an XY plot with 1 curve.

  draw(wks)
  frame(wks)
 

end