undef("get_mode_int")
function get_mode_int(x[*]:integer)
;
; No formulas ... brute force counting.
;
local typx, nx, xi, kuxi, xu, xmode, kmode, k, ku, ku_max, kumode, kind
begin
  typx = typeof(x)
  if (.not.(typx.eq."integer")) then
      print("get_mode_int: argument must be type integer: typeof(x)="+typx)
      exit
  end if
 
  nx   = dimsizes(x) 
  xi   = x
  qsort(xi)                       ; convenience: want mode(s) returned in ascending order

;---Basic stuff
  kuxi = count_unique_values(xi)  ; count unique values: scalar [1]
  xu   = get_unique_values(xi)    ; [*]

  if (kuxi.eq.nx) then
      xmode = xi                   
      xmode@information = "All input  values are unique" 
      kmode = 1                   ; one occurence of each
      return(([/xmode,kmode/]))
  end if

;---Count the number of occurrences of each unique integer value
  ku = conform(xu, 0, -1)         ; initialize counter to 0
  do k=0,kuxi-1
     ku(k) = num(xi.eq.xu(k))
  end do

  ku_max  = max(ku)                ; max kount
  ku_mode = ind(ku.eq.ku_max)      ; index/indices of max mode count
  kind    = dimsizes(ku_mode)

  if (kind.eq.1) then     ; only one mode
      xmode = xu(ku_mode)
      kmode = ku_max
  else                    ; multiple modes
      xmode = new( kind, typx, "No_FillValue")
      kmode = ku_max
      do nk=0,kind-1
         xmode(nk) = xu(ku_mode(nk))
      end do
  end if

  xmode!0 = "mode"
  if (isatt(x,"long_name")) then
      xmode@long_name = "MODE: "+x@long_name
  end if
  if (isatt(x,"units")) then
      xmode@units = x@units
  end if
  xmode@NCL_tag = "get_mode => get_mode_int"

  return(([/xmode,kmode/]))
end
;===========
undef("get_mode_round")
function get_mode_round(x[*]:numeric)
;
; This was writen for float/double using rounded values
; No formulas ... brute force counting.
;
; It groups values +/- of value
; eg: 82.5 ==> round(82.5,0) ==> 83.0
; eg: 83.4 ==> round(83.4,0) ==> 83.0
; eg: 82.4 ==> round(82.4,0) ==> 82.0
; eg: 81.5 ==> round(81.5,0) ==> 82.0
;
local typx, nx, xi, kuxr, xu, xmode, kmode, k, ku, ku_max, ku_mode, kind
begin
  typx = typeof(x)
  if (.not.(typx.eq."float" .or. typx.eq."double")) then
      print("get_mode_round: argument must be type float or double: typeof(x)="+typx)
      exit
  end if
 
  nx   = dimsizes(x) 
  xr   = round(x,0)
  qsort(xr)                       ; convenience: want mode(s) returned in ascending order

;---Basic stuff
  kuxr = count_unique_values(xr)  ; count unique rounded values [1]
  xu   = get_unique_values(xr)    ; [*]

  if (kuxr.eq.nx) then
     ;xmode = x                   
     ;xmode@information = "All input  values are unique" 
      xmode = xr                  
      xmode@information = "All rounded values are unique" 
      kmode = 1                   ; one occurence of each
      return(([/xmode,kmode/]))
  end if

;---Count the number of occurrences of each (rounded) unique value
  ku = conform(xu, 0, -1)         ; initialize counter to 0
  do k=0,kuxr-1
     ku(k) = num(xr.eq.xu(k))
  end do

  ku_max  = max(ku)                ; max kount
  ku_mode = ind(ku.eq.ku_max)      ; index/indices of max mode count
  kind    = dimsizes(ku_mode)

  if (kind.eq.1) then     ; only one mode
      xmode = xu(ku_mode)
      kmode = ku_max
  else
      xmode = new( kind, typs, "No_FillValue")
      kmode = kind
      do nk=0,kind-1
         xmode(nk) = xu(ku_mode(nk))
      end do
  end if

  xmode!0 = "mode"
  if (isatt(x,"long_name")) then
      xmode@long_name = "MODE: "+x@long_name
  end if
  if (isatt(x,"units")) then
      xmode@units = x@units
  end if
  xmode@NCL_tag = "get_mode => get_mode_round"

  return(([/xmode,kmode/]))
end
;---
; User interface
;---
undef("get_discrete_mode")
function get_discrete_mode(x[*]:numeric)
local typx
begin
  typx = typeof(x)
  if (typx.eq."integer") then
      return( get_mode_int(x) )
  elseif (typx.eq."float" .or. typx.eq."double") then
      return( get_mode_round(x) )    ; <== uses rounded values
  else
      print("get_discrete_mode: Only type integer/float/double allowed")
      exit
  end if
end