;----------------------------------------------------------
undef("calculate_hourly_values")
function calculate_hourly_values (x[*]:numeric, arith:string, nDim[1]:integer, opt[1]:logical)
;                  ************
; Untested
;                  ************
; Note that NCL's 'dimension reduction' which, in general, I view as a feature,
;      introduces some processing issues in the 'corner case' of one value.
;      This requires extra (nuisance) steps and leads to 'code bloat'.  :-(
;                  ************
local dnamx, TIME, TNAME, dimx, rankx, NTIM, utc_date, year, month, day, hour, minute \
     ,ntim, work, yyyy, mm, dd, dymon, it, nt, nmos, xSta, nval_crit, yrStrt, yrLast
begin

    if (nDim.ne.0) then
        print("calculate_hourly_values: unsupported or illegal: nDim="+nDim)
        exit
    end if
 
    dimx     = dimsizes(x)
    rankx    = dimsizes(dimx)                      ; # of dimensions

    if (.not.(arith.eq."ave" .or. arith.eq."avg" .or. arith.eq."sum"  \
                             .or. arith.eq."min" .or. arith.eq."max") ) then
        print("calculate_hourly_values: unrecognizezed 'arith' argument="+arith)
        exit
    end if

; an hourly mean will be calculated only if the # of values is >= nval_crit

    if (opt .and. isatt(opt,"nval_crit") .and. typeof(opt@nval_crit).eq."integer") then
        nval_crit = opt@nval_crit
    else
        nval_crit = 1
    end if

    dnamx    = getvardims(x)                       ; dimension names
    TNAME    = dnamx(nDim)                         ; typically TNAME="time"
    ntim     = dimx(nDim)                          ; size of input time dimension

;--- Create current yyyy/mm/dd/hh/mn
;   'cd_calendar' recognizes any 'calendar' attribute associated with ' x&$TNAME$ '.
;    Hence, the elements below (yyyy, month, ..) reflect that calendar.
;    No calendar attribute defaults to gregorian [standard].

    utc_date = cd_calendar(x&$TNAME$, 0)           ; (x&$TNAME$, 0) ; typically   x&time
    year     = floattointeger(utc_date(:,0))
    month    = floattointeger(utc_date(:,1))
    day      = floattointeger(utc_date(:,2))
    hour     = floattointeger(utc_date(:,3))
    minute   = floattointeger(utc_date(:,4))
   ;second   = utc_date(:,5)

    yrStrt   = min(year)
    yrLast   = max(year)

;--- Create hourly coordinate variable. Size must be determined because
;    each hour *may* have different number of time steps. 
;    An example might be a field program.

    nt       =  -1                                 ; number of individual hours
    work     = new (ntim, typeof(x&$TNAME$), getVarFillValue(x&$TNAME$)) 
    do yyyy=yrStrt,yrLast
      do mm=min(month),max(month)
         dymon  = days_in_month ( yyyy, mm) 
        do dd=1,dymon          
          do hh=0,23
             it  := ind(yyyy.eq.year .and. mm.eq.month .and. dd.eq.day .and. hh.eq.hour)
             nit  = dimsizes(it)
             if (.not.ismissing(it(0)) .and. nit.ge.nval_crit) then
                 nt = nt+1
                 work(nt) = (/ x&$TNAME$(it(0)) /)   ; use 1st time of current hour  
             end if
          end do   ; hh
        end do     ; dd
      end do       ; mm
    end do         ; yyyy
    NTIM = nt+1    ; get 1-based total

    TIME   = work(0:nt)
    TIME!0 = "TIME"
    if (isatt(x&$TNAME$,"long_name")) then
        TIME@long_name = x&$TNAME$@long_name 
    else
        TIME@long_name = "time corresponding to 1st time of current day"
    end if
    if (isatt(x&$TNAME$,"units")) then
        TIME@units     = x&$TNAME$@units 
    end if
    if (isatt(x&$TNAME$,"calendar")) then
        TIME@calendar  = x&$TNAME$@calendar
    end if

    TIME&TIME = TIME               ; create coordinate variable 
;;YMDH = cd_calendar(TIME,-3)
;;print(YMDH)

    delete(work)                   ; no longer needed

    nt   = -1
    nmos = 12

    if (rankx.eq.1) then
        xStat  = new ( (/NTIM/), typeof(x), getVarFillValue(x))

        do yyyy=min(year),max(year)
          do mm=min(month),max(month)
             dymon = days_in_month ( yyyy, mm) 
            do dd=1,dymon          ; min(day),max(day)
              do hh=0,23
                 it := ind(yyyy.eq.year .and. mm.eq.month .and. dd.eq.day .and. hh.eq.hour) 
                 nit = dimsizes(it)
                 if (.not.ismissing(it(0)) .and. nit.ge.nval_crit) then
                     nt = nt+1
                     if (arith.eq."avg" .or. arith.eq."ave") then
                         if (nit.gt.1) then ; work around for NCL's dimension reduction
                             xStat(nt)  = dim_avg_n( x(it), nDim)
                         else               ; one 'it'
                             xStat(nt)  = x(it(0))
                         end if
                     end if
                     if (arith.eq."sum") then
                         if (nit.gt.1) then ; work around for NCL's dimension reduction
                             xStat(nt)  = dim_sum_n( x(it), nDim)
                         else               ; one 'it'
                             xStat(nt)  = x(it(0))
                         end if
                     end if
                     if (arith.eq."min") then
                         if (nit.gt.1) then ; work around for NCL's dimension reduction
                             xStat(nt)  = dim_min_n( x(it), nDim)
                         else               ; one 'it'
                             xStat(nt)  = x(it(0))
                         end if
                     end if
                     if (arith.eq."max") then
                         if (nit.gt.1) then ; work around for NCL's dimension reduction
                             xStat(nt)  = dim_max_n( x(it), nDim)
                         else               ; one 'it'
                             xStat(nt)  = x(it(0))
                         end if
                     end if
                           
                     if (nval_crit.gt.1) then  ; make sure # of valid values => nit_valid
                         nit_valid  = dim_num_n(.not.ismissing(x(it)),0)
                         xStat(nt)  = where( nit_valid.ge.nval_crit,  xStat(nt),  xStat@_FillValue)
                     end if  ; nval_crit
  
                 end if  ; .not.ismissing(it(0)) .and. nit.ge.nval_crit
              end do   ; hh
            end do     ; dd
          end do       ; mm
        end do         ; yyyy

        copy_VarAtts(x, xStat)
    end if             ; rankx.eq.1 

;======
; All higher dimension stuff deleted
;=====

    xStat!nDim    = TNAME 
    xStat&$TNAME$ = TIME

    xStat@NCL_tag = "calculate_hourly_values: arith="+arith
    return(xStat)
end