Index: trunk/SRC/Grid/computegrid.pro
===================================================================
--- trunk/SRC/Grid/computegrid.pro	(revision 500)
+++ trunk/SRC/Grid/computegrid.pro	(revision 501)
@@ -264,16 +264,16 @@
 ;-
 PRO computegrid, startx, starty, stepxin, stepyin, nxin, nyin $
-                 , XAXIS=xaxis, YAXIS=yaxis, ZAXIS=zaxis, ZZAXIS=zzaxis $
-                 , MASK=mask, GLAMBOUNDARY=glamboundary $
-                 , XMINMESH=xminmesh, XMAXMESH=xmaxmesh $
-                 , YMINMESH=yminmesh, YMAXMESH=ymaxmesh $
-                 , ZMINMESH=zminmesh, ZMAXMESH=zmaxmesh $
-                 , ONEARTH=onearth, PERIODIC=periodic $
-                 , PLAIN=plain, SHIFT=shift, STRIDE=stride $
-                 , YREVERSE=yreverse, ZREVERSE=zreverse  $
-                 , FULLCGRID=fullcgrid, XYINDEX=xyindex $
-                 , UBASE2TBASE=ubase2tbase, VBASE2TBASE=vbase2tbase $
-                 , FBASE2TBASE=fbase2tbase, FFBASE2TBASE=ffbase2tbase $
-                 , STRCALLING=strcalling, ROMSH=romsh, _EXTRA=ex
+                 , XAXIS = xaxis, YAXIS = yaxis, ZAXIS = zaxis, ZZAXIS = zzaxis $
+                 , MASK = mask, GLAMBOUNDARY = glamboundary $
+                 , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
+                 , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
+                 , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
+                 , ONEARTH = onearth, PERIODIC = periodic $
+                 , PLAIN = plain, SHIFT = shift, STRIDE = stride $
+                 , YREVERSE = yreverse, ZREVERSE = zreverse  $
+                 , FULLCGRID = fullcgrid, XYINDEX = xyindex $
+                 , UBASE2TBASE = ubase2tbase, VBASE2TBASE = vbase2tbase $
+                 , FBASE2TBASE = fbase2tbase, FFBASE2TBASE = ffbase2tbase $
+                 , STRCALLING = strcalling, ROMSH = romsh, _EXTRA = ex
 ;
   compile_opt idl2, strictarrsubs
@@ -432,5 +432,5 @@
     nzz = jpk
   ENDIF ELSE BEGIN
-   jpkglo = long(nz)
+    jpkglo = long(nz)
     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
@@ -570,9 +570,9 @@
 ; check the periodicity if iyminmesh and iymaxmesh have the default definitions...
   IF NOT keyword_set(plain) AND key_onearth EQ 1 $
-    AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 AND jpj GE 3 AND jpi GE 2 THEN BEGIN
+     AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 AND jpj GE 3 AND jpi GE 2 THEN BEGIN
 
     CASE 1 OF
       ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
-        AND array_equal(gphit[1:*, jpj-1], reverse(gphit[1:*, jpj-3])) EQ 1:BEGIN
+         AND array_equal(gphit[1:*, jpj-1], reverse(gphit[1:*, jpj-3])) EQ 1:BEGIN
 ; T pivot
         ymaxmesh = -1
@@ -586,5 +586,5 @@
       END
       ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
-       AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN
+         AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN
 ; F pivot
         ymaxmesh = -1
@@ -604,5 +604,5 @@
 ; check x periodicity...
 ;====================================================
-IF n_elements(periodic) NE 0 THEN forcenoperio = 1 - keyword_set(periodic)
+  IF n_elements(periodic) NE 0 THEN forcenoperio = 1 - keyword_set(periodic)
 ; check the periodicity if ixminmesh and ixmaxmesh have the default definitions...
   IF NOT keyword_set(plain) AND NOT keyword_set(forcenoperio) AND key_onearth EQ 1 $
@@ -610,5 +610,5 @@
     CASE 0 OF
       total((glamt[0, *] - glamt[jpi-2, *]) MOD 360) $
-      + total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN
+         + total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN
         xminmesh = 1
         xmaxmesh = -1
@@ -685,6 +685,6 @@
       IF array_equal(sort(xtest), lindgen(jpi)) EQ 0 THEN BEGIN
         ras = report(['WARNING: we cannot sort the xaxis with a simple shift...', $
-        'we force key_periodic = 0 and key_shift = 0', $
-        'only horizontal plot may be ok...'])
+                      'we force key_periodic = 0 and key_shift = 0', $
+                      'only horizontal plot may be ok...'])
         key_periodic = 0
         xnotsorted = 1
@@ -751,5 +751,5 @@
       ELSE:BEGIN
         if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $
-          OR (keyword_set(key_periodic) AND key_irregular) then BEGIN
+           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)
@@ -758,9 +758,9 @@
           stepxf = min(abs(stepxf), dimension = 3)
           IF NOT keyword_set(key_periodic) THEN $
-            stepxf[jpi-1, *] = stepxf[jpi-2, *]
+             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, *] $
+             stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $
           ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *]
         ENDELSE
@@ -797,5 +797,5 @@
         IF jpi GT 1 THEN BEGIN
           if NOT keyword_set(key_periodic) THEN $
-            stepyf[jpi-1, *] = stepyf[jpi-2, *]
+             stepyf[jpi-1, *] = stepyf[jpi-2, *]
           stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2]
         ENDIF
@@ -831,5 +831,5 @@
       e1t = 0.5*(stepxu+shift(stepxu, 1, 0))
       IF NOT keyword_set(key_periodic) THEN $
-        e1t[0, *] = e1t[1, *]
+         e1t[0, *] = e1t[1, *]
     ENDIF ELSE e1t = replicate(stepx, jpi, jpj)
   ENDIF ELSE e1t = replicate(1b, jpi, jpj)
@@ -857,66 +857,66 @@
   IF jpj EQ 1 THEN e2t = reform(e2t, jpi, jpj, /over)
 ;
+;====================================================
+; 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
+;
+  IF jpj EQ 1 THEN glamu = reform(glamu, jpi, jpj, /over)
+;
+;====================================================
+; 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
+    ENDIF ELSE gphiu = gphit
+  ENDIF ELSE gphiu = gphit
+  IF key_onearth THEN gphiu = -90. > gphiu < 90.
+;
+  IF jpj EQ 1 THEN gphiu = reform(gphiu, jpi, jpj, /over)
+;
+;====================================================
+; 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.
+;
+  IF jpj EQ 1 THEN gphiv = reform(gphiv, jpi, jpj, /over)
+
 ;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  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
-;
-    IF jpj EQ 1 THEN glamu = reform(glamu, jpi, jpj, /over)
-;
-;====================================================
-; 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
-      ENDIF ELSE gphiu = gphit
-    ENDIF ELSE gphiu = gphit
-  IF key_onearth THEN gphiu = -90. > gphiu < 90.
-;
-  IF jpj EQ 1 THEN gphiu = reform(gphiu, jpi, jpj, /over)
-;
-;====================================================
-; 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.
-;
-    IF jpj EQ 1 THEN gphiv = reform(gphiv, jpi, jpj, /over)
+  IF keyword_set(fullcgrid) THEN BEGIN    
 ;
 ;====================================================
@@ -925,5 +925,5 @@
 ;
     IF jpi GT 1 AND n_elements(stepx) NE 1 THEN $
-      e1u = stepxu ELSE e1u = e1t
+       e1u = stepxu ELSE e1u = e1t
 ;
 ;====================================================
@@ -1004,10 +1004,44 @@
 ;
   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
+    IF keyword_set(key_irregular) THEN BEGIN
+; T
+      e1t = map_npoints(glamu, gphiu, shift(glamu, 1, 0), shift(gphiu, 1, 0), /two_by_two)
+      e1t = reform(e1t, jpi, jpj, /over)
+      IF NOT keyword_set(key_periodic) THEN e1t[0, *] = e1t[1, *]
+      e2t = map_npoints(glamv, gphiv, shift(glamv, 0, 1), shift(gphiv, 0, 1), /two_by_two)
+      e2t = reform(e2t, jpi, jpj, /over)
+      e2t[*, 0] = e2t[*, 1]
+      IF keyword_set(fullcgrid) THEN BEGIN
+; U
+        e1u = map_npoints(glamt, gphit, shift(glamt, -1, 0), shift(gphit, -1, 0), /two_by_two)
+        e1u = reform(e1u, jpi, jpj, /over)
+        IF NOT keyword_set(key_periodic) THEN e1u[jpi-1, *] = e1u[jpi-2, *]
+        e2u = map_npoints(glamf, gphif, shift(glamf, 0, 1), shift(gphif, 0, 1), /two_by_two)
+        e2u = reform(e2u, jpi, jpj, /over)
+        e2u[*, 0] = e2u[*, 1]
+; V
+        e1v = map_npoints(glamf, gphif, shift(glamf, 1, 0), shift(gphif, 1, 0), /two_by_two)
+        e1v = reform(e1v, jpi, jpj, /over)
+        IF NOT keyword_set(key_periodic) THEN e1v[0, *] = e1v[1, *]
+        e2v = map_npoints(glamt, gphit, shift(glamt, 0, -1), shift(gphit, 0, -1), /two_by_two)
+        e2v = reform(e2v, jpi, jpj, /over)
+        e2v[*, jpj-1] = e2v[*, jpj-2]
+; F
+        e1f = map_npoints(glamv, gphiv, shift(glamv, -1, 0), shift(gphiv, -1, 0), /two_by_two)
+        e1f = reform(e1f, jpi, jpj, /over)
+        IF NOT keyword_set(key_periodic) THEN e1f[jpi-1, *] = e1f[jpi-2, *]
+        e2f = map_npoints(glamu, gphiu, shift(glamu, 0, -1), shift(gphiu, 0, -1), /two_by_two)
+        e2f = reform(e2f, jpi, jpj, /over)
+        e2f[*, jpj-1] = e2f[*, jpj-2]
+      ENDIF
+
+    ENDIF ELSE 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
+    ENDELSE
   ENDIF
 ;
@@ -1042,33 +1076,33 @@
 ; 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[zmin:zmax]
-        IF n_elements(zreverse) EQ 0 THEN BEGIN
-          IF jpk GT 1 THEN BEGIN
-            if gdept[0] GT gdept[1] then key_zreverse = 1 ELSE key_zreverse = 0
-          ENDIF ELSE key_zreverse = 0
-        ENDIF ELSE key_zreverse = zreverse
-        IF keyword_set(key_zreverse) THEN gdept = reverse(gdept)
-      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
-    IF keyword_set(romsh) THEN gdepw = gdept
+  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[zmin:zmax]
+      IF n_elements(zreverse) EQ 0 THEN BEGIN
+        IF jpk GT 1 THEN BEGIN
+          if gdept[0] GT gdept[1] then key_zreverse = 1 ELSE key_zreverse = 0
+        ENDIF ELSE key_zreverse = 0
+      ENDIF ELSE key_zreverse = zreverse
+      IF keyword_set(key_zreverse) THEN gdept = reverse(gdept)
+    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
+  IF keyword_set(romsh) THEN gdepw = gdept
 ;
 ;====================================================
@@ -1076,9 +1110,9 @@
 ;====================================================
 ;
-    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
+  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
 ;
 ;====================================================
@@ -1228,5 +1262,5 @@
     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
+         OR (keyword_set(key_periodic) AND key_irregular) then BEGIN
         stepxf = (glamt + 720) MOD 360
         stepxf = shift(stepxf, -1, -1) - stepxf
@@ -1234,10 +1268,10 @@
         stepxf = min(abs(stepxf), dimension = 3)
         IF NOT keyword_set(key_periodic) THEN $
-          stepxf[jpi-1, *] = stepxf[jpi-2, *]
+           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, *]
+           stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $
+        ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *]
       ENDELSE
       IF jpj GT 1 THEN BEGIN
@@ -1259,5 +1293,5 @@
       IF jpi GT 1 THEN BEGIN
         if NOT keyword_set(key_periodic) THEN $
-          stepyf[jpi-1, *] = stepyf[jpi-2, *]
+           stepyf[jpi-1, *] = stepyf[jpi-2, *]
         stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2]
       ENDIF
@@ -1316,5 +1350,5 @@
 ;
   IF total(tmask) EQ jpi*jpj*jpk $
-    AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $
+     AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $
   ELSE BEGIN
 ; are we using ORCA2 ?
@@ -1357,24 +1391,24 @@
 
   ccmeshparameters = {filename:strcalling  $
-          , glaminfo:float(string(glaminfo, format = '(E11.4)')) $
-          , gphiinfo:float(string(gphiinfo, format = '(E11.4)')) $
-          , jpiglo:jpiglo, jpjglo:jpjglo, jpkglo:jpkglo $
-          , jpi:jpi, jpj:jpj, jpk:jpk $
-          , ixminmesh:ixminmesh, ixmaxmesh:ixmaxmesh $
-          , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $
-          , izminmesh:izminmesh, izmaxmesh:izmaxmesh $
-          , key_shift:key_shift, key_periodic:key_periodic $
-          , key_stride:key_stride, key_gridtype:key_gridtype $
-          , key_yreverse:key_yreverse, key_zreverse:key_zreverse $
-          , key_partialstep:key_partialstep, key_onearth:key_onearth}
+                      , glaminfo:float(string(glaminfo, format = '(E11.4)')) $
+                      , gphiinfo:float(string(gphiinfo, format = '(E11.4)')) $
+                      , jpiglo:jpiglo, jpjglo:jpjglo, jpkglo:jpkglo $
+                      , jpi:jpi, jpj:jpj, jpk:jpk $
+                      , ixminmesh:ixminmesh, ixmaxmesh:ixmaxmesh $
+                      , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $
+                      , izminmesh:izminmesh, izmaxmesh:izmaxmesh $
+                      , key_shift:key_shift, key_periodic:key_periodic $
+                      , key_stride:key_stride, key_gridtype:key_gridtype $
+                      , key_yreverse:key_yreverse, key_zreverse:key_zreverse $
+                      , key_partialstep:key_partialstep, key_onearth:key_onearth}
 
   ccreadparameters = {funclec_name:'read_ncdf' $
-          , jpidta:jpidta, jpjdta:jpjdta, jpkdta:jpkdta $
-          , ixmindta:ixmindta, ixmaxdta:ixmaxdta $
-          , iymindta:iymindta, iymaxdta:iymaxdta $
-          , izmindta:izmindta, izmaxdta:izmaxdta}
+                      , jpidta:jpidta, jpjdta:jpjdta, jpkdta:jpkdta $
+                      , ixmindta:ixmindta, ixmaxdta:ixmaxdta $
+                      , iymindta:iymindta, iymaxdta:iymaxdta $
+                      , izmindta:izmindta, izmaxdta:izmaxdta}
 ;------------------------------------------------------------
   IF keyword_set(key_performance) EQ 1 THEN $
-    print, 'time computegrid', systime(1)-time1
+     print, 'time computegrid', systime(1)-time1
 ;------------------------------------------------------------
   return
Index: trunk/SRC/Grid/ncdf_meshroms.pro
===================================================================
--- trunk/SRC/Grid/ncdf_meshroms.pro	(revision 500)
+++ trunk/SRC/Grid/ncdf_meshroms.pro	(revision 501)
@@ -273,6 +273,6 @@
     nothing = execute(command)
   ENDFOR
-  d1t = 1.e3*(shift(d1t, -1, 0) - d1t)
-  d2t = 1.e3*(shift(d2t, 0, -1) - d2t)
+  d1t = shift(d1t, -1, 0) - d1t
+  d2t = shift(d2t, 0, -1) - d2t
   for i = 0, n_elements(namebase2)-1 do begin
     command = namebase2[i]+'t = '+namebase2[i]+'t[0:jpi-2, 0:jpj-2]'
@@ -306,7 +306,7 @@
     nothing = execute(command)
   ENDFOR
-  tmpsave = 2. * 1.e3 * d1u[0, 0:jpj-2]
-  d1u = 1.e3*(shift(d1u, -1, 0) - d1u)
-  d2u = 1.e3*(shift(d2u, 0, -1) - d2u)
+  tmpsave = 2. * d1u[0, 0:jpj-2]
+  d1u = shift(d1u, -1, 0) - d1u
+  d2u = shift(d2u, 0, -1) - d2u
   for i = 0, n_elements(namebase2)-1 do begin
     command = namebase2[i]+'u = '+namebase2[i]+'u[*, 0:jpj-2]'
@@ -342,7 +342,7 @@
     nothing = execute(command)
   ENDFOR
-  d1v = 1.e3*(shift(d1v, -1, 0) - d1v)
-  tmpsave = 2. * 1.e3 * d2v[0:jpi-2, 0]
-  d2v = 1.e3*(shift(d2v, 0, -1) - d2v)
+  d1v = shift(d1v, -1, 0) - d1v
+  tmpsave = 2. * d2v[0:jpi-2, 0]
+  d2v = shift(d2v, 0, -1) - d2v
   for i = 0, n_elements(namebase2)-1 do begin
     command = namebase2[i]+'v = '+namebase2[i]+'v[0:jpi-2, *]'
@@ -372,8 +372,8 @@
     nothing = execute(command)
   ENDFOR
-  tmpsave1 = 2. * 1.e3 * d1f[0, *]
-  d1f = 1.e3*(shift(d1f, -1, 0) - d1f)
-  tmpsave2 = 2. * 1.e3 * d2f[*, 0]
-  d2f = 1.e3*(shift(d2f, 0, -1) - d2f)
+  tmpsave1 = 2. * d1f[0, *]
+  d1f = shift(d1f, -1, 0) - d1f
+  tmpsave2 = 2. * d2f[*, 0]
+  d2f = shift(d2f, 0, -1) - d2f
   fmaskredy = byte(maskf[jpi-1, *])
   IF jpk GT 1 THEN fmaskredy = reform(fmaskredy[*]#replicate(1b, jpk), 1, jpj, jpk, /overwrite)
Index: trunk/SRC/Grid/romsdepth.pro
===================================================================
--- trunk/SRC/Grid/romsdepth.pro	(revision 500)
+++ trunk/SRC/Grid/romsdepth.pro	(revision 501)
@@ -25,5 +25,5 @@
 ;
 ;-
-FUNCTION romsdepth
+FUNCTION romsdepth, new_s_coord = new_s_coord
 ;
   compile_opt idl2, strictarrsubs
@@ -46,8 +46,5 @@
   grille, -1, -1, -1, -1, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz
   hroms = hroms[firstx:lastx, firsty:lasty]
-  nt = n_elements(zeta)/nx/ny
-;
-  cff1 = 1./sinh(theta_s)
-  cff2 = 0.5/tanh(0.5*theta_s)
+  nt = n_elements(zeta[firstx:lastx, firsty:lasty, *])/nx/ny
 ;
   IF type EQ 'W' THEN BEGIN
@@ -59,16 +56,42 @@
   ENDELSE
 ;
-  cs = (1.-theta_b)*cff1*sinh(theta_s*sc)+theta_b*(cff2*tanh(theta_s*(sc+0.5))-0.5)
-  cff = hc*(sc-cs)
-  cff1 = cs
+  IF keyword_set(new_s_coord) THEN BEGIN
+
+    IF theta_s GT 0. THEN BEGIN
+      csrf = ( 1. - cosh(theta_s*sc) ) / ( cosh(theta_s) - 1. )
+    ENDIF ELSE BEGIN
+      csrf = - sc*sc
+    ENDELSE 
+    if theta_b GT 0. THEN BEGIN 
+      cs = ( exp(theta_b*csrf) - 1. ) / ( 1. - exp(-theta_b) )  
+    ENDIF ELSE BEGIN
+      cs = csrf
+    ENDELSE
+    
+    hh = (hroms[*]#replicate(1., jpk))
+    z0 = hh * ( replicate(1., nx*ny)#(hc*sc)[*] + hh * replicate(1., nx*ny)#cs[*] ) / ( hc + hh )
+    
+    
+  ENDIF ELSE BEGIN
+    
+    cff1 = 1./sinh(theta_s)
+    cff2 = 0.5/tanh(0.5*theta_s)
+
+    cs = (1.-theta_b)*cff1*sinh(theta_s*sc)+theta_b*(cff2*tanh(theta_s*(sc+0.5))-0.5)
+    cff = hc*(sc-cs)
+    cff1 = cs
+
+    z0 = replicate(1., nx*ny)#cff[*] + (replicate(1., nx*ny)#cff1[*]) * (hroms[*]#replicate(1., jpk))
+
+  ENDELSE
+
 ;
   hinv = 1./hroms
   hinv = hinv[*]#replicate(1., jpk)
 ; put a z dimension to zeta
-  zeta = transpose(temporary(zeta))
+  zeta = transpose((temporary(zeta))[firstx:lastx, firsty:lasty, *])
   zeta = reform((temporary(zeta))[*]#replicate(1., jpk), nt, ny, nx, jpk, /overwrite)
   zeta = transpose(temporary(zeta), [2, 1, 3, 0])
 
-  z0 = replicate(1., nx*ny)#cff[*] + (replicate(1., nx*ny)#cff1[*]) * (hroms[*]#replicate(1., jpk))
   z = (z0[*]#replicate(1., nt)) + temporary(zeta) * ((1.+z0*hinv)[*]#replicate(1., nt))
   z = reform(z, nx, ny, jpk, nt, /overwrite)
Index: trunk/SRC/Interpolation/compute_fromreg_imoms3_weigaddr.pro
===================================================================
--- trunk/SRC/Interpolation/compute_fromreg_imoms3_weigaddr.pro	(revision 500)
+++ trunk/SRC/Interpolation/compute_fromreg_imoms3_weigaddr.pro	(revision 501)
@@ -8,8 +8,8 @@
 ; Interpolation
 ;
-; @param alonin {in}{required}{type=2d array}
+; @param alonin {in}{required}{type=1d array}
 ; longitude of the input data
 ;
-; @param alatin {in}{required}{type=2d array}
+; @param alatin {in}{required}{type=1d array}
 ; latitude of the input data
 ;
@@ -133,5 +133,5 @@
   IF total(alon[indexlon] GT olon) NE 0 THEN stop
   IF total(alon[indexlon + 1L] LE olon) NE 0 THEN stop
-  IF (where(indexlon LE 1L     ))[0] NE -1 THEN stop
+  IF (where(indexlon LT 1L     ))[0] NE -1 THEN stop
   IF (where(indexlon GE jpia-3L))[0] NE -1 THEN stop
   indexlat = value_locate(alat, olat)
Index: trunk/SRC/Interpolation/neighbor.pro
===================================================================
--- trunk/SRC/Interpolation/neighbor.pro	(revision 500)
+++ trunk/SRC/Interpolation/neighbor.pro	(revision 501)
@@ -69,4 +69,5 @@
     IF arg_present(distance) THEN distance = sqrt(distance)
   ENDELSE
-  RETURN, where(distance EQ min(distance))
+  minval = min(distance, indmin)
+  RETURN, indmin
 END
Index: trunk/SRC/ReadWrite/ncdf_getatt.pro
===================================================================
--- trunk/SRC/ReadWrite/ncdf_getatt.pro	(revision 500)
+++ trunk/SRC/ReadWrite/ncdf_getatt.pro	(revision 501)
@@ -58,4 +58,9 @@
 ; if this attribute does not exist.
 ;
+; @keyword interval_operation
+; Set this keyword to a named variable in which the value of
+; interval_operation is returned. Return -1 if this attribute
+; does not exist.
+;
 ; @keyword _EXTRA
 ; defined only to be able to call <pro>ncdf_getatt</pro> with the
@@ -76,5 +81,5 @@
                  , MISSING_VALUE = missing_value, UNITS = units, CALENDAR = calendar $
                  , LONG_NAME = long_name, STANDARD_NAME = standard_name, DOUBLE = double $
-                 , _extra = ex
+                 , INTERVAL_OPERATION = interval_operation, _extra = ex
 ;
   compile_opt idl2, strictarrsubs
@@ -94,4 +99,5 @@
   long_name = ''
   standard_name = ''
+  interval_operation = -1
 ;
   varinq = ncdf_varinq(cdfid, varid)
@@ -136,4 +142,9 @@
         '_fillvalue':ncdf_attget, cdfid, varid, attname, missing_value
         'missing_value':ncdf_attget, cdfid, varid, attname, missing_value
+        'interval_operation':BEGIN
+          ncdf_attget, cdfid, varid, attname, tmp
+          tmp = string(tmp)
+          interval_operation = fix(strmid(tmp,0,strpos(tmp,' ')))
+        END
         ELSE:
       ENDCASE
@@ -148,6 +159,7 @@
       missing_value = tmp
     ENDIF ELSE BEGIN
-      ras = report('Warning: missing value is not a number: ' + missing_value)
-      missing_value = 'no'
+      ;; ras = report('Warning: missing value is not a number: ' + missing_value)
+      ;; missing_value = 'no'
+      missing_value = (byte(strlowcase(a)))[0]
     ENDELSE
   ENDIF
Index: trunk/SRC/ToBeReviewed/GRILLE/grille.pro
===================================================================
--- trunk/SRC/ToBeReviewed/GRILLE/grille.pro	(revision 500)
+++ trunk/SRC/ToBeReviewed/GRILLE/grille.pro	(revision 501)
@@ -396,5 +396,5 @@
 ; for the vertical sections with roms
   IF keyword_set(type) AND n_elements(romszinfos) NE 0 THEN BEGIN
-    romsdp = romsdepth()
+    romsdp = romsdepth(/new_s_coord)
     IF romsdp[0] NE -1 THEN BEGIN
       IF jpt EQ 1 THEN BEGIN
Index: trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro
===================================================================
--- trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro	(revision 500)
+++ trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro	(revision 501)
@@ -120,5 +120,5 @@
                     , GRID = grid, CALLITSELF = callitself, DIREC = direc $
                     , ZETAFILENAME = zetafilename, ZETAZERO = zetazero $
-                    , ZINVAR = zinvar, _EXTRA = ex
+                    , ZINVAR = zinvar, ISROMS = isroms, _EXTRA = ex
 ;
   compile_opt idl2, strictarrsubs
@@ -164,17 +164,4 @@
     dimnames[i] = tmp
   ENDFOR
-;------------------------------------------------------------
-; we check if the variable has a vertical dimension and/or
-; a record dimension. This is useful to define boxzoom
-; keyword in domdef
-;------------------------------------------------------------
-  dummy = where(varinq.dim EQ inq.recdim, tinvar)
-; check the presence of a vertical dimension according to the
-; number of dimensions and the presence of a record dimension
-  zinvar = (varinq.ndims EQ 3 AND tinvar NE 1) OR varinq.ndims EQ 4
-;------------------------------------------------------------
-; shall we redefine the grid parameters
-;------------------------------------------------------------
-  IF keyword_set(init) THEN initncdf, filename, _extra = ex
 ;------------------------------------------------------------
 ; check the time axis and the debut and ending dates
@@ -248,4 +235,18 @@
     END
   ENDCASE
+  IF inq.recdim EQ -1 AND jpt NE 1 THEN inq.recdim = varinq.dim[varinq.ndims-1]
+;------------------------------------------------------------
+; we check if the variable has a vertical dimension and/or
+; a record dimension. This is useful to define boxzoom
+; keyword in domdef
+;------------------------------------------------------------
+  dummy = where(varinq.dim EQ inq.recdim, tinvar)
+; check the presence of a vertical dimension according to the
+; number of dimensions and the presence of a record dimension
+  zinvar = (varinq.ndims EQ 3 AND tinvar NE 1) OR varinq.ndims EQ 4
+;------------------------------------------------------------
+; shall we redefine the grid parameters
+;------------------------------------------------------------
+  IF keyword_set(init) THEN initncdf, filename, _extra = ex
 ;------------------------------------------------------------
 ; Name of the grid on which the field refer to.
@@ -405,5 +406,5 @@
 ; if it is roms outputs, we need to get additional infos...
   IF NOT keyword_set(callitself) THEN BEGIN
-    IF strmid(dimnames[0], 0, 3) EQ 'xi_' AND strmid(dimnames[1], 0, 4) EQ 'eta_' THEN BEGIN
+    IF (strmid(dimnames[0], 0, 3) EQ 'xi_' AND strmid(dimnames[1], 0, 4) EQ 'eta_') OR keyword_set(isroms) THEN BEGIN
       ncdf_attget, cdfid, 'theta_s', theta_s, /global
       ncdf_attget, cdfid, 'theta_b', theta_b, /global
Index: trunk/SRC/ToBeReviewed/PLOTS/DESSINE/sbar_plot.pro
===================================================================
--- trunk/SRC/ToBeReviewed/PLOTS/DESSINE/sbar_plot.pro	(revision 500)
+++ trunk/SRC/ToBeReviewed/PLOTS/DESSINE/sbar_plot.pro	(revision 501)
@@ -53,6 +53,6 @@
    if NOT keyword_set(NOREINITPLT) then reinitplt, _extra = ex
 ; 2) I place the drawing on the screen like on the postscript
-   IF chkstru(ex, 'overplot') EQ 0 THEN placedessin, 'autre', _extra = ex $
-   ELSE IF ex.overplot EQ 0 THEN placedessin, 'autre', _extra = ex
+   IF chkstru(ex, 'overplot') EQ 0 THEN placedessin, 'autre', /nocolorbar, _extra = ex $
+   ELSE IF ex.overplot EQ 0 THEN placedessin, 'autre', /nocolorbar, _extra = ex
 ; 3) Drawing
    if n_elements(COLORS) NE 0 then BEGIN
Index: trunk/SRC/ToBeReviewed/PLOTS/DESSINE/splot.pro
===================================================================
--- trunk/SRC/ToBeReviewed/PLOTS/DESSINE/splot.pro	(revision 500)
+++ trunk/SRC/ToBeReviewed/PLOTS/DESSINE/splot.pro	(revision 501)
@@ -45,5 +45,5 @@
    if NOT keyword_set(NOREINITPLT) then reinitplt, _extra = ex
 ; 2) i put the drawing on the screen like on the postscript
-   placedessin, 'autre', _extra = ex
+   placedessin, 'autre', /nocolorbar, _extra = ex
 ; 3) Drawing
    if n_elements(y) EQ 0 then plot,  x, xstyle = 1, ystyle = 1, _EXTRA = ex $
