;---------------------------------------------------------------------
undef("calculate_segment_values")
function calculate_segment_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'.  :-(
;                  ************
local dnamx, TIME, TNAME, dimx, rankx, NTIM, utc_date, year, month, day \
     ,ntim, work, yyyy, mm, dd, dymon, it, nt, nmos, xSta, nval_crit, lseg, nseg
begin

    if (nDim.ne.0) then
        print("calculate_segment_values: only nDim=0 allowed: nDim="+nDim)
        exit
    end if
 
    dimx     = dimsizes(x)
    rankx    = dimsizes(dimx)                      ; # of dimensions
    if (rankx.gt.4) then
        print("calculate_segment_values does not support arrays of rank > 4: 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_segment_values: unrecognizezed 'arith' argument="+arith)
        exit
    end if

; define segment length in terms of 'day length' ... default is 7

    lseg = 7                                       ; default [7 days]
    if (opt .and. isatt(opt,"segment_length")) then
        lseg = toint(opt@segment_length)
    end if
    lseg1 = lseg-1

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

    nval_crit = 1                                  ; default [1 day]
    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"
    TIME     = x&$TNAME$                           ; extract for convenience
    ntim     = dimx(nDim)                          ; size of input time dimension

; cd_convert has bug recognizing calendar attribute
; This function currently works with the 'standard' and 'gregorian' calendars. 
; No calendar attribute defaults to 'gregorian.'

    if (isatt(TIME,"calendar") .and. \
        .not.(TIME@calendar.eq."standard" .or. TIME@calendar.eq."gregorian")) then
        print("calculate_segment_values: does not work with calendar="+TIME@calendar)
        exit
    end if 

;--- Create current yyyy/mm/dd

    utc_date = cd_calendar(TIME, 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)

;--- Trick: create a 'time' variable in units of 'days since' the start of the input array.

    timeDay  = cd_convert( TIME, "days since "+year(0)+"-"+month(0)+"-"+day(0)+" 00:00" )
    timeDay@info = "indicates the start time of the current segment"
  ;;print(timeDay)  ; for (non-standard; non-gregorian) this will have episodic 1 day gaps

    nday_span= round((timeDay(ntim-1)-timeDay(0)), 3)  ; last_day - start_day; start_day=0  
    nseg     = nday_span/lseg                          ; eg: lseg=7;  nseg= # weeks

    nt       = -1
    ndyStrt  = 0
    ndyLast  = lseg1
    timeSeg  = new( nseg, typeof(timeDay)) 

   ;print(year+"  "+month+"  "+day+"  "+timeDay)
   ;print("nday_span="+nday_span+"; nseg="+nseg)
   ;print("rankx="+rankx+"  ndyStrt="+ndyStrt+"  ndyLast="+ndyLast)

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

        do ns=0,nseg-1
           it := ind(timeDay.ge.ndyStrt .and. timeDay.le.ndyLast)
           if (.not.all(ismissing(it))) then
               nit = dimsizes(it)
               nt  = nt+1

               if (.not.ismissing(it(0))) then
                   timeSeg(nt) = timeDay(it(0))
    
                   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 if                     ; .not.all(ismissing(it))

           ndyStrt = ndyLast+1          ; update segment bounds
           ndyLast = ndyStrt+lseg1
        end do         ; ns (segment loop)         

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

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

        do ns=0,nseg-1
           it  := ind(timeDay.ge.ndyStrt .and. timeDay.le.ndyLast)
           if (.not.all(ismissing(it))) then
               nit  = dimsizes(it)
               nt   = nt+1

               if (.not.ismissing(it(0))) then
                   timeSeg(nt) = timeDay(it(0))
    
                   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 if                         ; .not.all(ismissing(it))

           ndyStrt = ndyLast+1            ; update segment bounds
           ndyLast = ndyStrt+lseg1
        end do                            ; ns (segment loop)         

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

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

        do ns=0,nseg-1
           it  := ind(timeDay.ge.ndyStrt .and. timeDay.le.ndyLast)
           if (.not.all(ismissing(it))) then
                nit  = dimsizes(it)
                nt   = nt+1

               if (.not.ismissing(it(0))) then
                   timeSeg(nt) = timeDay(it(0))
 
                   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 if                       ; .not.all(ismissing(it))

           ndyStrt = ndyLast+1          ; update segment bounds
           ndyLast = ndyStrt+lseg1
        end do         ; ns (segment loop)         

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

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

        do ns=0,nseg-1
           it  := ind(timeDay.ge.ndyStrt .and. timeDay.le.ndyLast)
           if (.not.all(ismissing(it))) then
                nit  = dimsizes(it)
                nt   = nt+1

                if (.not.ismissing(it(0))) then
                    timeSeg(nt) = timeDay(it(0))
     
                    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 if                       ; .not.all(ismissing(it))

           ndyStrt = ndyLast+1          ; update segment bounds
           ndyLast = ndyStrt+lseg1
        end do         ; ns (segment loop)         

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

    xStat!nDim    = TNAME 
    xStat&$TNAME$ = timeSeg

    xStat@NCL_tag = "calculate_segment_values: lseg="+lseg+" days: arith="+arith
    return(xStat)
end