pro SUB_accumulate_step5_tims, $
    var_combine_all, var_mean_all, var_std_all, $
    var_combine,     var_mean,     var_std, $
    nobs_all=var_nobs_all, nobs_this=nobs_this, $
    time_all = var_time_all, time_this = time_this

if n_elements(var_combine) eq 0 then begin ; if var_combine is empty, then insert dummy column; required by json format

  if n_elements(var_mean_all) gt 0 then begin 
    size_this = size(var_mean_all,/dim)
    var_mean_all= [var_mean_all, fltarr([1,size_this[1:-1]])+!values.f_nan]
    var_std_all= [var_std_all, fltarr([1,size_this[1:-1]])+!values.f_nan]
    if ARG_PRESENT(nobs_this) then begin
      var_nobs_all = [var_nobs_all, fltarr([1,size_this[1:-1]])+!values.f_nan]
    endif
    if ARG_PRESENT(time_this) then begin
      var_time_all= [var_time_all, time_this]
    endif
  endif

endif else begin

  if n_elements(var_mean_all) eq 0 then begin ; the first element, simply give the current value

    var_combine_all = var_combine
    var_mean_all =  reform(var_mean, [1, size(var_mean, /dim)])
    var_std_all  =  reform(var_std,  [1, size(var_std, /dim)])
    if ARG_PRESENT(nobs_this) then begin
      var_nobs_all = reform(nobs_this,  [1, size(nobs_this, /dim)])
    endif
    if ARG_PRESENT(time_this) then begin
      var_time_all = time_this
    endif

    ;print,'var_combine_all:', var_combine_all
    ;print,'var_nobs_all:', var_nobs_all

  endif else begin

    existing = 0
    ;print, 'var_combine:',size(var_combine), var_combine
    for combine_loop = 0, n_elements(var_combine)-1 do begin

      this_id = where(var_combine_all eq var_combine[combine_loop], this_num)

      if this_num eq 0 then begin ; new combine, rebuild the whole line for the new combine

;        print,'new combine'

        var_combine_all = [var_combine_all, var_combine[combine_loop]]
        size_this = size(var_mean_all,/dim)

        ; [time, n_combine, n_level, n_lat_list]
        ; switch dim 0 (time) and dim 1 (combine), to add new combine
        var_mean_all= [transpose(var_mean_all, [1,0,2,3]), fltarr([1, size_this[0], size_this[2:-1]])+!values.f_nan]
        var_mean_all = transpose(var_mean_all, [1,0,2,3])

        var_std_all= [transpose(var_std_all, [1,0,2,3]), fltarr([1, size_this[0], size_this[2:-1]])+!values.f_nan]
        var_std_all = transpose(var_std_all, [1,0,2,3])

        if ARG_PRESENT(nobs_this) then begin
          var_nobs_all= [transpose(var_nobs_all, [1,0,2,3]), fltarr([1, size_this[0], size_this[2:-1]])+!values.f_nan]
          var_nobs_all = transpose(var_nobs_all, [1,0,2,3])
        endif

        ;; if the new combine is the first in current combine list, then we need to additionally append column (along time) 
        ;; if not the first, not need 
        if combine_loop eq 0 then begin ; only append once
          size_this = size(var_mean_all,/dim)
          var_mean_all= [var_mean_all, fltarr([1, size_this[1], size_this[2:-1]])+!values.f_nan]
          var_std_all= [var_std_all, fltarr([1, size_this[1], size_this[2:-1]])+!values.f_nan]
          if ARG_PRESENT(nobs_this) then begin
            var_nobs_all= [var_nobs_all, fltarr([1, size_this[1], size_this[2:-1]])+!values.f_nan]
          endif
          if ARG_PRESENT(time_this) then begin
            var_time_all= [var_time_all, time_this]
          endif
          existing=1
        endif 

        var_mean_all[-1, -1, *, *] = var_mean[combine_loop, *, *]
        var_std_all[-1, -1, *, *] = var_std[combine_loop, *, *]
        if ARG_PRESENT(nobs_this) then begin
          var_nobs_all[-1, -1, *, *] = nobs_this[combine_loop, *, *]
        endif

        ;stop

      endif else begin  ; existing combine, append to column

;        print,'existing combine', existing
;        print, 'combine_loop',combine_loop

        if existing eq 0 then begin ; only append once
;        save, /variable, filename='test.sav'
;        stop

          size_this = size(var_mean_all,/dim)

;    print,'size(var_mean_all,/dim) before:', size(var_mean_all,/dim)
;    print, 'size(var_mean,/dim) before:',size(var_mean,/dim)

          var_mean_all= [var_mean_all, fltarr([1, size_this[1], size_this[2:-1]])+!values.f_nan]
          var_std_all= [var_std_all, fltarr([1,size_this[1], size_this[2:-1]])+!values.f_nan]
          if ARG_PRESENT(nobs_this) then begin
            var_nobs_all= [var_nobs_all, fltarr([1,size_this[1], size_this[2:-1]])+!values.f_nan]
          endif
          if ARG_PRESENT(time_this) then begin
            var_time_all= [var_time_all, time_this]
          endif
          existing = 1
        endif 

        var_mean_all[-1, this_id,*, *] = var_mean[combine_loop, *, *]
        var_std_all[-1, this_id,*, *] = var_std[combine_loop, *, *]

        if ARG_PRESENT(nobs_this) then begin
          var_nobs_all[-1, this_id,*, *] = nobs_this[combine_loop, *, *]
        endif

;    print,'size(var_mean_all,/dim) after:', size(var_mean_all,/dim)
;    print, 'size(var_mean,/dim) after:',size(var_mean,/dim)

      endelse
      ;print,'-----------'

    endfor ; combine_loop

  endelse  ; n_elements(var_mean_all) eq 0 
endelse ; n_elements(var_combine) eq 0

end