Index: trunk/SRC/Computation/grad.pro
===================================================================
--- trunk/SRC/ToBeReviewed/CALCULS/grad.pro	(revision 163)
+++ trunk/SRC/Computation/grad.pro	(revision 167)
@@ -4,26 +4,37 @@
 ;+
 ; @file_comments
-;
+; compute the gradient of a variable
 ;
 ; @categories
-;
+; Calculation
 ;
 ; @param FIELD
-;
-;
-; @param DIREC
-;
-;
-; @returns
-;
+; The field for which we want to compute the gradient.  A 2D (xy),
+; 3D (xyz or yt) or 4D (xyzt) array or a structure readable by litchamp
+; and containing a 2D (xy), 3D (xyz or yt) or 4D (xyzt) array.
+; note that the dimension of the arry must suit the domain dimension.
+;
+; @param DIREC {type=scalar string}
+; the gradient direction: 'x', 'y', 'z'
+;
+; @returns RES {type=2D, 3D or 4D array}
+; the gradient of the input data (with the same size)
 ;
 ; @uses
-;
+; cm_4cal, cm_4data, cm_4mmesh 
 ;
 ; @restrictions
-;
+; - Works only for Arakawa C-grid. 
+; - When computing the gradient, the result is not on the same grid point
+;   than the input data. In consequence, we update, vargrid and the grid position
+;   parameters (firstx[tuvf], lastx[tuvf], nx[tuvf], firsty[tuvf], lasty[tuvf],
+;   ny[tuvf], firstz[tw], lastz[tw], nz[tw]).
+; - points that cannot be computed (domain bondaries, coastline) are set to NaN
+; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases.
 ;
 ; @examples
-;
+; IDL> @tst_initorca2
+; IDL> plt, grad({arr:gphit,g:'T'}, 'x')
+; IDL> plt, grad({arr:gphit,g:'T'}, 'y')
 ;
 ; @history
@@ -33,6 +44,4 @@
 ; $Id$
 ;
-; @todo 
-; seb: remplir les truc!
 ;-
 ;------------------------------------------------------------
@@ -43,118 +52,123 @@
   compile_opt idl2, strictarrsubs
 ;
-@common
-;------------------------------------------------------------
-;
-   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
+@cm_4cal   ; for jpt
+@cm_4data  ; for varname, vargrid, vardate, varunit, valmask
+@cm_4mesh 
+;------------------------------------------------------------
+;
+  IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
      return, report(['This version of grad is based on Arakawa C-grid.' $
                      , 'U and V grids must therefore be defined'])
 ;
-   res = litchamp(field)
-   taille = size(res)
-   grille, mask, glam, gphi, gdep, nx, ny, nz $
-           , firstx, firsty, firstz, lastx, lasty, lastz     
-   if n_elements(valmask) EQ 0 then valmask = 1e20
-   case strupcase(vargrid) of
-      'T':BEGIN
-         case direc of
-            'x':BEGIN 
-               divi = e1u[firstx:lastx, firsty:lasty]
-               newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
-               vargrid = 'U'
-               firstxu = firstxt & lastxu = lastxt & nxu = nxt 
-               firstyu = firstyt & lastyu = lastyt & nyu = nyt
-            END
-            'y':BEGIN
-               divi = e2v[firstx:lastx, firsty:lasty]
-               newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
-               vargrid = 'V'
-               firstxv = firstxt & lastxv = lastxt & nxv = nxt 
-               firstyv = firstyt & lastyv = lastyt & nyv = nyt
-            END
-            'z':BEGIN
-               divi = e3w[firstz:lastz]
-               newmask = mask
-               vargrid = 'W'
-               firstzw = firstzt & lastzw = lastzt & nzw = nzt 
-            END
-            ELSE:return, report('Bad definition of direction argument')
-         ENDCASE
-      END
-      'W':BEGIN
-         case direc of
-            'x':BEGIN 
-               divi = e1u[firstx:lastx, firsty:lasty]
-               newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
-               vargrid = 'U'
-               firstxu = firstxt & lastxu = lastxt & nxu = nxt 
-               firstyu = firstyt & lastyu = lastyt & nyu = nyt
-            END
-            'y':BEGIN
-               divi = e2v[firstx:lastx, firsty:lasty]
-               newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
-               vargrid = 'V'
-               firstxv = firstxt & lastxv = lastxt & nxv = nxt 
-               firstyv = firstyt & lastyv = lastyt & nyv = nyt
-            END
-            'z':BEGIN
-               divi = e3t[firstz:lastz]
-               newmask = mask
-               vargrid = 'T'
-               firstzt = firstzw & lastzt = lastzw & nzt = nzw 
-            END
-            ELSE:return, report('Bad definition of direction argument')
-         endcase
-      END
-      'U':BEGIN
-         case direc of
-            'x':BEGIN
-               divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty]
-               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
-               vargrid = 'T'
-               firstxt = firstxu & lastxt = lastxu & nxt = nxu 
-               firstyt = firstyu & lastyt = lastyu & nyt = nyu
-            END
-            'y':BEGIN
-               divi = e2f[firstx:lastx, firsty:lasty]
-               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
-               vargrid = 'F'
-               firstxf = firstxu & lastxf = lastxu & nxf = nxu 
-               firstyf = firstyu & lastyf = lastyu & nyf = nyu
-            END
-            'z':BEGIN
-               divi = e3w[firstz:lastz]
-               newmask = mask
-               vargrid = 'W'
-               firstzw = firstzt & lastzw = lastzt & nzw = nzt 
-            END
-            ELSE:return, report('Bad definition of direction argument')
-         endcase
-      END
-      'V':BEGIN
-         case direc of
-            'x':BEGIN
-               divi = e1f[firstx:lastx, firsty:lasty]
-               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
-               vargrid = 'F'
-               firstxf = firstxv & lastxf = lastxv & nxf = nxv
-               firstyf = firstyv & lastyf = lastyv & nyf = nyv
-            END
-            'y':BEGIN
-               divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
-               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
-               vargrid = 'T'
-               firstxt = firstxv & lastxt = lastxv & nxt = nxv
-               firstyt = firstyv & lastyt = lastyv & nyt = nyv
-            END
-            'z':BEGIN
-               divi = e3w[firstz:lastz]
-               newmask = mask
-               vargrid = 'W'
-               firstzw = firstzt & lastzw = lastzt & nzw = nzt 
-            END
-            ELSE:return, report('Bad definition of direction argument')
-         endcase
-      END
-;       'F':BEGIN
+  res = litchamp(field)
+  szres = size(res)
+  grille, mask, glam, gphi, gdep, nx, ny, nz $
+          , firstx, firsty, firstz, lastx, lasty, lastz     
+;
+  if n_elements(valmask) EQ 0 then valmask = 1.e20
+  varname = 'grad of '+varname
+  varunit = varunit+'/m'
+  case strupcase(vargrid) of
+    'T':BEGIN
+      case direc of
+        'x':BEGIN 
+          divi = e1u[firstx:lastx, firsty:lasty]
+          newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
+          vargrid = 'U'
+          firstxu = firstxt & lastxu = lastxt & nxu = nxt 
+          firstyu = firstyt & lastyu = lastyt & nyu = nyt
+        END
+        'y':BEGIN
+          divi = e2v[firstx:lastx, firsty:lasty]
+          newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
+          vargrid = 'V'
+          firstxv = firstxt & lastxv = lastxt & nxv = nxt 
+          firstyv = firstyt & lastyv = lastyt & nyv = nyt
+        END
+        'z':BEGIN
+          divi = e3w[firstz:lastz]
+          newmask = mask
+          vargrid = 'W'
+          firstzw = firstzt & lastzw = lastzt & nzw = nzt 
+        END
+        ELSE:return, report('Bad definition of direction argument')
+      ENDCASE
+    END
+    'W':BEGIN
+      case direc of
+        'x':BEGIN 
+          divi = e1u[firstx:lastx, firsty:lasty]
+          newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
+          vargrid = 'U'
+          firstxu = firstxt & lastxu = lastxt & nxu = nxt 
+          firstyu = firstyt & lastyu = lastyt & nyu = nyt
+        END
+        'y':BEGIN
+          divi = e2v[firstx:lastx, firsty:lasty]
+          newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
+          vargrid = 'V'
+          firstxv = firstxt & lastxv = lastxt & nxv = nxt 
+          firstyv = firstyt & lastyv = lastyt & nyv = nyt
+        END
+        'z':BEGIN
+          divi = e3t[firstz:lastz]
+          newmask = mask
+          vargrid = 'T'
+          firstzt = firstzw & lastzt = lastzw & nzt = nzw 
+        END
+        ELSE:return, report('Bad definition of direction argument')
+      endcase
+    END
+    'U':BEGIN
+      case direc of
+        'x':BEGIN
+          divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty]
+          newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
+          vargrid = 'T'
+          firstxt = firstxu & lastxt = lastxu & nxt = nxu 
+          firstyt = firstyu & lastyt = lastyu & nyt = nyu
+        END
+        'y':BEGIN
+          divi = e2f[firstx:lastx, firsty:lasty]
+          newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
+          vargrid = 'F'
+          firstxf = firstxu & lastxf = lastxu & nxf = nxu 
+          firstyf = firstyu & lastyf = lastyu & nyf = nyu
+        END
+        'z':BEGIN
+          divi = e3w[firstz:lastz]
+          newmask = mask
+          vargrid = 'W'
+          firstzw = firstzt & lastzw = lastzt & nzw = nzt 
+        END
+        ELSE:return, report('Bad definition of direction argument')
+      endcase
+    END
+    'V':BEGIN
+      case direc of
+        'x':BEGIN
+          divi = e1f[firstx:lastx, firsty:lasty]
+          newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
+          vargrid = 'F'
+          firstxf = firstxv & lastxf = lastxv & nxf = nxv
+          firstyf = firstyv & lastyf = lastyv & nyf = nyv
+        END
+        'y':BEGIN
+          divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
+          newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
+          vargrid = 'T'
+          firstxt = firstxv & lastxt = lastxv & nxt = nxv
+          firstyt = firstyv & lastyt = lastyv & nyt = nyv
+        END
+        'z':BEGIN
+          divi = e3w[firstz:lastz]
+          newmask = mask
+          vargrid = 'W'
+          firstzw = firstzt & lastzw = lastzt & nzw = nzt 
+        END
+        ELSE:return, report('Bad definition of direction argument')
+      endcase
+    END
+    'F':BEGIN
 ;          case direc of
 ;             'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty]
@@ -163,140 +177,151 @@
 ;             ELSE:return, report('Bad definition of direction argument')
 ;          endcase
-;       END
-      ELSE:return, report('Bad definition of vargrid')
-   ENDCASE
-   res = fitintobox(res)
-   IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res
-   case 1 of
+      return, report('F grid: case not coded, please contact SAXO team...')
+    END
+    ELSE:return, report('Bad definition of vargrid')
+  ENDCASE
+  res = fitintobox(temporary(res))
+  IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res
+  case 1 of
 ;----------------------------------------------------------------------------
 ;xy
 ;----------------------------------------------------------------------------
-      taille[0] EQ 2:BEGIN
-         earth = where(mask[*, *, firstz] EQ 0)
-         if earth[0] NE -1 then res[earth] = !values.f_nan
-         case direc of
-            'x':BEGIN 
-               res = (shift(res, -1, 0)-res)/divi
-               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan
-               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0)
-            END
-            'y':BEGIN 
-               res = (shift(res, 0, -1)-res)/divi
-               res[*, ny-1] = !values.f_nan
-               if vargrid EQ 'T' OR vargrid EQ 'U' then res =   shift(res, 0, 1)           
-            END
-            ELSE:return,  report('Bad definition of direction argument for the type of array')
-         ENDCASE
-         earth = where(newmask[*, *, firstz] EQ 0)
-         if earth[0] NE -1 then res[earth] = valmask
-      END
+    szres[0] EQ 2:BEGIN
+      land = where((temporary(mask))[*, *, firstz] EQ 0)
+      if land[0] NE -1 then res[temporary(land)] = !values.f_nan
+      case direc of
+        'x':BEGIN 
+          res = (shift(res, -1, 0)-res)/temporary(divi)
+          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan
+          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0)
+        END
+        'y':BEGIN 
+          res = (shift(res, 0, -1)-res)/temporary(divi)
+          res[*, ny-1] = !values.f_nan
+          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1)           
+        END
+        ELSE:return,  report('Bad definition of direction argument for the type of array')
+      ENDCASE
+      land = where((temporary(newmask))[*, *, firstz] EQ 0)
+      if land[0] NE -1 then res[temporary(land)] = valmask
+    END
 ;----------------------------------------------------------------------------
 ;xyt
 ;----------------------------------------------------------------------------
-      taille[0] EQ 3 AND jpt NE 1:BEGIN
-         earth = where(mask[*, *, firstz] EQ 0)
-         if earth[0] NE -1 then BEGIN
-            earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt))
-            res[earth[*]] = !values.f_nan
-         ENDIF
-         divi = divi[*]#replicate(1, jpt)
-         case direc of
-            'x':BEGIN 
-               res = (shift(res, -1, 0, 0)-res)/divi
-               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
-               if vargrid EQ 'T' OR vargrid EQ 'V'  then res = shift(res, 1, 0, 0)
-            END
-            'y':BEGIN 
-               res = (shift(res, 0, -1, 0)-res)/divi
-               res[*, ny-1, *] = !values.f_nan
-               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)        
-            END
-            ELSE:return,  report('Bad definition of direction argument for the type of array')
-         ENDCASE
-         earth = where(newmask[*, *, firstz] EQ 0)
-         if earth[0] NE -1 THEN res[earth] = valmask
-      END
+    szres[0] EQ 3 AND jpt NE 1:BEGIN
+      land = where((temporary(mask))[*, *, firstz] EQ 0, cnt)
+      if land[0] NE -1 then BEGIN
+        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
+        res[temporary(land)] = !values.f_nan
+      ENDIF
+      divi = (temporary(divi))[*]#replicate(1., jpt)
+      case direc of
+        'x':BEGIN 
+          res = (shift(res, -1, 0, 0)-res)/temporary(divi)
+          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
+          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0)
+        END
+        'y':BEGIN 
+          res = (shift(res, 0, -1, 0)-res)/temporary(divi)
+          res[*, ny-1, *] = !values.f_nan
+          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)        
+        END
+        ELSE:return,  report('Bad definition of direction argument for the type of array')
+      ENDCASE
+      land = where((temporary(newmask))[*, *, firstz] EQ 0, cnt)
+      if land[0] NE -1 then BEGIN
+        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
+        res[temporary(land)] = valmask
+      ENDIF
+    END
 ;----------------------------------------------------------------------------
 ;xyz
 ;----------------------------------------------------------------------------
-      taille[0] EQ 3 AND jpt EQ 1:BEGIN
-         earth = where(mask EQ 0)
-         if earth[0] NE -1 then res[earth] = !values.f_nan
-         case direc OF
-            'x':BEGIN 
-               divi = divi[*]#replicate(1, nz)
-               res = (shift(res, -1, 0, 0)-res)/divi
-               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
-               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0)
-            END
-            'y':BEGIN 
-               divi = divi[*]#replicate(1, nz)
-               res = (shift(res, 0, -1, 0)-res)/divi
-               res[*, ny-1, *] = !values.f_nan
-               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)        
-            END
-            'z':BEGIN
-               divi = reform(replicate(1, nx*ny)#divi, nx, ny, nz)
-               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz)
-               if vargrid EQ 'W' THEN BEGIN 
-                  res = (shift(res, 0, 0, 1)-res)/divi
-                  res[*, *, 0] = !values.f_nan
-               ENDIF ELSE BEGIN
-                  res = (res-shift(res, 0, 0, -1))/divi
-                  res[*, *, nz-1] = !values.f_nan
-               ENDELSE
-               if earth[0] NE -1 then res[earth] = valmask
-            END
-         ENDCASE
-      END
-;------------------------------------------------------------
+    szres[0] EQ 3 AND jpt EQ 1:BEGIN
+      land = where(mask EQ 0)
+      if land[0] NE -1 then res[temporary(land)] = !values.f_nan
+      case direc OF
+        'x':BEGIN 
+          divi = (temporary(divi))[*]#replicate(1., nz)
+          res = (shift(res, -1, 0, 0)-res)/temporary(divi)
+          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
+          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0)
+        END
+        'y':BEGIN 
+          divi = (temporary(divi))[*]#replicate(1., nz)
+          res = (shift(res, 0, -1, 0)-res)/temporary(divi)
+          res[*, ny-1, *] = !values.f_nan
+          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)        
+        END
+        'z':BEGIN
+          divi = replicate(1., nx*ny)#(temporary(divi))[*]
+          if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, /overwrite)
+          if vargrid EQ 'W' THEN BEGIN 
+            res = (shift(res, 0, 0, 1)-res)/temporary(divi)
+            res[*, *, 0] = !values.f_nan
+          ENDIF ELSE BEGIN
+            res = (res-shift(res, 0, 0, -1))/temporary(divi)
+            res[*, *, nz-1] = !values.f_nan
+          ENDELSE
+        END
+      ENDCASE
+      land = where(temporary(newmask) EQ 0)
+      if land[0] NE -1 then res[temporary(land)] = valmask
+    END
 ;----------------------------------------------------------------------------
 ;xyzt
 ;----------------------------------------------------------------------------
-      taille[0] EQ 4:BEGIN
-         earth = where((mask[*]#replicate(1, jpt)) EQ 0)
-         if earth[0] NE -1 then res[earth] = !values.f_nan
-         case direc OF
-            'x':BEGIN 
-               divi = divi[*]#replicate(1, nz*jpt)
-               res = (shift(res, -1, 0, 0, 0)-res)/divi
-               if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan
-               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0, 0)
-            END
-            'y':BEGIN 
-               divi = divi[*]#replicate(1, nz*jpt)
-               res = (shift(res, 0, -1, 0, 0)-res)/divi
-               res[*, ny-1, *, *] = !values.f_nan
-               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0, 0)        
-            END
-            'z':BEGIN
-               divi = replicate(1, nx*ny)#divi
-               divi = reform(divi[*]#replicate(1, jpt), nx, ny, nz, jpt, /over)
-               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt)
-               if vargrid EQ 'W' THEN BEGIN 
-                  res = (shift(res, 0, 0, 1, 0)-res)/divi
-                  res[*, *, 0, *] = !values.f_nan
-               ENDIF ELSE BEGIN
-                  res = (res-shift(res, 0, 0, -1, 0))/divi
-                  res[*, *, nz-1, *] = !values.f_nan
-               ENDELSE
-            END
-         ENDCASE
-         if earth[0] NE -1 then res[earth] = valmask
-      END
-;------------------------------------------------------------
-;------------------------------------------------------------
-   endcase
-   varname = 'grad of '+varname
-   varunit = varunit+'/m'
-
-
-;------------------------------------------------------------
-   return, res
-end
-
-
-
-
-
-
+    szres[0] EQ 4:BEGIN
+      land = where(temporary(mask) EQ 0, cnt)
+      if land[0] NE -1 then BEGIN
+        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt))
+        res[temporary(land)] = !values.f_nan
+      ENDIF
+      case direc OF
+        'x':BEGIN 
+          divi = (temporary(divi))[*]#replicate(1., nz*jpt)
+          res = (shift(res, -1, 0, 0, 0)-res)/temporary(divi)
+          if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan
+          if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0, 0)
+        END
+        'y':BEGIN 
+          divi = (temporary(divi))[*]#replicate(1., nz*jpt)
+          res = (shift(res, 0, -1, 0, 0)-res)/temporary(divi)
+          res[*, ny-1, *, *] = !values.f_nan
+          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0, 0)        
+        END
+        'z':BEGIN
+          divi = replicate(1., nx*ny)#(temporary(divi))[*]
+          divi = (temporary(divi))[*]#replicate(1L, jpt)
+          if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt, /overwrite)
+          if vargrid EQ 'W' THEN BEGIN 
+            res = (shift(res, 0, 0, 1, 0)-res)/temporary(divi)
+            res[*, *, 0, *] = !values.f_nan
+          ENDIF ELSE BEGIN
+            res = (res-shift(res, 0, 0, -1, 0))/temporary(divi)
+            res[*, *, nz-1, *] = !values.f_nan
+          ENDELSE
+        END
+      ENDCASE
+      land = where(newmask EQ 0, cnt)
+      if land[0] NE -1 then BEGIN
+        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt))
+        res[temporary(land)] = valmask
+      ENDIF
+    END
+;------------------------------------------------------------
+;------------------------------------------------------------
+    ELSE:return, report('input array must have 2, 3 or 4 dimensions')
+  ENDCASE
+;------------------------------------------------------------
+
+
+;------------------------------------------------------------
+  return, res
+END
+
+
+
+
+
+
