;+
;
; @file_comments
; extrapolate data (zinput) where maskinput equal 0 by filling step by
; step the coastline points with the mean value of the 8 neighbors
; (weighted by their mask values).
;
; @categories
; Interpolation
;
; @param zinput {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
;
; @param nb_iteration {in}{optional}{type=integer}{default=large enough to fill everything}
; Maximum number of iterations done in the extrapolation process. If there
; is no more masked values we exit extrapolate before reaching nb_iteration
;
; @keyword X_PERIODIC {type=scalar, 0 or 1}{default=0}
; put 1 to specify that the data are periodic along x axis
;
; @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.
;
; @keyword 
; _EXTRA to be able to call extrapolate with _extra keyword
;
; @returns
; the extrapolated 2d array
;
; @examples
;
;   IDL> a=extrapolate(dist(jpi,jpj),tmask[*,*,0],/x_periodic)
;   IDL> tvplus, a
;   IDL> tvplus, a*(1-tmask[*,*,0])
;
; get the coastline:
;
;   IDL> a=extrapolate(tmask[*,*,0],tmask[*,*,0],1,/x_periodic)
;   IDL> tvplus, a-tmask[*,*,0]
;
; @history
;  Originaly written by G. Roullet
;  Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
FUNCTION extrapolate, zinput, maskinput, nb_iteration, X_PERIODIC=x_periodic $
                      , MINVAL=minval, MAXVAL=maxval, GE0=ge0, _EXTRA=ex
;
  compile_opt idl2, strictarrsubs
;
; check the number of iteration used in the extrapolation.
  szin = size(zinput)
  IF szin[0] NE 2 THEN return, -1. ELSE szin = szin[1:2]
  nx = szin[0]
  ny = szin[1]
  IF n_elements(nb_iteration) EQ 0 THEN nb_iteration = max(szin)
  IF nb_iteration EQ 0 THEN return, zinput
; take care of the boundary conditions...
;
; for the x direction, we put 2 additional columns at the left and
; right side of the array.
; for the y direction, we put 2 additional lines at the bottom and
; top side of the array.
; These changes allow us to use shift function without taking care of
; the x and y periodicity.
;
  ztmp = bytarr(nx+2, ny+2)
  IF n_elements(maskinput) EQ 1 AND maskinput[0] EQ -1 THEN maskinput = replicate(1b, nx, ny)
  IF n_elements(maskinput) NE nx*ny THEN BEGIN
    ras = report('input grid mask do not have the good size')
    return, -1
  ENDIF
  ztmp[1:nx, 1:ny] = byte(maskinput)
  msk = temporary(ztmp)
;
  ztmp = replicate(1.e20, nx+2, ny+2)
  ztmp[1:nx, 1:ny] = zinput
  if keyword_set(x_periodic) then begin
    ztmp[0, 1:ny] = zinput[nx-1, *]
    ztmp[nx+1, 1:ny] = zinput[0, *]
  ENDIF
; remove NaN points if there is some...
  nan = where(finite(ztmp) EQ 0, cnt_nan)
  IF cnt_nan NE 0 THEN ztmp[temporary(nan)] = 1.e20
  z = temporary(ztmp)
  nx2 = nx+2
  ny2 = ny+2
;---------------------------------------------------------------
;---------------------------------------------------------------
; extrapolation
;---------------------------------------------------------------
  sqrtinv = 1./sqrt(2)
;
  cnt = 1
; When we look for the coastline, we don't want to select the
; borderlines of the array. -> we force the value of the mask for
; those lines.
  msk[0, *] = 1b
  msk[nx+1, *] = 1b
  msk[*, 0] = 1b
  msk[*, ny+1] = 1b
; find the land points
  land = where(msk EQ 0, cnt_land)
;---------------------------------------------------------------
  WHILE cnt LE nb_iteration AND cnt_land NE 0 DO BEGIN
;---------------------------------------------------------------
; find the coastline points...
;---------------------------------------------------------------
; Once the land points list has been found, we change back the
; mask values for the boundary conditions.
    msk[0, *] = 0b
    msk[nx+1, *] = 0b
    msk[*, 0] = 0b
    msk[*, ny+1] = 0b
    if keyword_set(x_periodic) then begin
      msk[0, *] = msk[nx, *]
      msk[nx+1, *] = msk[1, *]
    endif
;
; we compute the weighted number of sea neighbors.
; those 4 neighbors have a weight of 1:
;    *
;   *+*
;    *
;
; those 4 neighbors have a weight of 1/sqrt(2):
;
;    * *
;     +
;    * *
;
; As we make sure that none of the land points are located on the
; border of the array, we can compute the weight without shift
; (faster).
;
    weight = msk[land+1]+msk[land-1]+msk[land+nx2]+msk[land-nx2] $
             +sqrtinv*(msk[land+nx2+1]+msk[land+nx2-1] $
                       +msk[land-nx2+1]+msk[land-nx2-1])
; list all the points that have sea neighbors
    ok = where(weight GT 0)
; the coastline points
    coast = land[ok]
; their weighted number of sea neighbors.
    weight = weight[temporary(ok)]
;---------------------------------------------------------------
; fill the coastline points
;---------------------------------------------------------------
    z = temporary(z)*msk
;
    zcoast = z[1+coast]+z[-1+coast]+z[nx2+coast]+z[-nx2+coast] $
             +1./sqrt(2)*(z[nx2+1+coast]+z[nx2-1+coast] $
                          +z[-nx2+1+coast]+z[-nx2-1+coast])
;
    IF keyword_set(ge0) THEN zcoast = 0. > temporary(zcoast)
    IF n_elements(minval) NE 0 THEN zcoast = minval > temporary(zcoast)
    IF n_elements(maxval) NE 0 THEN zcoast = temporary(zcoast) < maxval
    z[coast] = temporary(zcoast)/temporary(weight)
; we update the boundary conditions of z
    if keyword_set(x_periodic) then begin
      z[0, *] = z[nx, *]
      z[nx+1, *] = z[1, *]
    endif
;---------------------------------------------------------------
; we update the mask
;---------------------------------------------------------------
    msk[temporary(coast)] = 1
;
    cnt = cnt + 1
; When we look for the coast line, we don't want to select the
; borderlines of the array. -> we force the value of the mask for
; those lines.
    msk[0, *] = 1b
    msk[nx+1, *] = 1b
    msk[*, 0] = 1b
    msk[*, ny+1] = 1b
; find the land points
    land = where(msk EQ 0, cnt_land)
;
  ENDWHILE
;---------------------------------------------------------------
; we return the original size of the array
;---------------------------------------------------------------
;
  return, z[1:nx, 1:ny]
END