;----------------------------------------------------------
undef("calculate_monthly_values")
function calculate_monthly_values (x:numeric, arith:string, nDim[1]:integer, opt[1]:logical)
;
; calculate monthly values [avg, sum, min, max]
; x: numeric array of 5D or less [eg: time,lev,lat,lon]
;    *must* have time coordinate recognized by cd_calendar
;    if 5D [case,time,lev,lat,lon]
; arith: "avg" [also, "ave"], "sum","min","max"     others may be added later
; nDim : scalar integer that specifies the 'time' dimension [generally, 0]
; opt  : only used to eliminate a warning message
;        opt= 0 ... default ... print Warning
;        opt=-1 ... do not print Warning
;
; Sample usage:  x(time,lat,lon) where time are n-hrly or daily values.
;        xMonthAvg = calculate_monthly_values(x, "avg", 0, False)
;        xMonthSum = calculate_monthly_values(x, "sum", 0, False)
;        xMonthMin = calculate_monthly_values(x, "min", 0, False)
;        xMonthMax = calculate_monthly_values(x, "max", 0, False)
;        xMonthVar = calculate_monthly_values(x, "var", 0, False)
;        xMonthStd = calculate_monthly_values(x, "std", 0, False)
;                  ************
; 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 monmean foo_hourly_or_daily.nc  foo_monthly_mean.nc
;      cdo monmin  foo_hourly_or_daily.nc  foo_monthly_min.nc
;      cdo monmax  foo_hourly_or_daily.nc  foo_monthly_max.nc
;      cdo monsum  foo_hourly_or_daily.nc  foo_monthly_sum.nc
;                  ************
;
local dnamx, TNAME, TIME, dimx, rankx, NTIM, utc_date, year, month, ntim, work
     ,yyyy, mm, nt, it, nmos, xStat, nval_crit, yrStrt, yrLast
begin

    if (nDim.lt.0 .or. nDim.gt.1) then
        print("calculate_monthly_values: unsupported or illegal: nDim="+nDim)
        exit
    end if
 
    dimx     = dimsizes(x)
    rankx    = dimsizes(dimx)                      ; # of dimensions
    if (rankx.gt.5) then
        print("calculate_monthly_values does not support arrays of rank > 5: rank="+rankx)
        exit
    end if

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

; a monthly mean will be calculated only if the # of values is >= nval_crit

    nval_crit = 1                                  ; default
    if (opt .and. isatt(opt,"nval_crit")) then
        nval_crit = toint(opt@nval_crit)
    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
;   'cd_calendar' recognizes any 'calendar' attribute associated with ' x&$TNAME$ '.
;    Hence, the elements below (yyyy, month, ..) reflect that calendar.
;    No calendar attributes 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 monthly coordinate variable. Size must be determined because
;    each month *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$), getVarFillValue(x&$TNAME$)) 

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

    do yyyy=yrStrt,yrLast
      do mm=min(month),max(month)
         it := ind(yyyy.eq.year .and. mm.eq.month)
         if (.not.all(ismissing(it))) then
             nt   = nt+1
             work(nt) = (/ x&$TNAME$(it(0)) /)     ; use 1st encountered time  
             nit  = dimsizes(it)                   ; for print only
           ;;print("yyyy="+yyyy+"  mm="+mm +"  nit="+nit+"  it(0)="+it(0)+" nt="+nt)
         end if    ; .not.all(ismissing(it))
      end do       ; mm
    end do         ; yyyy

    if (nt.lt.0) then
        print("WARNING: calculate_monthly_values: no values for period: "+ \
               yrStrt+"-"+yrLast)
        xStat      = getVarFillValue(x)
        xStat!0    = "time"
        xStat&time = getVarFillValue(x)
        xStat@long_name = "All Missing Values"
        return(xStat)
    end if
    NTIM = nt+1           ; number of months

    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 month"
        TIME@long_name = "time corresponding to 1st encountered time of current month"
    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), getVarFillValue(x))

        do yyyy=yrStrt,yrLast
          do mm=min(month),max(month)
             it  := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.all(ismissing(it))) then
                 nit  = dimsizes(it)
                 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 (arith.eq."var") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt)  = dim_variance_n( x(it), nDim)
                     else               ; one 'it'; variance set to zero
                         xStat(nt)  = 0.0       
                     end if
                 end if
                 if (arith.eq."std") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt)  = dim_stddev_n( x(it), nDim)
                     else               ; one 'it'; std. dev set to zero
                         xStat(nt)  = 0.0       
                     end if
                 end if
                         
                 if (nval_crit.gt.1) then
                     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.all(ismissing(it))
          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), getVarFillValue(x))

        do yyyy=yrStrt,yrLast
          do mm=min(month),max(month)
             it  := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.all(ismissing(it))) then
                 nit = dimsizes(it)
                 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 (arith.eq."var") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt,:)  = dim_variance_n( x(it,:), nDim)
                     else               ; one 'it'; variance set to zero
                         xStat(nt,:)  = 0.0       
                     end if
                 end if
                 if (arith.eq."std") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt,:)  = dim_stddev_n( x(it,:), nDim)
                     else               ; one 'it'; std. dev set to zero
                         xStat(nt,:)  = 0.0       
                     end if
                 end if
                         
                 if (nval_crit.gt.1) then
                     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.all(ismissing(it))
          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), getVarFillValue(x))

        do yyyy=yrStrt,yrLast
          do mm=min(month),max(month)
             it  := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.all(ismissing(it))) then
                 nt  = nt+1
                 nit = dimsizes(it)
                 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 (arith.eq."var") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt,:,:)  = dim_variance_n( x(it,:,:), nDim)
                     else               ; one 'it'; variance set to zero
                         xStat(nt,:,:)  = 0.0       
                     end if
                 end if
                 if (arith.eq."std") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt,:,:)  = dim_stddev_n( x(it,:,:), nDim)
                     else               ; one 'it'; std. dev set to zero
                         xStat(nt,:,:)  = 0.0       
                     end if
                 end if
                         
                 if (nval_crit.gt.1) then
                     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.all(ismissing(it))
          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), getVarFillValue(x))

        do yyyy=yrStrt,yrLast
          do mm=min(month),max(month)
             it  := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.all(ismissing(it))) then
                 nit = dimsizes(it)
                 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 (arith.eq."var") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt,:,:,:)  = dim_variance_n( x(it,:,:,:), nDim)
                     else               ; one 'it'; variance set to zero
                         xStat(nt,:,:,:)  = 0.0       
                     end if
                 end if
                 if (arith.eq."std") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt,:,:,:)  = dim_stddev_n( x(it,:,:,:), nDim)
                     else               ; one 'it'; std. dev set to zero
                         xStat(nt,:,:,:)  = 0.0       
                     end if
                 end if
                         
                 if (nval_crit.gt.1) then
                     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.all(ismissing(it))
          end do       ; mm
        end do         ; yyyy

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

    if (rankx.eq.5 .and. nDim.eq.1) then      ; special case (nmember,time,lev,lat,lon)
        xStat  = new ( (/dimx(0),NTIM,dimx(2),dimx(3),dimx(4)/), typeof(x), getVarFillValue(x))

        do yyyy=yrStrt,yrLast
          do mm=min(month),max(month)
             it  := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.all(ismissing(it))) then
                 nit  = dimsizes(it)
                 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 (arith.eq."var") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(:,nt,:,:,:)  = dim_variance_n( x(:,it,:,:,:), nDim)
                     else               ; one 'it'; set variance to zero
                         xStat(:,nt,:,:,:)  = 0.0
                     end if
                 end if
                 if (arith.eq."std") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(:,nt,:,:,:)  = dim_stddev_n( x(:,it,:,:,:), nDim)
                     else               ; one 'it'; set set. dev. to zero
                         xStat(:,nt,:,:,:)  = 0.0
                     end if
                 end if

                 xStat!0 = "case"
                         
                 if (nval_crit.gt.1) then
                     nit_valid         = dim_num_n(.not.ismissing(x(:,it,:,:,:)),nDim)
                     xStat(:,nt,:,:,:) = where( nit_valid.ge.nval_crit,  xStat(:,nt,:,:,:),  xStat@_FillValue)
                 end if  ; nval_crit

                 copy_VarMeta(x(:,0,:,:,:), xStat(:,0,:,:,:))
                 xStat!1 = "case"
             end if    ; .not.all(ismissing(it))
          end do       ; mm
        end do         ; yyyy
    end if            ; rankx=5 and nDim=1

    if (rankx.eq.5 .and. nDim.eq.0) then    ; special case (time,?nmember?,lev,lat,lon)
        xStat  = new ( (/NTIM,dimx(1),dimx(2),dimx(3),dimx(4)/), typeof(x), getVarFillValue(x))

        do yyyy=yrStrt,yrLast
          do mm=min(month),max(month)
             it  := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.all(ismissing(it))) then
                 nit  = dimsizes(it)
                 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 (arith.eq."var") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt,:,:,:,:)  = dim_variance_n( x(it,:,:,:,:), nDim)
                     else               ; one 'it'; set variance to zero
                         xStat(nt,:,:,:,:)  = 0.0
                     end if
                 end if
                 if (arith.eq."std") then
                     if (nit.gt.1) then ; work around for NCL's dimension reduction
                         xStat(nt,:,:,:,:)  = dim_stddev_n( x(it,:,:,:,:), nDim)
                     else               ; one 'it'; set set. dev. to zero
                         xStat(nt,:,:,:,:)  = 0.0
                     end if
                 end if
                         
                 if (nval_crit.gt.1) then
                     nit_valid         = dim_num_n(.not.ismissing(x(it,:,:,:,:)),nDim)
                     xStat(nt,:,:,:,:) = where( nit_valid.ge.nval_crit,  xStat(nt,:,:,:,:),  xStat@_FillValue)
                 end if  ; nval_crit

                 copy_VarMeta(x(0,:,:,:,:), xStat(0,:,:,:,:))
                 xStat!1 = "case"
             end if    ; .not.all(ismissing(it))
          end do       ; mm
        end do         ; yyyy
     end if            ; rankx=5 and nDim=1

    xStat!nDim    = TNAME 
    xStat&$TNAME$ = TIME
    xStat@NCL_tag = "calculate_monthly_values: arith="+arith
    return(xStat)
end
;----------------------------------------------------------
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, yrStrt, yrLast
begin

    if (nDim.lt.0 .or. nDim.gt.1) then
        print("calculate_daily_values: unsupported or illegal: nDim="+nDim)
        exit
    end if
 
    dimx     = dimsizes(x)
    rankx    = dimsizes(dimx)                      ; # of dimensions
    if (rankx.gt.5) then
        print("calculate_daily_values does not support arrays of rank > 5: rank="+rankx)
        exit
    end if

    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"
    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)

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

;--- 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$), getVarFillValue(x&$TNAME$)) 
    do yyyy=yrStrt,yrLast
      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), 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)
               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
                         
                   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     ; 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), 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)
               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
                        
                   if (nval_crit.gt.1) then
                       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     ; 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), getVarFillValue(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
                         
                   if (nval_crit.gt.1) then
                       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     ; 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), getVarFillValue(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
                        
                   if (nval_crit.gt.1) then
                       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))
            end do     ; dd
          end do       ; mm
        end do         ; yyyy

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

    if (rankx.eq.5 .and. nDim.eq.1) then     ; (ncase,time,level,lat,lon); (0,1,2,3,4)
        xStat  = new ( (/dimx(0),NTIM,dimx(2),dimx(3),dimx(4)/), typeof(x), getVarFillValue(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
                         
                   if (nval_crit.gt.1) then
                       nit_valid         = dim_num_n(.not.ismissing(x(:,it,:,:,:)),nDim)
                       xStat(:,nt,:,:,:) = where( nit_valid.ge.nval_crit,  xStat(:,nt,:,:,:),  xStat@_FillValue)
                   end if  ; nval_crit

               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 .and. nDim.eq.1 

    if (rankx.eq.5 .and. nDim.eq.0) then      ; (time,ncase,level,lat,lon); (0,1,2,3,4)
        xStat  = new ( (/NTIM,dimx(1),dimx(2),dimx(3),dimx(4)/), typeof(x), getVarFillValue(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
                        
                   if (nval_crit.gt.1) then
                       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))
            end do     ; dd
          end do       ; mm
        end do         ; yyyy

    end if             ; rankx.eq.5 


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

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