Changeset 167 for trunk/SRC/Computation/curl.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)
• curl.pro (moved) (moved from trunk/SRC/ToBeReviewed/CALCULS/curl.pro) (5 diffs)

trunk/SRC/Computation/curl.pro

@@ -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
+; with a mean done in the DIREC direction
 ;
 ; @returns RES
-; A 2d matrix
+; the vertical component of the curl 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 ae 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 (firstxf, lastxf, nxf, firstyf, lastyf, nyf).
+; - points that cannot be computed (domain bondaries, coastline) are set to NaN
+;
+; @examples
+; IDL> @tst_initorca2
+; IDL> plt, curl(dist(jpi,jpj), dist(jpi,jpj))
 ;
 ; @history 
 ; Guillaume Roullet (grlod\@ipsl.jussieu.fr)
-; 
 ; Sebastien Masson (smasson\@lodyc.jussieu.fr)
 ; adaptation to work with a reduce domain
@@ -48 +59 @@
 ; $Id$
 ;
+; @todo 
+; code the 4d case
 ;-
 ;------------------------------------------------------------
 ;------------------------------------------------------------
 ;------------------------------------------------------------
-FUNCTION curl, uu, vv
+FUNCTION curl, uu, vv, DIREC = DIREC
 ;
   compile_opt idl2, strictarrsubs
 ;
-@common
-   tempsun = systime(1)         ; To key_performance
-;
-   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)          ; To key_performance
+;
+  IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
      return, report(['This version of curl 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)
-   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 = 'F'
+  varname = 'vorticity'
+  varunits = 's-1'
+  if n_elements(valmask) EQ 0 THEN valmask = 1e20
+  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 
-         indice2d = lindgen(jpi, jpj)
-         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+    szu[0] EQ 3 AND jpt EQ 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 curl
-;------------------------------------------------------------
-         coefu = (e1u[indice2d])[*]#replicate(1, nzt)
-         coefu = reform(coefu, nx, ny, nzt, /over)
-         coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
-         terreu = where(coefu EQ 0)
-         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
-         
-         coefv = (e2v[indice2d])[*]#replicate(1, nzt)
-         coefv = reform(coefv, nx, ny, nzt, /over)
-         coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
-         terrev = where(coefv EQ 0)
-         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
-
-         tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
-         div = (e1f[indice2d]*e2f[indice2d])[*]#replicate(1, nzt)
-         div = reform(div, nx, ny, nzt, /over)
-         tabf = tabf/div
-;
-         zu = u*temporary(coefu)
-         zv = v*temporary(coefv)
-
-         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
-         psi = tabf*psi
+      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
+;------------------------------------------------------------
+; curl computation
+;------------------------------------------------------------
+      coefu = ((e1u[indice2d])[*]#replicate(1., nzt)) $
+              * (umask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
+      landu = where(coefu EQ 0)
+      if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan
+      
+      coefv = ((e2v[indice2d])[*]#replicate(1., nzt)) $
+              *(vmask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
+      landv = where(coefv EQ 0)
+      if landv[0] NE -1 then coefv[temporary(landv)] = !values.f_nan
+
+      tabf = (fmask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] $
+             / ((e1f[indice2d]*e2f[indice2d])[*]#replicate(1., nzt))
+      landf =  where(tabf EQ 0)
+;
+      zu = temporary(u) * temporary(coefu)
+      zv = temporary(v) * temporary(coefv)
+
+      psi = (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
+      psi = temporary(tabf) * temporary(psi)
 ;------------------------------------------------------------
 ; Edging put at !values.f_nan
 ;------------------------------------------------------------
-         if NOT keyword_set(key_periodic)  OR nx NE jpi then begin
-            psi[0, *, *] = !values.f_nan
-            psi[nx-1, *, *] = !values.f_nan
-         endif
-         psi[*, 0, *] = !values.f_nan
-         psi[*, ny-1, *] = !values.f_nan
-;
-         if n_elements(valmask) EQ 0 THEN valmask = 1e20
-         terref =  where(tabf EQ 0)
-         if terref[0] NE -1 then psi[temporary(terref)] = valmask
-;------------------------------------------------------------
-; For the graphic drawing
-;------------------------------------------------------------
-         domdef, (glagmt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
-         if keyword_set(direc) then psi = moyenne(psi,direc,/nan)
-
-;----------------------------------------------------------------------------
-      END
+      if NOT keyword_set(key_periodic)  OR nx NE jpi then begin
+        psi[0, *, *] = !values.f_nan
+        psi[nx-1, *, *] = !values.f_nan
+      endif
+      psi[*, 0, *] = !values.f_nan
+      psi[*, ny-1, *] = !values.f_nan
+;
+      if landf[0] NE -1 then psi[temporary(landf)] = valmask
+      if keyword_set(direc) then psi = moyenne(psi, direc, /nan)
+    END
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
@@ -162 +175 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 
-         indice2d = lindgen(jpi, jpj)
-         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
+    szu[0] EQ 3 AND jpt GT 1:BEGIN 
 ;------------------------------------------------------------
 ; extraction of U and V on the appropriated domain
 ;------------------------------------------------------------
-         case 1 of
-            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
-             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
-               if nxu NE nx then $
-                if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
-               IF nxv NE nx THEN $
-                if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
-               IF nyu NE ny THEN $
-                if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
-               IF nyv NE ny THEN $
-                if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
-            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 
-               print, 'problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs'
-               return, -1
-            end
-         endcase
-;----------------------------------------------------------------------------
-; Calculation of the curl
-;----------------------------------------------------------------------------
-         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
-         terreu = where(coefu EQ 0)
-         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
-         coefu = temporary(coefu[*])#replicate(1, jpt)
-         coefu = reform(coefu, nx, ny, jpt, /over)
-;
-         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
-         terrev = where(coefv EQ 0)
-         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
-         coefv = temporary(coefv[*])#replicate(1, jpt)
-         coefv = reform(coefv, nx, ny, jpt, /over)
-;
-         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
-         tabf = temporary(tabf[*])#replicate(1, jpt)
-         tabf = reform(tabf, nx, ny, jpt, /over)
-;------------------------------------------------------------
-; Calculation of the curl
-;------------------------------------------------------------
-         zu = u*temporary(coefu)
-         zv = v*temporary(coefv)
-;
-         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
-         psi = tabf*psi
+      case 1 of
+        szu[1] EQ nxu AND szu[2] EQ nyu AND $
+           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 
+          if nxu NE nx then $
+             if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
+          IF nxv NE nx THEN $
+             if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
+          IF nyu NE ny THEN $
+             if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
+          IF nyv NE ny THEN $
+             if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
+        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:BEGIN 
+          print, 'problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs'
+          return, -1
+        end
+      endcase
+;----------------------------------------------------------------------------
+; curl computation
+;----------------------------------------------------------------------------
+      coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
+      landu = where(coefu EQ 0)
+      if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan
+      coefu = temporary(coefu[*])#replicate(1., jpt)
+;
+      coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
+      landv = where(coefv EQ 0)
+      if landv[0] NE -1 then coefv[temporary(landv)] = !values.f_nan
+      coefv = temporary(coefv[*])#replicate(1., jpt)
+;
+      tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
+      tabf = temporary(tabf[*])#replicate(1., jpt)
+      landf = where(tabf EQ 0)
+;
+      zu = temporary(u) * temporary(coefu)
+      zv = temporary(v) * temporary(coefv)
+;
+      psi = (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
+      psi = temporary(tabf) * temporary(psi)
 ;------------------------------------------------------------
 ; extraction of U and V on the appropriated domain
 ;------------------------------------------------------------
-         if NOT keyword_set(key_periodic) OR nx NE jpi then begin
-            psi[0, *, *] = !values.f_nan
-            psi[nx-1, *, *] = !values.f_nan
-         endif
-         psi[*, 0, *] = !values.f_nan
-         psi[*, ny-1, *] = !values.f_nan
-         if n_elements(valmask) EQ 0 THEN valmask = 1e20
-         terref =  where(tabf EQ 0)
-         if terref[0] NE -1 then psi[temporary(terref)] = valmask
-;----------------------------------------------------------------------------
-         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
-         if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan)
-;----------------------------------------------------------------------------
-      END
+      if NOT keyword_set(key_periodic) OR nx NE jpi then begin
+        psi[0, *, *] = !values.f_nan
+        psi[nx-1, *, *] = !values.f_nan
+      endif
+      psi[*, 0, *] = !values.f_nan
+      psi[*, ny-1, *] = !values.f_nan
+      if landf[0] NE -1 then psi[temporary(landf)] = valmask
+      if keyword_set(direc) then psi = grossemoyenne(psi, direc, /nan)
+    END
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
@@ -238 +240 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      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 +248 @@
 ;----------------------------------------------------------------------------
 ;----------------------------------------------------------------------------
-      ELSE:BEGIN                ;xy
-         indice2d = lindgen(jpi, jpj)
-         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
-;------------------------------------------------------------
-;------------------------------------------------------------
-         case 1 of
-            (size(u))[0] NE 2 OR (size(v))[0] NE 2: 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 
-               if nxu NE nx then $
-                if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
-               IF nxv NE nx THEN $
-                if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
-               IF nyu NE ny THEN $
-                if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
-               IF nyv NE ny THEN $
-                if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
-            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 curl
-;------------------------------------------------------------
-         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
-         terreu = where(coefu EQ 0)
-         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
-         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
-         terrev = where(coefv EQ 0)
-         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
-         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
-;
-         zu = u*temporary(coefu)
-         zv = v*temporary(coefv)
-
-         psi =  (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1))
-         psi = tabf*psi
-
+    szu[0] EQ 2:BEGIN 
+;------------------------------------------------------------
+;------------------------------------------------------------
+      case 1 of
+        szu[1] EQ nxu AND szu[2] EQ nyu AND $
+           szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 
+          if nxu NE nx then $
+             if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
+          IF nxv NE nx THEN $
+             if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
+          IF nyu NE ny THEN $
+             if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
+          IF nyv NE ny THEN $
+             if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
+        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
+;------------------------------------------------------------
+; curl computation
+;------------------------------------------------------------
+      coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
+      landu = where(coefu EQ 0)
+      if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan
+      coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
+      landv = where(coefv EQ 0)
+      if landv[0] NE -1 then coefv[temporary(landv)] = !values.f_nan
+      tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
+      landf =  where(tabf EQ 0)
+;
+      zu = temporary(u) * temporary(coefu)
+      zv = temporary(v) * temporary(coefv)
+
+      psi = (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1))
+      psi = temporary(tabf) * temporary(psi)
 ;------------------------------------------------------------
 ; Edging put at !values.f_nan
 ;------------------------------------------------------------
-         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
-            psi[0, *] = !values.f_nan
-            psi[nx-1, *] = !values.f_nan
-         endif
-         psi[*, 0] = !values.f_nan
-         psi[*, ny-1] = !values.f_nan
-;
-         if n_elements(valmask) EQ 0 THEN valmask = 1e20
-         terref =  where(tabf EQ 0)
-         if terref[0] NE -1 then psi[temporary(terref)] = valmask
-;------------------------------------------------------------
-; for the graphic drawing
-;------------------------------------------------------------
-         domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
-         if keyword_set(direc) then psi = moyenne(psi,direc,/nan)
-
-      END
-;----------------------------------------------------------------------------
-   endcase
-;------------------------------------------------------------
-   if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun 
-;------------------------------------------------------------
-   vargrid = 'F'
-   varname = 'vorticity'
-   return, psi
+      if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
+        psi[0, *] = !values.f_nan
+        psi[nx-1, *] = !values.f_nan
+      endif
+      psi[*, 0] = !values.f_nan
+      psi[*, ny-1] = !values.f_nan
+;
+      if landf[0] NE -1 then psi[temporary(landf)] = valmask
+      if keyword_set(direc) then psi = moyenne(psi, 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 curl', systime(1)-tempsun 
+
+  return, psi
 end