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)
;
; It is assumed that there will be multiple elements for the dim_avg_n
; calculation.
;
local dnamx, TNAME, TIME, dimx, rankx, NTIM, utc_date, year, month, ntim, work
     ,yyyy, mm, nt, it, nmos, xStat
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_monthly_values: unrecognizezed 'arith' argument="+arith)
        exit
    end if

    dnamx    = getvardims(x)                       ; dimension names
    TNAME    = dnamx(nDim)                         ; typically TNAME="time"
 
    dimx     = dimsizes(x)                         ; dimensions sizes
    rankx    = dimsizes(dimx)                      ; # of dimensions
    ntim     = dimx(nDim)                          ; size of input time dimension

;--- Create current yyyy/mm

    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$), getFillValue(x&$TNAME$)) 
    do yyyy=min(year),max(year)
      do mm=min(month),max(month)
         it := ind(yyyy.eq.year .and. mm.eq.month) 
         if (.not.ismissing(it(0))) then
             nt = nt+1
             work(nt) = (/ x&$TNAME$(it(0)) /)     ; use 1st time of current month  
         end if
      end do       ; mm
    end do         ; yyyy
    NTIM = nt+1

    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"
    end if
    if (isatt(x&$TNAME$,"units")) then
        TIME@units     = x&$TNAME$@units 
    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)
             it := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.ismissing(it(0))) then
                 nt = nt+1
                 if (arith.eq."avg" .or. arith.eq."ave") then
                     xStat(nt) = dim_avg_n( x(it), nDim)
                 end if
                 if (arith.eq."sum") then
                     xStat(nt) = dim_sum_n( x(it), nDim)
                 end if
                 if (arith.eq."min") then
                     xStat(nt) = dim_min_n( x(it), nDim)
                 end if
                 if (arith.eq."max") then
                     xStat(nt) = dim_max_n( x(it), nDim)
                 end if
             end if    ; .not.ismissing(it(0))
          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)
             it := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.ismissing(it(0))) then
                 nt = nt+1
                 if (arith.eq."avg" .or. arith.eq."ave") then
                     xStat(nt,:) = dim_avg_n( x(it,:), nDim)
                 end if
                 if (arith.eq."sum") then
                     xStat(nt,:) = dim_sum_n( x(it,:), nDim)
                 end if
                 if (arith.eq."min") then
                     xStat(nt,:) = dim_min_n( x(it,:), nDim)
                 end if
                 if (arith.eq."max") then
                     xStat(nt,:) = dim_max_n( x(it,:), nDim)
                 end if
             end if  ; .not.ismissing(it(0))
          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)
             it := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.ismissing(it(0))) then
                 nt = nt+1
                 if (arith.eq."avg" .or. arith.eq."ave") then
                     xStat(nt,:,:) = dim_avg_n( x(it,:,:), nDim)
                 end if
                 if (arith.eq."sum") then
                     xStat(nt,:,:) = dim_sum_n( x(it,:,:), nDim)
                 end if
                 if (arith.eq."min") then
                     xStat(nt,:,:) = dim_min_n( x(it,:,:), nDim)
                 end if
                 if (arith.eq."max") then
                     xStat(nt,:,:) = dim_max_n( x(it,:,:), nDim)
                 end if
             end if    ; .not.ismissing(it(0))
          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)
             it := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.ismissing(it(0))) then
                 nt = nt+1
                 if (arith.eq."avg" .or. arith.eq."ave") then
                     xStat(nt,:,:,:) = dim_avg_n( x(it,:,:,:), nDim)
                 end if
                 if (arith.eq."sum") then
                     xStat(nt,:,:,:) = dim_sum_n( x(it,:,:,:), nDim)
                 end if
                 if (arith.eq."min") then
                     xStat(nt,:,:,:) = dim_min_n( x(it,:,:,:), nDim)
                 end if
                 if (arith.eq."max") then
                     xStat(nt,:,:,:) = dim_max_n( x(it,:,:,:), nDim)
                 end if
             end if  ; .not.ismissing(it(0))
          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)
             it := ind(yyyy.eq.year .and. mm.eq.month) 
             if (.not.ismissing(it(0))) then
                 nt = nt+1
                 if (arith.eq."avg" .or. arith.eq."ave") then
                     xStat(:,nt,:,:,:) = dim_avg_n( x(:,it,:,:,:), nDim)
                 end if
                 if (arith.eq."sum") then
                     xStat(:,nt,:,:,:) = dim_sum_n( x(:,it,:,:,:), nDim)
                 end if
                 if (arith.eq."min") then
                     xStat(:,nt,:,:,:) = dim_min_n( x(:,it,:,:,:), nDim)
                 end if
                 if (arith.eq."max") then
                     xStat(:,nt,:,:,:) = dim_max_n( x(:,it,:,:,:), nDim)
                 end if
             end if  ; .not.ismissing(it(0))
          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_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)
local dnamx, TIME, TNAME, dimx, rankx, NTIM, utc_date, year, month, day \
     ,ntim, work, yyyy, mm, dd, dymon, it, nt, nmos, xStat
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

    dnamx    = getvardims(x)                       ; dimension names
    TNAME    = dnamx(nDim)                         ; typically TNAME="time"
 
    dimx     = dimsizes(x)                         ; dimensions sizes
    rankx    = dimsizes(dimx)                      ; # of dimensions
    ntim     = dimx(nDim)                          ; size of input time dimension

;--- Create current yyyy/mm/dd

    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) 
           if (.not.ismissing(it(0))) 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

    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

    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) 
               if (.not.ismissing(it(0))) then
                   nt = nt+1
                   if (arith.eq."avg" .or. arith.eq."ave") then
                       xStat(nt) = dim_avg_n( x(it), nDim)
                   end if
                   if (arith.eq."sum") then
                       xStat(nt) = dim_sum_n( x(it), nDim)
                   end if
                   if (arith.eq."min") then
                       xStat(nt) = dim_min_n( x(it), nDim)
                   end if
                   if (arith.eq."max") then
                       xStat(nt) = dim_max_n( x(it), nDim)
                   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) 
               if (.not.ismissing(it(0))) then
                   nt = nt+1
                   if (arith.eq."avg" .or. arith.eq."ave") then
                       xStat(nt,:) = dim_avg_n( x(it,:), nDim)
                   end if
                   if (arith.eq."sum") then
                       xStat(nt,:) = dim_sum_n( x(it,:), nDim)
                   end if
                   if (arith.eq."min") then
                       xStat(nt,:) = dim_min_n( x(it,:), nDim)
                   end if
                   if (arith.eq."max") then
                       xStat(nt,:) = dim_max_n( x(it,:), nDim)
                   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) 
               if (.not.ismissing(it(0))) then
                   nt = nt+1
                   if (arith.eq."avg" .or. arith.eq."ave") then
                       xStat(nt,:,:) = dim_avg_n( x(it,:,:), nDim)
                   end if
                   if (arith.eq."sum") then
                       xStat(nt,:,:) = dim_sum_n( x(it,:,:), nDim)
                   end if
                   if (arith.eq."min") then
                       xStat(nt,:,:) = dim_min_n( x(it,:,:), nDim)
                   end if
                   if (arith.eq."max") then
                       xStat(nt,:,:) = dim_max_n( x(it,:,:), nDim)
                   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) 
               if (.not.ismissing(it(0))) then
                   nt = nt+1
                   if (arith.eq."avg" .or. arith.eq."ave") then
                       xStat(nt,:,:,:) = dim_avg_n( x(it,:,:,:), nDim)
                   end if
                   if (arith.eq."sum") then
                       xStat(nt,:,:,:) = dim_sum_n( x(it,:,:,:), nDim)
                   end if
                   if (arith.eq."min") then
                       xStat(nt,:,:,:) = dim_min_n( x(it,:,:,:), nDim)
                   end if
                   if (arith.eq."max") then
                       xStat(nt,:,:,:) = dim_max_n( x(it,:,:,:), nDim)
                   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 

    xStat!0       = TNAME 
    xStat&$TNAME$ = TIME

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