Changeset 13

Timestamp:

04/28/06 14:18:03 (18 years ago)

Author:

pinsard

Message:

upgrade of GRILLE/Utilities according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files

Location:

trunk

Files:

• GRILLE (deleted)
• Grid/computegrid.pro (copied) (copied from trunk/GRILLE/computegrid.pro) (2 diffs)
• Grid/micromeshmask.pro (added)
• Grid/n128gaussian.pro (added)
• Grid/n160gaussian.pro (added)
• Grid/n256gaussian.pro (added)
• Grid/n48gaussian.pro (added)
• Grid/n80gaussian.pro (added)
• Grid/ncdf_meshread.pro (added)
• Grid/restoreboxparam.pro (added)
• Grid/saveboxparam.pro (added)
• Grid/smallmeshmask.pro (added)
• Obsolete/meshlec.pro (copied) (copied from trunk/GRILLE/meshlec.pro) (1 diff)
• Obsolete/ncdf_meshlec.pro (copied) (copied from trunk/GRILLE/ncdf_meshlec.pro) (1 diff)
• ToBeReviewed/GRILLE/changegrid.pro (copied) (copied from trunk/GRILLE/changegrid.pro) (2 diffs)
• ToBeReviewed/GRILLE/cmpgrid.pro (copied) (copied from trunk/GRILLE/cmpgrid.pro) (1 diff)
• ToBeReviewed/GRILLE/domdef.pro (copied) (copied from trunk/GRILLE/domdef.pro) (4 diffs)
• ToBeReviewed/GRILLE/f2v.pro (copied) (copied from trunk/GRILLE/f2v.pro) (6 diffs)
• ToBeReviewed/GRILLE/fmask.pro (copied) (copied from trunk/GRILLE/fmask.pro) (1 diff)
• ToBeReviewed/GRILLE/grille.pro (copied) (copied from trunk/GRILLE/grille.pro) (5 diffs)
• ToBeReviewed/GRILLE/t2v.pro (copied) (copied from trunk/GRILLE/t2v.pro) (7 diffs)
• ToBeReviewed/GRILLE/tracegrille.pro (copied) (copied from trunk/GRILLE/tracegrille.pro) (5 diffs)
• ToBeReviewed/GRILLE/u2t.pro (copied) (copied from trunk/GRILLE/u2t.pro) (6 diffs)
• ToBeReviewed/GRILLE/umask.pro (copied) (copied from trunk/GRILLE/umask.pro) (1 diff)
• ToBeReviewed/GRILLE/v2t.pro (copied) (copied from trunk/GRILLE/v2t.pro) (7 diffs)
• ToBeReviewed/GRILLE/vmask.pro (copied) (copied from trunk/GRILLE/vmask.pro) (1 diff)
• ToBeReviewed/GRILLE/whichgrid.pro (copied) (copied from trunk/GRILLE/whichgrid.pro) (12 diffs)

trunk/Grid/computegrid.pro

@@ -5 +5 @@
 ; NAME:computegrid
 ;
-; PURPOSE:calcule les parametres (necessaires aux traces) d''une
-; grille reguliere avec 
-;          - les points T, U et V sont confondus
-;          - un seul niveau vertical.
-;
-; CATEGORY:grille
-;
-; CALLING SEQUENCE:computegrid, departx, departy, pasx, pasy, nx, ny
+; PURPOSE:compute the grid parameters from cm_4mesh common:
+;
+; horizontal parameters:
+;       glam[tf], gphi[tf], e1t and e2t
+;    and if FULLCGRID keyword is defined:
+;       glam[uv], gphi[uv], e1[uvf] and e2[uvf]
+;      
+; verticals parameters:
+;       gdep[tw], e3[tw]
+;
+; masks: tmask
+;        and if FULLCGRID keyword is defined:[uv]maskred fmaskred[xy]
+;
+; triangulation: triangles_list
+;
+; key_ parameters:
+;       key_shift, key_periodic, key_zreverse, key_yreverse,
+;       key_stride, key_onearth, key_partialstep
+;
+; CATEGORY:grid
+;
+; CALLING SEQUENCE:
+;
+;       computegrid, startx, starty, stepx, stepy, nx, ny
+;       computegrid, startx, starty, stepx, stepy
+;       computegrid, xaxis = xaxis, yaxis = yaxis
+;       or a suitable mix...
 ; 
-; INPUTS:des saclaires.
-;       departx:la longitude minimum
-;       departy:la latitude minimum
-;       pasx:la pas suivant la longitude
-;       pasy:la pas suivant la latitude
-;       nx:le nombre de points en longitude
-;       ny:le nombre de points en latitude
+; INPUTS: 
+;       startx:scalar, x starting point
+;       starty:scalar, y starting point
+;       stepx:scalar or vector: x direction step, must be > 0
+;             if vector nx is not used
+;       stepy:scalar or vector: y direction step, 
+;             could be > 0 (south to north) or < 0 (north to south)
+;             if vector ny is not used
+;       nx:scalar, number of points in x direction 
+;       ny:scalar, number of points in y direction
 ;
 ; KEYWORD PARAMETERS:
 ;
+;       /FULLCGRID: activate to specify that you want to compute
+;       all the paremeters of a C grid. Computation of glam[uv],
+;       gphi[uv], e1[uvf], e2[uvf], [uv]maskred and fmaskred[xy]
+;       will be add to the default computations
+;
+;       GLAMBOUNDARY: a 2 elements vector, [lon1,lon2], the longitute
+;       boundaries that should be used to visualize the data.
+;       we must have lon2 > lon1 and lon2 - lon1 le 360
+;       key_shift will be defined automaticaly computed according to
+;       glamboundary by using the FIRST LINE of glamt but 
+;       key_shift will /= 0 only if key_periodic = 1 
+;
+;       MASK: to specify the mask with a 2 or 3 dimension array
+;
+;       ONEARTH = 0 or 1: to force the manual definition of
+;       key_onearth (to specify if the data are on earth -> use longitude
+;       /latitude etc...). By default, key_onearth = 1.
+;       note that ONEARTH = 0 forces PERIODIC = 0, SHIFT = 0, 
+;       and is cancelling GLAMBOUNDARY
+;
+;       PERIODIC = 0 or 1: to force the manual definition of
+;       key_periodic. By default, key_periodic is automaticaly
+;       computed by using the first line of glamt. 
+;
+;       SHIFT = scalar to force the manual definition of key_shift. By
+;       debault, key_shift is automaticaly computed according to
+;       glamboundary (when defined) by using the FIRST LINE of glamt. if
+;       key_periodic=0 then in any case key_shift = 0. 
+;
+;       STRIDE = : a 3 elements vector to specify the stride in x, y, z
+;       direction. Default definition is [1, 1, 1]. The resulting value
+;       will be stored in the common (cm_4mesh) variable key_stride
+;
+;       XAXIS: to specify longitude1 with a 1 or 2 dimension array, in 
+;       this case startx, stepx and nx are not used but could be
+;       necessary if the y axis is not defined with yaxis. It must be
+;       possible to sort the first line of xaxis in the increasing
+;       order by shifting its elements.
+;
+;       YAXIS: to specify latitudes with a 1 or 2 dimension array, in 
+;       this case starty, stepy and ny are not used but starty and
+;       stepy could be necessary if the x axis is not defined with xaxis.
+;       It must be sorted in the increasing or deceasing order
+;       (along each column if 2d array).
+;
+;       /XYINDEX: activate to specify that the horizontal grid should 
+;       be simply defined by using the index of the points
+;          (xaxis = findgen(nx) and yaxis = findgen(ny))
+;       using this keyword forces key_onearth=0
+;
+;       [XYZ]MINMESH: to define the common variables i[xyz]minmesh
+;       used to define the grid only in a zoomed part of the original
+;       grid. Defaut values are 0L, max value is [XYZ]MAXMESH
+;
+;       [XYZ]MAXMESH: to define the common variables i[xyz]maxmesh
+;       used to define the grid only in a zoomed part of the original
+;       grid. Defaut values are jp[ijk]glo-1, max value is
+;       jp[ijk]glo-1. if [XYZ]MAXMESH is negative, then we define
+;       i[xyz]maxmesh as jp[ijk]glo - 1 + [XYZ]MAXMESH instead of
+;       [XYZ]MAXMESH    
+;
+;       ZAXIS: to specify the vertical axis with a 1 dimension
+;       array. Must be sorted in the increasing or deceasing order
+;
 ; OUTPUTS:
 ;
-; COMMON BLOCKS:common.pro
-;
-; SIDE EFFECTS:definit aussi les parametres pour les grilles U et V
-; meme si on n''en a pas besoin.
-;
-; RESTRICTIONS:programme vite fait, avec les latitudes qui doivent
-; aller du sud vers le nord.
+; COMMON BLOCKS: cm_4mesh cm_4data cm_4cal
+;
+; SIDE EFFECTS:
+;
+; RESTRICTIONS:FUV points definition...
 ;
 ; EXAMPLE:
@@ -38 +122 @@
 ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
 ;                      2000-04-20
+;  Sept 2004, several bug fixs to suit C grid type...
+;  Aug 2005, rewritte almost everything...
 ;-
 ;------------------------------------------------------------
 ;------------------------------------------------------------
 ;------------------------------------------------------------
-PRO computegrid, departx, departy, pasx, pasy, nx, ny
-@common
+PRO computegrid, startx, starty, stepxin, stepyin, nxin, nyin $
+                 , XAXIS = xaxis, YAXIS = yaxis, ZAXIS = zaxis $
+                 , MASK = mask, GLAMBOUNDARY = glamboundary $
+                 , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
+                 , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
+                 , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
+                 , ONEARTH = onearth, PERIODIC = periodic $
+                 , SHIFT = shift, STRIDE = stride $
+                 , FULLCGRID = fullcgrid, XYINDEX = xyindex $
+                 , FBASE2TBASE = fbase2tbase, _extra = ex 
+;---------------------------------------------------------
+@cm_4mesh
+@cm_4data
+@cm_4cal
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+@updatekwd
+  ENDIF
+;---------------------------------------------------------
 ;------------------------------------------------------------
-   jpiglo = long(nx)
-   jpjglo = long(ny)
-   jpkglo = 1l
-;
-   jpi = jpiglo
-   jpj = jpjglo
-   jpk = jpkglo
-;
-   ixminmesh  =0l
-   ixmaxmesh  =long(jpi-1)
-   iyminmesh  =0l
-   iymaxmesh  =long(jpj-1)
-   izminmesh  =0l
-   izmaxmesh  =long(jpk-1)
-;
-   glamt = departx+findgen(jpi)*pasx
-   glamt = glamt#replicate(1, jpj)
-   glamf = glamt+pasx/2.
-   gphit = departy+findgen(jpj)*pasy
-   gphit = replicate(1, jpi)#gphit
-   gphif = gphit+pasy/2.
-   glamu = glamf
-   glamv = glamt
-   gphiu = gphit
-   gphiv = gphif
-   gdept = 5.
-   gdepw = 5.
-;
-   pasx = abs(pasx)
-   pasy = abs(pasy)
-; facteurs d'echelle:
-   r = 6371000. 
-   rprojete = r*cos(gphit[0,*]*!pi/180.)
-   e1t = rprojete*pasx*!pi/180.
-   e1t = reform(replicate(1, jpi)#e1t[*], jpi, jpj)
-   e2t = replicate(r*pasy*!pi/180., jpi, jpj)
-;
-   e1u = e1t
-   e2u = e2t
-;
-   rprojete = r*cos(gphiv[0,*]*!pi/180.)
-   e1v = rprojete*pasx*!pi/180.
-   e1v = reform(replicate(1, jpi)#e1v[*], jpi, jpj)
-   e2v = e2t
-;
-   e1f = e1v
-   e2f = e2t
-;
-   e3t = 10.
-   e3w = 10.
-; mask par defauyt a 1
-   tmask = replicate(1, jpi, jpj)
-   umaskred = replicate(1, jpj)
-   vmaskred = replicate(1, jpi)
-   fmaskredy = replicate(1, jpj)
-   fmaskredx = replicate(1, jpi)
+  time1 = systime(1)            ; for key_performance
+;------------------------------------------------------------
+;
+;====================================================
+; Check input parameters
+;====================================================
+;
+; xaxis related parameters
+;
+  if n_elements(xaxis) NE 0 then BEGIN 
+    CASE (size(xaxis))[0] OF
+      0:nx = 1L
+      1:nx = (size(xaxis))[1]
+      2:nx = (size(xaxis))[1]
+    ENDCASE
+  ENDIF ELSE BEGIN 
+    IF n_elements(startx) EQ 0 THEN BEGIN 
+      dummy = report('If xaxis is not given, startx must be defined')
+      return
+    ENDIF
+    CASE n_elements(stepxin) OF
+      0:BEGIN 
+        dummy = report('If xaxis is not given, stepxin must be defined')
+        return
+      END
+      1:BEGIN 
+        IF n_elements(nxin) EQ 0 THEN BEGIN 
+          dummy = report('If xaxis is not given and stepxin has only one element, nx must be defined')
+          return
+        ENDIF ELSE nx = nxin
+      END
+      ELSE:nx = n_elements(stepxin)    
+    ENDCASE
+  ENDELSE 
+;
+; yaxis related parameters
+;
+  if n_elements(yaxis) NE 0 then BEGIN 
+    CASE (size(yaxis))[0] OF
+      0:ny = 1L
+      1:ny = (size(yaxis))[1]
+      2:ny = (size(yaxis))[2]
+    ENDCASE
+  ENDIF ELSE BEGIN 
+    IF n_elements(starty) EQ 0 THEN BEGIN 
+      dummy = report('If yaxis is not given, starty must be defined')
+      return
+    ENDIF
+    CASE n_elements(stepyin) OF
+      0:BEGIN 
+        dummy = report('If yaxis is not given, stepyin must be defined')
+        return
+      END
+      1:BEGIN 
+        IF n_elements(nyin) EQ 0 THEN BEGIN 
+          dummy = report('If yaxis is not given and stepyin has only one element, ny must be defined')
+          return
+        ENDIF ELSE ny = nyin
+      END
+      ELSE:ny = n_elements(stepyin)    
+    ENDCASE
+  ENDELSE
+;
+; zaxis related parameters
+;
+  if n_elements(zaxis) NE 0 then BEGIN 
+    CASE (size(zaxis))[0] OF
+      0:nz = 1L
+      1:nz = (size(zaxis))[1]
+      ELSE:BEGIN
+        print, 'not coded'
+        stop
+      END
+    ENDCASE
+  ENDIF ELSE nz = 1L
+;
+;====================================================
+; Others automatic definitions...
+;====================================================
+;
+  jpiglo = long(nx)
+  jpjglo = long(ny)
+  jpkglo = long(nz)
+;
+  IF n_elements(xminmesh) NE 0 THEN ixminmesh = long(xminmesh[0]) ELSE ixminmesh  = 0l
+  IF n_elements(xmaxmesh) NE 0 THEN ixmaxmesh = long(xmaxmesh[0]) ELSE ixmaxmesh  = jpiglo-1
+  IF n_elements(yminmesh) NE 0 THEN iyminmesh = long(yminmesh[0]) ELSE iyminmesh  = 0l
+  IF n_elements(ymaxmesh) NE 0 THEN iymaxmesh = long(ymaxmesh[0]) ELSE iymaxmesh  = jpjglo-1
+  IF n_elements(zminmesh) NE 0 THEN izminmesh = long(zminmesh[0]) ELSE izminmesh  = 0l
+  IF n_elements(zmaxmesh) NE 0 THEN izmaxmesh = long(zmaxmesh[0]) ELSE izmaxmesh  = jpkglo-1
+;
+  iymaxmesh = iymaxmesh-keyword_set(fbase2tbase) 
+;
+  IF ixmaxmesh LT 0 THEN ixmaxmesh = jpiglo -1 + ixmaxmesh
+  IF iymaxmesh LT 0 THEN iymaxmesh = jpjglo -1 + iymaxmesh
+  IF izmaxmesh LT 0 THEN izmaxmesh = jpkglo -1 + izmaxmesh
+; avoid basics errors...
+  ixmaxmesh = 0 > ixmaxmesh < (jpiglo-1)
+  ixminmesh = 0 > ixminmesh < ixmaxmesh
+  iymaxmesh = 0 > iymaxmesh < (jpjglo-1)
+  iyminmesh = 0 > iyminmesh < iymaxmesh
+  izmaxmesh = 0 > izmaxmesh < (jpkglo-1)
+  izminmesh = 0 > izminmesh < izmaxmesh
+;
+  jpi = ixmaxmesh-ixminmesh+1
+  jpj = iymaxmesh-iyminmesh+1
+  jpk = izmaxmesh-izminmesh+1
+;
+  jpidta = jpiglo
+  jpjdta = jpjglo
+  jpkdta = jpkglo
+  ixmindta = 0
+  ixmaxdta = jpidta-1
+  iymindta = 0
+  iymaxdta = jpjdta-1
+  izmindta = 0
+  izmaxdta = jpkdta-1
+;
+  key_partialstep = 0
+  if n_elements(stride) eq 3 then key_stride = stride $
+  ELSE key_stride = [1, 1, 1]
+  key_gridtype = 'c'
+;
+; check xyindex and its consequences
+;
+  if keyword_set(xyindex) then onearth = 0
+;
+; check onearth and its consequences
+; 
+  IF n_elements(onearth) EQ 0 THEN key_onearth = 1b $
+  ELSE key_onearth = keyword_set(onearth) 
+  IF NOT key_onearth THEN BEGIN
+    periodic = 0
+    shift = 0
+  ENDIF;
 
+  r = 6371000.
+;
+;====================================================
+; X direction : glamt
+;====================================================
+;
+; def of glamt
+;
+  if n_elements(xaxis) NE 0 then BEGIN 
+    if keyword_set(xyindex) THEN glamt = findgen(jpiglo) ELSE glamt = xaxis
+  ENDIF ELSE BEGIN
+    if keyword_set(xyindex) THEN stepx = 1. ELSE stepx = stepxin
+    CASE 1 OF
+      n_elements(stepx):glamt = startx + findgen(jpiglo)*stepx
+      size(stepx, /n_dimensions):glamt = startx + total(stepx, /cumulative)
+      ELSE:BEGIN 
+        dummy = report('Wrong definition of stepx...')
+        return
+      END 
+    ENDCASE
+  ENDELSE
+;
+; apply glamboundary
+;
+  IF keyword_set(glamboundary) AND key_onearth THEN BEGIN
+    IF glamboundary[0] GE glamboundary[1] THEN stop
+    IF glamboundary[1]-glamboundary[0] GT 360 THEN stop
+    glamt = glamt MOD 360
+    smaller = where(glamt LT glamboundary[0])
+    if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360.
+    bigger = where(glamt GE glamboundary[1])
+    if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360.
+  ENDIF
+;
+; force glamt to have 2 dimensions
+;
+  CASE size(reform(glamt), /n_dimensions) OF
+    0:glamt = replicate(glamt, jpi, jpj)
+    1:glamt = glamt[ixminmesh:ixmaxmesh]#replicate(1, jpj)
+    2:glamt = glamt[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh]
+  ENDCASE
+; keep 2d array even with degenereted dimension
+  IF jpj EQ 1 THEN glamt = reform(glamt, jpi, jpj, /over)
+;
+;====================================================
+; Y direction : gphit
+;====================================================
+;
+; def of gphit
+;
+  if n_elements(yaxis) NE 0 THEN BEGIN
+    if keyword_set(xyindex) THEN gphit = findgen(jpjglo) ELSE gphit = yaxis
+  ENDIF ELSE BEGIN
+    if keyword_set(xyindex) THEN stepy = 1. ELSE stepy = stepyin
+    CASE 1 OF
+      n_elements(stepy):gphit = starty + findgen(jpjglo)*stepy
+      size(stepy, /n_dimensions):gphit = starty + total(stepy, /cumulative)
+      ELSE:BEGIN 
+        dummy = report('Wrong definition of stepy...')
+        return
+      END 
+    ENDCASE
+  ENDELSE
+;
+; force gphit to have 2 dimensions
+;
+  CASE size(reform(gphit), /n_dimensions) OF
+    0:gphit = replicate(gphit, jpi, jpj)
+    1:gphit = replicate(1, jpi)#gphit[iyminmesh:iymaxmesh]
+    2:gphit = gphit[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh]
+  ENDCASE
+; keep 2d array even with degenereted dimension
+  IF jpj EQ 1 THEN gphit = reform(gphit, jpi, jpj, /over) 
+;
+;====================================================
+; def of key_shift
+;====================================================
+;
+; definition of key_shift by shiftting the array to have the min
+; values of glamt[*, 0] in glamt[0, 0]
+;
+  IF n_elements(shift) EQ 0 THEN BEGIN
+    IF jpi GT 1 then BEGIN
+      xtest = glamt[*, 0]
+      key_shift = (where(xtest EQ min(xtest)))[0]
+      IF key_shift NE 0 THEN key_shift = jpi - key_shift
+    ENDIF ELSE key_shift = 0
+  ENDIF ELSE key_shift = shift
+;
+;====================================================
+; def of key_periodic
+;====================================================
+;
+  IF n_elements(periodic) EQ 0 THEN BEGIN
+    IF jpi GT 1 THEN BEGIN
+      xtest = shift(glamt[*, 0], key_shift)
+; check that xtest is now sorted in the increasing order
+      IF array_equal(sort(xtest), lindgen(jpi)) EQ 0 THEN BEGIN
+        print, 'WARNING: we cannot sort the xaxis with a simple shift...'
+        print, 'we force key_periodic = 0 and key_shift = 0'
+        print, 'only horizontal plot may be ok...' 
+        key_periodic = 0
+        xnotsorted = 1
+      ENDIF ELSE BEGIN 
+        key_periodic = (xtest[jpi-1]+2*(xtest[jpi-1]-xtest[jpi-2])) $
+                       GE (xtest[0]+360)
+      ENDELSE
+    ENDIF ELSE key_periodic = 0
+  ENDIF ELSE key_periodic = keyword_set(periodic) 
+;
+; update key_shift
+;
+  key_shift = key_shift * (key_periodic EQ 1)
+;
+  IF NOT keyword_set(key_periodic) AND keyword_set(fbase2tbase) THEN BEGIN
+    ixmaxmesh = ixmaxmesh-1
+    jpi = jpi-1
+  ENDIF
+;
+;====================================================
+; apply key_shift
+;====================================================
+;
+  if keyword_set(key_shift) then BEGIN 
+    glamt = shift(glamt, key_shift, 0)
+    gphit = shift(gphit, key_shift, 0)
+    IF jpj EQ 1 THEN BEGIN 
+      glamt = reform(glamt, jpi, jpj, /over)
+      gphit = reform(gphit, jpi, jpj, /over)
+    ENDIF
+  ENDIF
+;
+;====================================================
+; def key_yreverse
+;====================================================
+;
+  IF jpj GT 1 THEN BEGIN
+    if gphit[0, 1] LT gphit[0, 0] then begin
+      key_yreverse = 1
+      gphit = reverse(gphit, 2)
+      glamt = reverse(glamt, 2)
+    ENDIF ELSE key_yreverse = 0
+  ENDIF ELSE key_yreverse = 0
+;
+;====================================================
+; Are we using a "regular" grid (that can be described
+; with x vector and y vector)?
+;====================================================
+;
+; to get faster, we first test the most basic cases before
+; testing the full array.
+;
+  CASE 1 OF
+    keyword_set(xyindex):key_irregular = 0b
+    jpi EQ 1 OR jpj EQ 1:key_irregular = 0b
+    n_elements(xaxis) EQ 0 AND n_elements(yaxis) EQ 0:key_irregular = 0b
+    size(reform(xaxis), /n_dimensions) EQ 1 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b
+    n_elements(xaxis) EQ 0 AND size(reform(yaxis), /n_dimensions) EQ 1:key_irregular = 0b
+    n_elements(yaxis) EQ 0 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b
+    array_equal(glamt[*, 0], glamt[*, jpj-1]) EQ 0:key_irregular = 1b
+    array_equal(gphit[0, *], gphit[jpi-1, *]) EQ 0:key_irregular = 1b
+    array_equal(glamt, glamt[*, 0]#replicate(1, jpj)) EQ 0:key_irregular = 1b
+    array_equal(gphit, replicate(1, jpi)#(gphit[0, *])[*]) EQ 0:key_irregular = 1b
+    ELSE:key_irregular = 0b
+  ENDCASE
+;
+;====================================================
+; def of glamf: defined as the middle of T(i,j) T(i+1,j+1)
+;====================================================
+;
+  IF jpi GT 1 THEN BEGIN 
+; we must compute stepxf: x distance between T(i,j) T(i+1,j+1)
+    CASE 1 OF
+      n_elements(stepx):stepxf = stepx
+      size(stepx, /n_dimensions):stepxf = stepx#replicate(1, jpj)
+      ELSE:BEGIN
+        if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $
+          OR (keyword_set(key_periodic) AND key_irregular) then BEGIN
+          stepxf = (glamt + 720) MOD 360
+          IF jpj EQ 1 THEN stepxf = reform(stepxf, jpi, jpj, /over)
+          stepxf = shift(stepxf, -1, -1) - stepxf
+          stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ]
+          stepxf = min(abs(stepxf), dimension = 3)
+          IF NOT keyword_set(key_periodic) THEN $
+            stepxf[jpi-1, *] = stepxf[jpi-2, *]
+        ENDIF ELSE BEGIN
+          stepxf = shift(glamt, -1, -1) - glamt
+          IF keyword_set(key_periodic) THEN $
+            stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $
+          ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *]
+        ENDELSE
+        IF jpj GT 1 THEN BEGIN 
+          stepxf[*, jpj-1] = stepxf[*, jpj-2]
+          stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2]
+        ENDIF
+      END
+    ENDCASE
+    glamf = glamt + 0.5 * stepxf
+    IF jpj EQ 1 THEN glamf = reform(glamf, jpi, jpj, /over)
+  ENDIF ELSE glamf = glamt + 0.5
+  IF keyword_set(glamboundary) AND key_onearth THEN $
+    glamf = glamboundary[0] > temporary(glamf) < glamboundary[1]
+;
+;====================================================
+; def of gphif: defined as the middle of T(i,j) T(i+1,j+1)
+;====================================================
+;
+  IF jpj GT 1 THEN BEGIN 
+; we must compute stepyf: y distance between T(i,j) T(i+1,j+1)
+    CASE 1 OF
+      n_elements(stepy):stepyf = stepy
+      size(stepy, /n_dimensions):stepyf = replicate(1, jpi)#stepy
+      ELSE:BEGIN
+        stepyf = shift(gphit, -1, -1) - gphit
+        stepyf[*, jpj-1] = stepyf[*, jpj-2]
+        IF jpi GT 1 THEN BEGIN
+          if NOT keyword_set(key_periodic) THEN $
+            stepyf[jpi-1, *] = stepyf[jpi-2, *]
+          stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2]
+        ENDIF 
+      END 
+    ENDCASE
+    gphif = gphit + 0.5 * stepyf
+  ENDIF ELSE gphif = gphit + 0.5
+  IF key_onearth THEN gphif = -90. > gphif < 90.
+;
+;====================================================
+; e1t: x distance between U(i-1,j) and U(i,j)
+;====================================================
+;
+; *-|-*---|---*---|
+;
+  IF jpi GT 1 THEN BEGIN
+    IF n_elements(stepx) NE 1 THEN BEGIN
+      IF keyword_set(irregular) THEN BEGIN
+; we must compute stepxu: x distance between T(i,j) T(i+1,j)
+        IF keyword_set(key_periodic) THEN BEGIN
+          stepxu = (glamt + 720) MOD 360
+          stepxu = shift(stepxu, -1, 0) - stepxu
+          stepxu = [ [[stepxu]], [[stepxu + 360]], [[stepxu - 360]] ]
+          stepxu = min(abs(stepxu), dimension = 3)
+        ENDIF ELSE BEGIN 
+          stepxu = shift(glamt, -1, 0) - glamt
+          stepxu[jpi-1, *] = stepxf[jpi-2, *]
+        ENDELSE 
+      ENDIF ELSE stepxu = stepxf
+      IF jpj EQ 1 THEN stepxu = reform(stepxu, jpi, jpj, /over)
+      e1t = 0.5*(stepxu+shift(stepxu, 1, 0))
+      IF NOT keyword_set(key_periodic) THEN $
+        e1t[0, *] = e1t[1, *] 
+      IF jpj EQ 1 THEN e1t = reform(e1t, jpi, jpj, /over)
+    ENDIF ELSE e1t = replicate(stepx, jpi, jpj)
+  ENDIF ELSE e1t = replicate(1b, jpi, jpj)
+;
+;====================================================
+; e2t: y distance between V(i,j-1) and V(i,j)
+;====================================================
+;
+  IF jpj GT 1 THEN BEGIN
+; we must compute stepyv: y distance between T(i,j) T(i,j+1)
+    IF n_elements(stepy) NE 1 THEN BEGIN
+      IF keyword_set(key_irregular) THEN BEGIN
+        stepyv = shift(gphit, 0, -1) - gphit
+        stepyv[*, jpj-1] = stepyv[*, jpj-2]
+      ENDIF ELSE stepyv = stepyf
+      e2t = 0.5*(stepyv+shift(stepyv, 0, 1))
+      e2t[*, 0] = e2t[*, 1]
+    ENDIF ELSE e2t = replicate(stepy, jpi, jpj)
+  ENDIF ELSE e2t = replicate(1b, jpi, jpj)
+;
+  IF key_onearth THEN e2t = r * !pi/180. * temporary(e2t)
+;
+;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  IF keyword_set(fullcgrid) THEN BEGIN
+;
+;====================================================
+; def of glamu: defined as the middle of T(i,j) T(i+1,j)
+;====================================================
+;
+    IF keyword_set(irregular) THEN BEGIN
+      glamu = glamt + 0.5 * stepxu
+      IF keyword_set(glamboundary) AND key_onearth THEN $
+        glamu = glamboundary[0] > temporary(glamu) < glamboundary[1]
+    ENDIF ELSE glamu = glamf
+;
+;====================================================
+; def of gphiu: defined as the middle of T(i,j) T(i+1,j)
+;====================================================
+;
+    IF jpi GT 1 THEN BEGIN
+ ; we must compute stepyu: y distance between T(i+1,j) T(i,j)
+      IF keyword_set(key_irregular) THEN BEGIN
+       stepyu = shift(gphit, -1, 0) - gphit
+        IF NOT keyword_set(key_periodic) THEN $
+          stepyu[jpi-1, *] = stepyu[jpi-2, *]
+        gphiu = gphit + 0.5 * stepyu
+        IF jpj EQ 1 THEN gphiu = reform(gphiu, jpi, jpj, /over)
+      ENDIF ELSE gphiu = gphit
+    ENDIF ELSE gphiu = gphit
+  IF key_onearth THEN gphiu = -90. > gphiu < 90.
+;
+;====================================================
+; def of glamv: defined as the middle of T(i,j) T(i,j+1)
+;====================================================
+;
+    IF jpj GT 1 THEN BEGIN
+ ; we must compute stepxv: x distance between T(i,j) T(i,j+1)
+      IF keyword_set(irregular) THEN BEGIN
+        IF keyword_set(key_periodic) THEN BEGIN
+          stepxv = (glamt + 720) MOD 360
+          stepxv = shift(stepxv, 0, -1) - stepxv
+          stepxv = [ [[stepxv]], [[stepxv + 360]], [[stepxv - 360]] ]
+          stepxv = min(abs(stepxv), dimension = 3)
+        ENDIF ELSE stepxv = shift(glamt, 0, -1) - glamt
+        stepxv[*, jpj-1] = stepxv[*, jpj-2]
+        glamv = glamt + 0.5 * stepxv
+        IF keyword_set(glamboundary) AND key_onearth THEN $
+          glamv = glamboundary[0] > temporary(glamv) < glamboundary[1]
+      ENDIF ELSE glamv = glamt
+    ENDIF ELSE glamv = glamt
+;
+;====================================================
+; def of gphiv: defined as the middle of T(i,j) T(i,j+1)
+;====================================================
+;
+    IF keyword_set(key_irregular) THEN $        
+      gphiv = gphit + 0.5 * stepyv $
+    ELSE gphiv = gphif
+    IF key_onearth THEN gphiv = -90. > gphiv < 90.
+;
+;====================================================
+; e1u: x distance between T(i,j) and T(i+1,j)
+;====================================================
+;
+    IF jpi GT 1 AND n_elements(stepx) NE 1 THEN $
+      e1u = stepxu ELSE e1u = e1t
+;
+;====================================================
+; e2u: y distance between F(i,j-1) and F(i,j)
+;====================================================
+;
+    IF keyword_set(key_irregular) THEN BEGIN
+      e2u = gphif - shift(gphif, 0, 1) 
+      e2u[*, 0] = e2u[*, 1]
+      IF key_onearth THEN e2u = r * !pi/180. * temporary(e2u)
+    ENDIF ELSE e2u = e2t
+;
+;====================================================
+; e1v: x distance between F(i-1,j) and F(i,j)
+;====================================================
+;
+    IF keyword_set(irregular) THEN BEGIN
+      IF keyword_set(key_periodic) THEN BEGIN
+        e1v = (glamf + 720) MOD 360
+        e1v = e1v - shift(e1v, 1, 0) 
+        e1v = [ [[e1v]], [[e1v + 360]], [[e1v - 360]] ]
+        e1v = min(abs(e1v), dimension = 3)
+      ENDIF ELSE BEGIN 
+        e1v = glamf - shift(glamf, 1, 0)
+        e1v[0, *] = stepxf[1, *]
+      ENDELSE 
+    ENDIF ELSE e1v = e1t 
+;
+;====================================================
+; e2v: y distance between T(i,j) and T(i+1,j)
+;====================================================
+;
+    IF jpj GT 1 and n_elements(stepy) NE 1 THEN BEGIN
+      e2v = stepyv
+      IF key_onearth THEN e2v = r * !pi/180. * temporary(e2v)
+    ENDIF ELSE e2v = e2t
+;
+;====================================================
+; e1f: x distance between V(i,j) and V(i+1,j)
+;====================================================
+;
+    IF keyword_set(irregular) THEN BEGIN
+      IF keyword_set(key_periodic) THEN BEGIN
+        e1f = (glamv + 720) MOD 360
+        e1f = shift(e1f, -1, 0) - e1f
+        e1f = [ [[e1f]], [[e1f + 360]], [[e1f - 360]] ]
+        e1f = min(abs(e1f), dimension = 3)
+      ENDIF ELSE BEGIN 
+        e1f = shift(glamv, -1, 0) - glamt
+        e1f[jpi-1, *] = stepxf[jpi-2, *]
+      ENDELSE 
+    ENDIF ELSE e1f = e1u
+;
+;====================================================
+; e2f: y distance between U(i,j) and U(i,j+1)
+;====================================================
+;
+    IF keyword_set(key_irregular) THEN BEGIN
+      e2f = shift(gphiu, 0, -1) - gphiu
+      e2f[*, jpj-1] = e2f[*, jpj-2]
+      IF key_onearth THEN e2f = r * !pi/180. * temporary(e2f)
+    ENDIF ELSE e2f = e2v
+;
+  ENDIF
+;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+;
+;
+;====================================================
+; e1[tuvf] from degree to meters
+;====================================================
+;
+  IF keyword_set(key_onearth)  THEN BEGIN
+    e1t = r * !pi/180. * temporary(e1t) * cos(!pi/180.*gphit)
+    IF keyword_set(fullcgrid) THEN BEGIN
+      e1u = r * !pi/180. * temporary(e1u) * cos(!pi/180.*gphiu)
+      e1v = r * !pi/180. * temporary(e1v) * cos(!pi/180.*gphiv)
+      e1f = r * !pi/180. * temporary(e1f) * cos(!pi/180.*gphif)
+    ENDIF
+  ENDIF
+;
+;====================================================
+; if not fullcgrid: make sure we don't use glam[uv], gphi[uv], e[12][uvf]
+;====================================================
+;
+  IF NOT keyword_set(fullcgrid) THEN BEGIN
+    glamu = !values.f_nan & glamv = !values.f_nan
+    gphiu = !values.f_nan & gphiv = !values.f_nan
+    e1u = !values.f_nan & e1v = !values.f_nan & e1f = !values.f_nan
+    e2u = !values.f_nan & e2v = !values.f_nan & e2f = !values.f_nan
+    firstxu = !values.f_nan & lastxu = !values.f_nan & nxu = !values.f_nan  
+    firstyu = !values.f_nan & lastyu = !values.f_nan & nyu = !values.f_nan  
+    firstxv = !values.f_nan & lastxv = !values.f_nan & nxv = !values.f_nan  
+    firstyv = !values.f_nan & lastyv = !values.f_nan & nyv = !values.f_nan  
+  ENDIF
+;
+;====================================================
+; Z direction
+;====================================================
+;
+; z axis   
+;
+    CASE n_elements(zaxis) OF
+      0:BEGIN
+        gdept = 0.
+        key_zreverse = 0
+      END
+      1:BEGIN
+        gdept = zaxis
+        key_zreverse = 0
+      END
+      ELSE:BEGIN
+        gdept = zaxis[izminmesh:izmaxmesh] 
+        IF jpk GT 1 THEN BEGIN
+          if gdept[0] GT gdept[1] then begin
+            gdept = reverse(gdept)
+            key_zreverse = 1
+          ENDIF ELSE key_zreverse = 0
+        ENDIF ELSE key_zreverse = 0
+      END
+    ENDCASE
+;
+    if n_elements(gdept) GT 1 then BEGIN
+      stepz = shift(gdept, -1)-gdept
+      stepz[jpk-1] = stepz[jpk-2]
+      gdepw = 0. > (gdept-stepz/2.)
+    ENDIF ELSE BEGIN
+      stepz = 1.
+      gdepw = gdept
+    ENDELSE
+;
+;====================================================
+; e3[tw]:
+;====================================================
+;
+    e3t = stepz
+    IF n_elements(stepz) GT 1 THEN BEGIN
+      e3w = 0.5*(stepz+shift(stepz, 1))
+      e3w[0] = 0.5*e3t[0]
+    ENDIF ELSE e3w = e3t
+;
+;====================================================
+; Mask
+;====================================================
+;
+; defaut mask eq 1
+  if NOT keyword_set(mask) then mask = -1
+;
+  if mask[0] NE -1 then BEGIN
+    tmask = byte(mask[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh, izminmesh:izmaxmesh])
+    tmask = reform(tmask, jpi, jpj, jpk, /over)
+    if key_shift NE 0 then tmask = shift(tmask, key_shift, 0, 0)
+; because tmask = reverse(tmask, 2) is not working if the 3rd
+; dimension of tmask = 1, we call reform.
+    IF jpk EQ 1 THEN tmask = reform(tmask, /over)
+    IF key_yreverse EQ 1 THEN tmask = reverse(tmask, 2)
+    IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
+    IF key_zreverse EQ 1 THEN tmask = reverse(tmask, 3)
+    IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
+    IF keyword_set(fullcgrid) THEN BEGIN
+      IF keyword_set(key_periodic) THEN BEGIN
+        msk = tmask*shift(tmask, -1, 0, 0)
+        umaskred = msk[jpi-1, *, *]
+      ENDIF ELSE umaskred = tmask[jpi-1, *, *]
+      vmaskred = tmask[*, jpj-1, *]
+      fmaskredy = tmask[jpi-1, *, *]
+      fmaskredx = tmask[*, jpj-1, *]
+    ENDIF 
+  ENDIF ELSE BEGIN 
+    tmask = replicate(1b, jpi, jpj, jpk)
+    IF keyword_set(fullcgrid) THEN BEGIN
+      umaskred  = replicate(1b, jpj, jpk)
+      vmaskred  = replicate(1b, jpi, jpk)
+      fmaskredy = replicate(1b, jpj, jpk)
+      fmaskredx = replicate(1b, jpi, jpk)
+    ENDIF 
+  ENDELSE
+;
+  IF NOT keyword_set(fullcgrid) THEN BEGIN
+    umaskred = !values.f_nan
+    vmaskred = !values.f_nan
+    fmaskredy = !values.f_nan
+    fmaskredx = !values.f_nan
+  ENDIF
+;
+;====================================================
+; stride...
+;====================================================
+;
+  IF total(key_stride) GT 3 THEN BEGIN
+    IF key_shift NE 0 THEN BEGIN
+; for explanation, see header of read_ncdf_varget.pro
+      jpiright = key_shift
+      jpileft = jpi - key_shift - ( (key_stride[0]-1)-((key_shift-1) MOD key_stride[0]) )
+      jpi = ((jpiright-1)/key_stride[0]+1) + ((jpileft-1)/key_stride[0]+1)
+    ENDIF ELSE jpi = (jpi-1)/key_stride[0]+1
+    jpj = (jpj-1)/key_stride[1]+1
+    jpk = (jpk-1)/key_stride[2]+1
+;
+    glamt = (temporary(glamt))[0:*:stride[0], 0:*:stride[1]]
+    gphit = (temporary(gphit))[0:*:stride[0], 0:*:stride[1]]
+    e1t = (temporary(e1t))[0:*:stride[0], 0:*:stride[1]]
+    e2t = (temporary(e2t))[0:*:stride[0], 0:*:stride[1]]
+    tmask = (temporary(tmask))[0:*:stride[0], 0:*:stride[1], 0:*:stride[2]]
+    gdept = gdept[0:*:stride[2]]
+    gdepw = gdepw[0:*:stride[2]]
+    e3t = e3t[0:*:stride[2]]
+    e3w = e3w[0:*:stride[2]]
+; we must recompute glamf and gphif...
+   IF jpi GT 1 THEN BEGIN 
+     if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $
+       OR (keyword_set(key_periodic) AND key_irregular) then BEGIN
+       stepxf = (glamt + 720) MOD 360
+       stepxf = shift(stepxf, -1, -1) - stepxf
+       stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ]
+       stepxf = min(abs(stepxf), dimension = 3)
+       IF NOT keyword_set(key_periodic) THEN $
+         stepxf[jpi-1, *] = stepxf[jpi-2, *]
+     ENDIF ELSE BEGIN
+       stepxf = shift(glamt, -1, -1) - glamt
+       IF keyword_set(key_periodic) THEN $
+         stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $
+       ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *]
+     ENDELSE
+     IF jpj GT 1 THEN BEGIN 
+       stepxf[*, jpj-1] = stepxf[*, jpj-2]
+       stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2]
+     ENDIF
+     glamf = glamt + 0.5 * stepxf
+     IF jpj EQ 1 THEN glamf = reform(glamf, jpi, jpj, /over)
+   ENDIF ELSE glamf = glamt + 0.5
+   IF jpj GT 1 THEN BEGIN 
+; we must compute stepyf: y distance between T(i,j) T(i+1,j+1)
+     stepyf = shift(gphit, -1, -1) - gphit
+     stepyf[*, jpj-1] = stepyf[*, jpj-2]
+     IF jpi GT 1 THEN BEGIN
+       if NOT keyword_set(key_periodic) THEN $
+         stepyf[jpi-1, *] = stepyf[jpi-2, *]
+       stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2]
+     ENDIF 
+     gphif = gphit + 0.5 * stepyf
+   ENDIF ELSE gphif = gphit + 0.5
+;
+    IF keyword_set(fullcgrid) THEN BEGIN
+      glamu = (temporary(glamu))[0:*:stride[0], 0:*:stride[1]]
+      gphiu = (temporary(gphiu))[0:*:stride[0], 0:*:stride[1]]
+      e1u = (temporary(e1u))[0:*:stride[0], 0:*:stride[1]]
+      e2u = (temporary(e2u))[0:*:stride[0], 0:*:stride[1]]
+      glamv = (temporary(glamv))[0:*:stride[0], 0:*:stride[1]]
+      gphiv = (temporary(gphiv))[0:*:stride[0], 0:*:stride[1]]
+      e1v = (temporary(e1v))[0:*:stride[0], 0:*:stride[1]]
+      e2v = (temporary(e2v))[0:*:stride[0], 0:*:stride[1]]
+      e1f = (temporary(e1f))[0:*:stride[0], 0:*:stride[1]]
+      e2f = (temporary(e2f))[0:*:stride[0], 0:*:stride[1]]
+      umaskred = (temporary(umaskred))[0, 0:*:stride[1], 0:*:stride[2]]
+      vmaskred = (temporary(vmaskred))[0:*:stride[0], 0, 0:*:stride[2]]
+      fmaskredy = (temporary(fmaskredy))[0, 0:*:stride[1], 0:*:stride[2]]
+      fmaskredx = (temporary(fmaskredx))[0:*:stride[0], 0, 0:*:stride[2]]
+    ENDIF
+  ENDIF
+;
+;====================================================
+; apply all the grid parameters
+;====================================================
+;
+  @updateold
+  domdef
+;
+;====================================================
+; Triangulation
+;====================================================
+;
+  IF total(mask) EQ jpi*jpj*jpk $
+    AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $
+  ELSE triangles_list = triangule(/keep_cont)
+;
+;====================================================
+; time axis (default definition)
+;====================================================
+;
+  jpt = 1
+  time = 0
+;
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updateold
+  ENDIF
 ;------------------------------------------------------------
-   return
+  IF keyword_set(key_performance) EQ 1 THEN $
+    print, 'time computegrid', systime(1)-time1
+;------------------------------------------------------------
+  return
 end
 
-

trunk/Obsolete/meshlec.pro

@@ -330 +330 @@
    endif
 ;-------------------------------------------------------
+   key_yreverse = 0
+   key_zreverse = 0
+   key_partialstep = 0
+   key_stride = [1, 1, 1]
+   key_gridtype = 'c'
+;
    if not keyword_set(pasblabla) then print,'lecture '+nomfich+' finie'
    widget_control, noticebase, bad_id = toto, /destroy

trunk/Obsolete/ncdf_meshlec.pro

@@ -5 +5 @@
 ; NAME:ncdf_meshlec
 ;
-; PURPOSE:lit le fichier meshmask au format netcdf cree par OPA ds sa
-; version tout NETCDF!
+; PURPOSE:obsolete, use ncdf_meshread instead...
 ;
-; CATEGORY:lecture de grille
+; MODIFICATION HISTORY:
 ;
-; CALLING SEQUENCE:ncdf_meshlec [,' nomfich']
-; 
-; INPUTS:
-;
-;       nomfich: un string definissant le nom du fichier a lire. Si
-; ce string ne contient pas l''adresse entiere du fichier nomfich est
-; complete par iodir. Par defaut nomfich=iodir+'meshmask.nc'
-;
-; KEYWORD PARAMETERS:
-;
-;    GLAMBOUNDARY: un vecteur de 2 elements specifaint le min et le
-;    max qui doivent etre imposes en longitude (obligatoire si le
-;    tableau depasse 360 degres)
-;
-;    /CHECKDAT: cherche si il existe un ficher du meme nom que le
-;    fichier meshmask mais terminant par .dat au lieu de .dat. 
-;             Si ce fichier n''existe pas, ncdf_meshlec lit le fichier
-;         NetCdf et sauve tous les tableaux lus ds le fichier
-;         ...dat. Rq: Ce fichier est bcp plus petit que le fichier
-;         NetCdf d''origine car par ex on ne sauve pas les tableaux
-;         umask, vmask, fmask mais simplement certains de leurs bords
-;         (cf. umask.pro vmask.pro et fmask.pro) et ces masks sont
-;         sauves en bytes et toues les autres tableaux en
-;         float. Cependant ce fichier .dat est cree en fonction des
-;         parametres rentres ds le init... (ixminmesh.....) et Je ne
-;         sais rien sur la portabilite du fichier .dat.
-;            Si ce fichier .dat existe on le lit avec restore.
-;
-; OUTPUTS:mise a jours des variables de common glam, common
-; gphi,common e12,common mask,common tab (cf. common.pro)
-;
-; COMMON BLOCKS:common.pro
-;
-; SIDE EFFECTS:prend en compte (si il sont definits) les
-; ixminmesh,...et definit jpiglo,jpi,....
-;
-; RESTRICTIONS:
-;
-;  La definition de ixminmesh,ixmaxmesh,iyminmesh,iymaxmesh
-;  ,izminmesh,izmaxmesh doit etre faite avant l''entree dans cette
-;  routine. pour attribuer automatiquement ces valeurs au maximum
-;  possible les mettre toutes a -1 et ncdf_meshlec les calculera.
-;
-; EXAMPLE:
-;
-; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
-;                      12/1999
+; Aug. 2005, Sebastien Masson: switch to ncdf_meshread
 ;-
 ;------------------------------------------------------------
 ;------------------------------------------------------------
 ;------------------------------------------------------------
-PRO ncdf_meshlec, nomfich, GLAMBOUNDARY = glamboundary, GETDIMENSIONS = getdimensions, CHECKDAT = checkdat
-@common
-   tempsun = systime(1)         ; pour key_performance
-;----------------------------------------------------------
-;       constitution de l''adresse s_fichier et
-;     ouverture du fichier a l''adresse s_fichier
-;-------------------------------------------------------
-; def de nomfich par defaut
-   IF n_params() EQ 0 then nomfich = 'meshmask.nc'
-;------------------
-   if keyword_set(CHECKDAT) then begin
-      nomfichdat = strmid(nomfich, 0, strlen(nomfich)-2)+'dat'
-      thisOS = strupcase(strmid(!version.os_family, 0, 3))
-      CASE thisOS OF 
-         'MAC':sep = ':'
-         'WIN':sep = '\'
-         ELSE:sep = '/'
-      ENDCASE
-      if strpos(nomfichdat, sep) EQ -1 then begin
-; si iodir ne finit pas par sep on le rajoute
-         if rstrpos(iodir, sep) NE strlen(iodir)-1 then iodir = iodir+sep
-         nomfichdat = iodir+nomfichdat
-      endif
-      test = findfile(nomfichdat)  ; le nom cherche correspond bien a un fichier?
-      if test[0] NE '' then begin
-         noticebase=xnotice('Lecture du fichier !C '+nomfichdat+'!C ...')
-         restore, nomfichdat
-         widget_control, noticebase, bad_id = nothing, /destroy
-         if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun 
-         return
-      ENDIF
-   endif
-;------------------
-   s_fichier = isafile(file = nomfich, iodir = iodir)
+PRO ncdf_meshlec, filename, _EXTRA = ex
 ;
-   noticebase=xnotice('Lecture du fichier !C '+s_fichier+'!C ...')
+  CASE n_params() OF
+    0:ncdf_meshread, _EXTRA = ex
+    1:ncdf_meshread, filename, _EXTRA = ex
+  ENDCASE
 ;
-   cdfid=ncdf_open(s_fichier)
-   contient=ncdf_inquire(cdfid)
-
-;------------------------------------------------------------
-; differentes dimensions
-;------------------------------------------------------------
-   ncdf_diminq,cdfid,'x',name,jpiglo
-   ncdf_diminq,cdfid,'y',name,jpjglo
-   ncdf_diminq,cdfid,'z',name,jpkglo
-;
-   if keyword_set(getdimensions) then begin
-      ncdf_close,  cdfid
-      return
-   endif
-;-------------------------------------------------------
-;-------------------------------------------------------
-   if n_elements(ixminmesh) EQ 0 THEN ixminmesh = 0
-   if n_elements(ixmaxmesh) EQ 0 then ixmaxmesh = jpiglo-1
-   if ixminmesh EQ -1 THEN ixminmesh = 0
-   IF ixmaxmesh EQ -1 then ixmaxmesh = jpiglo-1
-   if n_elements(iyminmesh) EQ 0 THEN iyminmesh = 0
-   IF n_elements(iymaxmesh) EQ 0 then iymaxmesh = jpjglo-1
-   if iyminmesh EQ -1 THEN iyminmesh = 0
-   IF iymaxmesh EQ -1 then iymaxmesh = jpjglo-1
-   if n_elements(izminmesh) EQ 0 THEN izminmesh = 0
-   IF n_elements(izmaxmesh) EQ 0 then izmaxmesh = jpkglo-1
-   if izminmesh EQ -1 THEN izminmesh = 0
-   IF izmaxmesh EQ -1 then izmaxmesh = jpkglo-1
-;
-   ixmindtasauve = testvar(var = ixmindta)
-   iymindtasauve = testvar(var = iymindta)
-   izmindtasauve = testvar(var = izmindta)
-;
-   ixmindta = 0l
-   iymindta = 0l
-   izmindta = 0l
-;
-   jpi    = long(ixmaxmesh-ixminmesh+1)
-   jpj    = long(iymaxmesh-iyminmesh+1)
-   jpk    = long(izmaxmesh-izminmesh+1)
-;
-   if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1]
-   key_stride = 1l > long(key_stride)
-;
-   jpitotal = jpi
-   key_shift = long(testvar(var = key_shift))
-   key = long(key_shift MOD jpitotal)
-   if key LT 0 then key = key+jpitotal
-   ixmin = ixminmesh-ixmindta
-;
-   jpi = jpi/key_stride[0]
-   jpj = jpj/key_stride[1]
-   jpk = jpk/key_stride[2]
-;
-   jpt = 1
-   premiertps = 0
-;
-   premierx = 0
-   dernierx = jpi-1
-   premiery = 0
-   derniery = jpj-1
-   premierz = 0
-   dernierz = jpk-1
-   nx = jpi
-   ny = jpj
-   nz = 1
-   izminmeshsauve = izminmesh
-   izminmesh = 0
-;
-;-------------------------------------------------------
-; tableaux 2d:
-;-------------------------------------------------------
-   if total(key_stride) EQ 3 then $
-    namevar = ['glamt','glamu','glamv','glamf','gphit','gphiu','gphiv','gphif','e1t','e1u','e1v','e1f','e2t','e2u','e2v','e2f'] $
-   ELSE namevar = ['glamt','glamu','glamv','gphit','gphiu','gphiv']
-;
-   for i = 0, n_elements(namevar)-1 do begin
-      varcontient=ncdf_varinq(cdfid, namevar[i])
-      name = varcontient.name
-@read_ncdf_varget
-      commande = namevar[i]+'=float(res)'
-      rien = execute(commande)
-   ENDFOR
-;
-   if total(key_stride) NE 3 then BEGIN
-      glamf = glamt
-      glamf = [glamf, reform(2.*reform(glamt[jpi-1, *])-reform(glamt[jpi-2, *]), 1, jpj)]
-      glamf = [ [glamf], [reform(2*glamf[*, jpj-1]-glamf[*, jpj-2], jpi+1, 1)] ]
-      glamf = glamf+shift(glamf, -1, -1)
-      glamf = .5*glamf[0:jpi-1, 0:jpj-1]
-;
-      gphif = gphit
-      gphif = [gphif, reform(2.*reform(gphit[jpi-1, *])-reform(gphit[jpi-2, *]), 1, jpj)]
-      gphif = [ [gphif], [reform(2*gphif[*, jpj-1]-gphif[*, jpj-2], jpi+1, 1)] ]
-      gphif = gphif+shift(gphif, -1, -1)
-      gphif = .5*gphif[0:jpi-1, 0:jpj-1]
-;
-      e1t = replicate(1, jpi, jpj)
-      e2t = e1t
-      e1u = e1t
-      e2u = e1t
-      e1v = e1t
-      e2v = e1t
-      e1f = e1t
-      e2f = e1t
-   ENDIF
-;-------------------------------------------------------
-; tableaux 3d:
-;-------------------------------------------------------
-   nz = jpk
-   izminmesh = izminmeshsauve
-;
-   varcontient=ncdf_varinq(cdfid,'tmask')
-   name = varcontient.name
-@read_ncdf_varget
-   tmask = byte(res)
-;
-   varcontient=ncdf_varinq(cdfid,'umask')
-   name = varcontient.name
-   nx = 1
-   premierx = jpi-1
-@read_ncdf_varget
-   umaskred = byte(res)
-   umaskred = reform(umaskred, /over)
-;
-   varcontient=ncdf_varinq(cdfid,'vmask')
-   name = varcontient.name
-   nx = jpi
-   premierx = 0
-   ny = 1
-   premiery = jpj-1
-@read_ncdf_varget
-   vmaskred = byte(res)
-   vmaskred = reform(vmaskred, /over)
-;
-   varcontient=ncdf_varinq(cdfid,'fmask')
-   name = varcontient.name
-   nx = 1
-   premierx = jpi-1
-   ny = jpj
-   premiery = 0
-@read_ncdf_varget
-   fmaskredy = byte(res)
-   coast = where(fmaskredy NE 0 and fmaskredy NE 1)
-   IF coast[0] NE -1 THEN fmaskredy[coast] = 0b
-   fmaskredy= reform(fmaskredy, /over)
-   nx = jpi
-   premierx = 0
-   ny = 1
-   premiery = jpj-1
-@read_ncdf_varget
-   fmaskredx = byte(res)
-   coast = where(fmaskredx NE 0 and fmaskredx NE 1)
-   IF coast[0] NE -1 THEN fmaskredx[coast] = 0b
-   fmaskredx = reform(fmaskredx, /over)
-;-------------------------------------------------------
-; tableaux 1d:
-;-------------------------------------------------------
-   namevar = ['e3t','e3w','gdept','gdepw']
-   for i = 0, n_elements(namevar)-1 do begin
-      varcontient=ncdf_varinq(cdfid, namevar[i])
-      if n_elements(varcontient.dim) EQ 4 THEN BEGIN
-         commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $
-          +',offset = [0,0,izminmesh,0], count = [1,1,jpk,1]'
-         if key_stride[2] NE 1 then commande = commande+', stride=[1,1,key_stride[2],1]'
-      ENDIF ELSE BEGIN
-         commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $
-          +',offset = [izminmesh], count = [jpk]'
-         if key_stride[2] NE 1 then commande = commande+', stride=key_stride[2]'
-      ENDELSE
-      rien = execute(commande)
-      commande = namevar[i]+'=float('+namevar[i]+')'
-      rien = execute(commande)
-      commande = namevar[i]+' = reform('+namevar[i]+', /over)'
-      rien = execute(commande)
-   ENDFOR
-;-------------------------------------------------------
-   ncdf_close,  cdfid
-;-------------------------------------------------------
-; bornes de glam qui ne doivent pas depasser 360 degres... 
-;-------------------------------------------------------
-;    minglam = min(glamt, max = maxglam)
-;    if maxglam-minglam GE 360 AND NOT keyword_set(glamboundary) then $
-;     nothing = execute('glamboundary = '+xquestion('What are the longitudes boundary?', '[-180,180]',/chkwidget))
-;-------------------------------------------------------
-   if keyword_set(glamboundary) then BEGIN
-      if glamboundary[0] NE glamboundary[1] then BEGIN
-         glamt = glamt MOD 360
-         smaller = where(glamt LT glamboundary[0])
-         if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360
-         bigger = where(glamt GE glamboundary[1])
-         if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360
-         glamu = glamu MOD 360
-         smaller = where(glamu LT glamboundary[0])
-         if smaller[0] NE -1 then glamu[smaller] = glamu[smaller]+360
-         bigger = where(glamu GE glamboundary[1])
-         if bigger[0] NE -1 then glamu[bigger] = glamu[bigger]-360
-         glamv = glamv MOD 360
-         smaller = where(glamv LT glamboundary[0])
-         if smaller[0] NE -1 then glamv[smaller] = glamv[smaller]+360
-         bigger = where(glamv GE glamboundary[1])
-         if bigger[0] NE -1 then glamv[bigger] = glamv[bigger]-360
-         glamf = glamf MOD 360
-         smaller = where(glamf LT glamboundary[0])
-         if smaller[0] NE -1 then glamf[smaller] = glamf[smaller]+360
-         bigger = where(glamf GE glamboundary[1])
-         if bigger[0] NE -1 then glamf[bigger] = glamf[bigger]-360
-      endif
-   endif
-;-------------------------------------------------------
-   ixmindta = ixmindtasauve
-   iymindta = iymindtasauve
-   izmindta = izmindtasauve
-;-------------------------------------------------------
-   widget_control, noticebase, bad_id = nothing, /destroy
-;
-   if keyword_set(checkdat) then BEGIN
-   noticebase=xnotice('Ecriture du fichier !C '+nomfichdat+'!C ...')
-;
-      save, jpi, jpj, jpk, ixminmesh, ixmaxmesh, iyminmesh, iymaxmesh, izminmesh, izmaxmesh, glamt,glamu,glamv,glamf,gphit,gphiu,gphiv,gphif,e1t,e1u,e1v,e1f,e2t,e2u,e2v,e2f, tmask, umaskred, vmaskred, fmaskredx, fmaskredy, gdept, gdepw, e3t, e3w, filename = nomfichdat
-;
-      widget_control, noticebase, bad_id = nothing, /destroy
-   endif
-;
-   if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun 
-
-;-------------------------------------------------------
-   return
-end
+  return
+END

trunk/ToBeReviewed/GRILLE/changegrid.pro

@@ -28 +28 @@
    if newgrid.filetype EQ 'batch file' THEN BEGIN 
       createpro, '@'+strmid(newgrid.filename[0], 0, strlen(newgrid.filename)-4) $
-          , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $
-          +'for_createpro.pro'
+          , filename = myuniquetmpdir +'for_createpro.pro'
       return, 1
    ENDIF ELSE BEGIN 
@@ -36 +35 @@
 ;
 ;
-   key_periodique = newgrid.key_periodique
+   key_periodic = newgrid.key_periodic
+;   
+   @updateold
    domdef
-   if newgrid.triangulation EQ 1 then triangles = triangule() ELSE triangles = -1
+   if newgrid.triangulation EQ 1 then triangles_list = triangule() $
+   ELSE triangles_list = -1
+;
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+   @updateold
+  ENDIF
 ;
    return, 1

trunk/ToBeReviewed/GRILLE/cmpgrid.pro

@@ -14 +14 @@
 ; we compare the structure which caracterise the grid whith
 ; ccmeshparameters
+;
 ;
    case 1 of

trunk/ToBeReviewed/GRILLE/domdef.pro

@@ -10 +10 @@
 ; CATEGORY:
 ;
-; CALLING SEQUENCE:domdef [,lon1, lon2, lat1, lat2[,prof1,prof2]] ou
+; CALLING SEQUENCE:domdef [,lon1, lon2, lat1, lat2[,vert1,vert2]] ou
 ;             bien domdef,vecteur
 ;
 ; INPUTS:(facultatif), [vecteur a] 2, 4 ou 6 elements:;
 ; sans l''utilisation des mots cles index,xindex,yindex,zindex:
-; *prof1, prof2: pour un domaine 3D dont la partie horizontale couvre tout
+; *vert1, vert2: pour un domaine 3D dont la partie horizontale couvre tout
 ;  glam et gphi
 ; *lon1, lon2, lat1, lat2:
 ; definissant les longitudes min. max et les latitudes min, max du  domaine a
 ; etudier (tous les niveaux sont selectiones)
-; *lon1,lon2,lat1,lat2,prof1,prof2 pour specifier les profondeurs.
+; *lon1,lon2,lat1,lat2,vert1,vert2 pour specifier les profondeurs.
 ;
 ; KEYWORD PARAMETERS:
+;
+;       ENDPOINTS: a four elements vector [x1,y1,x2,y2] used to specify
+;       that domdef must define the box used to make a plot (pltz, pltt,
+;       plt1d) done strictly along the line (that can have any direction)
+;       starting at (x1, y1) ending at (x2, y2). When defining endpoints,
+;       you must also define TYPE which define the type of plots
+;       ('pltz', 'xt', 'yt', 'zt', 'x', 'y', 'z', 't') will used
+;       ENDPOINTS keywords
 ;
 ;       FINDALWAYS:oblige a redefinir une boite meme qd auqun point
@@ -28 +36 @@
 ;       grille.
 ;
-;       GRILLE:un string ou un vecteur de string contennant le nom des
+;       GRIDTYPE:un string ou un vecteur de string contennant le nom des
 ;       grilles (determinees uniquement par : 'T','U','V','W','F') pour
 ;       lesquelles le calcul doit etre fait. par
@@ -61 +69 @@
 ;   -nzt,w:entier qui contient le nombre de pts en profondeur de
 ;    la grille reduite au domaine 3D
-;   -premierxt,u,v,f: le premier indice qui delimite le sous domaine
+;   -firstxt,u,v,f: le first indice qui delimite le sous domaine
 ;   suivant x
-;   -premieryt,u,v,f: le premier indice qui delimite le sous domaine
+;   -firstyt,u,v,f: le first indice qui delimite le sous domaine
 ;   suivant y
-;   -premierzt,w: le premier indice qui delimite le sous domaine
+;   -firstzt,w: le first indice qui delimite le sous domaine
 ;   suivant z
-;   -dernierxt,u,v,f: le dernier indice qui delimite le sous domaine
+;   -lastxt,u,v,f: le last indice qui delimite le sous domaine
 ;   suivant x
-;   -dernieryt,u,v,f: le dernier indice qui delimite le sous domaine
+;   -lastyt,u,v,f: le last indice qui delimite le sous domaine
 ;   suivant y
-;   -dernierzt,w: le dernier indice qui delimite le sous domaine
+;   -lastzt,w: le last indice qui delimite le sous domaine
 ;   suivant z
 ;
@@ -85 +93 @@
 ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
 ;                       8/2/98
+; rewrite everything, debug and spee-up Sebastien Masson April 2005
 ;-
 ;------------------------------------------------------------
 ;------------------------------------------------------------
 ;------------------------------------------------------------
-pro domdef,x1,x2,y1,y2,z1,z2, FINDALWAYS = findalways, GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
-@common
-   tempsun = systime(1)         ; pour key_performance
-;-------------------------------------------------------------------
-; determination en fonction du nombre d''arguments de
-; lon1,lon2,lat1,lat2,prof1,prof2
-;-------------------------------------------------------------------
-   if not keyword_set(grille) then grille=['T', 'U', 'V', 'W', 'F'] ELSE $ 
-    for i = 0, n_elements(grille)-1 do grille[i] = strupcase(grille[i])
-   if keyword_set(memeindices) then grille=['T', grille]
-   CASE N_PARAMS() OF
-      0: BEGIN
+pro domdef, x1, x2, y1, y2, z1, z2, FINDALWAYS = findalways $
+            , GRIDTYPE = gridtype, MEMEINDICES = memeindices $
+            , XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex $
+            , ENDPOINTS = endpoints, TYPE = type $
+            , INDEX = index, _extra = ex
+;------------------------------------------------------------
+; include commons
+@cm_4mesh
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+@updatekwd
+  ENDIF
+;---------------------
+  tempsun = systime(1)          ; pour key_performance
+;
+  CASE N_PARAMS() OF
+    0:
+    1:
+    2:
+    4:
+    6:
+    ELSE:BEGIN
+      ras = report('Bad number of parameter in the call of domdef.')
+      RETURN
+    END
+  ENDCASE
+;
+  IF keyword_set(endpoints) THEN BEGIN
+    IF NOT keyword_set(type) THEN BEGIN 
+      dummy = report(['If domdef is used do find the box associated' $
+                      , 'to endpoints, you must also specify type keyword'])
+      return
+    ENDIF 
+    CASE N_PARAMS() OF
+      0:
+      1:boxzoom = [x1]
+      2:boxzoom = [x1, x2]
+      4:boxzoom = [x1, x2, y1, y2]
+      6:boxzoom = [x1, x2, y1, y2, z1, z2]
+    ENDCASE
+    section, BOXZOOM = boxzoom, ENDPOINTS = endpoints, TYPE = type, /ONLYBOX
+    return
+  ENDIF   
+
+;-------------------------------------------------------------------
+; recall domdef when there is only one input parameter
+;-------------------------------------------------------------------
+;
+  IF N_PARAMS() EQ 1 THEN BEGIN
+    CASE n_elements(x1) OF
+      2:domdef, x1[0], x1[1], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
+      4:domdef, x1[0], x1[1], x1[2], x1[3], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
+      6:domdef, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
+      ELSE:BEGIN
+        ras = report('Bad number of elements in x1')
+        RETURN
+      END
+    ENDCASE
+    RETURN  
+  ENDIF
+;-------------------------------------------------------------------
+; default definitions and checks
+;-------------------------------------------------------------------
+  IF NOT keyword_set(gridtype) THEN gridtype = ['T', 'U', 'V', 'W', 'F'] $
+  ELSE gridtype = strupcase(gridtype)
+  IF keyword_set(memeindices) THEN gridtype = ['T', gridtype]
+; default definitions
+  lon1t = 99999. & lon2t = -99999. & lat1t = 99999. & lat2t = -99999.
+  lon1u = 99999. & lon2u = -99999. & lat1u = 99999. & lat2u = -99999.
+  lon1v = 99999. & lon2v = -99999. & lat1v = 99999. & lat2v = -99999.
+  lon1f = 99999. & lon2f = -99999. & lat1f = 99999. & lat2f = -99999.
+  vert1t = 99999. & vert2t = -99999. & vert1w = 99999. & vert2w = -99999.
+;
+  IF N_PARAMS() EQ 2 THEN GOTO, vertical
+;
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+; define all horizontal parameters ...
 ; lon1 et lon2
-         test = 1
-         IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, glamt[*]]
-         IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $
-          test = [test, glamt[*]]
-         IF (where(grille eq 'U'))[0] NE -1 AND n_elements(glamu) NE 0 THEN test = [test, glamu[*]]
-         IF (where(grille eq 'V'))[0] NE -1 AND n_elements(glamv) NE 0 THEN test = [test, glamv[*]]
-         IF (where(grille eq 'F'))[0] NE -1 AND n_elements(glamf) NE 0 THEN test = [test, glamf[*]]
-         test = test[1:n_elements(test)-1] 
-         lon1=min([test], max = max)
-         lon2=max
 ; lat1 et lat2
-         test = 1
-         IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, gphit[*]]
-         IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $
-          test = [test, gphit[*]]
-         IF (where(grille eq 'U'))[0] NE -1 AND n_elements(gphiu) NE 0 THEN test = [test, gphiu[*]]
-         IF (where(grille eq 'V'))[0] NE -1 AND n_elements(gphiv) NE 0 THEN test = [test, gphiv[*]]
-         IF (where(grille eq 'F'))[0] NE -1 AND n_elements(gphif) NE 0 THEN test = [test, gphif[*]]
-         test = test[1:n_elements(test)-1] 
-         lat1=min([test], max = max)
-         lat2=max
-; prof1 et prof2
-         prof1=0
-         test = 1
-         IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN test = [test, gdept[jpk-1]]
-         IF (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN test = [test, gdepw[jpk-1]]
-         test = test[1:n_elements(test)-1] 
-         prof2=max([test])
-      end
-      1: BEGIN
-         xx1 = reform(x1)
-         taille = size(xx1)
-         CASE taille[1] OF
-            2:domdef, xx1[0], xx1[1]                                , FINDALWAYS = findalways,GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
-            4:domdef, xx1[0], xx1[1], xx1[2], xx1[3]                , FINDALWAYS = findalways,GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
-            6:domdef, xx1[0], xx1[1], xx1[2], xx1[3], xx1[4], xx1[5], FINDALWAYS = findalways,GRILLE=grille, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
-            else: begin
-               ras = report('Mauvais nombre de parametre dans l''appel de DOMDEF')
-               return
-            end
-         ENDCASE
-         return
-      end
-      2: begin
-; lon1 et lon2
-         IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN
-            ouca = where(gdept ge prof1 and gdept le prof2,nzt)
-            if ouca[0] NE -1 then zt=gdept[ouca] ELSE zt = -1
-         ENDIF
-         if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then BEGIN
-            ouca = where(gdepw ge prof1 and gdepw le prof2,nzw)
-            if ouca[0] NE -1 then zw=gdepw[ouca] ELSE zw = -1
-         ENDIF
-         test = 1
-         IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, glamt[*]]
-         IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $
-          test = [test, glamt[*]]
-         IF (where(grille eq 'U'))[0] NE -1 AND n_elements(glamu) NE 0 THEN test = [test, glamu[*]]
-         IF (where(grille eq 'V'))[0] NE -1 AND n_elements(glamv) NE 0 THEN test = [test, glamv[*]]
-         IF (where(grille eq 'F'))[0] NE -1 AND n_elements(glamf) NE 0 THEN test = [test, glamf[*]]
-         test = test[1:n_elements(test)-1] 
-         lon1=min([test], max = max)
-         lon2=max
-; lat1 et lat2
-         test = 1
-         IF (where(grille eq 'T'))[0] NE -1 THEN test = [test, gphit[*]]
-         IF (where(grille eq 'W'))[0] NE -1 AND (where(grille eq 'T'))[0] EQ -1 THEN $
-          test = [test, gphit[*]]
-         IF (where(grille eq 'U'))[0] NE -1 AND n_elements(gphiu) NE 0 THEN test = [test, gphiu[*]]
-         IF (where(grille eq 'V'))[0] NE -1 AND n_elements(gphiv) NE 0 THEN test = [test, gphiv[*]]
-         IF (where(grille eq 'F'))[0] NE -1 AND n_elements(gphif) NE 0 THEN test = [test, gphif[*]]
-         test = test[1:n_elements(test)-1] 
-         lat1=min([test], max = max)
-         lat2=max
-; prof1 et prof2
-         if keyword_set(zindex) OR keyword_set(index) then BEGIN
-            z1 = x1
-            z2 = x2
-            IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN prof1 = gdept[z1] ELSE prof1=gdepw[z1]
-            if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then prof2=gdepw[z2] ELSE prof2=gdept[z2]
-         ENDIF ELSE BEGIN
-            prof1=x1
-            prof2=x2
-         ENDELSE 
-      end
-      4: begin 
-; lon1 et lon2
-         if keyword_set(xindex) OR keyword_set(index) then BEGIN
-            if keyword_set(yindex) OR keyword_set(index) then BEGIN
-               if n_elements(glamf) NE 0 then BEGIN
-               lon1 = min([glamt[x1:x2, y1:y2], glamf[x1:x2, y1:y2]], max = lon2)
-               ENDIF ELSE BEGIN 
-               lon1 = min(glamt[x1:x2, y1:y2], max = lon2)
-               ENDELSE
-            ENDIF ELSE BEGIN 
-               lon1 = glamt[x1, 0]
-               if n_elements(glamf) NE 0 then lon2 = glamf[x2, 0] $
-               ELSE lon2 = glamt[x2, 0]
+; firstx[tuvf], lastx[tuvf], nx[tuvf]
+;-------------------------------------------------------------------
+; check if the grid is defined for U and V points. If not, take care
+; of the cases gridtype eq 'U' or 'V'
+;-------------------------------------------------------------------
+  errstatus = 0
+  IF (finite(glamu[0]*gphiu[0]) EQ 0 OR n_elements(glamu) EQ 0 OR n_elements(gphiu) EQ 0) AND (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
+    firstxu = !values.f_nan
+    lastxu = !values.f_nan
+    nxu = !values.f_nan
+    okgrid = where(gridtype NE 'U', count)
+    IF count NE 0 THEN gridtype = gridtype[okgrid] $
+    ELSE errstatus = report('U grid is undefined. Impossible to call domdef with vargid = ''U''')
+  ENDIF
+;
+  IF (finite(glamv[0]*gphiv[0]) EQ 0 OR n_elements(glamv) EQ 0 OR n_elements(gphiv) EQ 0) AND (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
+    firstxv = !values.f_nan
+    lastxv = !values.f_nan
+    nxv = !values.f_nan
+    okgrid = where(gridtype NE 'V', count)
+    IF count NE 0 THEN gridtype = gridtype[okgrid] $
+    ELSE errstatus = report('V grid is undefined. Impossible to call domdef with vargid = ''V''')
+  ENDIF
+  IF errstatus EQ -1 THEN return
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+; horizontal domain defined with lon1, lon2, lat1 and lat2
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+  IF N_PARAMS() EQ 0 $
+    OR ( (N_PARAMS() EQ 4 OR N_PARAMS() EQ 6) $
+         AND NOT keyword_set(xindex) AND NOT keyword_set(yindex) AND NOT keyword_set(index) ) THEN BEGIN
+    IF N_PARAMS() EQ 0 THEN BEGIN
+; find lon1 and lon2 the longitudinal boudaries of the full domain
+      IF (where(gridtype eq 'T'))[0] NE -1 THEN lon1t = min(glamt, max = lon2t)
+      IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lon1t = min(glamt, max = lon2t)
+      IF (where(gridtype eq 'U'))[0] NE -1 THEN lon1u = min(glamu, max = lon2u)
+      IF (where(gridtype eq 'V'))[0] NE -1 THEN lon1v = min(glamv, max = lon2v)
+      IF (where(gridtype eq 'F'))[0] NE -1 THEN lon1f = min(glamf, max = lon2f)
+      lon1 = min([lon1t, lon1u, lon1v, lon1f])
+      lon2 = max([lon2t, lon2u, lon2v, lon2f])
+; find lat1 and lat2 the latitudinal boudaries of the full domain
+      IF (where(gridtype eq 'T'))[0] NE -1 THEN lat1t = min(gphit, max = lat2t)
+      IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lat1t = min(gphit, max = lat2t)
+      IF (where(gridtype eq 'U'))[0] NE -1 THEN lat1u = min(gphiu, max = lat2u)
+      IF (where(gridtype eq 'V'))[0] NE -1 THEN lat1v = min(gphiv, max = lat2v)
+      IF (where(gridtype eq 'F'))[0] NE -1 THEN lat1f = min(gphif, max = lat2f)
+      lat1 = min([lat1t, lat1u, lat1v, lat1f])
+      lat2 = max([lat2t, lat2u, lat2v, lat2f])
+    ENDIF ELSE BEGIN 
+      lon1 = min([x1, x2], max = lon2)
+      lat1 = min([y1, y2], max = lat2)
+    ENDELSE
+; find firstxt, firstxt, nxt and nyt according to lon1, lon2, lat1 and lat2
+    IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
+      dom = where( (glamt GE lon1) AND (glamt LE lon2) $
+                   AND (gphit GE lat1) AND (gphit LE lat2) )
+      IF (dom[0] EQ -1) THEN BEGIN
+        IF keyword_set(findalways) THEN BEGIN
+          domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+          IF N_PARAMS() EQ 6 THEN $
+            domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
+          RETURN
+        ENDIF ELSE BEGIN
+          ras = report('WARNING, The box does not contain any T points...')
+          firstxt = -1 & lastxt = -1 & nxt = 0
+          firstyt = -1 & lastyt = -1 & nyt = 0
+        ENDELSE
+      ENDIF ELSE BEGIN
+        jyt = dom / jpi
+        ixt = temporary(dom) MOD jpi
+        firstxt = min(temporary(ixt), max = lastxt)
+        firstyt = min(temporary(jyt), max = lastyt)
+        nxt = lastxt - firstxt + 1
+        nyt = lastyt - firstyt + 1
+      ENDELSE
+    ENDIF
+; find firstxu, firstxu, firstyu, firstyu, nxu and nyu
+; according to lon1, lon2, lat1 and lat2
+    IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
+      IF keyword_set(memeindices) THEN BEGIN
+        firstxu = firstxt & lastxu = lastxt & nxu = nxt
+        firstyu = firstyt & lastyu = lastyt & nyu = nyt
+      ENDIF ELSE BEGIN
+        dom = where( (glamu GE lon1) AND (glamu LE lon2) $
+                     AND (gphiu GE lat1) AND (gphiu LE lat2) )
+        IF (dom[0] EQ -1) THEN BEGIN
+          IF keyword_set(findalways) THEN BEGIN
+            domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+            IF N_PARAMS() EQ 6 THEN $
+              domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
+            RETURN
+          ENDIF ELSE BEGIN
+            ras = report('WARNING, The box does not contain any U points...')
+            firstxu = -1 & lastxu = -1 & nxu = 0
+            firstyu = -1 & lastyu = -1 & nyu = 0
+          ENDELSE
+        ENDIF ELSE BEGIN
+          jyu = dom / jpi
+          ixu = temporary(dom) MOD jpi
+          firstxu = min(temporary(ixu), max = lastxu)
+          firstyu = min(temporary(jyu), max = lastyu)
+          nxu = lastxu - firstxu + 1
+          nyu = lastyu - firstyu + 1
+        ENDELSE
+      ENDELSE
+    ENDIF
+; find firstxv, firstxv, firstyv, firstyv, nxv and nyv
+;  according to lon1, lon2, lat1 and lat2
+    IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
+      IF keyword_set(memeindices) THEN BEGIN
+        firstxv = firstxt & lastxv = lastxt & nxv = nxt
+        firstyv = firstyt & lastyv = lastyt & nyv = nyt
+      ENDIF ELSE BEGIN
+        dom = where( (glamv GE lon1) AND (glamv LE lon2) $
+                     AND (gphiv GE lat1) AND (gphiv LE lat2) )
+        IF (dom[0] EQ -1) THEN BEGIN
+          IF keyword_set(findalways) THEN BEGIN
+            domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+            IF N_PARAMS() EQ 6 THEN $
+              domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
+            RETURN
+          ENDIF ELSE BEGIN
+            ras = report('WARNING, The box does not contain any V points...')
+            firstxv = -1 & lastxv = -1 & nxv = 0
+            firstyv = -1 & lastyv = -1 & nyv = 0
+          ENDELSE
+        ENDIF ELSE BEGIN
+          jyv = dom / jpi
+          ixv = temporary(dom) MOD jpi
+          firstxv = min(temporary(ixv), max = lastxv)
+          firstyv = min(temporary(jyv), max = lastyv)
+          nxv = lastxv - firstxv + 1
+          nyv = lastyv - firstyv + 1
+        ENDELSE
+      ENDELSE
+    ENDIF
+; find firstxf, firstxf, firstyf, firstyf, nxf and nyf
+; according to lon1, lon2, lat1 and lat2
+    IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN
+      IF keyword_set(memeindices) THEN BEGIN
+        firstxf = firstxt & lastxf = lastxt & nxf = nxt
+        firstyf = firstyt & lastyf = lastyt & nyf = nyt
+      ENDIF ELSE BEGIN
+        dom = where( (glamf GE lon1) AND (glamf LE lon2) $
+                     AND (gphif GE lat1) AND (gphif LE lat2) )
+        IF (dom[0] EQ -1) THEN BEGIN
+          IF keyword_set(findalways) THEN BEGIN
+            domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+            IF N_PARAMS() EQ 6 THEN $
+              domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
+            RETURN
+          ENDIF ELSE BEGIN
+            ras = report('WARNING, The box does not contain any F points...')
+            firstxf = -1 & lastxf = -1 & nxf = 0
+            firstyf = -1 & lastyf = -1 & nyf = 0
+          ENDELSE
+        ENDIF ELSE BEGIN
+          jyf = dom / jpi
+          ixf = temporary(dom) MOD jpi
+          firstxf = min(temporary(ixf), max = lastxf)
+          firstyf = min(temporary(jyf), max = lastyf)
+          nxf = lastxf - firstxf + 1
+          nyf = lastyf - firstyf + 1
+        ENDELSE
+      ENDELSE
+    ENDIF
+;
+  ENDIF ELSE BEGIN
+    CASE 1 OF
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+; horizontal domain defined with the X and Y indexes
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+      (keyword_set(xindex) AND keyword_set(yindex)) OR keyword_set(index):BEGIN
+        fstx = min([x1, x2], max = lstx)
+        fsty = min([y1, y2], max = lsty)
+        IF fstx LT 0 OR lstx GE jpi THEN BEGIN
+          ras = report('Bad definition of X1 or X2')
+          return
+        ENDIF
+        IF fsty LT 0 OR lsty GE jpj THEN BEGIN
+          ras = report('Bad definition of Y1 or Y2')
+          return
+        ENDIF
+        nx = lstx - fstx + 1
+        ny = lsty - fsty + 1
+; find lon1t, lon2t, lat1t, lat2t, firstxt, firstxt, nxt and nyt
+; according to x1, x2, y1, y2
+        IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1 THEN BEGIN
+          lon1t = min(glamt[fstx:lstx, fsty:lsty], max = lon2t)
+          lat1t = min(gphit[fstx:lstx, fsty:lsty], max = lat2t)
+          firstxt = fstx & lastxt = lstx
+          firstyt = fsty & lastyt = lsty
+          nxt = nx & nyt = ny
+        ENDIF
+; find lon1u, lon2u, lat1u, lat2u, firstxu, firstxu, nxu and nyu
+; according to x1, x2, y1, y2
+        IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
+          lon1u = min(glamu[fstx:lstx, fsty:lsty], max = lon2u)
+          lat1u = min(gphiu[fstx:lstx, fsty:lsty], max = lat2u)
+          firstxu = fstx & lastxu = lstx
+          firstyu = fsty & lastyu = lsty
+          nxu = nx & nyu = ny
+        ENDIF
+; find lon1v, lon2v, lat1v, lat2v, firstxv, firstxv, nxv and nyv
+; according to x1, x2, y1, y2
+        IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
+          lon1v = min(glamv[fstx:lstx, fsty:lsty], max = lon2v)
+          lat1v = min(gphiv[fstx:lstx, fsty:lsty], max = lat2v)
+          firstxv = fstx & lastxv = lstx
+          firstyv = fsty & lastyv = lsty
+          nxv = nx & nyv = ny
+        ENDIF         
+; find lon1f, lon2f, lat1f, lat2f, firstxf, firstxf, nxf and nyf
+; according to x1, x2, y1, y2
+        IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN
+          lon1f = min(glamf[fstx:lstx, fsty:lsty], max = lon2f)
+          lat1f = min(gphif[fstx:lstx, fsty:lsty], max = lat2f)
+          firstxf = fstx & lastxf = lstx
+          firstyf = fsty & lastyf = lsty
+          nxf = nx & nyf = ny
+        ENDIF
+        lon1 = min([lon1t, lon1u, lon1v, lon1f])
+        lon2 = max([lon2t, lon2u, lon2v, lon2f])
+        lat1 = min([lat1t, lat1u, lat1v, lat1f])
+        lat2 = max([lat2t, lat2u, lat2v, lat2f])
+      END
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+; horizontal domain defined with the X index and lat1/lat2
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+      keyword_set(xindex):BEGIN
+        fstx = min([x1, x2], max = lstx)
+        IF fstx LT 0 OR lstx GE jpi THEN BEGIN
+          ras = report('Bad definition of X1 or X2')
+          return
+        ENDIF
+        nx = lstx - fstx + 1
+        lat1 = min([y1, y2], max = lat2)
+; find lon1t, lon2t, firstxt, firstxt, firstyt, firstyt, nxt
+; and nyt according to x1, x2, lat1 and lat2
+        IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
+          firstxt = fstx & lastxt = lstx & nxt = nx
+          dom = where( (gphit[fstx:lstx, *] GE lat1) AND (gphit[fstx:lstx, *] LE lat2) )
+          IF (dom[0] EQ -1) THEN BEGIN
+            IF keyword_set(findalways) THEN BEGIN
+              domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+              IF N_PARAMS() EQ 6 THEN $
+                domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
+              RETURN
+            ENDIF ELSE BEGIN
+              ras = report('WARNING, The box does not contain any T points...')
+              firstyt = -1 & lastyt = -1 & nyt = 0
             ENDELSE
-         ENDIF ELSE BEGIN
-            lon1=x1
-            lon2=x2
-         ENDELSE 
-; lat1 et lat2
-         if keyword_set(yindex) OR keyword_set(index) then BEGIN
-            if keyword_set(xindex) OR keyword_set(index) then begin
-               lat1 = min([gphit[x1, y1], gphit[x2, y1]])
-               if n_elements(gphif) NE 0 then lat2 = max([gphif[x1, y2], gphif[x2, y2]]) $
-               ELSE lat2 = max([gphit[x1, y2], gphit[x2, y2]])
-            ENDIF ELSE BEGIN 
-               lat1 = gphit[0, y1]
-               if n_elements(gphif) NE 0 then lat2 = gphif[0, y2] $
-               ELSE lat2 = gphit[0, y2]
+          ENDIF ELSE BEGIN
+            jyt = temporary(dom) / nx
+            firstyt = min(temporary(jyt), max = lastyt)
+            nyt = lastyt - firstyt + 1
+          ENDELSE
+          IF nyt NE 0 THEN lon1t = min(glamt[firstxt:lastxt, firstyt:lastyt], max = lon2t)
+        ENDIF
+; find lon1u, lon2u, firstxu, firstxu, firstyu, firstyu, nxu
+; and nyu according to x1, x2, lat1 and lat2
+        IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
+          firstxu = fstx & lastxu = lstx & nxu = nx
+          IF keyword_set(memeindices) THEN BEGIN
+            firstyu = firstyt & lastyu = lastyt & nyu = nyt
+          ENDIF ELSE BEGIN
+            dom = where( (gphiu[fstx:lstx, *] GE lat1) AND (gphiu[fstx:lstx, *] LE lat2) )
+            IF (dom[0] EQ -1) THEN BEGIN
+              IF keyword_set(findalways) THEN BEGIN
+                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+                IF N_PARAMS() EQ 6 THEN $
+                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
+                RETURN
+              ENDIF ELSE BEGIN
+                ras = report('WARNING, The box does not contain any U points...')
+                firstyu = -1 & lastyu = -1 & nyu = 0
+              ENDELSE
+            ENDIF ELSE BEGIN
+              jyu = temporary(dom) / nx
+              firstyu = min(temporary(jyu), max = lastyu)
+              nyu = lastyu - firstyu + 1
             ENDELSE
-         ENDIF ELSE BEGIN
-            lat1=y1
-            lat2=y2
-         ENDELSE 
-; prof1 et prof2
-         prof1=0
-         test = 1
-         IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN test = [test, gdept[jpk-1]]
-         IF (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN test = [test, gdepw[jpk-1]]
-         test = test[1:n_elements(test)-1] 
-         prof2=max([test])
-      end
-      6: begin 
-; lon1 et lon2
-         if keyword_set(xindex) OR keyword_set(index) then BEGIN
-            if keyword_set(yindex) OR keyword_set(index) then begin
-               lon1 = min([glamt[x1, y1], glamt[x1, y2]])
-               if n_elements(glamf) NE 0 then lon2 = max([glamf[x2, y1], glamf[x2, y2]]) $
-               ELSE lon2 = max([glamt[x2, y1], glamt[x2, y2]])
-            ENDIF ELSE BEGIN 
-               lon1 = glamt[x1, 0]
-               if n_elements(glamf) NE 0 then lon2 = glamf[x2, 0] $
-               ELSE lon2 = glamt[x2, 0]
+          ENDELSE
+          IF nyu NE 0 THEN lon1u = min(glamu[firstxu:lastxu, firstyu:lastyu], max = lon2u)
+        ENDIF
+; find lon1v, lon2v, firstxv, firstxv, firstyv, firstyv, nxv
+; and nyv according to x1, x2, lat1 and lat2
+        IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
+          firstxv = fstx & lastxv = lstx & nxv = nx
+          IF keyword_set(memeindices) THEN BEGIN
+            firstyv = firstyt & lastyv = lastyt & nyv = nyt
+          ENDIF ELSE BEGIN
+            dom = where( (gphiv[fstx:lstx, *] GE lat1) AND (gphiv[fstx:lstx, *] LE lat2) )
+            IF (dom[0] EQ -1) THEN BEGIN
+              IF keyword_set(findalways) THEN BEGIN
+                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+                IF N_PARAMS() EQ 6 THEN $
+                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
+                RETURN
+              ENDIF ELSE BEGIN
+                ras = report('WARNING, The box does not contain any V points...')
+                firstyv = -1 & lastyv = -1 & nyv = 0
+              ENDELSE
+            ENDIF ELSE BEGIN
+              jyv = temporary(dom) / nx
+              firstyv = min(temporary(jyv), max = lastyv)
+              nyv = lastyv - firstyv + 1
             ENDELSE
-         ENDIF ELSE BEGIN
-            lon1=x1
-            lon2=x2
-         ENDELSE 
-; lat1 et lat2
-         if keyword_set(yindex) OR keyword_set(index) then BEGIN
-            if keyword_set(xindex) OR keyword_set(index) then begin
-               lat1 = min([gphit[x1, y1], gphit[x2, y1]])
-               if n_elements(gphif) NE 0 then lat2 = max([gphif[x1, y2], gphif[x2, y2]]) $
-               ELSE lat2 = max([gphit[x1, y2], gphit[x2, y2]])
-            ENDIF ELSE BEGIN 
-               lat1 = gphit[0, y1]
-               if n_elements(gphif) NE 0 then lat2 = gphif[0, y2] $
-               ELSE lat2 = gphit[0, y2]
+          ENDELSE
+          IF nyv NE 0 THEN lon1v = min(glamv[firstxv:lastxv, firstyv:lastyv], max = lon2v)
+        ENDIF
+; find lon1f, lon2f, firstxf, firstxf, firstyf, firstyf, nxf
+; and nyf according to x1, x2, lat1 and lat2
+        IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN
+          firstxf = fstx & lastxf = lstx & nxf = nx
+          IF keyword_set(memeindices) THEN BEGIN
+            firstyf = firstyt & lastyf = lastyt & nyf = nyt
+          ENDIF ELSE BEGIN
+            dom = where( (gphif[fstx:lstx, *] GE lat1) AND (gphif[fstx:lstx, *] LE lat2) )
+            IF (dom[0] EQ -1) THEN BEGIN
+              IF keyword_set(findalways) THEN BEGIN
+                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+                IF N_PARAMS() EQ 6 THEN $
+                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
+                RETURN
+              ENDIF ELSE BEGIN
+                ras = report('WARNING, The box does not contain any F points...')
+                firstyf = -1 & lastyf = -1 & nyf = 0
+              ENDELSE
+            ENDIF ELSE BEGIN
+              jyf = temporary(dom) / nx
+              firstyf = min(temporary(jyf), max = lastyf)
+              nyf = lastyf - firstyf + 1
             ENDELSE
-         ENDIF ELSE BEGIN
-            lat1=y1
-            lat2=y2
-         ENDELSE 
-; prof1 et prof2
-         if keyword_set(zindex) OR keyword_set(index) then begin
-            IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN prof1=gdept[z1] ELSE prof1=gdepw[z1]
-            if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then prof2=gdepw[z2] ELSE prof2=gdept[z2]
-         ENDIF ELSE BEGIN
-            prof1=z1
-            prof2=z2
-         ENDELSE 
-      end
-      else: begin
-         ras = report('Mauvais nombre de parametre dans l''appel de DOMDEF')
-         return
-      end
-   endcase  
-;----------------------------------------------------------------------       
-; determination de nzt et nwt
-;----------------------------------------------------------------------       
-   IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN
-      if keyword_set(zindex) OR keyword_set(index) then nzt = z2-z1+1 $
-      ELSE BEGIN 
-         ouca = where(gdept ge prof1 and gdept le prof2,nzt)
-         if ouca[0] NE -1 then zt=gdept[ouca] ELSE zt = -1
-      ENDELSE
-   ENDIF
-   if (where(grille eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 then BEGIN
-      if keyword_set(zindex) OR keyword_set(index) then nzw = z2-z1+1 $
-      ELSE BEGIN 
-         ouca = where(gdepw ge prof1 and gdepw le prof2,nzw)
-         if ouca[0] NE -1 then zw=gdepw[ouca] ELSE zw = -1
-      ENDELSE
-   ENDIF
-;----------------------------------------------------------------------       
-; determination de premierzt, dernierzt, premierzw, dernierzw
-;----------------------------------------------------------------------       
-   IF (inter(byte(grille), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN
-      if keyword_set(zindex) OR keyword_set(index) then BEGIN 
-         premierzt = z1
-         dernierzt = z2
+          ENDELSE
+          IF nyf NE 0 THEN lon1f = min(glamf[firstxf:lastxf, firstyf:lastyf], max = lon2f)
+        ENDIF
+        lon1 = min([lon1t, lon1u, lon1v, lon1f])
+        lon2 = max([lon2t, lon2u, lon2v, lon2f])
+      END
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+; horizontal domain defined with the Y index and lon1/lon2
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+      keyword_set(yindex):BEGIN
+        fsty = min([y1, y2], max = lsty)
+        IF fsty LT 0 OR lsty GE jpj THEN BEGIN
+          ras = report('Bad definition of Y1 or Y2')
+          return
+        ENDIF
+        ny = lsty - fsty + 1
+        lon1 = min([x1, x2], max = lon2)
+; find lat1t, lat2t, firstxt, firstxt, firstyt, firstyt, nxt
+; and nyt according to x1, x2, lon1 and lon2
+        IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
+          firstyt = fsty & lastyt = lsty & nyt = ny
+          dom = where( (glamt[*, fsty:lsty] GE lon1) AND (glamt[*, fsty:lsty] LE lon2) )
+          IF (dom[0] EQ -1) THEN BEGIN
+            IF keyword_set(findalways) THEN BEGIN
+              domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+              IF N_PARAMS() EQ 6 THEN $
+                domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
+              RETURN
+            ENDIF ELSE BEGIN
+              ras = report('WARNING, The box does not contain any T points...')
+              firstxt = -1 & lastxt = -1 & nxt = 0
+            ENDELSE
+          ENDIF ELSE BEGIN
+            jxt = temporary(dom) MOD jpi
+            firstxt = min(temporary(jxt), max = lastxt)
+            nxt = lastxt - firstxt + 1
+          ENDELSE
+          IF nxt NE 0 THEN lat1t = min(gphit[firstxt:lastxt, firstyt:lastyt], max = lat2t)
+        ENDIF
+; find lat1u, lat2u, firstxu, firstxu, firstyu, firstyu, nxu
+; and nyu according to x1, x2, lon1 and lon2
+        IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
+          firstyu = fsty & lastyu = lsty & nyu = ny
+          IF keyword_set(memeindices) THEN BEGIN
+            firstxu = firstyt & lastxu = lastyt & nxu = nxt
+          ENDIF ELSE BEGIN
+            dom = where( (glamu[*, fsty:lsty] GE lon1) AND (glamu[*, fsty:lsty] LE lon2) )
+            IF (dom[0] EQ -1) THEN BEGIN
+              IF keyword_set(findalways) THEN BEGIN
+                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+                IF N_PARAMS() EQ 6 THEN $
+                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
+                RETURN
+              ENDIF ELSE BEGIN
+                ras = report('WARNING, The box does not contain any U points...')
+                firstxu = -1 & lastxu = -1 & nxu = 0
+              ENDELSE
+            ENDIF ELSE BEGIN
+              jxu = temporary(dom) MOD jpi
+              firstxu = min(temporary(jxu), max = lastxu)
+              nxu = lastxu - firstxu + 1
+            ENDELSE
+          ENDELSE
+          IF nxu NE 0 THEN lat1u = min(gphiu[firstxu:lastxu, firstyu:lastyu], max = lat2u)
+        ENDIF
+; find lat1v, lat2v, firstxv, firstxv, firstyv, firstyv, nxv
+; and nyv according to x1, x2, lon1 and lon2
+        IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
+          firstyv = fsty & lastyv = lsty & nyv = ny
+          IF keyword_set(memeindices) THEN BEGIN
+            firstxv = firstyt & lastxv = lastyt & nxv = nxt
+          ENDIF ELSE BEGIN
+            dom = where( (glamv[*, fsty:lsty] GE lon1) AND (glamv[*, fsty:lsty] LE lon2) )
+            IF (dom[0] EQ -1) THEN BEGIN
+              IF keyword_set(findalways) THEN BEGIN
+                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+                IF N_PARAMS() EQ 6 THEN $
+                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
+                RETURN
+              ENDIF ELSE BEGIN
+                ras = report('WARNING, The box does not contain any V points...')
+                firstxv = -1 & lastxv = -1 & nxv = 0
+              ENDELSE
+            ENDIF ELSE BEGIN
+              jxv = temporary(dom) MOD jpi
+              firstxv = min(temporary(jxv), max = lastxv)
+              nxv = lastxv - firstxv + 1
+            ENDELSE
+          ENDELSE
+          IF nxv NE 0 THEN lat1v = min(gphiv[firstxv:lastxv, firstyv:lastyv], max = lat2v)
+        ENDIF
+; find lat1f, lat2f, firstxf, firstxf, firstyf, firstyf, nxf
+; and nyf according to x1, x2, lon1 and lon2
+        IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN
+          firstyf = fsty & lastyf = lsty & nyf = ny
+          IF keyword_set(memeindices) THEN BEGIN
+            firstxf = firstyt & lastxf = lastyt & nxf = nxt
+          ENDIF ELSE BEGIN
+            dom = where( (glamf[*, fsty:lsty] GE lon1) AND (glamf[*, fsty:lsty] LE lon2) )
+            IF (dom[0] EQ -1) THEN BEGIN
+              IF keyword_set(findalways) THEN BEGIN
+                domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
+                IF N_PARAMS() EQ 6 THEN $
+                  domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex 
+                RETURN
+              ENDIF ELSE BEGIN
+                ras = report('WARNING, The box does not contain any F points...')
+                firstxf = -1 & lastyf = -1 & nxf = 0
+              ENDELSE
+            ENDIF ELSE BEGIN
+              jxf = temporary(dom) MOD jpi
+              firstxf = min(temporary(jxf), max = lastxf)
+              nxf = lastxf - firstxf + 1
+            ENDELSE
+          ENDELSE
+          IF nxf NE 0 THEN lat1f = min(gphif[firstxf:lastxf, firstyf:lastyf], max = lat2f)
+        ENDIF
+        lat1 = min([lat1t, lat1u, lat1v, lat1f])
+        lat2 = max([lat2t, lat2u, lat2v, lat2f])
+      END
+    ENDCASE
+  ENDELSE 
+;-------------------------------------------------------------------
+; The extracted domain is it regular or not? 
+;-------------------------------------------------------------------
+  CASE 1 OF
+    ((where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1) AND nxt NE 0 AND nyt NE 0:BEGIN
+; to get faster, we first test the most basic cases befor testing the
+; full array.
+      CASE 0 OF
+        array_equal(glamt[firstxt:lastxt, firstyt], glamt[firstxt:lastxt, lastyt]):key_irregular = 1b
+        array_equal(gphit[firstxt, firstyt:lastyt], gphit[lastxt, firstyt:lastyt]):key_irregular = 1b
+        array_equal(glamt[firstxt:lastxt, firstyt:lastyt] $
+                    , glamt[firstxt:lastxt, firstyt]#replicate(1, nyt)):key_irregular = 1b
+        array_equal(gphit[firstxt:lastxt, firstyt:lastyt] $
+                    , replicate(1, nxt)#(gphit[firstxt, firstyt:lastyt])[*]):key_irregular = 1b
+        ELSE:key_irregular = 0b
+      ENDCASE
+    END
+    (where(gridtype eq 'U'))[0] NE -1 AND nxu NE 0 AND nyu NE 0:BEGIN
+      CASE 0 OF
+        array_equal(glamu[firstxu:lastxu, firstyu], glamu[firstxu:lastxu, lastyu]):key_irregular = 1b
+        array_equal(gphiu[firstxu, firstyu:lastyu], gphiu[lastxu, firstyu:lastyu]):key_irregular = 1b
+        array_equal(glamu[firstxu:lastxu, firstyu:lastyu] $
+                    , glamu[firstxu:lastxu, firstyu]#replicate(1, nyu)):key_irregular = 1b
+        array_equal(gphiu[firstxu:lastxu, firstyu:lastyu] $
+                    , replicate(1, nxu)#(gphiu[firstxu, firstyu:lastyu])[*]):key_irregular = 1b
+        ELSE:key_irregular = 0b
+      ENDCASE
+    END
+    (where(gridtype eq 'V'))[0] NE -1 AND nxv NE 0 AND nyv NE 0:BEGIN
+      CASE 0 OF
+        array_equal(glamv[firstxv:lastxv, firstyv], glamv[firstxv:lastxv, lastyv]):key_irregular = 1b
+        array_equal(gphiv[firstxv, firstyv:lastyv], gphiv[lastxv, firstyv:lastyv]):key_irregular = 1b
+        array_equal(glamv[firstxv:lastxv, firstyv:lastyv] $
+                    , glamv[firstxv:lastxv, firstyv]#replicate(1, nyv)):key_irregular = 1b
+        array_equal(gphiv[firstxv:lastxv, firstyv:lastyv] $
+                    , replicate(1, nxv)#(gphiv[firstxv, firstyv:lastyv])[*]):key_irregular = 1b
+        ELSE:key_irregular = 0b
+      ENDCASE
+    END
+    (where(gridtype eq 'F'))[0] NE -1 AND nxf NE 0 AND nyf NE 0:BEGIN
+      CASE 0 OF
+        array_equal(glamf[firstxf:lastxf, firstyf], glamf[firstxf:lastxf, lastyf]):key_irregular = 1b
+        array_equal(gphif[firstxf, firstyf:lastyf], gphif[lastxf, firstyf:lastyf]):key_irregular = 1b
+        array_equal(glamf[firstxf:lastxf, firstyf:lastyf] $
+                    , glamf[firstxf:lastxf, firstyf]#replicate(1, nyf)):key_irregular = 1b
+        array_equal(gphif[firstxf:lastxf, firstyf:lastyf] $
+                    , replicate(1, nxf)#(gphif[firstxf, firstyf:lastyf])[*]):key_irregular = 1b
+        ELSE:key_irregular = 0b
+      ENDCASE
+    END
+    ELSE:
+  ENDCASE
+;
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+; define all vertical parameters ...
+; vert1, vert2
+; firstz[tw], lastz[tw], nz[tw]
+;-------------------------------------------------------------------
+;-------------------------------------------------------------------
+;
+  vertical:
+;
+;-------------------------------------------------------------------
+; vertical domain defined with vert1, vert2
+;-------------------------------------------------------------------
+  IF NOT (keyword_set(zindex) OR keyword_set(index)) THEN BEGIN
+; define vert1 et vert2
+    CASE N_PARAMS() OF
+      2:vert1 = min([x1, x2], max = vert2)
+      6:vert1 = min([z1, z2], max = vert2)
+      ELSE:BEGIN
+        IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN $
+          vert1t = min(gdept, max = vert2t)
+        IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN $
+          vert1w = min(gdepw, max = vert2w)
+        vert1 = min([vert1t, vert1w])
+        vert2 = max([vert2t, vert2w])
+      END
+    ENDCASE
+; define firstzt, firstzt, nzt
+    IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN
+      domz = where(gdept ge vert1 and gdept le vert2, nzt)
+      IF nzt NE 0 THEN BEGIN
+        firstzt = domz[0]
+        lastzt = domz[nzt-1]
       ENDIF ELSE BEGIN 
-         if zt[0] NE -1 then begin
-            premierzt = (where(gdept eq zt[0]))[0]
-            dernierzt = (where(gdept eq zt[nzt-1]))[0]
-         ENDIF ELSE BEGIN 
-            premierzt = -1
-            dernierzt = -1
-         ENDELSE
-      ENDELSE
-   endif
-   if (where(grille eq 'W'))[0] NE -1  AND n_elements(gdepw) NE 0 then begin
-      if keyword_set(zindex) OR keyword_set(index) then BEGIN 
-         premierzw = z1
-         dernierzw = z2
+        ras = report('WARNING, The box does not contain any T level...')
+        firstzt = -1
+        lastzt = -1
+      ENDELSE 
+    ENDIF
+; define firstzw, firstzw, nzw
+    IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
+      IF keyword_set(memeindices) THEN BEGIN
+        firstzw = firstzt
+        lastzw = lastzt
+        nzw = nzt
       ENDIF ELSE BEGIN 
-         if zw[0] NE -1 then begin
-            premierzw = (where(gdepw eq zw[0]))[0]
-            dernierzw = (where(gdepw eq zw[nzw-1]))[0]
-         ENDIF ELSE BEGIN 
-            premierzw = -1
-            dernierzw = -1
-         ENDELSE
-      ENDELSE
-   endif
-;-------------------------------------------------------------------
-;-------------------------------------------------------------------
-; calcul pour les pts t et w
-;-------------------------------------------------------------------
-;-------------------------------------------------------------------
-   tempdeux = systime(1)        ; pour key_performance =2
-   if (where(grille eq 'T'))[0] NE -1 or (where(grille EQ 'W'))[0] NE -1 then begin
-      if (keyword_set(xindex) AND keyword_set(yindex)) OR keyword_set(index) then begin
-         premierxt = x1
-         premieryt = y1
-         dernierxt = x2
-         dernieryt = y2
-      ENDIF ELSE BEGIN
-         dom =   where( (glamt ge lon1) $
-                        and (glamt le lon2) $
-                        and (gphit ge lat1) $
-                        and (gphit le lat2) )
-         if (dom[0] eq -1) then BEGIN
-            if keyword_set(findalways) then BEGIN
-               domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex
-               return
-            ENDIF ELSE BEGIN
-               ras = report('ATTENTION, la boite ne contient pas de points t')
-               goto, ptu
-            ENDELSE
-         endif 
-         IF testvar(var = key_performance) EQ 2 THEN $
-          print, 'temps domdef: trouver le sous domaine', systime(1)-tempdeux
-;-------------------------------------------------------------------
-; on recupere les abscisses et les ordonees des points selectionnes
-;-------------------------------------------------------------------
-         tempdeux = systime(1)  ; pour key_performance =2
-         ordonnee = dom/jpi
-         abscisse = dom-ordonnee*jpi
-;-------------------------------------------------------------------
-; on les ordonnes et on elimine les elements en double
-;-------------------------------------------------------------------
-         ordonnee = ordonnee
-         abscisse = abscisse
-         ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee)))
-         abscisse = abscisse(uniq(abscisse, sort(abscisse)))
-;-------------------------------------------------------------------
-; indices des bornes du domaine
-;-------------------------------------------------------------------
-         premierxt = abscisse[0]
-         premieryt = ordonnee[0]
-         dernierxt = abscisse[n_elements(abscisse) -1]
-         dernieryt = ordonnee[n_elements(ordonnee) -1]
+        domz = where(gdepw ge vert1 and gdepw le vert2, nzw)
+        IF nzw NE 0 THEN BEGIN
+          firstzw = domz[0]
+          lastzw = domz[nzw-1]
+        ENDIF ELSE BEGIN 
+          ras = report('WARNING, The box does not contain any W level...')
+          firstzw = -1
+          lastzw = -1
+        ENDELSE 
       ENDELSE 
-;-------------------------------------------------------------------
-; nombre de points selectionnes
-;-------------------------------------------------------------------
-      nxt = dernierxt-premierxt+1
-      nyt = dernieryt-premieryt+1
-      IF testvar(var = key_performance) EQ 2 THEN $
-       print, 'temps domdef: calculs de premier, dernier,...', systime(1)-tempdeux
-;-------------------------------------------------------------------
-; la grille est reguliere ou non ??
-;-------------------------------------------------------------------
-      tempdeux = systime(1)     ; pour key_performance =2
-      if (total(glamt[premierxt:dernierxt,premieryt:dernieryt] $
-                NE glamt[premierxt:dernierxt,0]#replicate(1, nyt)) eq 0 $
-          AND total(gphit[premierxt:dernierxt,premieryt:dernieryt] $
-                    NE replicate(1, nxt)#(gphit[0,premieryt:dernieryt])[*]) eq 0 ) $
-       then key_irregular = 0 ELSE key_irregular = 1
-      IF testvar(var = key_performance) EQ 2 THEN $
-       print, 'temps domdef: la grille est reguliere ou non ', systime(1)-tempdeux
-;-------------------------------------------------------------------
-      if keyword_set(memeindices) then begin
-         premierxu = premierxt & premierxv = premierxt & premierxf = premierxt 
-         dernierxu = dernierxt & dernierxv = dernierxt & dernierxf = dernierxt
-         premieryu = premieryt & premieryv = premieryt & premieryf = premieryt 
-         dernieryu = dernieryt & dernieryv = dernieryt & dernieryf = dernieryt
-         premierzw = premierzt 
-         dernierzw = dernierzt 
-         nxu = nxt & nxv = nxt & nxf = nxt
-         nyu = nyt & nyv = nyt & nyf = nyt
-         nzw = nzt
-         lon1 = glamt[premierxt, 0]
-         lon2 = glamf[dernierxf, 0]
-         lat1 = gphit[0, premieryt]
-         lat2 = gphif[0, dernieryf]
-         return
-      endif
-ptu:
-   endif
-;-------------------------------------------------------------------
-; calcul pour les pts u
-;-------------------------------------------------------------------
-   if (where(grille eq 'U'))[0] NE -1 then begin
-      if (keyword_set(xindex) AND keyword_set(xindex)) OR keyword_set(index) then begin
-         premierxu = x1
-         premieryu = y1
-         dernierxu = x2
-         dernieryu = y2
-      ENDIF ELSE BEGIN
-         dom =   where( (glamu ge lon1) $
-                        and (glamu le lon2) $
-                        and (gphiu ge lat1) $
-                        and (gphiu le lat2) )
-         if (dom[0] eq -1) then begin
-            if keyword_set(findalways) then begin
-               domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex
-               return
-            ENDIF ELSE BEGIN
-               ras = report('ATTENTION, la boite ne contient pas de points u')
-               goto, ptv
-            ENDELSE
-         ENDIF
-         ordonnee = dom/jpi
-         abscisse = dom-ordonnee*jpi
-         ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee)))
-         abscisse = abscisse(uniq(abscisse, sort(abscisse)))
-         premierxu = abscisse[0]
-         premieryu = ordonnee[0]
-         dernierxu = abscisse[n_elements(abscisse) -1]
-         dernieryu = ordonnee[n_elements(ordonnee) -1]
-      ENDELSE 
-      nxu = dernierxu-premierxu+1
-      nyu = dernieryu-premieryu+1
-      if (total(glamu[premierxu:dernierxu,premieryu:dernieryu] $
-                NE glamu[premierxu:dernierxu,0]#replicate(1, nyu)) eq 0 $
-          AND total(gphiu[premierxu:dernierxu,premieryu:dernieryu] $
-                    NE replicate(1, nxu)#(gphiu[0,premieryu:dernieryu])[*]) eq 0 ) $
-       then key_irregular = 0 ELSE key_irregular = 1
-ptv:
-   endif
-;-------------------------------------------------------------------
-; calcul pour les pts v
-;-------------------------------------------------------------------
-   if (where(grille eq 'V'))[0] NE -1 then begin
-      if (keyword_set(xindex) AND keyword_set(xindex)) OR keyword_set(index) then begin
-         premierxv = x1
-         premieryv = y1
-         dernierxv = x2
-         dernieryv = y2
-      ENDIF ELSE BEGIN
-         dom =   where( (glamv ge lon1) $
-                        and (glamv le lon2) $
-                        and (gphiv ge lat1) $
-                        and (gphiv le lat2) )
-         if (dom[0] eq -1) then begin
-            if keyword_set(findalways) then begin
-               domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex
-               return
-            ENDIF ELSE BEGIN
-               ras = report('ATTENTION, la boite ne contient pas de points v')
-               goto, ptf
-            ENDELSE
-         ENDIF
-         ordonnee = dom/jpi
-         abscisse = dom-ordonnee*jpi
-         ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee)))
-         abscisse = abscisse(uniq(abscisse, sort(abscisse)))
-         premierxv = abscisse[0]
-         premieryv = ordonnee[0]
-         dernierxv = abscisse[n_elements(abscisse) -1]
-         dernieryv = ordonnee[n_elements(ordonnee) -1]
-      ENDELSE 
-      nxv = dernierxv-premierxv+1
-      nyv = dernieryv-premieryv+1
-      if (total(glamv[premierxv:dernierxv,premieryv:dernieryv] $
-                NE glamv[premierxv:dernierxv,0]#replicate(1, nyv)) eq 0 $
-          AND total(gphiv[premierxv:dernierxv,premieryv:dernieryv] $
-                    NE replicate(1, nxv)#(gphiv[0,premieryv:dernieryv])[*]) eq 0 ) $
-       then key_irregular = 0 ELSE key_irregular = 1
-ptf:
-   endif
-;-------------------------------------------------------------------
-; calcul pour les pts f
-;-------------------------------------------------------------------
-   if (where(grille eq 'F'))[0] NE -1 then begin
-      if (keyword_set(xindex) AND keyword_set(xindex)) OR keyword_set(index) then begin
-         premierxf = x1
-         premieryf = y1
-         dernierxf = x2
-         dernieryf = y2
-      ENDIF ELSE BEGIN
-         dom =   where( (glamf ge lon1) $
-                        and (glamf le lon2) $
-                        and (gphif ge lat1) $
-                        and (gphif le lat2) )
-         if (dom[0] eq -1) then begin
-            if keyword_set(findalways) then begin
-               domdef, prof1,prof2, GRILLE=grille, MEMEINDICES = memeindices, _extra = ex
-               return
-            ENDIF ELSE BEGIN
-               ras = report('ATTENTION, la boite ne contient pas de points f')
-               return
-            ENDELSE
-         ENDIF
-         ordonnee = dom/jpi
-         abscisse = dom-ordonnee*jpi
-         ordonnee = ordonnee(uniq(ordonnee, sort(ordonnee)))
-         abscisse = abscisse(uniq(abscisse, sort(abscisse)))
-         premierxf = abscisse[0]
-         premieryf = ordonnee[0]
-         dernierxf = abscisse[n_elements(abscisse) -1]
-         dernieryf = ordonnee[n_elements(ordonnee) -1]
-      ENDELSE
-      nxf = dernierxf-premierxf+1
-      nyf = dernieryf-premieryf+1
-      if (total(glamf[premierxf:dernierxf,premieryf:dernieryf] $
-                NE glamf[premierxf:dernierxf,0]#replicate(1, nyf)) eq 0 $
-          AND total(gphif[premierxf:dernierxf,premieryf:dernieryf] $
-                    NE replicate(1, nxf)#(gphif[0,premieryf:dernieryf])[*]) eq 0 ) $
-       then key_irregular = 0 ELSE key_irregular = 1
-   endif
-;-------------------------------------------------------------------
-   if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun 
+    ENDIF
+;-------------------------------------------------------------------
+; vertical domain defined with the Z index
+;-------------------------------------------------------------------
+  ENDIF ELSE BEGIN 
+    CASE N_PARAMS() OF
+      2:fstz = min([x1, x2], max = lstz)
+      4:return
+      6:fstz = min([z1, z2], max = lstz)
+    ENDCASE
+    IF fstz LT 0 OR lstz GE jpk THEN BEGIN
+      ras = report('Bad definition of X1, X2, Z1 or Z2')
+      return
+    ENDIF
+    nz = lstz - fstz + 1
+; find vert1t, vert2t, firstzt, firstzt, nzt
+; according to (x1, x2) or (z1, z2)
+    IF (where(gridtype eq 'T'))[0] NE -1 THEN BEGIN
+      vert1t = min(gdept[fstz:lstz], max = vert2t)
+      firstzt = fstz & lastzt = lstz & nzt = nz
+    ENDIF      
+; find vert1w, vert2w, firstzw, firstzw, nzw
+; according to (x1, x2) or (z1, z2)
+    IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
+      vert1w = min(gdepw[fstz:lstz], max = vert2w)
+      firstzw = fstz & lastzw = lstz & nzw = nz
+    ENDIF      
+    vert1 = min([vert1t, vert1w])
+    vert2 = max([vert2t, vert2w])
+  ENDELSE
+;
+;-------------------------------------------------------------------
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updateold
+  ENDIF 
+;-------------------------------------------------------------------
+  if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun 
 ;-------------------------------------------------------------------
 
 ;-------------------------------------------------------------------
-   return
+  return
 end

trunk/ToBeReviewed/GRILLE/f2v.pro

@@ -38 +38 @@
 ;------------------------------------------------------------
 FUNCTION f2v, temp
-@common
+;---------------------------------------------------------
+@cm_4mesh
+@cm_4data
+@cm_4cal
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
    res = temp
 ;on force nxt=nxf, etc ...
-   premierxv = premierxf
-   dernierxv = dernierxf
-   premieryv = premieryf
-   dernieryv = dernieryf
+   firstxv = firstxf
+   lastxv = lastxf
+   firstyv = firstyf
+   lastyv = lastyf
    nxv = nxf
    nyv = nyf
    vargrid = 'V'
    if NOT keyword_set(valmask) then valmask = 1e20
-   lon1 = glamv[premierxv, 0]
-   lon2 = glamf[dernierxf, 0]
+   lon1 = glamv[firstxv, 0]
+   lon2 = glamf[lastxf, 0]
 
 ; cas sur la taille du tableau et application 
@@ -60 +67 @@
             taille[1] eq nxf and taille[2] eq nyf:
             taille[1] eq jpi and taille[2] eq jpj: $
-             res=res[premierxf:dernierxf, premieryf:dernieryf]
+             res=res[firstxf:lastxf, firstyf:lastyf]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          endcase
-         mask = (fmask())[premierxf:dernierxf, premieryf:dernieryf, 0]
+         mask = (fmask())[firstxf:lastxf, firstyf:lastyf, 0]
          terre = where(mask EQ 0) 
          IF terre[0] NE -1 THEN res[terre] = !values.f_nan
          res = 0.5*(res + shift(res, 1, 0))
-         if NOT (keyword_set(key_periodique) AND nxf EQ jpi) then res[0, *] = !values.f_nan
-         mask = (vmask())[premierxf:dernierxf, premieryf:dernieryf, 0]
+         if NOT (keyword_set(key_periodic) AND nxf EQ jpi) then res[0, *] = !values.f_nan
+         mask = (vmask())[firstxf:lastxf, firstyf:lastyf, 0]
          terre = where(mask EQ 0) 
          IF terre[0] NE -1 THEN res[terre] = valmask
@@ -77 +84 @@
             taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ nzt:
             taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ jpk: $
-             res=res[*, *, premierzt:dernierzt]
+             res=res[*, *, firstzt:lastzt]
             taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ jpt:
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $
-             res=res[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]
+             res=res[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt]
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $
-             res=res[premierxf:dernierxf, premieryf:dernieryf, *]
+             res=res[firstxf:lastxf, firstyf:lastyf, *]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          ENDCASE
          if taille[3] EQ jpt then begin
-            mask = (fmask())[premierxf:dernierxf, premieryf:dernieryf, dernierzt*(nzt NE jpk)]
+            mask = (fmask())[firstxf:lastxf, firstyf:lastyf, lastzt*(nzt NE jpk)]
             mask = temporary(mask[*])#replicate(1, jpt)
             mask = reform(mask, nxf, nyf, jpt, /over)
-         ENDIF ELSE mask = (fmask())[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]
+         ENDIF ELSE mask = (fmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt]
          terre = where(temporary(mask) EQ 0) 
          IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan
          res = 0.5*(res + shift(res, 1, 0, 0))
-         if NOT (keyword_set(key_periodique) AND nxf EQ jpi) then res[0, *, *] = !values.f_nan
+         if NOT (keyword_set(key_periodic) AND nxf EQ jpi) then res[0, *, *] = !values.f_nan
          if taille[3] EQ jpt then BEGIN
-            mask = tmask[premierxf:dernierxf, premieryf:dernieryf, dernierzt*(nzt NE jpk)]
+            mask = tmask[firstxf:lastxf, firstyf:lastyf, lastzt*(nzt NE jpk)]
             mask = temporary(mask[*])#replicate(1, jpt)
             mask = reform(mask, nxf, nyf, jpt, /over)
-         ENDIF ELSE mask = (vmask())[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]
+         ENDIF ELSE mask = (vmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt]
          terre = where(temporary(mask) EQ 0) 
          IF terre[0] NE -1 THEN res[temporary(terre)] = valmask
@@ -107 +114 @@
             taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ nzt AND taille[4] EQ jpt:
             taille[1] eq nxf and taille[2] eq nyf AND taille[3] EQ jpk AND taille[4] EQ jpt: $
-             res=res[*, *, premierzt:dernierzt, *]
+             res=res[*, *, firstzt:lastzt, *]
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $
-             res=res[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt, *]
+             res=res[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt, *]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          ENDCASE
-         mask = (fmask())[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]
+         mask = (fmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt]
          mask = temporary(mask[*])#replicate(1, jpt)
          mask = reform(mask, nxf, nyf, nzt, jpt, /over)
@@ -119 +126 @@
          IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan
          res = 0.5*(res + shift(res, 1, 0, 0, 0))
-         if NOT (keyword_set(key_periodique) AND nxf EQ jpi) then res[0, *, *, *] = !values.f_nan
-         mask = (vmask())[premierxf:dernierxf, premieryf:dernieryf, premierzt:dernierzt]
+         if NOT (keyword_set(key_periodic) AND nxf EQ jpi) then res[0, *, *, *] = !values.f_nan
+         mask = (vmask())[firstxf:lastxf, firstyf:lastyf, firstzt:lastzt]
          mask = temporary(mask[*])#replicate(1, jpt)
          mask = reform(mask, nxf, nyf, nzt, jpt, /over)
@@ -128 +135 @@
    endcase
 
+   IF NOT keyword_set(key_forgetold) THEN BEGIN
+   @updateold
+   ENDIF 
+   
    return, res
 END

trunk/ToBeReviewed/GRILLE/fmask.pro

@@ -25 +25 @@
 ;------------------------------------------------------------
 FUNCTION fmask
-@common
-   tempsun = systime(1)         ; pour key_performance
-   if jpk EQ 1 then begin
-      res = tmask*shift(tmask, -1, 0)*shift(tmask, 0, -1)*shift(tmask, -1, -1)
-      if NOT keyword_set(key_periodique) then res[jpi-1, *] = fmaskredy
-      res[*, jpj-1] = fmaskredx
-   ENDIF ELSE BEGIN
-      res = tmask*shift(tmask, -1, 0, 0)*shift(tmask, 0, -1, 0)*shift(tmask, -1, -1, 0)
-      if NOT keyword_set(key_periodique) then res[jpi-1, *, *] = fmaskredy
-      res[*, jpj-1, *] = fmaskredx
-   ENDELSE
-   if keyword_set(key_performance) THEN print, 'temps fmask', systime(1)-tempsun 
-
-   return, res
+;---------------------------------------------------------
+@cm_4mesh
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
+  tempsun = systime(1)          ; pour key_performance
+; 
+  CASE size(tmask, /n_dimensions) OF
+    2:res = tmask*shift(tmask, -1, 0)*shift(tmask, 0, -1)*shift(tmask, -1, -1)
+    3:res = tmask*shift(tmask, -1, 0, 0)*shift(tmask, 0, -1, 0)*shift(tmask, -1, -1, 0)
+  ENDCASE
+;
+  if NOT keyword_set(key_periodic) then res[jpi-1, *, *] = fmaskredy
+  res[*, jpj-1, *] = fmaskredx
+;
+  if keyword_set(key_performance) THEN print, 'temps fmask', systime(1)-tempsun 
+  
+  return, res
 end

trunk/ToBeReviewed/GRILLE/grille.pro

@@ -13 +13 @@
 ;
 ; CALLING SEQUENCE:
-;  grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz,e1,e2,e3
+;  grille,mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz,e1,e2,e3
 ;
 ; INPUTS:rien. ATTENTION les choix de la grille se fait a partir de la
@@ -32 +32 @@
 ;         /NOTRI: utile seulement qd TRI est active. dans ce cas
 ;         grille retourne -1 ds la variable tri meme si la variable du
-;         common triangles est definie et differente de -1
-;
-; OUTPUTS:mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,
-;         dernierx,derniery,dernierz,e1,e2,e3
+;         common triangles_list est definie et differente de -1
+;
+;         /WDEPTH: to specify that the field is at W depth instad of T 
+;         depth (automatically activated if vargrid eq 'W')
+;
+; OUTPUTS:mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz,
+;         lastx,lasty,lastz,e1,e2,e3
 ;
 ;         pour leur definition cf domdef et la gestion des sous
@@ -43 +46 @@
 ;         mask, glam et gphi il suffit de taper grille, mask, glam, gphi
 ;
-; COMMON BLOCKS:
-;       common.pro      congridseb.pro
+; COMMON BLOCKS: cm_4mesh and cm_4data
 ;
 ; SIDE EFFECTS: utilise la variable globale vargird
@@ -59 +61 @@
 ;------------------------------------------------------------
 ;------------------------------------------------------------
-pro grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz, e1,e2,e3, TRI = tri, NOTRI = notri, TOUT = tout, FORPLT = forplt, _EXTRA = ex
-@common
-   tempsun = systime(1)         ; pour key_performance
-;------------------------------------------------------------
-   if keyword_set(tout) then begin
-      oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
-      domdef, grille = vargrid, _EXTRA = ex
-   endif
-   tempdeux = systime(1)        ; pour key_performance =2
-   CASE 1 OF 
-;------------------------------------------------------------
-; grille T
-;------------------------------------------------------------
-      vargrid eq 'OPAPTDHT' or vargrid eq 'OPAPT3DT' $
-       or vargrid eq 'T': begin
+pro grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, TRI = tri, NOTRI = notri, TOUT = tout, FORPLT = forplt, IFPLTZ = ifpltz, WDEPTH = wdepth, _EXTRA = ex
+;------------------------------------------------------------
+; include commons
+@cm_4mesh
+@cm_4data
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------
+  tempsun = systime(1)          ; pour key_performance
+;------------------------------------------------------------
+  vargrid = strupcase(strmid(vargrid,0,/reverse_offset))
+;
+  if vargrid eq 'W' then wdepth = 1
+  if keyword_set(tout) then begin
+    savedbox = 1b
+    saveboxparam, 'boxparam4grille.dat'
+    domdef, gridtype = vargrid, _EXTRA = ex
+  endif
+  tempdeux = systime(1)         ; pour key_performance =2
+;------------------------------------------------------------
+;------------------------------------------------------------
+  IF keyword_set(wdepth) THEN BEGIN
+    firstz = firstzw 
+    lastz = lastzw
+    nz = nzw
+  ENDIF ELSE BEGIN
+    firstz = firstzt
+    lastz = lastzt
+    nz = nzt
+  ENDELSE
+;------------------------------------------------------------
+;------------------------------------------------------------
+  CASE 1 OF 
+;------------------------------------------------------------
+; grille T and W
+;------------------------------------------------------------
+    vargrid eq 'T' OR vargrid eq 'W' : begin
 ;scalaires
-         nx=nxt
-         ny=nyt
-         nz=nzt
-         premierx = premierxt
-         premiery = premieryt
-         premierz = premierzt
-         dernierx = dernierxt
-         derniery = dernieryt
-         dernierz = dernierzt
+      nx = nxt
+      ny = nyt
+      firstx = firstxt
+      firsty = firstyt
+      lastx = lastxt
+      lasty = lastyt
 ;vecteurs 2d
-         glam=glamt[premierx:dernierx, premiery:derniery]
-         gphi=gphit[premierx:dernierx, premiery:derniery]
-         e1  =e1t[premierx:dernierx, premiery:derniery]
-         e2  =e2t[premierx:dernierx, premiery:derniery]
+      IF n_elements(glam) NE 1 THEN glam = glamt[firstx:lastx, firsty:lasty]
+      IF n_elements(gphi) NE 1 THEN gphi = gphit[firstx:lastx, firsty:lasty]
+      IF n_elements(e1) NE 1 THEN e1 = e1t[firstx:lastx, firsty:lasty]
+      IF n_elements(e2) NE 1 THEN e2 = e2t[firstx:lastx, firsty:lasty]
 ;vecteurs 3d
-         mask=tmask[premierx:dernierx, premiery:derniery, premierz:dernierz]
-      end
-;------------------------------------------------------------
-; grille W
-;------------------------------------------------------------
-      vargrid eq 'OPAPT3DW' or vargrid eq 'W': begin
+      IF keyword_set(forplt) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz] $
+      ELSE IF n_elements(mask) NE 1 THEN mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
+    end
+;------------------------------------------------------------
+; grille U
+;------------------------------------------------------------
+    vargrid eq 'U': begin
 ;scalaires
-         nx=nxt
-         ny=nyt
-         nz=nzw
-         premierx = premierxt
-         premiery = premieryt
-         premierz = premierzw
-         dernierx = dernierxt
-         derniery = dernieryt
-         dernierz = dernierzw
+      nx = nxu
+      ny = nyu
+      firstx = firstxu
+      firsty = firstyu
+      lastx = lastxu
+      lasty = lastyu
 ;vecteurs 2d
-         terre = where(tmask[*, *, 0] EQ 0)
-         glam=glamt[premierx:dernierx, premiery:derniery]
-         gphi=gphit[premierx:dernierx, premiery:derniery]
-         e1  =e1t[premierx:dernierx, premiery:derniery]
-         e2  =e2t[premierx:dernierx, premiery:derniery]
+      IF n_elements(glam) NE 1 THEN glam = glamu[firstx:lastx, firsty:lasty]
+      IF n_elements(gphi) NE 1 THEN gphi = gphiu[firstx:lastx, firsty:lasty]
+      if keyword_set(forplt) then BEGIN
+        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
+        eastboarder = mask-shift(mask, 1, 0)*mask
+        westboarder = mask-shift(mask, -1, 0)*mask
+        if key_periodic NE 1 OR nx NE jpi then westboarder[nx-1, *] = 0b
+        tmp1 = shift(eastboarder, 0, 1)
+        tmp1[*, 0] = 0b
+        tmp2 = shift(eastboarder, 0, -1)
+        tmp2[*, ny-1] = 0b
+        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
+        eastboarder = temporary(eastboarder)+temporary(add)
+        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
+        tmp1[*, ny-1] = 1b
+        tmp1[*, 0] = 1b
+        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
+        if key_periodic NE 1 OR nx NE jpi then begin
+          tmp2[nx-1, *] = 1b
+          tmp2[0, *] = 0b
+        endif
+        no1 = temporary(tmp1)*temporary(tmp2)
+        tmp = temporary(eastboarder)*temporary(no1)*mask
+        mask[0:nx-2, *] = 0b 
+        tmp = temporary(tmp)+temporary(mask)
+        tmp = where(tmp GE 1)
+        if tmp[0] NE -1 then begin
+          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
+          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
+        endif
+      ENDIF
+      IF n_elements(e1) NE 1 THEN e1  = e1u[firstx:lastx, firsty:lasty]
+      IF n_elements(e2) NE 1 THEN e2  = e2u[firstx:lastx, firsty:lasty]
 ;vecteurs 3d
-         mask=tmask[premierx:dernierx, premiery:derniery, premierz:dernierz]
-      end
-;------------------------------------------------------------
-; grille U
-;------------------------------------------------------------
-      vargrid eq 'OPAPTDHU' or vargrid eq 'OPAPT3DU' $
-       or vargrid eq 'U': begin
+      IF keyword_set(forplt) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz] $
+      ELSE IF n_elements(mask) NE 1 THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
+    end
+;------------------------------------------------------------
+; grille V
+;------------------------------------------------------------
+    vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $
+      or vargrid eq 'V': begin
 ;scalaires
-         nx=nxu
-         ny=nyu
-         nz=nzt
-         premierx = premierxu
-         premiery = premieryu
-         premierz = premierzt
-         dernierx = dernierxu
-         derniery = dernieryu
-         dernierz = dernierzt
+      nx = nxv
+      ny = nyv
+      firstx = firstxv
+      firsty = firstyv
+      lastx = lastxv
+      lasty = lastyv
 ;vecteurs 2d
-         glam=glamu[premierx:dernierx, premiery:derniery]
-         gphi=gphiu[premierx:dernierx, premiery:derniery]
-         if keyword_set(forplt) then BEGIN 
-            mask = tmask[premierx:dernierx, premiery:derniery, niveau-1]
-            terre = where(mask EQ 0)
-            if terre[0] NE -1 then begin
-               glam[terre] = (glamt[premierx:dernierx, premiery:derniery])[terre]
-               gphi[terre] = (gphit[premierx:dernierx, premiery:derniery])[terre]
-            endif
-         ENDIF
-         e1  =e1u[premierx:dernierx, premiery:derniery]
-         e2  =e2u[premierx:dernierx, premiery:derniery]
+      IF n_elements(glam) NE 1 THEN glam = glamv[firstx:lastx, firsty:lasty]
+      IF n_elements(gphi) NE 1 THEN gphi = gphiv[firstx:lastx, firsty:lasty]
+      if keyword_set(forplt) then BEGIN 
+        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
+        northboarder = mask-shift(mask, 0, 1)*mask
+        southboarder = mask-shift(mask, 0, -1)*mask
+        southboarder[*, ny-1] = 0b
+        tmp1 = shift(northboarder, -1, 0)
+        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
+        tmp2 = shift(northboarder, 1, 0)
+        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
+        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
+        northboarder = temporary(northboarder)+temporary(add)
+        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
+        tmp1[*, ny-1] = 1b
+        tmp1[*, 0] = 0b
+        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
+        if key_periodic NE 1 OR nx NE jpi then begin
+          tmp2[nx-1, *] = 1b
+          tmp2[0, *] = 1b
+        endif
+        no1 = temporary(tmp1)*temporary(tmp2)
+        tmp = temporary(northboarder)*mask*temporary(no1)
+        mask[*, 0:ny-2] = 0b
+        tmp = temporary(tmp)+temporary(mask)
+        tmp = where(tmp GE 1)
+        if tmp[0] NE -1 then begin
+          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
+          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
+        endif
+      ENDIF
+      IF n_elements(e1) NE 1 THEN e1  = e1v[firstx:lastx, firsty:lasty]
+      IF n_elements(e2) NE 1 THEN e2  = e2v[firstx:lastx, firsty:lasty]
 ;vecteurs 3d
-         mask=(umask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
-      end
-;------------------------------------------------------------
-; grille V
-;------------------------------------------------------------
-      vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $
-       or vargrid eq 'V': begin
+      IF keyword_set(forplt) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz] $
+      ELSE IF n_elements(mask) NE 1 THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
+    end
+;------------------------------------------------------------
+; grille F
+;------------------------------------------------------------
+    vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $ 
+      or vargrid eq 'F': begin
 ;scalaires
-         nx=nxv
-         ny=nyv
-         nz=nzt
-         premierx = premierxv
-         premiery = premieryv
-         premierz = premierzt
-         dernierx = dernierxv
-         derniery = dernieryv
-         dernierz = dernierzt
+      nx = nxf
+      ny = nyf
+      firstx = firstxf
+      firsty = firstyf
+      lastx = lastxf
+      lasty = lastyf
 ;vecteurs 2d
-         glam=glamv[premierx:dernierx, premiery:derniery]
-         gphi=gphiv[premierx:dernierx, premiery:derniery]
-         if keyword_set(forplt) then BEGIN 
-            mask = tmask[premierx:dernierx, premiery:derniery, niveau-1]
-            terre = where(mask EQ 0)
-            if terre[0] NE -1 then begin
-               glam[terre] = (glamt[premierx:dernierx, premiery:derniery])[terre]
-               gphi[terre] = (gphit[premierx:dernierx, premiery:derniery])[terre]
-            endif
-         ENDIF
-         e1  =e1v[premierx:dernierx, premiery:derniery]
-         e2  =e2v[premierx:dernierx, premiery:derniery]
+      IF n_elements(glam) NE 1 THEN glam = glamf[firstx:lastx, firsty:lasty]
+      IF n_elements(gphi) NE 1 THEN gphi = gphif[firstx:lastx, firsty:lasty]
+      if keyword_set(forplt) then BEGIN 
+        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
+        eastboarder = mask-shift(mask, 1, 0)*mask
+        westboarder = mask-shift(mask, -1, 0)*mask
+        westboarder[nx-1, *] = 0b
+        northboarder = mask-shift(mask, 0, 1)*mask
+        southboarder = mask-shift(mask, 0, -1)*mask
+        southboarder[*, ny-1] = 0b
+        tmp1 = shift(northboarder, -1, 0)
+        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
+        tmp2 = shift(northboarder, 1, 0)
+        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
+        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
+        northboarder = temporary(northboarder)+temporary(add)
+        tmp1 = shift(eastboarder, 0, 1)
+        tmp1[*, 0] = 0b
+        tmp2 = shift(eastboarder, 0, -1)
+        tmp2[*, ny-1] = 0b
+        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
+        eastboarder = temporary(eastboarder)+temporary(add)
+        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
+        tmp1[*, ny-1] = 1b
+        tmp1[*, 0] = 1b
+        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
+        if key_periodic NE 1 OR nx NE jpi then begin
+          tmp2[nx-1, *] = 1b
+          tmp2[0, *] = 1b
+        endif
+        no1 = temporary(tmp1)*temporary(tmp2)
+        tmp = (temporary(northboarder)+temporary(eastboarder))*mask*temporary(no1)
+        mask[0:nx-2, *] = 0b 
+        mask[*, 0:ny-2] = 0b 
+        tmp = temporary(tmp)+temporary(mask)
+        tmp = where(tmp GE 1)
+        if tmp[0] NE -1 then begin
+          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
+          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
+        endif
+      ENDIF
+      IF n_elements(e1) NE 1 THEN e1  = e1f[firstx:lastx, firsty:lasty]
+      IF n_elements(e2) NE 1 THEN e2  = e2f[firstx:lastx, firsty:lasty]
 ;vecteurs 3d
-         mask=(vmask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
-      end
-;------------------------------------------------------------
-; grille F
-;------------------------------------------------------------
-      vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $ 
-       or vargrid eq 'F': begin
-;scalaires
-         nx=nxf
-         ny=nyf
-         nz=nzt
-         premierx = premierxf
-         premiery = premieryf
-         premierz = premierzt
-         dernierx = dernierxf
-         derniery = dernieryf
-         dernierz = dernierzt
-;vecteurs 2d
-         glam=glamf[premierx:dernierx, premiery:derniery]
-         gphi=gphif[premierx:dernierx, premiery:derniery]
-         if keyword_set(forplt) then BEGIN 
-            mask = tmask[premierx:dernierx, premiery:derniery, niveau-1]
-            terre = where(mask EQ 0)
-            if terre[0] NE -1 then begin
-               glam[terre] = (glamt[premierx:dernierx, premiery:derniery])[terre]
-               gphi[terre] = (gphit[premierx:dernierx, premiery:derniery])[terre]
-            endif
-         ENDIF
-         e1  =e1f[premierx:dernierx, premiery:derniery]
-         e2  =e2f[premierx:dernierx, premiery:derniery]
-;vecteurs 3d
-         mask=(fmask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
+      IF keyword_set(forplt) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz] $
+      ELSE IF n_elements(mask) NE 1 THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
+    END
+;------------------------------------------------------------
+    ELSE:BEGIN 
+      ras = report('Wrong definition of vargrid = '+vargrid+'. Only T, U, V, W or F are acceptable')
+      stop
+    END
+  ENDCASE
+  IF testvar(var = key_performance) EQ 2 THEN $
+    print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux
+;
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+; Variables se rapportant a la dimension verticale
+;------------------------------------------------------------
+;------------------------------------------------------------
+;------------------------------------------------------------
+;
+;
+  tempdeux = systime(1)         ; pour key_performance =2
+  if keyword_set(wdepth) then begin
+    gdep = gdepw[firstz:lastz]
+    e3 = e3w[firstz:lastz]
+  endif else begin
+    gdep = gdept[firstz:lastz]
+    e3 = e3t[firstz:lastz]
+  ENDELSE
+; for the vertical sections with partial steps
+  IF keyword_set(ifpltz) AND keyword_set(key_partialstep) THEN BEGIN
+    CASE 1 OF
+      ifpltz EQ 'xz' AND ny EQ 1:BEGIN
+        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
+        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
+        bottom = lindgen(nx)+(bottom-1l+keyword_set(wdepth))*nx
+        IF good[0] NE -1 THEN BEGIN
+          bottom = bottom[good]
+          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
+          gdep = replicate(1, nx)#gdep
+          if keyword_set(wdepth) THEN $
+            truegdep = hdepw[firstx:lastx, firsty:lasty] $
+          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
+          gdep[bottom] = truegdep[good]
+        ENDIF
       END
-;------------------------------------------------------------
-      ELSE:BEGIN 
-         ras = report('Vauvaise definition de la variable vargrid: '+vargrid+'ceete variable doit etre T, U, V, W ou F')
-         stop
+      ifpltz EQ 'yz' AND nx EQ 1:BEGIN
+        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
+        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
+        bottom = lindgen(ny)+(bottom-1l+keyword_set(wdepth))*ny
+        IF good[0] NE -1 THEN BEGIN
+          bottom = bottom[good]
+          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
+          gdep = replicate(1, ny)#gdep
+          if keyword_set(wdepth) THEN $
+            truegdep = hdepw[firstx:lastx, firsty:lasty] $
+          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
+          gdep[bottom] = truegdep[good]
+        ENDIF
       END
-   ENDCASE
-   IF testvar(var = key_performance) EQ 2 THEN $
-    print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux
-;------------------------------------------------------------
-; Variables se rapportant a la dimension verticale
-;------------------------------------------------------------
-   tempdeux = systime(1)        ; pour key_performance =2
-   if vargrid eq 'OPAPT3DW' or vargrid eq 'W' then begin
-      gdep = gdepw[premierz:dernierz]
-      e3=e3w[premierz:dernierz]
-   endif else begin
-      gdep = gdept[premierz:dernierz]
-      e3=e3t[premierz:dernierz]
-   endelse
-   IF testvar(var = key_performance) EQ 2 THEN $
+      ELSE:
+    ENDCASE
+  ENDIF
+  IF testvar(var = key_performance) EQ 2 THEN $
     print, 'temps grille: Variables se rapportant a la dimension verticale ', systime(1)-tempdeux
 ;------------------------------------------------------------
 ; vecteur triangulation Qd TRI est active
 ;------------------------------------------------------------
-   if arg_present(TRI) then $
-    if triangles[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN 
-      tempdeux = systime(1)     ; pour key_performance =2
-      msk = bytarr(jpi, jpj)
-      msk[premierx:dernierx,premiery:derniery] = 1
-      ind = where( msk[triangles[0, *]]*msk[triangles[1, *]]*msk[triangles[2, *]] EQ 1 )
-      tri =triangles[*, ind]-(premierx+premiery*jpi)
-      y = tri/jpi
-      x = tri-y*jpi
-      tri = x+y*nx
-      IF testvar(var = key_performance) EQ 2 THEN $
-       print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux
-   ENDELSE
+  if arg_present(TRI) then $
+    if triangles_list[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN 
+    tempdeux = systime(1)       ; pour key_performance =2
+    msk = bytarr(jpi, jpj)
+    msk[firstx:lastx, firsty:lasty] = 1
+    ind = where( msk[triangles_list[0, *]]*msk[triangles_list[1, *]]*msk[triangles_list[2, *]] EQ 1 )
+    tri = triangles_list[*, ind]-(firstx+firsty*jpi)
+    y = tri/jpi
+    x = tri-y*jpi
+    tri = x+y*nx
+    IF testvar(var = key_performance) EQ 2 THEN $
+      print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux
+  ENDELSE
 ;------------------------------------------------------------------
 ; pour s'assurer qu'il n'y a pas de dimension degenerees (=1)
@@ -256 +363 @@
 ;    e3=reform(e3, /over)
 
-   if keyword_set(tout) then domdef, oldboite, grille = vargrid
-   if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun 
-
-   return
+  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grille.dat'
+  if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun 
+
+;------------------------------------------------------------
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updateold
+  ENDIF 
+;---------------------
+  return
 
 end

trunk/ToBeReviewed/GRILLE/t2v.pro

@@ -38 +38 @@
 ;------------------------------------------------------------
 FUNCTION t2v, temp
-@common
+;---------------------------------------------------------
+@cm_4mesh
+@cm_4data
+@cm_4cal
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
    res = temp
 
 ;on force nxt=nxv, etc ...
-   premierxv = premierxt
-   dernierxv = dernierxt
-   premieryv = premieryt
-   dernieryv = dernieryt
+   firstxv = firstxt
+   lastxv = lastxt
+   firstyv = firstyt
+   lastyv = lastyt
    nxv = nxt
    nyv = nyt
    vargrid = 'V'
    if NOT keyword_set(valmask) then valmask = 1e20
-   lat1 = gphit[0, premieryt]
-   lat2 = gphiv[0, dernieryv]
+   lat1 = gphit[0, firstyt]
+   lat2 = gphiv[0, lastyv]
 
 ; cas sur la taille du tableau et application 
@@ -61 +68 @@
             taille[1] eq nxt and taille[2] eq nyt:
             taille[1] eq jpi and taille[2] eq jpj: $
-             res=res[premierxt:dernierxt, premieryt:dernieryt]
+             res=res[firstxt:lastxt, firstyt:lastyt]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          endcase
-         mask = tmask[premierxt:dernierxt, premieryt:dernieryt, 0]
+         mask = tmask[firstxt:lastxt, firstyt:lastyt, 0]
          terre = where(mask EQ 0) 
          IF terre[0] NE -1 THEN res[terre] = !values.f_nan
          res = 0.5*(res + shift(res, 0, -1))
          res[*, nyt-1] = !values.f_nan
-         mask = (vmask())[premierxt:dernierxt, premieryt:dernieryt, 0]
+         mask = (vmask())[firstxt:lastxt, firstyt:lastyt, 0]
          terre = where(mask EQ 0) 
          IF terre[0] NE -1 THEN res[terre] = valmask
@@ -78 +85 @@
             taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ nzt:
             taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ jpk: $
-             res=res[*, *, premierzt:dernierzt]
+             res=res[*, *, firstzt:lastzt]
             taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ jpt:
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $
-             res=res[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]
+             res=res[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt]
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $
-             res=res[premierxt:dernierxt, premieryt:dernieryt, *]
+             res=res[firstxt:lastxt, firstyt:lastyt, *]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          ENDCASE
          if taille[3] EQ jpt then begin
-            mask = tmask[premierxt:dernierxt, premieryt:dernieryt, dernierzt*(nzt NE jpk)]
+            mask = tmask[firstxt:lastxt, firstyt:lastyt, lastzt*(nzt NE jpk)]
             mask = temporary(mask[*])#replicate(1, jpt)
             mask = reform(mask, nxt, nyt, jpt, /over)
-         ENDIF ELSE mask = tmask[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]
+         ENDIF ELSE mask = tmask[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt]
          terre = where(temporary(mask) EQ 0) 
          IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan
@@ -97 +104 @@
          res[*, nyt-1, *] = !values.f_nan
          if taille[3] EQ jpt then BEGIN
-            mask = (vmask())[premierxt:dernierxt, premieryt:dernieryt, dernierzt*(nzt NE jpk)]
+            mask = (vmask())[firstxt:lastxt, firstyt:lastyt, lastzt*(nzt NE jpk)]
             mask = temporary(mask[*])#replicate(1, jpt)
             mask = reform(mask, nxt, nyt, jpt, /over)
-         ENDIF ELSE mask = (vmask())[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]
+         ENDIF ELSE mask = (vmask())[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt]
          terre = where(temporary(mask) EQ 0) 
          IF terre[0] NE -1 THEN res[temporary(terre)] = valmask
@@ -108 +115 @@
             taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ nzt AND taille[4] EQ jpt:
             taille[1] eq nxt and taille[2] eq nyt AND taille[3] EQ jpk AND taille[4] EQ jpt: $
-             res=res[*, *, premierzt:dernierzt, *]
+             res=res[*, *, firstzt:lastzt, *]
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $
-             res=res[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt, *]
+             res=res[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt, *]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          ENDCASE
-         mask = tmask[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]
+         mask = tmask[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt]
          mask = temporary(mask[*])#replicate(1, jpt)
          mask = reform(mask, nxt, nyt, nzt, jpt, /over)
@@ -121 +128 @@
          res = 0.5*(res + shift(res, 0, -1, 0, 0))
          res[*, nyt-1, *, *] = !values.f_nan
-         mask = (vmask())[premierxt:dernierxt, premieryt:dernieryt, premierzt:dernierzt]
+         mask = (vmask())[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt]
          mask = temporary(mask[*])#replicate(1, jpt)
          mask = reform(mask, nxt, nyt, nzt, jpt, /over)
@@ -127 +134 @@
          IF terre[0] NE -1 THEN res[temporary(terre)] = valmask
       END
-   endcase
+    ENDCASE
+
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+   @updateold
+  ENDIF 
 
    return, res

trunk/ToBeReviewed/GRILLE/tracegrille.pro

@@ -11 +11 @@
 ; CALLING SEQUENCE:tracegrille
 ; 
-; INPUTS:none
+; INPUTS:glam et gphi, les tableaux 1d ou 2d des position en
+; longitude/latitude des points de la grille. Si glam et gphi ne sont
+; pas specifies, trace la grille specifiee par vargrid, sur le domaine
+; definit par le dernier domdef.
 ;
 ; KEYWORD PARAMETERS:
@@ -21 +24 @@
 ;        qu''une ligne de j constant tout les ystride points
 ;
-;        OCEAN: pour ne tracer la grille que sur les points oceans
+;        /OCEAN: pour ne tracer la grille que sur les points oceans
 ; 
-;        EARTH: pour ne tracer la grille que sur les points terre
+;        /EARTH: pour ne tracer la grille que sur les points terre
+;
+;        /RMOUT:select to remove all cell having one corner out of the
+;        plot boundaries (!x.range, !y.range)
 ;
 ;        + tous les mots clefs de la procedure PLOTS
@@ -31 +37 @@
 ; COMMON BLOCKS:common.pro
 ;
-; SIDE EFFECTS:trace la grille specifiee par vargrid, sur le domaine
-; definit par le dernier domdef.
+; SIDE EFFECTS:
 ;
 ; RESTRICTIONS:
@@ -38 +43 @@
 ; EXAMPLE:
 ;
-;     IDL> plt,indgen(jpi,jpj),/nocontour,/nocouleur
+;     IDL> plt,indgen(jpi,jpj),/nocontour,/nofill
+;     IDL> vargrid='T'
 ;     IDL> tracegrille,/ocean,color=20
 ;     IDL> tracegrille,/earth,color=80
@@ -49 +55 @@
 ;------------------------------------------------------------
 ;------------------------------------------------------------
-PRO tracegrille, OCEAN = ocean, EARTH = earth, XSTRIDE = xstride, YSTRIDE = ystride, _extra = extra
-@common
-   tempsun = systime(1)         ; pour key_performance
-   if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c'
+PRO tracegrille, glamin, gphiin, OCEAN = ocean, EARTH = earth $
+                 , XSTRIDE = xstride, YSTRIDE = ystride, RMOUT = rmout $
+                 , _extra = extra
+;---------------------------------------------------------
+@cm_4mesh
+@cm_4data
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
+  tempsun = systime(1)          ; pour key_performance
+; to avoid warning message
+  oldexcept = !except
+  !except = 0
+  if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c'
+;
+  if n_elements(glamin) * n_elements(gphiin) EQ 0 then BEGIN
+    grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz
+    IF keyword_set(ocean) AND key_gridtype EQ 'c' THEN BEGIN
+; we reduce the mask to take into account the point located ON the coastline.
+      CASE vargrid OF
+        'U':BEGIN
+          mask = tmask[firstx:lastx, firsty:lasty]
+          IF NOT keyword_set(key_periodic) OR nx NE jpi $
+            THEN tmpx = mask[nx-1, *]
+          mask = (mask+shift(mask, -1, 0)) < 1
+          IF NOT keyword_set(key_periodic) OR nx NE jpi $
+            THEN mask[nx-1, *] = temporary(tmpx)
+        END
+        'V':BEGIN
+          mask = tmask[firstx:lastx, firsty:lasty]
+          tmpy = mask[*, ny-1]
+          mask = (mask+shift(mask, 0, -1)) < 1
+          mask[*, ny-1] = temporary(tmpy)
+        END
+        'F':BEGIN
+          mask = tmask[firstx:lastx, firsty:lasty]
+          IF NOT keyword_set(key_periodic) OR nx NE jpi $
+            THEN tmpx = mask[nx-1, *]
+          tmpy = mask[*, ny-1]
+          mask = (mask+shift(mask, -1, 0)+shift(mask, 0, -1)+shift(mask, -1, -1)) < 1
+          mask[*, ny-1] = temporary(tmpy)
+          IF NOT keyword_set(key_periodic) OR nx NE jpi $
+            THEN mask[nx-1, *] = temporary(tmpx)
+        END
+        ELSE:
+      ENDCASE
+    ENDIF
+  ENDIF ELSE BEGIN
+    glam = glamin
+    gphi = gphiin
+    IF (size(glam))[0] EQ 1 AND (size(gphi))[0] EQ 1 THEN BEGIN
+      nx = n_elements(glam) 
+      ny = n_elements(gphi) 
+      glam = glam#replicate(1, ny)
+      gphi = replicate(1, nx)#gphi
+    ENDIF ELSE BEGIN
+      nx = (size(glam))[1]
+      ny = (size(glam))[2]
+    ENDELSE
+  ENDELSE
+  if n_elements(mask) EQ 0 then mask = replicate(1b, nx, ny)
+  if (size(mask))[0] EQ 3 then mask = mask[*, *, 0]
+;
+  IF keyword_set(RMOUT) THEN BEGIN
+    out = where(glam GT max(!x.range) OR glam LT min(!x.range) $
+                OR gphi GT max(!y.range) OR gphi LT min(!y.range))
+    IF out[0] NE -1 THEN BEGIN
+      glam[out] = !values.f_nan
+      gphi[out] = !values.f_nan
+    ENDIF
+  ENDIF
+;
+  IF keyword_set(ocean) then BEGIN
+    earth = where(mask EQ 0)
+    if earth[0] NE -1 then begin
+      glam[earth] = !values.f_nan
+      gphi[earth] = !values.f_nan
+    ENDIF
+    earth = 0
+  ENDIF
+;
+  IF keyword_set(earth) THEN BEGIN
+    ocean = where(mask EQ 1)
+    if ocean[0] NE -1 then begin
+      glam[ocean] = !values.f_nan
+      gphi[ocean] = !values.f_nan
+    ENDIF
+    ocean = 0
+  ENDIF
+;
+  if NOT keyword_set(xstride) then xstride = 1
+  if NOT keyword_set(ystride) then ystride = 1
+  case key_gridtype of
+    'c':BEGIN
+      for i = 0, ny-1, ystride do begin
+        plots,  glam[*, i], gphi[*, i], _extra = extra
+      endfor
+      for i = 0, nx-1, xstride do begin
+        plots,  glam[i, *], gphi[i, *], _extra = extra
+      endfor
+    END
+    'e':BEGIN
+      shifted = glam[0, 0] LT glam[0, 1]
+      glam2 = glam+(glam[1]-glam[0])/2.
+      if shifted then begin
+        for i = 0, ny-2 do BEGIN
+          xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*]
+          yy = (transpose([[gphi[*, i]], [gphi[*, i+1]]]))[*]
+          plots, xx[0:2*nx-2], yy[0:2*nx-2], _extra = extra
+        ENDFOR
+      ENDIF ELSE BEGIN
+        for i = 1, ny-1 do BEGIN
+          xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*]
+          yy = (transpose([[gphi[*, i]], [gphi[*, i-1]]]))[*]
+          plots, xx[0:2*nx-2], yy[0:2*nx-2], _extra = extra
+        ENDFOR
+      ENDELSE
+      for i = 1, (ny-1)/2 do $
+        plots, [glam[0, 2*i-1], glam[0, 2*i]] $
+        , [gphi[0, 2*i-1], gphi[0, 2*i]], _extra = extra
+      for i = 0, (ny-2)/2 do $
+        plots, [glam[nx-1, 2*i], glam[nx-1, 2*i+1]] $
+        , [gphi[nx-1, 2*i], gphi[nx-1, 2*i+1]], _extra = extra
+    END
+  endcase
 
-   grille, mask, glam, gphi, gdep, nx, ny
-   if (size(mask))[0] EQ 3 then mask = mask[*, *, 0]
-   case 1 of
-      keyword_set(ocean):BEGIN
-         earth = where(mask EQ 0)
-         if earth[0] NE -1 then begin
-            glam[earth] = !values.f_nan
-            gphi[earth] = !values.f_nan
-         endif
-      END
-      keyword_set(earth):BEGIN
-         ocean = where(mask EQ 1)
-         if ocean[0] NE -1 then begin
-            glam[ocean] = !values.f_nan
-            gphi[ocean] = !values.f_nan
-         endif
-      END
-      ELSE:
-   endcase
-   if NOT keyword_set(xstride) then xstride = 1
-   if NOT keyword_set(ystride) then ystride = 1
-   case key_gridtype of
-      'c':BEGIN
-         for i = 0, ny-1, ystride do begin
-            plots,  glam[*, i], gphi[*, i], color = c_cote, _extra = extra
-         endfor
-         for i = 0, nx-1, xstride do begin
-            plots,  glam[i, *], gphi[i, *], color = c_cote, _extra = extra
-         endfor
-      END
-      'e':BEGIN
-         shifted = glam[0, 0] LT glam[0, 1]
-         glam2 = glam+(glam[1]-glam[0])/2.
-         if shifted then begin
-            for i = 0, ny-2 do BEGIN
-               xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*]
-               yy = (transpose([[gphi[*, i]], [gphi[*, i+1]]]))[*]
-               plots, xx[0:2*nx-2], yy[0:2*nx-2], color = c_cote, _extra = extra
-            ENDFOR
-         ENDIF ELSE BEGIN
-            for i = 1, ny-1 do BEGIN
-               xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*]
-               yy = (transpose([[gphi[*, i]], [gphi[*, i-1]]]))[*]
-               plots, xx[0:2*nx-2], yy[0:2*nx-2], color = c_cote, _extra = extra
-            ENDFOR
-         ENDELSE
-         for i = 1, (ny-1)/2 do $
-          plots, [glam[0, 2*i-1], glam[0, 2*i]] $
-          , [gphi[0, 2*i-1], gphi[0, 2*i]], color = c_cote, _extra = extra
-         for i = 0, (ny-2)/2 do $
-          plots, [glam[nx-1, 2*i], glam[nx-1, 2*i+1]] $
-          , [gphi[nx-1, 2*i], gphi[nx-1, 2*i+1]], color = c_cote, _extra = extra
-      END
-   endcase
+  if keyword_set(key_performance) THEN print, 'temps trace grille', systime(1)-tempsun 
+  !except = oldexcept
 
-   if keyword_set(key_performance) THEN print, 'temps trace grille', systime(1)-tempsun 
-   return
+  return
 end

trunk/ToBeReviewed/GRILLE/u2t.pro

@@ -38 +38 @@
 ;------------------------------------------------------------
 FUNCTION u2t, temp
-@common
+;---------------------------------------------------------
+@cm_4mesh
+@cm_4data
+@cm_4cal
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
    res = temp
 ;on force nxt=nxu, etc ...
-   premierxt = premierxu
-   dernierxt = dernierxu
-   premieryt = premieryu
-   dernieryt = dernieryu
+   firstxt = firstxu
+   lastxt = lastxu
+   firstyt = firstyu
+   lastyt = lastyu
    nxt = nxu
    nyt = nyu
    vargrid = 'T'
    if NOT keyword_set(valmask) then valmask = 1e20
-   lon1 = glamt[premierxt, 0]
-   lon2 = glamu[dernierxu, 0]
+   lon1 = glamt[firstxt, 0]
+   lon2 = glamu[lastxu, 0]
 ;
 ; cas sur la taille du tableau et application 
@@ -60 +67 @@
             taille[1] eq nxu and taille[2] eq nyu:
             taille[1] eq jpi and taille[2] eq jpj: $
-             res=res[premierxu:dernierxu, premieryu:dernieryu]
+             res=res[firstxu:lastxu, firstyu:lastyu]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          endcase
-         mask = (umask())[premierxu:dernierxu, premieryu:dernieryu, 0]
+         mask = (umask())[firstxu:lastxu, firstyu:lastyu, 0]
          terre = where(mask EQ 0) 
          IF terre[0] NE -1 THEN res[terre] = !values.f_nan
          res = 0.5*(res + shift(res, 1, 0))
-         if NOT (keyword_set(key_periodique) AND nxu EQ jpi) then res[0, *] = !values.f_nan
-         mask = tmask[premierxu:dernierxu, premieryu:dernieryu, 0]
+         if NOT (keyword_set(key_periodic) AND nxu EQ jpi) then res[0, *] = !values.f_nan
+         mask = tmask[firstxu:lastxu, firstyu:lastyu, 0]
          terre = where(mask EQ 0) 
          IF terre[0] NE -1 THEN res[terre] = valmask
@@ -77 +84 @@
             taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ nzt:
             taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ jpk: $
-             res=res[*, *, premierzt:dernierzt]
+             res=res[*, *, firstzt:lastzt]
             taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ jpt:
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $
-             res=res[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]
+             res=res[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt]
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $
-             res=res[premierxu:dernierxu, premieryu:dernieryu, *]
+             res=res[firstxu:lastxu, firstyu:lastyu, *]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          ENDCASE
          if taille[3] EQ jpt then begin
-            mask = (umask())[premierxu:dernierxu, premieryu:dernieryu, dernierzt*(nzt NE jpk)]
+            mask = (umask())[firstxu:lastxu, firstyu:lastyu, lastzt*(nzt NE jpk)]
             mask = temporary(mask[*])#replicate(1, jpt)
             mask = reform(mask, nxu, nyu, jpt, /over)
-         ENDIF ELSE mask = (umask())[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]
+         ENDIF ELSE mask = (umask())[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt]
          terre = where(temporary(mask) EQ 0) 
          IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan
          res = 0.5*(res + shift(res, 1, 0, 0))
-         if NOT (keyword_set(key_periodique) AND nxu EQ jpi) then res[0, *, *] = !values.f_nan
+         if NOT (keyword_set(key_periodic) AND nxu EQ jpi) then res[0, *, *] = !values.f_nan
          if taille[3] EQ jpt then BEGIN
-            mask = tmask[premierxu:dernierxu, premieryu:dernieryu, dernierzt*(nzt NE jpk)]
+            mask = tmask[firstxu:lastxu, firstyu:lastyu, lastzt*(nzt NE jpk)]
             mask = temporary(mask[*])#replicate(1, jpt)
             mask = reform(mask, nxu, nyu, jpt, /over)
-         ENDIF ELSE mask = tmask[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]
+         ENDIF ELSE mask = tmask[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt]
          terre = where(temporary(mask) EQ 0) 
          IF terre[0] NE -1 THEN res[temporary(terre)] = valmask
@@ -107 +114 @@
             taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ nzt AND taille[4] EQ jpt:
             taille[1] eq nxu and taille[2] eq nyu AND taille[3] EQ jpk AND taille[4] EQ jpt: $
-             res=res[*, *, premierzt:dernierzt, *]
+             res=res[*, *, firstzt:lastzt, *]
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $
-             res=res[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt, *]
+             res=res[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt, *]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          ENDCASE
-         mask = (umask())[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]
+         mask = (umask())[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt]
          mask = temporary(mask[*])#replicate(1, jpt)
          mask = reform(mask, nxu, nyu, nzt, jpt, /over)
@@ -119 +126 @@
          IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan
          res = 0.5*(res + shift(res, 1, 0, 0, 0))
-         if NOT (keyword_set(key_periodique) AND nxu EQ jpi) then res[0, *, *, *] = !values.f_nan
-         mask = tmask[premierxu:dernierxu, premieryu:dernieryu, premierzt:dernierzt]
+         if NOT (keyword_set(key_periodic) AND nxu EQ jpi) then res[0, *, *, *] = !values.f_nan
+         mask = tmask[firstxu:lastxu, firstyu:lastyu, firstzt:lastzt]
          mask = temporary(mask[*])#replicate(1, jpt)
          mask = reform(mask, nxu, nyu, nzt, jpt, /over)
@@ -128 +135 @@
    endcase
 
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+   @updateold
+  ENDIF 
+  
    return, res
 END

trunk/ToBeReviewed/GRILLE/umask.pro

@@ -39 +39 @@
 ;------------------------------------------------------------
 FUNCTION umask
-   tempsun = systime(1)         ; pour key_performance
-@common
-   if jpk EQ 1 then begin
-      res = tmask*shift(tmask, -1, 0)
-      if NOT keyword_set(key_periodique) then res[jpi-1, *] = umaskred
-   ENDIF ELSE BEGIN 
-      res = tmask*shift(tmask, -1, 0, 0)
-      if NOT keyword_set(key_periodique) then res[jpi-1, *, *] = umaskred
-   ENDELSE
+;---------------------------------------------------------
+@cm_4mesh
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
+  tempsun = systime(1)          ; pour key_performance
 ;
-   if keyword_set(key_performance) THEN print, 'temps umask', systime(1)-tempsun 
+  CASE size(tmask, /n_dimensions) OF
+    2:res = tmask*shift(tmask, -1, 0)
+    3:res = tmask*shift(tmask, -1, 0, 0)
+  ENDCASE
+;   
+  if NOT keyword_set(key_periodic) then res[jpi-1, *, *] = umaskred
+  if keyword_set(key_performance) THEN print, 'temps umask', systime(1)-tempsun 
 ;
-   return, res
+  return, res
 end

trunk/ToBeReviewed/GRILLE/v2t.pro

@@ -38 +38 @@
 ;------------------------------------------------------------
 FUNCTION v2t, temp
-@common
+;---------------------------------------------------------
+@cm_4mesh
+@cm_4data
+@cm_4cal
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
    res = temp
 ;on force nxt=nxv, etc ...
-   premierxt = premierxv
-   dernierxt = dernierxv
-   premieryt = premieryv
-   dernieryt = dernieryv
+   firstxt = firstxv
+   lastxt = lastxv
+   firstyt = firstyv
+   lastyt = lastyv
    nxt = nxv
    nyt = nyv
    vargrid = 'T'
    if NOT keyword_set(valmask) then valmask = 1e20
-   lat1 = gphit[0, premieryt]
-   lat2 = gphiv[0, dernieryv]
+   lat1 = gphit[0, firstyt]
+   lat2 = gphiv[0, lastyv]
 
 ; cas sur la taille du tableau et application 
@@ -60 +67 @@
             taille[1] eq nxv and taille[2] eq nyv:
             taille[1] eq jpi and taille[2] eq jpj: $
-             res=res[premierxv:dernierxv, premieryv:dernieryv]
+             res=res[firstxv:lastxv, firstyv:lastyv]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          endcase
-         mask = (vmask())[premierxv:dernierxv, premieryv:dernieryv, 0]
+         mask = (vmask())[firstxv:lastxv, firstyv:lastyv, 0]
          terre = where(mask EQ 0) 
          IF terre[0] NE -1 THEN res[terre] = !values.f_nan
          res = 0.5*(res + shift(res, 0, +1))
          res[*, 0] = !values.f_nan
-         mask = tmask[premierxv:dernierxv, premieryv:dernieryv, 0]
+         mask = tmask[firstxv:lastxv, firstyv:lastyv, 0]
          terre = where(mask EQ 0) 
          IF terre[0] NE -1 THEN res[terre] = valmask
@@ -77 +84 @@
             taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ nzt:
             taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ jpk: $
-             res=res[*, *, premierzt:dernierzt]
+             res=res[*, *, firstzt:lastzt]
             taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ jpt:
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk: $
-             res=res[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]
+             res=res[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt]
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpt: $
-             res=res[premierxv:dernierxv, premieryv:dernieryv, *]
+             res=res[firstxv:lastxv, firstyv:lastyv, *]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          ENDCASE
          if taille[3] EQ jpt then begin
-            mask = (vmask())[premierxv:dernierxv, premieryv:dernieryv, dernierzt*(nzt NE jpk)]
+            mask = (vmask())[firstxv:lastxv, firstyv:lastyv, lastzt*(nzt NE jpk)]
             mask = temporary(mask[*])#replicate(1, jpt)
             mask = reform(mask, nxv, nyv, jpt, /over)
-         ENDIF ELSE mask = (vmask())[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]
+         ENDIF ELSE mask = (vmask())[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt]
          terre = where(temporary(mask) EQ 0) 
          IF terre[0] NE -1 THEN res[temporary(terre)] = !values.f_nan
@@ -96 +103 @@
          res[*, 0, *] = !values.f_nan
          if taille[3] EQ jpt then BEGIN
-            mask = tmask[premierxv:dernierxv, premieryv:dernieryv, dernierzt*(nzt NE jpk)]
+            mask = tmask[firstxv:lastxv, firstyv:lastyv, lastzt*(nzt NE jpk)]
             mask = temporary(mask[*])#replicate(1, jpt)
             mask = reform(mask, nxv, nyv, jpt, /over)
-         ENDIF ELSE mask = tmask[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]
+         ENDIF ELSE mask = tmask[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt]
          terre = where(temporary(mask) EQ 0) 
          IF terre[0] NE -1 THEN res[temporary(terre)] = valmask
@@ -107 +114 @@
             taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ nzt AND taille[4] EQ jpt:
             taille[1] eq nxv and taille[2] eq nyv AND taille[3] EQ jpk AND taille[4] EQ jpt: $
-             res=res[*, *, premierzt:dernierzt, *]
+             res=res[*, *, firstzt:lastzt, *]
             taille[1] eq jpi and taille[2] eq jpj AND taille[3] EQ jpk AND taille[4] EQ jpt: $
-             res=res[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt, *]
+             res=res[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt, *]
             else: $
              return, report('Probleme d''adequation entre les tailles du domaine et de la boite.')
          ENDCASE
-         mask = (vmask())[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]
+         mask = (vmask())[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt]
          mask = temporary(mask[*])#replicate(1, jpt)
          mask = reform(mask, nxv, nyv, nzt, jpt, /over)
@@ -120 +127 @@
          res = 0.5*(res + shift(res, 0, +1, 0, 0))
          res[*, 0, *, *] = !values.f_nan
-         mask = tmask[premierxv:dernierxv, premieryv:dernieryv, premierzt:dernierzt]
+         mask = tmask[firstxv:lastxv, firstyv:lastyv, firstzt:lastzt]
          mask = temporary(mask[*])#replicate(1, jpt)
          mask = reform(mask, nxv, nyv, nzt, jpt, /over)
@@ -128 +135 @@
    endcase
 
-   return, res
+   IF NOT keyword_set(key_forgetold) THEN BEGIN
+   @updateold
+   ENDIF 
+
+    return, res
 END
 

trunk/ToBeReviewed/GRILLE/vmask.pro

@@ -27 +27 @@
 FUNCTION vmask
 @common
-   tempsun = systime(1)         ; pour key_performance
-   if jpk EQ 1 then begin
-      res = tmask*shift(tmask, 0, -1)
-      res[*, jpj-1] = vmaskred
-   ENDIF ELSE BEGIN 
-      res = tmask*shift(tmask, 0, -1, 0)
-      res[*, jpj-1, *] = vmaskred
-   ENDELSE 
-   if keyword_set(key_performance) THEN print, 'temps vmask', systime(1)-tempsun 
+  tempsun = systime(1)          ; pour key_performance
 ;
-   return, res
+  CASE size(tmask, /n_dimensions) OF
+    2:res = tmask*shift(tmask, 0, -1)
+    3:res = tmask*shift(tmask, 0, -1, 0)
+  ENDCASE
+;
+  res[*, jpj-1, *] = vmaskred
+  if keyword_set(key_performance) THEN print, 'temps vmask', systime(1)-tempsun 
+;
+  return, res
 end

trunk/ToBeReviewed/GRILLE/whichgrid.pro

@@ -4 +4 @@
    case filetype of
       'OPA C-Grid.Binary IEEE: Meshmask':BEGIN
-         meshlec, filename, /pasblabla, _extra = ex
+         meshread, filename, /pasblabla, _extra = ex
       END
       'OPA C-Grid.Net Cdf: Meshmask':BEGIN
-         ncdf_meshlec, filename, _extra = ex
+         ncdf_meshread, filename, _extra = ex
       END
       'Regular 2D.Binary IEEE: mask':
@@ -28 +28 @@
 ;*********************************************************************
 FUNCTION getgridparameter, top
-@common
-
+;---------------------------------------------------------
+@cm_4mesh
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
    widget_control, widget_info(top, find_by_uname = 'xmesh'), get_value = answer
    jpiglo = long(answer[0])
@@ -48 +52 @@
     , get_value = answer
    key_shift = long(answer[0])
-   widget_control, widget_info(top, find_by_uname = 'key_periodique') $
+   widget_control, widget_info(top, find_by_uname = 'key_periodic') $
     , get_value = answer
-   key_periodique = long(answer[0])
+   key_periodic = long(answer[0])
    widget_control, widget_info(top, find_by_uname = 'triangulation') $
     , get_value = answer
@@ -61 +65 @@
           , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $
           , izminmesh:izminmesh, izmaxmesh:izmaxmesh $
-          , key_shift:key_shift, key_periodique:key_periodique $
+          , key_shift:key_shift, key_periodic:key_periodic $
           , triangulation:triangulation, boundary:boundary}
+
+   @updateold
 
    return, res
@@ -68 +74 @@
 ;*********************************************************************
 pro showgridparameter, basetop, EDITABLE = editable, _EXTRA = ex
-;
-@common
-;------------------------------------------------------------
-;
+;---------------------------------------------------------
+@cm_4mesh
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
    widget_control, basetop, update = 0
    base=widget_base(basetop, /COLUMN, /align_center, _EXTRA = ex)
@@ -82 +90 @@
    if NOT keyword_set(key_shift) then key_shift = 0
    nothing = widget_text(basea, value = strtrim(key_shift*(1-keyword_set(editable)),1), uname = 'key_shift', xsize = 4, EDITABLE = editable)
-   nothing = widget_label(basea, value = 'key_periodique')
-   if NOT keyword_set(key_periodique) then key_periodique = 0
-   nothing = widget_text(basea, value = strtrim(key_periodique*(1-keyword_set(editable)),1), uname = 'key_periodique', xsize = 4, EDITABLE = editable)
+   nothing = widget_label(basea, value = 'key_periodic')
+   if NOT keyword_set(key_periodic) then key_periodic = 0
+   nothing = widget_text(basea, value = strtrim(key_periodic*(1-keyword_set(editable)),1), uname = 'key_periodic', xsize = 4, EDITABLE = editable)
    baseb=widget_base(base, /row, /align_center)
    nothing = widget_label(baseb, value = 'use a triangulation (y/n) ?')
@@ -159 +167 @@
 ;*********************************************************************
 PRO whichgrid_event, event
-@common
+;---------------------------------------------------------
+@cm_4mesh
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
    widget_control, event.id,  get_uvalue = eventuvalue
    IF chkstru(eventuvalue,'name') EQ 0 THEN return
@@ -181 +194 @@
          if rstrpos(filename, '.pro') EQ strlen(filename)-4 then begin
             createpro, '@'+strmid(filename, 0, strlen(filename)-4) $
-          , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $
-          +'for_createpro.pro'
+          , filename = myuniquetmpdir +'for_createpro.pro'
             showgridparameter, event.handler, group_leader = event.handler,/frame, uname = 'showgridparameter'
          ENDIF ELSE BEGIN
@@ -204 +216 @@
          if rstrpos(filename, '.pro') EQ strlen(filename)-4 then begin
             createpro, '@'+strmid(filename, 0, strlen(filename)-4) $
-          , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $
-          +'for_createpro.pro'
+          , filename = myuniquetmpdir +'for_createpro.pro'
             showgridparameter, event.handler, group_leader = event.handler,/frame, uname = 'showgridparameter'
          ENDIF ELSE BEGIN
@@ -223 +234 @@
          case filetype of
             'OPA C-Grid.Binary IEEE: Meshmask':BEGIN
-               meshlec, filename, /pasblabla, /getdimensions
+               meshread, filename, /pasblabla, /getdimensions
                showgridparameter, event.handler, /editable, group_leader = event.handler,/frame, uname = 'showgridparameter'
             END
             'OPA C-Grid.Net Cdf: Meshmask':BEGIN
-               ncdf_meshlec, filename, /getdimensions
+               ncdf_meshread, filename, /getdimensions
                showgridparameter, event.handler, /editable, group_leader = event.handler,/frame, uname = 'showgridparameter'
             END
@@ -269 +280 @@
 END
 ;*********************************************************************
-FUNCTION whichgrid, name, IODIRECTORY = iodirectory, PARENT = parent, _EXTRA = ex
-@common
-;
+FUNCTION whichgrid, name, PARENT = parent, _EXTRA = ex
+;---------------------------------------------------------
+@cm_4mesh
+  IF NOT keyword_set(key_forgetold) THEN BEGIN
+@updatenew
+  ENDIF
+;---------------------------------------------------------
 ;
    if n_elements(name) NE 0 then begin
-      filename = isafile(filename = name)
+      filename = isafile(filename = name, IODIRECTORY = iodir, _extra = ex)
       if size(filename, /type) NE 7 then return, -1
    ENDIF ELSE filaname = 'no file'
@@ -301 +316 @@
       if rstrpos(filename, '.pro') EQ strlen(filename)-4 then begin
          createpro, '@'+strmid(filename, 0, strlen(filename)-4) $
-          , filename = isadirectory(io = homedir, title = 'Bad definition of Homedir') $
-          +'for_createpro.pro'
+          , filename = myuniquetmpdir +'for_createpro.pro'
          showgridparameter, base, group_leader = base,/frame, uname = 'showgridparameter'
       ENDIF ELSE BEGIN