;+ NAME: 
; xdeflag 
;PURPOSE:
; Replaces FLAGs in arrays with interpolated or other values
;CALLING SEQUENCE:
; xdeflag, method, t, y, flag=flag, _extra=_extra
;INPUT:
; method = set to "repeat", this will repeat the last good value.
;          set to "linear", then linear interpolation is used, but for
;          the edges, the closest value is used, there is no
;          extrapolation
; t = time array, in any useable tplot format
; y = the input array, n_elements(t) by n
;OUTPUT:
; y = either interpolated or repated, where the value is > 0.98*flag,
;     or NaN
;KEYWORDS:
; flag = the value that flagged data are set to, the default is
;           6.8792e28
; maxgap = the maximum number of rows that can be filled? the default
;           is n_elements(t)
;HISTORY:
; 2-feb-2007, jmm, jimm.ssl.berkeley.edu from Vassilis' clip_deflag.pro
;
;$LastChangedBy$
;$LastChangedDate$
;$LastChangedRevision$
;$URL$
;-
Pro xdeflag, method0, t0, y, flag=flag, maxgap = maxgap, _extra=_extra

  method = strtrim(strlowcase(method0), 2)
  If(keyword_set(flag)) Then big = flag Else big = 6.879e28
  big98 = 0.98*big
  t = time_double(t0)
  nrows = n_elements(t)
  if (keyword_set(maxgap)) then mxgp = maxgap $
  else mxgp = nrows
;
  If(n_elements(y[*, 0]) ne nrows) Then Begin
    message, /info, 'number of rows does not agree between time and yarray(s)'
    help, t0
    help, y
    Return
  Endif
  ndims = size(y, /dimensions)
  ncols = 0 & nlayers = 0
  if (n_elements(ndims) ge 2) then ncols = ndims[1]
  if (n_elements(ndims) eq 3) then nlayers = ndims[2]
  if (ncols eq 0) then nycols = 1 else nycols = ncols
  if (nlayers eq 0) then nylayers = 1 else nylayers = nlayers
  nycolayers = nycols*nylayers
  y = reform(y, nrows, nycolayers)
  for j=0,nycolayers-1 do begin
    jiwhere = where((y(*, j) gt BIG98) Or (finite(y[*, j]) Eq 0), jiany)
    if ((jiany gt 0) and (jiany lt nrows))  then begin
      if (jiany eq 1) then begin
        ngaps = 1l
        kbegin = lonarr(1)
        kbegin(0) = long(jiwhere(0))
        kend = kbegin
        ksize = kend-kbegin+1l
        kslope = dblarr(1)
        tk0 = dblarr(1)
        yk0 = dblarr(1)
        if (kbegin(0) eq 0) then begin
          tk0(0) = 0l
          yk0(0) = y(kend(0)+1, j)
          kslope(0) = 0l
        endif
        if (kend(0) eq (nrows-1)) then begin
          tk0(0) = t(kbegin(0)-1)
          yk0(0) = y(kbegin(0)-1, j)
          kslope(0) = 0l
        endif
        if ((kbegin(0) ne 0) and (kend(0) ne (nrows-1))) then begin
          tk0(0) = t(kbegin(0)-1)
          yk0(0) = y(kbegin(0)-1, j)
          kslope(0) = (y(kend(0)+1, j)-y(kbegin(0)-1, j))/(t(kend(0)+1)-tk0(0))
        endif
      endif else begin
        kany = make_array(jiany-1, /long, /index)
        difji = jiwhere(kany+1)-jiwhere(kany)
        kdif = where(difji gt 1, knum)
        ngaps = knum+1
        kbegin = lonarr(ngaps)
        kend = lonarr(ngaps)
        kbegin(0) = jiwhere(0)
        kend(ngaps-1) = jiwhere(jiany-1)
        if (ngaps gt 1) then begin
          kindex1 = make_array(ngaps-1, /long, /index)
          kend(kindex1) = jiwhere(kdif(kindex1))
          kbegin(kindex1+1) = jiwhere(kdif(kindex1)+1)
        endif
        ksize = kend-kbegin+1
        kslope = dblarr(ngaps)
        tk0 = dblarr(ngaps)
        yk0 = dblarr(ngaps)
        if (kbegin(0) eq 0) then begin
          tk0(0) = 0l
          yk0(0) = y(kend(0)+1, j)
          kslope(0) = 0l
        endif else begin
          if (kend(0) ne (nrows-1)) then begin
            tk0(0) = t(kbegin(0)-1)
            yk0(0) = y(kbegin(0)-1, j)
            kslope(0) = (y(kend(0)+1, j)-y(kbegin(0)-1, j))/(t(kend(0)+1)-tk0(0))
          endif
        endelse
        if (kend(ngaps-1) eq (nrows-1)) then begin
          tk0(ngaps-1) = t(kbegin(ngaps-1)-1)
          yk0(ngaps-1) = y(kbegin(ngaps-1)-1, j)
          kslope(ngaps-1) = 0l
        endif else begin
          if (kbegin(ngaps-1) ne 0) then begin
            tk0(ngaps-1) = t(kbegin(ngaps-1)-1)
            yk0(ngaps-1) = y(kbegin(ngaps-1)-1, j)
            kslope(ngaps-1) = (y(kend(ngaps-1)+1, j)-y(kbegin(ngaps-1)-1, j))/(t(kend(ngaps-1)+1)-tk0(ngaps-1))
          endif
        endelse
        if (ngaps gt 2) then begin
          kindex2 = make_array(ngaps-2, /long, /index)
          tk0(kindex2+1) = t(kbegin(kindex2+1)-1)
          yk0(kindex2+1) = y(kbegin(kindex2+1)-1, j)
          kslope(kindex2+1) = (y(kend(kindex2+1)+1, j)-y(kbegin(kindex2+1)-1, j))/(t(kend(kindex2+1)+1)-tk0(kindex2+1))
        endif
      endelse
      case method of
        "repeat":begin
          repeat_loop:
          kthgood = where(ksize le mxgp, ianygood)
          if (ianygood gt 0) then begin
;                print,'ianygood=',ianygood
            for kthgap = 0l, long(ianygood-1) do begin
              kindices = kbegin(kthgood(kthgap))+$
                indgen(ksize(kthgood(kthgap)))
              y(kindices, j) = yk0(kthgood(kthgap))
            endfor
          endif
        end
        "linear":begin
          kthgood = where(ksize le mxgp, ianygood)
;                print,'ianygood=',ianygood
          if (ianygood gt 0) then begin
            for kthgap = 0l, long(ianygood-1) do begin
              kindices = kbegin(kthgood(kthgap))+$
                indgen(ksize(kthgood(kthgap)))
              y(kindices, j) = yk0(kthgood(kthgap))+$
                (t(kindices)-tk0(kthgood(kthgap)))*kslope(kthgood(kthgap))
            endfor
          endif
        end
        else:Begin
          message, /info, 'Invalid Method input, Set to ''repeat'''
          goto, repeat_loop
        end
      endcase
    endif
    if (jiany eq nrows) then print,'One row is all FLAGs: left as is by deflag'
    if ((ncols gt 0) and (nlayers gt 0)) then y=reform(y,nrows,ncols,nlayers)
    if ((ncols gt 0) and (nlayers eq 0)) then y=reform(y,nrows,ncols)
    if ((ncols eq 0) and (nlayers eq 0)) then y=reform(y,nrows)
  endfor


  Return
End