;+
; NAME: extrap.pro (based on remplit.pro)
;
; PURPOSE: Extrapolates ocean data values on land points
;
; CATEGORY: Subroutine
;
; CALLING SEQUENCE: extrap, zdata, zmask, it, val
;
; INPUTS:
;          zdata : field to extrapolate
;          zmask : field's mask
;          it    : iteration
;          val   : field value on masked points
;
; KEYWORD PARAMETERS: None
;
; OUTPUTS:
;          zdata : extrapolated field
;          zmask : mask after extrapolation
;
; COMMON BLOCKS:
;          common_interp.pro
;
;
; MODIFICATION HISTORY: 19/11/99 Arnaud Jouzeau
;                       25/02/00 Sebastien Masson (remplit.pro)
;-
PRO extrap, z, mmmask, it, val
@common_interp
;
; 1. Looking for land points
; ==========================
;
; IF keyword_set(key_performance) THEN temp = systime(1)
; znbvoisins=shift(zmask,1,0)+shift(zmask,-1,0)+shift(zmask,0,1)+shift(zmask,0,-1) $
;           +shift(zmask,1,1)+shift(zmask,-1,1)+shift(zmask,1,-1)+shift(zmask,-1,-1)
; IF keyword_set(key_performance) THEN print, systime(1)-temp
;
; 2. Extrapolating when needed
; ============================
;
; IF keyword_set(key_performance) THEN temp = systime(1)
; zdfm = zdata*zmask
; zremplit =shift(zdfm,1,0)+shift(zdfm,-1,0)+shift(zdfm,0,1)+shift(zdfm, 0,-1) $
;          +shift(zdfm,1,1)+shift(zdfm,-1,1)+shift(zdfm,1,-1)+shift(zdfm, -1,-1)
; IF keyword_set(key_performance) THEN print, systime(1)-temp
;
;  ... coastal points
;
; IF keyword_set(key_performance) THEN temp = systime(1)
; ocean = where(zmask NE 0.)
; IF keyword_set(key_performance) THEN print, systime(1)-temp
; IF keyword_set(key_performance) THEN temp = systime(1)
; ind = where(znbvoisins NE 0. AND znbvoisins NE 8.)
; IF keyword_set(key_performance) THEN print, systime(1)-temp
; IF keyword_set(key_performance) THEN temp = systime(1)
; znbvoisins(ind) = 1.d/temporary(znbvoisins(ind))
; IF keyword_set(key_performance) THEN print, systime(1)-temp
; IF keyword_set(key_performance) THEN temp = systime(1)
; zremplit(ind) = temporary(zremplit(ind))*znbvoisins(ind)
; IF keyword_set(key_performance) THEN print, systime(1)-temp
; IF keyword_set(key_performance) THEN temp = systime(1)
; zmask(ind) = 1.d
; IF keyword_set(key_performance) THEN print, systime(1)-temp
;
;  ... inland points
;
; IF keyword_set(key_performance) THEN temp = systime(1)
; ocean2 = where(zmask NE 0.)
; IF keyword_set(key_performance) THEN print, systime(1)-temp
; IF keyword_set(key_performance) THEN temp = systime(1)
; IF n_elements(ocean2) NE jpiatm*(jpjatm+4) THEN BEGIN
; land = where(zmask EQ 0.)
; zremplit(land) = val
; ENDIF
; IF keyword_set(key_performance) THEN print, systime(1)-temp
;
;  ... filling with ocean values
;
; IF keyword_set(key_performance) THEN temp = systime(1)
; zremplit(ocean)=zdata(ocean)
; zdata = zremplit
; IF keyword_set(key_performance) THEN print, systime(1)-temp
; IF keyword_set(key_performance) THEN print, 'fin iteration', it
;
; on trouve les points cote
   tempdeux = systime(1)        ; pour key_performance =2
; les points du bord du cadre ne doivent pas etre dans la cote
   mmmask[0, *] = 1
   mmmask[jpiatm+1, *] = 1
   mmmask[*, 0] = 1
   mmmask[*, jpjatm+1] = 1
; liste des points terre restant
   terre = where(mmmask EQ 0)
   if terre[0] EQ -1 then GOTO, fini
; les points du bord du cadre doivent maintenant etre dans la terre
   mmmask[0, *] = 0
   mmmask[jpiatm+1, *] = 0
   mmmask[*, 0] = 0
   mmmask[*, jpjatm+1] = 0
; liste des voisins mer
   voisins = shift(mmmask,-1,0)+shift(mmmask,1,0)+shift(mmmask,0,-1)+shift(mmmask,0,1) $
    +1./sqrt(2)*(shift(mmmask,-1,-1)+shift(mmmask,1,1)+shift(mmmask,1,-1)+shift(mmmask,-1,1))
; liste des points cote
   cote = where(voisins[terre] GT 0)
   cote = terre[cote]
;
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps remplit: trouver la cote ', systime(1)-tempdeux
;---------------------------------------------------------------
; remplissage des points cote
;---------------------------------------------------------------
   tempdeux = systime(1)        ; pour key_performance =2
; on masque z
   z = z*mmmask
;
   zcote = z[1+cote]+z[-1+cote]+z[(jpiatm+2)+cote]+z[-(jpiatm+2)+cote] $
    +1./sqrt(2)*(z[(jpiatm+2)+1+cote]+z[(jpiatm+2)-1+cote]+z[-(jpiatm+2)+1+cote]+z[-(jpiatm+2)-1+cote])
   poids = voisins[cote]

   z[cote] = zcote/poids
;---------------------------------------------------------------
; IV) on reduit le masque
;---------------------------------------------------------------
   mmmask[cote] = 1
   z = z*mmmask + val*(1-mmmask)
;   print, max(z)
;
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps remplit: une iteration ', systime(1)-tempdeux
fini:

return
END