;----------------------------------------------------------
undef("calculate_daily_values")
function calculate_daily_values (x:numeric, arith:string, nDim[1]:integer, opt[1]:logical)
;                  ************
; 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'.  :-(
;                  ************
; If this function is slow, consider using the Climate Data Operator (CDO)
;      cdo daymean foo_hourly.nc  foo_daily_mean.nc
;      cdo daymin  foo_hourly.nc  foo_daily_min.nc
;      cdo daymax  foo_hourly.nc  foo_daily_max.nc
;      cdo daysum  foo_hourly.nc  foo_daily_sum.nc
;                  ************
local dnamx, TIME, TNAME, dimx, rankx, NTIM, utc_date, year, month, day \
     ,ntim, work, yyyy, mm, dd, dymon, it, nt, nmos, xSta, nval_crit
begin
    if (.not.(arith.eq."ave" .or. arith.eq."avg" .or. arith.eq."sum"  \
                             .or. arith.eq."min" .or. arith.eq."max") ) then
        print("calculate_daily_values: unrecognizezed 'arith' argument="+arith)
        exit
    end if

; a daily 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"
 
    dimx     = dimsizes(x)
    rankx    = dimsizes(dimx)                      ; # of dimensions
    ntim     = dimx(nDim)                          ; size of input time dimension

;--- Create current yyyy/mm/dd
;   '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)

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

    nt       =  -1                                 ; number of individual days
    work     = new (ntim, typeof(x&$TNAME$), getFillValue(x&$TNAME$)) 
    do yyyy=min(year),max(year)
      do mm=min(month),max(month)
         dymon  = days_in_month ( yyyy, mm) 
        do dd=1,dymon          
           it  := ind(yyyy.eq.year .and. mm.eq.month .and. dd.eq.day) 
           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 day  
           end if
        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 
    delete(work)                   ; no longer needed

    nt   = -1
    nmos = 12

    if (rankx.eq.1) then
        xStat  = new ( (/NTIM/), typeof(x), getFillValue(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)
               it := ind(yyyy.eq.year .and. mm.eq.month .and. dd.eq.day) 
               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
               end if  ; .not.ismissing(it(0))
            end do     ; dd
          end do       ; mm
        end do         ; yyyy

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

    if (rankx.eq.2) then
        xStat  = new ( (/NTIM,dimx(1)/), typeof(x), getFillValue(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)
               it := ind(yyyy.eq.year .and. mm.eq.month .and. dd.eq.day) 
               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
               end if  ; .not.ismissing(it(0))
            end do     ; dd
          end do       ; mm
        end do         ; yyyy

        copy_VarMeta(x(0,:), xStat(0,:))
    end if             ; rankx.eq.2 

    if (rankx.eq.3) then
        xStat  = new ( (/NTIM,dimx(1),dimx(2)/), typeof(x), getFillValue(x))

        do yyyy=min(year),max(year)
          do mm=min(month),max(month)
             dymon = days_in_month ( yyyy, mm) 
            do dd=min(day),max(day)
               it := ind(yyyy.eq.year .and. mm.eq.month .and. dd.eq.day) 
               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
               end if  ; .not.ismissing(it(0))
            end do     ; dd
          end do       ; mm
        end do         ; yyyy

        copy_VarMeta(x(0,:,:), xStat(0,:,:))
    end if             ; rankx.eq.3 

    if (rankx.eq.4) then
        xStat  = new ( (/NTIM,dimx(1),dimx(2),dimx(3)/), typeof(x), getFillValue(x))

        do yyyy=min(year),max(year)
          do mm=min(month),max(month)
             dymon = days_in_month ( yyyy, mm) 
            do dd=min(day),max(day)
               it := ind(yyyy.eq.year .and. mm.eq.month .and. dd.eq.day) 
               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
               end if  ; .not.ismissing(it(0))
            end do     ; dd
          end do       ; mm
        end do         ; yyyy

        copy_VarMeta(x(0,:,:,:), xStat(0,:,:,:))
    end if             ; rankx.eq.4 

    if (rankx.eq.5) then
        xStat  = new ( (/dimx(0),NTIM,dimx(2),dimx(3),dimx(4)/), typeof(x), getFillValue(x))
        do yyyy=min(year),max(year)
          do mm=min(month),max(month)
             dymon = days_in_month ( yyyy, mm) 
            do dd=min(day),max(day)
               it := ind(yyyy.eq.year .and. mm.eq.month .and. dd.eq.day) 
               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
               end if  ; .not.ismissing(it(0))
            end do     ; dd
          end do       ; mm
        end do         ; yyyy

        copy_VarMeta(x(:,0,:,:,:), xStat(:,0,:,:,:))
    end if             ; rankx.eq.5 

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

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