source: trunk/CALCULS/remplit.pro @ 2

FUNCTION remplit, zinput, NAN = nan, NITE = nite, BASIQUE = basique, mask = mask, _extra = ex
@common
   tempsun = systime(1)         ; pour key_performance
;+
   ;;
   ;; Extrapole zinout[jpi,jpj] sur les continents en utilisant les 4
   ;; plus proches voisins masques oceaniquement et construit un nouveau
;   masque
   ;; contenant l'ancien masque oceanique PLUSles points extrapoles.
   ;; Reitere le processus nite fois.
   ;; C'est pas clair, essayez !
   ;;
   ;;
;
;    /Nan: to fill the point which have the value
;    !values.f_nan. Whitout this keyword, these point are not filling
;    and stays at !values.f_nan.
;
;



;-
; les points non remplis sont masques a valmask
   z = zinput
   if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c'
   if keyword_set(basique) then begin
      oldkey_gridtype = key_gridtype
      key_gridtype = 'c'
   endif
   IF n_elements(nite) EQ 0 THEN nite = 1
   IF nite EQ 0 THEN return, zinput
   if keyword_set(basique)  then begin
      nx = (size(zinput))[1]
      ny = (size(zinput))[2]
      if NOT keyword_set(mask) then mmmask = basique ELSE mmmask = mask
      if  key_gridtype eq 'e' then begin
         case vargrid of
            'T':glam = glamt[premierxt:dernierxt, premieryt:dernieryt]
            'U':glam = glamu[premierxu:dernierxu, premieryu:dernieryu]
         endcase
      endif
   ENDIF ELSE grille, mmmask, glam, gphi, gdep, nx, ny, nz, _extra = ex
   if keyword_set(mask) then mmmask = mask
;---------------------------------------------------------------
   if (size(mmmask))[0] EQ 3 THEN BEGIN 
      if (size(mmmask))[3] EQ jpk THEN mmmask = mmmask[*, *, niveau-1] $
      ELSE  mmmask = mmmask[*, *, nz-1]
   ENDIF
;
   if n_elements(mmmask) EQ 1 then mmmask = replicate(1, nx, ny)
   if keyword_set(nan) then begin
      nanpoint = where(finite(z) EQ 0)
      if nanpoint[0] NE -1 then begin
         mmmask[nanpoint] = 0
         z[nanpoint] = 0
      endif
   endif
;---------------------------------------------------------------
; on ajoute un cadre de zero a z, mask, e1, e2
; comme ca apres on peut faire des shifts ds tous les sens sans se
; soucier des bords du domaine!
;---------------------------------------------------------------
   tempdeux = systime(1)        ; pour key_performance =2
   nx2 = nx+2
   case key_gridtype of
      'c':BEGIN
         zero = fltarr(nx+2, ny+2)
         zero[1:nx, 1:ny] = mmmask
         mmmask = zero
         zero = fltarr(nx+2, ny+2)
         zero[1:nx, 1:ny] = z
         if keyword_set(key_periodique) AND nx EQ jpi then begin
            zero[0, 1:ny] = z[jpi-1, *]
            zero[nx+1, 1:ny] = z[0, *]
         endif
         z = zero
      END
      'e':BEGIN
         zero = fltarr(nx+2, ny+4)
         zero[1:nx, 2:ny+1] = mmmask
         mmmask = zero
         zero = fltarr(nx+2, ny+4)
         zero[1:nx, 2:ny+1] = z
         if keyword_set(key_periodique) AND nx EQ jpi then begin
            zero[0, 2:ny+1] = z[jpi-1, *]
            zero[nx+1, 2:ny+1] = z[0, *]
         endif
         z = zero
      END
   endcase
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps remplit: on ajoute un cadre de zero ', systime(1)-tempdeux
;---------------------------------------------------------------
;---------------------------------------------------------------
; iteration 
;---------------------------------------------------------------
;---------------------------------------------------------------
   FOR n = 1, nite DO BEGIN 
; 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
      case key_gridtype of
         'c':BEGIN
            mmmask[0, *] = 1
            mmmask[nx+1, *] = 1
            mmmask[*, 0] = 1
            mmmask[*, ny+1] = 1
            if keyword_set(key_periodique) AND nx EQ jpi then begin
               mmmask[0,*] = mmmask[nx-2, *]
               mmmask[nx+1,*] = mmmask[1, *]
            endif
         END
         'e':BEGIN
            mmmask[0, *] = 1
            mmmask[nx+1, *] = 1
            mmmask[*, 0:1] = 1
            mmmask[*, ny+2:ny+3] = 1
            if keyword_set(key_periodique) AND nx EQ jpi then begin
               mmmask[0,*] = mmmask[nx-2, *]
               mmmask[nx+1,*] = mmmask[1, *]
            endif
         END
      endcase
; 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
      case key_gridtype of
         'c':BEGIN
            mmmask[0, *] = 0
            mmmask[nx+1, *] = 0
            mmmask[*, 0] = 0
            mmmask[*, ny+1] = 0
         END
         'e':BEGIN
            mmmask[0, *] = 0
            mmmask[nx+1, *] = 0
            mmmask[*, 0:1] = 0
            mmmask[*, ny+2:ny+3] = 0
         END
      endcase
; liste des voisins mer
      case key_gridtype of
         'c':BEGIN
            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))
            cote = where(voisins[terre] GT 0)
; liste des points cote
            cote = terre[cote]
            poids = voisins[cote]
         END
         'e':BEGIN
            shifted = glam[0, 0] LT glam[0, 1]
            oddeven = (terre/nx2+1-shifted) MOD 2
            voisins = mmmask[terre+1]+mmmask[terre-1] $
             +mmmask[2*nx2+terre]+mmmask[-2*nx2+terre] $
             +sqrt(2)*(mmmask[terre-nx2+oddeven]+mmmask[terre+nx2+oddeven] $
                       +mmmask[terre+(-nx2-1)+oddeven]+mmmask[terre+(nx2-1)+oddeven])
            cote = where(voisins GT 0)
            poids = voisins[cote]
            cote = terre[cote]
         END
      endcase
;
      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
;
      case key_gridtype of
         'c':BEGIN
            zcote = z[1+cote]+z[-1+cote]+z[nx2+cote]+z[-nx2+cote] $
             +1./sqrt(2)*(z[nx2+1+cote]+z[nx2-1+cote] $
                          +z[-nx2+1+cote]+z[-nx2-1+cote])
         END
         'e':BEGIN
            oddeven = (cote/nx2+1-shifted) MOD 2
            zcote = z[1+cote]+z[-1+cote]+z[2*nx2+cote]+z[-2*nx2+cote] $
             +sqrt(2)*(z[-nx2+oddeven+cote]+z[-nx2-1+oddeven+cote] $
                       +z[nx2+oddeven+cote]+z[nx2-1+oddeven+cote])
         END
      endcase
      
      z[cote] = zcote/poids
;---------------------------------------------------------------
; IV) on reduit le masque
;---------------------------------------------------------------
      mmmask[cote] = 1
;
      IF testvar(var = key_performance) EQ 2 THEN $
       print, 'temps remplit: une iteration ', systime(1)-tempdeux
   ENDFOR 
fini:
;
; on masque les valeurs sur les terres restantes
;
   IF n_elements(valmask) EQ 0 then valmask = 1e20
   z = z*mmmask+valmask*(1-mmmask)
;---------------------------------------------------------------
; on redecoupe le tableau pour retirer le cadre!
;---------------------------------------------------------------
   case key_gridtype of
      'c':BEGIN
         z = z[1:nx, 1:ny]
      END
      'e':BEGIN
         z = z[1:nx, 2:ny+1]
      END
   endcase
;
   if keyword_set(basique) then key_gridtype = oldkey_gridtype
;---------------------------------------------------------------
   if keyword_set(key_performance) THEN print, 'temps remplit', systime(1)-tempsun 
   return, z
END