Changeset 167 for trunk/SRC/Computation/div.pro

Timestamp:

09/05/06 14:24:07 (18 years ago)

Author:

smasson

Message:

clean div, curl, grad + minors bugfix

Location:

trunk/SRC/Computation

Files:

• . (added)
• div.pro (moved) (moved from trunk/SRC/ToBeReviewed/CALCULS/div.pro) (5 diffs)

trunk/SRC/Computation/div.pro

@@ -5 +5 @@
 ;
 ; @file_comments
-; calculation of the divergence of a 2d field
+; compute the horizontal divergence of a vectors field
 ;
 ; @categories
@@ -11 +11 @@
 ;
 ; @param UU 
-; Matrix representing coordinates of a field of vectors
+; Matrix representing the zonal coordinates (U point) of a field of vectors
+; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing
+; a 2D (xy), 3D (xyz or yt) array (4d case is not coded yet).
+; note that the dimension of the arry must suit the domain dimension.
 ;
 ; @param VV 
-; Matrix representing coordinates of a field of vectors
+; Matrix representing the meridional coordinates (V point) of a field of vectors
+; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing
+; a 2D (xy), 3D (xyz or yt) array (4d case is not coded yet).
+; note that the dimension of the arry must suit the domain dimension.
+;
+; @keyword DIREC {type=scalar string}
+; Use if you want to call moyenne or grossemoyenne after the div computation
+; (stupid ?) with a mean done in the DIREC direction
 ;
 ; @returns RES
-; A 2d matrix
+; the divergence of the input data (with the same size)
 ;
 ; @uses
-; common.pro
+; cm_4cal, cm_4data, cm_4mmesh 
 ;
 ; @restrictions
-; U and V matrices can be 2 or 4d.
-; Beware, to discern different configuration of U and V (xy, xyz, xyt, xyzt), 
-; we look at the variable of the common 
-;        -time which contain the calendar in IDL julian days to which U and 
-; V refered to, in the same way as the variable 
-;        -jpt which is the number of time's step to consider in time.
-; U and V arrays are cut in the same geographic domain. Because of the gap of 
-; T, U, V and F grids, it is possible that these two arrays has not the same 
-; size and refered to different indexes. In this case, arrays are re-cut on 
-; common indexes and the domain is redefined to match with these common 
-; indexes. To avoid these re-cuts, use the keyword /memeindice in domdef.pro
-;
-;
-; Points on the drawing edge are at !values.f_nan 
-;
+;
+; - Works only for Arakawa C-grid. 
+; - UU must be on U grid, VV must be on V grid
+; - 4d case is not coded yet
+; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases.
+; - U and V arrays are cut in the same geographic domain. Because of the shift between 
+;   T, U, V and F grids, it is possible that these two arrays do not have the same 
+;   size and refer to different indexes. In this case, arrays are re-cut on 
+;   common indexes. To avoid these re-cuts, use the keyword /memeindice in domdef.pro
+; - When computing the divergence, we update, vargrid, varname, varunits and the
+;   grid position parameters (firstxt, lastxt, nxt, firstyt, lastyt, nyt).
+; - points that cannot be computed (domain bondaries, coastline) are set to NaN
+;
+; @examples
+; IDL> @tst_initorca2
+; IDL> plt, div(dist(jpi,jpj), dist(jpi,jpj))
 ;
 ; @history
-; Guillaume Roullet (grlod\@ipsl.jussieu.fr)
-;                      Creation : printemps 1998
-;                      Sebastien Masson (smasson\@lodyc.jussieu.fr)
-;                      adaptation to work with a reduce domain
-;                      12/1/2000;
+; Guillaume Roullet (grlod\@ipsl.jussieu.fr): creation; spring 1998
+; Sebastien Masson (smasson\@lodyc.jussieu.fr)
+; adaptation to work with a reduce domain; 12/1/2000
 ;
 ; @version
 ; $Id$
 ;
+; @todo 
+; code the 4d case
 ;-
 ;------------------------------------------------------------
 ;------------------------------------------------------------
 ;------------------------------------------------------------
-FUNCTION div, uu, vv
+FUNCTION div, uu, vv, DIREC = DIREC
 ;
   compile_opt idl2, strictarrsubs
 ;
-   tempsun = systime(1)         ; For key_performance
-@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 
+;
+  tempsun = systime(1)          ; For key_performance
+;
+  IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
      return, report(['This version of div is based on Arakawa C-grid.' $
                      , 'U and V grids must therefore be defined'])
 ;
-   u = litchamp(uu)
-   v = litchamp(vv)
-;
-   date1 = time[0]
-   if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
-   if (size(u))[0] NE (size(v))[0] then return,  -1
+  u = litchamp(uu)
+  v = litchamp(vv)
+;
+  szu = size(u)
+  szv = size(v)
+  
+  if szu[0] NE szv[0] then return, report('U and V input data must have the same number of dimensions')
 
 ;------------------------------------------------------------
 ; We find common points between U and V
 ;------------------------------------------------------------
-   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
-   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
-   indicex = inter(indicexu, indicexv)
-   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
-   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
-   indicey = inter(indiceyu, indiceyv)
-   nx = n_elements(indicex) 
-   ny = n_elements(indicey)
-   indice2d = lindgen(jpi, jpj)
-   indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
-   case 1 of
-;----------------------------------------------------------------------------
+  indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
+  indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
+  indicex = inter(indicexu, indicexv)
+  indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
+  indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
+  indicey = inter(indiceyu, indiceyv)
+  nx = n_elements(indicex) 
+  ny = n_elements(indicey)
+  indice2d = lindgen(jpi, jpj)
+  indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1]
+;----------------------------------------------------------------------------
+  vargrid = 'T'
+  varname = 'div'
+  varunits = '1.e6*s-1'
+  if n_elements(valmask) EQ 0 THEN valmask = 1.e20
+  firstxt = indicex[0] & lastxt = indicex[0]+nx-1 & nxt = nx 
+  firstyt = indicey[0] & lastyt = indicey[0]+ny-1 & nyt = ny
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+  case 1 of
 ;----------------------------------------------------------------------------
 ;xyz
 ;----------------------------------------------------------------------------
-      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 
+    szu[0] EQ 3 AND jpt EQ 1:BEGIN 
 ;------------------------------------------------------------
 ; extraction of U and V on the appropriated domain
 ;------------------------------------------------------------
-         case 1 of
-            (size(v))[0] NE 3: return,  -1
-            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
-             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
-               case 1 of
-                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
-                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
-                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
-                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
-                  ELSE :
-               endcase
-            END
-            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
-             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
-               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
-               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
-            END
-            ELSE:BEGIN 
-               zdiv = -1
-               GOTO, sortie
-            end
-            
-         endcase
-;------------------------------------------------------------
-; calcul de la divergence
-;------------------------------------------------------------
-         zu = (e2u[indice2d])[*]#replicate(1, nzt)
-         zu = reform(zu, nx, ny, nzt, /over)
-         zu = temporary(u)*temporary(zu) $
-          *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
-         terreu = where(zu EQ 0)
-         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
-;
-         zv = (e1v[indice2d])[*]#replicate(1, nzt)
-         zv = reform(zv, nx, ny, nzt, /over)
-         zv = temporary(v)*temporary(zv) $
-          *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
-         terrev = where(zv EQ 0)
-         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
-;
-         zdiv = 1e6/(e1t[indice2d]*e2t[indice2d])
-         zdiv = (zdiv)[*]#replicate(1, nzt)
-         zdiv = reform(zdiv, nx, ny, nzt, /over)
-         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $
-          *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
+      case 1 of
+        szu[1] EQ nxu AND szu[2] EQ nyu AND $
+           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 
+          case 1 of
+            nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
+            nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
+            nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
+            nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
+            ELSE :
+          endcase
+        END
+        szu[1] EQ jpi AND szu[2] EQ jpj AND $
+           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN 
+          u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *]
+          v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *]
+        END
+        ELSE:return, -1
+      endcase
+;------------------------------------------------------------
+; divergence computation
+;------------------------------------------------------------
+      zu = (e2u[indice2d])[*]#replicate(1., nzt)
+      landu = where((umask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
+      if landu[0] NE -1 then zu[temporary(landu)] = !values.f_nan
+      zu = temporary(u) * temporary(zu)
+;
+      zv = (e1v[indice2d])[*]#replicate(1., nzt)
+      landv = where((vmask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
+      if landv[0] NE -1 then zv[temporary(landv)] = !values.f_nan
+      zv = temporary(v) * temporary(zv)
+;
+      zdiv = (e1t[indice2d]*e2t[indice2d])[*]#replicate(1.e6, nzt)
+      zdiv = ( zu - shift(zu, 1, 0, 0) + zv - shift(zv, 0, 1, 0) ) * temporary(zdiv)
 ;------------------------------------------------------------
 ; Edging put at !values.f_nan
 ;------------------------------------------------------------
-         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
-            zdiv[0, *, *] = !values.f_nan
-            zdiv[nx-1, *, *] = !values.f_nan
-         endif
-         zdiv[*, 0, *] = !values.f_nan
-         zdiv[*, ny-1, *] = !values.f_nan
-;
-         zdiv = temporary(zdiv)
-;
-         if n_elements(valmask) EQ 0 THEN valmask = 1e20
-         terre =  where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
-         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
-;------------------------------------------------------------
-; For the graphic drawing
-;------------------------------------------------------------
-         vargrid = 'T'
-         varname = 'div'
-         varunits = '1e6*s-1'
-         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
-         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan)
-      END
+      if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
+        zdiv[0, *, *] = !values.f_nan
+        zdiv[nx-1, *, *] = !values.f_nan
+      endif
+      zdiv[*, 0, *] = !values.f_nan
+      zdiv[*, ny-1, *] = !values.f_nan
+;
+      land = where(tmask[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
+      if land[0] NE -1 then zdiv[temporary(land)] = valmask
+      if keyword_set(direc) then  zdiv = moyenne(zdiv, direc, /nan)
+    END
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
@@ -167 +168 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 
+    szu[0] EQ 3 AND jpt GT 1:BEGIN 
 ;------------------------------------------------------------
 ; extraction of U and V on the appropriated domain
 ;------------------------------------------------------------
-         case 1 of
-            (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1
-            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
-             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
-               case 1 of
-                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
-                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
-                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
-                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
-                  ELSE :
-               endcase
-            END
-            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
-             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
-               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
-               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
-            END
-            ELSE:return,  -1
-         endcase
-;------------------------------------------------------------
-; Calculation of the divergence
-;------------------------------------------------------------
-         zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
-         terreu = where(zu EQ 0)
-         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
-         zu = (zu)[*]#replicate(1, jpt)
-         zu = reform(zu, nx, ny, jpt, /over)
-         zu = temporary(u)*temporary(zu)
-;
-         zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
-         terrev = where(zv EQ 0)
-         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
-         zv = (zv)[*]#replicate(1, jpt)
-         zv = reform(zv, nx, ny, jpt, /over)
-         zv = temporary(v)*temporary(zv)
-;
-         zdiv = 1e6*tmask[indice2d+jpi*jpj*firstzt]/(e1t[indice2d]*e2t[indice2d])
-         zdiv = (zdiv)[*]#replicate(1, jpt)
-         zdiv = reform(zdiv, nx, ny, jpt, /over)
-         terre =  where(zdiv EQ 0)
-         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0))
+      case 1 of
+        szu[1] EQ nxu AND szu[2] EQ nyu AND $
+           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 
+          case 1 of
+            nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
+            nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
+            nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
+            nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
+            ELSE :
+          endcase
+        END
+        szu[1] EQ jpi AND szu[2] EQ jpj AND $
+           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN 
+          u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *]
+          v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *]
+        END
+        ELSE:return, -1
+      endcase
+;------------------------------------------------------------
+; divergence computation
+;------------------------------------------------------------
+      zu = e2u[indice2d]
+      landu = where((umask())[indice2d+jpi*jpj*firstzt] EQ 0)
+      if landu[0] NE -1 then zu[temporary(landu)] = !values.f_nan
+      zu = (temporary(zu))[*]#replicate(1., jpt)
+      zu = temporary(u) * temporary(zu)
+;
+      zv = e1v[indice2d]
+      landv = where((vmask())[indice2d+jpi*jpj*firstzt] EQ 0)
+      if landv[0] NE -1 then zv[temporary(landv)] = !values.f_nan
+      zv = (temporary(zv))[*]#replicate(1., jpt)
+      zv = temporary(v) * temporary(zv)
+;
+      zdiv = (e1t[indice2d]*e2t[indice2d])[*]#replicate(1.e6, jpt)
+      zdiv = ( zu - shift(zu, 1, 0, 0) + zv - shift(zv, 0, 1, 0) ) * temporary(zdiv)
 ;------------------------------------------------------------
 ; Edging put at !values.f_nan
 ;------------------------------------------------------------
-         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
-            zdiv[0, *, *] = !values.f_nan
-            zdiv[nx-1, *, *] = !values.f_nan
-         endif
-         zdiv[*, 0, *] = !values.f_nan
-         zdiv[*, ny-1, *] = !values.f_nan
-;
-         if n_elements(valmask) EQ 0 THEN valmask = 1e20
-         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
-;------------------------------------------------------------
-; for the graphic drawing
-;------------------------------------------------------------
-         vargrid = 'T'
-         varname = 'div'
-         varunits = '1e6*s-1'
-         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
-         if keyword_set(direc) then  zdiv = grossemoyenne(zdiv,direc,/nan)
-      END
+      if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
+        zdiv[0, *, *] = !values.f_nan
+        zdiv[nx-1, *, *] = !values.f_nan
+      endif
+      zdiv[*, 0, *] = !values.f_nan
+      zdiv[*, ny-1, *] = !values.f_nan
+;
+      land = where(tmask[indice2d+jpi*jpj*firstzt] EQ 0, cnt)
+      if land[0] NE -1 then BEGIN 
+        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
+        zdiv[temporary(land)] = valmask
+      ENDIF
+      if keyword_set(direc) then  zdiv = grossemoyenne(zdiv, direc, /nan)
+    END
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
@@ -238 +229 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 
-         return, report('non code!')
-      END
+    szu[0] EQ 4:BEGIN 
+      return, report('Case not coded contact saxo team or make a do loop!')
+    END
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
@@ -246 +237 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      ELSE:BEGIN                ;xy
-         indice3d = lindgen(jpi, jpj, jpk)
-         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt]
+    szu[0] EQ 2:BEGIN 
 ;------------------------------------------------------------
 ; extraction of U and V on the appropriated domain
 ;------------------------------------------------------------
-         case 1 of
-            (size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN 
-               zdiv = -1
-               GOTO, sortie
-            end
-            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
-             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
-               case 1 of
-                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
-                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
-                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
-                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
-                  ELSE :
-               endcase
-            END
-            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
-             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
-               u = u[indice2d]
-               v = v[indice2d]
-            END
-            ELSE:return,  -1
-         endcase
-;------------------------------------------------------------
-; Calculation of the divergence
-;------------------------------------------------------------
-         zu = temporary(u)*e2u[indice2d]*(umask())[indice3d]
-         terreu = where(zu EQ 0)
-         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
-         zv = temporary(v)*e1v[indice2d]*(vmask())[indice3d]
-         terrev = where(zv EQ 0)
-         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
-         zdiv = zu-shift(zu, 1, 0)+zv-shift(zv, 0, 1)
-         zdiv = temporary(zdiv)*tmask[indice3d]/(e1t[indice2d]*e2t[indice2d])
+      case 1 of
+        szu[1] EQ nxu AND szu[2] EQ nyu AND $
+           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 
+          case 1 of
+            nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
+            nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
+            nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
+            nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
+            ELSE :
+          endcase
+        END
+        szu[1] EQ jpi AND szu[2] EQ jpj AND $
+           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN 
+          u = u[indice2d]
+          v = v[indice2d]
+        END
+        ELSE:return, -1
+      endcase
+;------------------------------------------------------------
+; divergence computation
+;------------------------------------------------------------
+      zu = e2u[indice2d]
+      landu = where((umask())[indice2d+jpi*jpj*firstzt] EQ 0)
+      if landu[0] NE -1 then zu[temporary(landu)] = !values.f_nan
+      zu = temporary(u) * temporary(zu)
+
+      zv = e1v[indice2d]
+      landv = where((vmask())[indice2d+jpi*jpj*firstzt] EQ 0)
+      if landv[0] NE -1 then zv[temporary(landv)] = !values.f_nan
+      zv = temporary(v) * temporary(zv)
+
+      zdiv = 1.e6 / (e1t[indice2d]*e2t[indice2d])
+      zdiv = ( zu - shift(zu, 1, 0) + zv - shift(zv, 0, 1) ) * temporary(zdiv)
 ;------------------------------------------------------------
 ; Edging put at !values.f_nan
 ;------------------------------------------------------------
-         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
-            zdiv[0, *] = !values.f_nan
-            zdiv[nx-1, *] = !values.f_nan
-         endif
-         zdiv[*, 0] = !values.f_nan
-         zdiv[*, ny-1] = !values.f_nan
-;
-         zdiv = temporary(zdiv)*1e6
-;
-         if n_elements(valmask) EQ 0 THEN valmask = 1e20
-         terre =  where(tmask[indice3d] EQ 0)
-         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
-;------------------------------------------------------------
-; for the graphic drawing
-;------------------------------------------------------------
-         vargrid = 'T'
-         varname = 'div'
-         varunits = '1e6*s-1'
-         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
-         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan)
-
-      END
-   endcase
-   
-;------------------------------------------------------------
-sortie:
-   if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun 
-   
-   return, zdiv
+      if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
+        zdiv[0, *] = !values.f_nan
+        zdiv[nx-1, *] = !values.f_nan
+      endif
+      zdiv[*, 0] = !values.f_nan
+      zdiv[*, ny-1] = !values.f_nan
+;
+      land =  where(tmask[indice2d+jpi*jpj*firstzt] EQ 0)
+      if land[0] NE -1 then zdiv[temporary(land)] = valmask
+      if keyword_set(direc) then zdiv = moyenne(zdiv, direc, /nan)
+    END
+;----------------------------------------------------------------------------
+;----------------------------------------------------------------------------
+    ELSE:return, report('U and V input arrays must have 2, 3 or 4 dimensions')
+  ENDCASE   
+;------------------------------------------------------------
+  if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun 
+  
+  return, zdiv
 end