;+
; @file_comments
; similar to extrapolate but could to the job in a better way because the
; extrapolated values are smoothed... takes more time than extrapolate.
; extrapolate data where maskinput eq 0 by filling
; step by step the coastline points with the mean value of the 8 neighbourgs.
;
;
; @categories
; Interpolation
;
; @param in {in}{required}{type=2d array}
; data to be extrapolate
;
; @param maskinput {in}{required}{type=2d array or -1}
; a 2D array, the land-sea mask of the output data (1 on ocean, 0 on land)
; put -1 if input data are not masked
;
; @keyword MINVAL {type=scalar}{default=not used}
; to specify a minimum value to the extrapolated values
;
; @keyword MAXVAL {type=scalar}{default=not used}
; to specify a maximum value to the extrapolated values
;
; @keyword GE0 {type=scalar 0 or 1}{default=0}
; put 1 to force the extrapolated values to be larger than 0, same as using minval=0.
;
; @returns
; the extrapolated array with no more masked values
;
; @restrictions
; You cannot specify the number of iterations done in the extrapolation process
;
; @examples
; IDL> a=extrapsmooth(dist(jpi,jpj),tmask[*,*,0],/x_periodic)
; IDL> tvplus, a
; compare to extrapolate result:
; IDL> b=extrapolate(dist(jpi,jpj),tmask[*,*,0],/x_periodic)
; IDL> tvplus, b, window = 1
;
; @history
;  January 2007: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;-
FUNCTION extrapsmooth, in, mskin, x_periodic = x_periodic, MINVAL = minval, MAXVAL = maxval, GE0 = ge0
;
  compile_opt strictarr, strictarrsubs
;
  sz = size(reform(in))
  IF sz[0] NE 2 THEN BEGIN
    print, 'Input arrays must have 2 dimensions'
    return, -1
  ENDIF
  nx = sz[1]
  ny = sz[2]
  IF n_elements(mskin) EQ 1 AND mskin[0] EQ -1 THEN mskin = replicate(1b, nx, ny)
  IF n_elements(mskin) NE nx*ny THEN BEGIN
    print, 'input grid mask do not have the good size'
    return, -1
  ENDIF
;
  out = reform(in)
  whmsk = where(mskin EQ 0, nbr)
  IF nbr NE 0 THEN out[temporary(whmsk)] = !values.f_nan
; add values on each side of the array to avoid boundary effects
  nx2 = nx/2
  ny2 = ny/2
  add = replicate(!values.f_nan, nx, ny2)
  out = [[add], [temporary(out)], [add]]
  IF keyword_set(x_periodic)  THEN BEGIN
    add1 = out[0:nx2, *]
    add2 = out[nx-nx2:nx-1, *]
    out = [add2, [temporary(out)], add1]
  ENDIF ELSE BEGIN
    add = replicate(!values.f_nan, nx2, ny+2*ny2)
    out = [add, [temporary(out)], add]
  ENDELSE
;
  msk0 = where(finite(out) EQ 0)
  nnan = total(finite(out, /nan))
  i = 1
  WHILE nnan NE 0 DO BEGIN
    tmp = smooth(out, 2*i + 1, /nan)
    ; find only the changed values that where on land
    new = where(finite(out) EQ 0 AND finite(tmp) EQ 1, nnew)
    IF nnew EQ 0 then nnan = 0 ELSE BEGIN
      IF keyword_set(ge0) THEN tmp = 0. > temporary(tmp)
      IF n_elements(minval) NE 0 THEN tmp = minval > temporary(tmp)
      IF n_elements(maxval) NE 0 THEN tmp = temporary(tmp) < maxval
      out[new] = tmp[new]
      i = i+1
      nnan = total(finite(out, /nan))
    ENDELSE
  ENDWHILE
; a final smooth is needed
  out = smooth(temporary(out), 5, /nan)
  IF keyword_set(ge0) THEN out = 0. > temporary(out)
  IF n_elements(minval) NE 0 THEN out = minval > temporary(out)
  IF n_elements(maxval) NE 0 THEN out = temporary(out) < maxval
; get back to the original domain
  out = (temporary(out))[nx2:nx+nx2-1, ny2:ny+ny2-1]
; put back the non-maskqed values
  whmsk = where(mskin NE 0)
  out[whmsk] = in[whmsk]
;
  return, out
END