;-------------------------------------------------------------
; ncl-talk: Melanie O' Hanoly [ 2/12/2020 ]
; How do I calculate the average in the same months, lat and lon for the variable (NPP)?
;-------------------------------------------------------------
;---netCDF file created via a different script
;-------------------------------------------------------------
  dirNC  = "./"
  filNC  = "NPP_EDITED.nc"
  pthNC  = dirNC+filNC

  f      = addfile(pthNC, "r")
  month  = f->MONTH
  lat    = f->LAT
  lon    = f->LON
  npp    = f->NPP

;-------------------------------------------------------------
;---Create unique location identification string
;-------------------------------------------------------------
  locid  = tostring(lat)+"_"+tostring(lon) ; create location id string
  locunq = get_unique_values(locid)        ; unique lat/lon locations
  nunq   = dimsizes(locunq)                ; number of unique locations

  print("nunq="+nunq)
  print("======")
;;print(locunq)                          ; not pretty
;;print("======")

;-------------------------------------------------------------
;---Monthly averages for each month at each unique location
;-------------------------------------------------------------
  nmos   = 12
  nppLat = new(  nunq       , typeof(lat), "No_FillValue")
  nppLon = new(  nunq       , typeof(lon), "No_FillValue")
  nppAvg = new((/nunq,nmos/), typeof(npp), getVarFillValue(npp))
  nppKnt = new((/nunq,nmos/), "integer"  , -99)   ; number of values used
                                                  ; to compute average
  do nu=0,nunq-1                         ; loop over each unique location
    do mon=1,nmos                        ; loop over each month
       locmon := ind(locid.eq.locunq(nu) .and. month.eq.mon)
       if (.not.ismissing(locmon(0))) then
           nppLat = tofloat(str_get_field(locunq(nu),1,"_"))
           nppLon = tofloat(str_get_field(locunq(nu),2,"_"))
           nppAvg(nu,mon-1) = avg(npp(locmon))
           nppKnt(nu,mon-1) = num(.not.ismissing(npp(locmon)))
       end if
    end do               ; mon
  end do                 ; nu: unique locations

;-------------------------------------------------------------
;---Create netCDF containing the statistics
;-------------------------------------------------------------
 
  nppAvg!0 = "location"
  nppAvg!1 = "month"
  nppAvg&month     =  ispan(1,nmos,1)
  nppAvg@long_name = "NPP average"
  nppAvg@units     =  npp@units 

  nppKnt!0 = "location"
  nppKnt!1 = "month"
  nppKnt&month=  ispan(1,nmos,1)
  nppKnt@long_name = "Number of values used to compute average"

  nppLat!0 = "location"
  nppLat@long_name = "Latitude"
  nppLat@units     = lat@units

  nppLon!0 = "location"
  nppLon@long_name = "Longitude"
  nppLon@units     = lon@units

  diro   = "./"
  filo   = "NPP_AVG.nc"
  ptho   = diro+filo
  
  system("/bin/rm -f "+ptho)    ; remove any pre-existing file
  ncdf = addfile(ptho  ,"c")    ; open output netCDF file

  fAtt               = True            ; assign file attributes
  fAtt@title         = "NPP csv file Monthly Averages"
  fAtt@source_file   = filNC  
  fAtt@ncltalk       = "Melanie O Hanoly: 2/12/2002"
  fAtt@Conventions   = "None"
  fAtt@creation_date = systemfunc ("date")
  fileattdef( ncdf, fAtt )            ; copy file attributes

  ncdf->LAT     = nppLat
  ncdf->LON     = nppLon       
  ncdf->NPP_AVG = nppAvg
  ncdf->NPP_KNT = nppKnt