Changeset 384 for trunk

Timestamp:

11/06/08 03:32:33 (16 years ago)

Author:

smasson

Message:

add fill(xy)dir keyword to extrapolate

File:

• trunk/SRC/Interpolation/extrapolate.pro (modified) (6 diffs)

trunk/SRC/Interpolation/extrapolate.pro

@@ -20 +20 @@
 ; is no more masked values we exit extrapolate before reaching nb_iteration
 ;
+; @keyword FILLXDIR {type=scalar, 0 or 1}{default=0}
+; put 1 to specify that filling of the data must be done only along x direction
+;
+; @keyword FILLYDIR {type=scalar, 0 or 1}{default=0}
+; put 1 to specify that filling of the data must be done only along y direction
+;
 ; @keyword X_PERIODIC {type=scalar, 0 or 1}{default=0}
 ; put 1 to specify that the data are periodic along x axis
@@ -57 +63 @@
 ;
 ;-
-FUNCTION extrapolate, zinput, maskinput, nb_iteration, X_PERIODIC=x_periodic $
-                      , MINVAL=minval, MAXVAL=maxval, GE0=ge0, _EXTRA=ex
+FUNCTION extrapolate, zinput, maskinput, nb_iteration, FILLXDIR = fillxdir, FILLYDIR = fillydir $
+                      , X_PERIODIC = x_periodic, MINVAL=minval, MAXVAL=maxval, GE0=ge0, _EXTRA=ex
 ;
   compile_opt idl2, strictarrsubs
@@ -118 +124 @@
 ; find the land points
   land = where(msk EQ 0, cnt_land)
+  IF keyword_set(fillxdir) THEN BEGIN 
+    tst = total(msk[1:nx, 1:ny], 1)
+    cnt_land = total(tst NE 0 AND tst NE nx)
+  ENDIF 
+  IF keyword_set(fillydir) THEN BEGIN 
+    tst = total(msk[1:nx, 1:ny], 2)
+    cnt_land = total(tst NE 0 AND tst NE ny)
+  ENDIF 
 ;---------------------------------------------------------------
   WHILE cnt LE nb_iteration AND cnt_land NE 0 DO BEGIN
@@ -149 +163 @@
 ; 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])
+; 
+    CASE 1 OF
+      keyword_set(fillxdir):weight = msk[land+1]+msk[land-1]
+      keyword_set(fillydir):weight = msk[land+nx2]+msk[land-nx2]
+      ELSE: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])
+    ENDCASE
 ; list all the points that have sea neighbors
     ok = where(weight GT 0)
@@ -164 +182 @@
     z = temporary(z)*msk
 ;
-    zcoast = z[1+coast]+z[-1+coast]+z[nx2+coast]+z[-nx2+coast] $
+    CASE 1 OF
+      keyword_set(fillxdir):zcoast = z[1+coast]+z[-1+coast]
+      keyword_set(fillydir):zcoast = z[nx2+coast]+z[-nx2+coast]
+      ELSE: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])
+    ENDCASE
 ;
     IF keyword_set(ge0) THEN zcoast = 0. > temporary(zcoast)
@@ -192 +214 @@
 ; find the land points
     land = where(msk EQ 0, cnt_land)
+    IF keyword_set(fillxdir) THEN BEGIN 
+      tst = total(msk[1:nx, 1:ny], 1)
+      cnt_land = total(tst NE 0 AND tst NE nx)
+    ENDIF 
+    IF keyword_set(fillydir) THEN BEGIN 
+      tst = total(msk[1:nx, 1:ny], 2)
+      cnt_land = total(tst NE 0 AND tst NE ny)
+    ENDIF 
 ;
   ENDWHILE