Index: trunk/SRC/Grid/computegrid.pro
===================================================================
--- trunk/SRC/Grid/computegrid.pro	(revision 163)
+++ trunk/SRC/Grid/computegrid.pro	(revision 172)
@@ -4,95 +4,114 @@
 ;+
 ;
-; @file_comments compute the grid parameters from cm_4mesh common:
-;
-;       computegrid, startx, starty, stepx, stepy, nx, ny
-;       computegrid, startx, starty, stepx, stepy
-;       computegrid, xaxis = xaxis, yaxis = yaxis
-;       or a suitable mix...
-;
-; glamt
-; glamf
-; gphit
-; gphit
-; e1t
-; e2t
-; horizontal parameters
-;
-; glamu {in}
-; glamv {in}
-; gphiu {in}
-; gphiv {in}
-; e1u {in}
-; e1v {in}
-; e1f {in}
-; e2u {in}
-; e2v {in}
-; e2f {in}
-; horizontal parameters if FULLCGRID keyword is defined
-;
-; gdept
-; gdepw
-; e3t
-; e3w
-; verticals parameters
-;
+; @file_comments 
+;
+; compute the grid parameters (cm_4mesh) common
+;
+; domains sizes:
+;---------------
+; jpi, jpj, jpk, jpiglo, jpjglo, jpkglo, jpidta, jpjdta, jpkdta
+;
+; domains positions regarding to the original grid:
+;--------------------------------------------------
+; ixminmesh, ixmaxmesh, iyminmesh, iymaxmesh, izminmesh, izmaxmesh
+; ixmindta, ixmaxdta, iymindta, iymaxdta, izmindta, izmaxdta
+;
+; horizontal parameters:
+; ----------------------
+; glamt, glamf, gphit, gphit, e1t, e2t
+;
+; additional horizontal parameters if FULLCGRID keyword is defined:
+; -----------------------------------------------------------------
+; glamu, glamv, gphiu, gphiv, e1u, e1v, e1f, e2u, e2v, e2f
+;
+; verticals parameters:
+; ---------------------
+; gdept, gdepw, e3t, e3w
+;
+; masks:
+; ------
 ; tmask
-; masks
-;
-; umaskred {in}
-; vmaskred {in}
-; fmaskredx {in}
-; fmaskredy {in}
-; masks if FULLCGRID keyword is defined
-;
-; triangles_list
+;
+; additional masks if FULLCGRID keyword is defined:
+; -------------------------------------------------
+; umaskred, vmaskred, fmaskredx, fmaskredy
+;
+; triangles_list:
+; ---------------
 ; triangulation
+;
+; key_* variables:
+; ----------------
+; key_onearth, key_periodic, key_shift, key_stride, key_partialstep,
+; key_yreverse, key_zreverse, key_gridtype
+;
+; xxx related variables:
+; ----------------------
+; ccmeshparameters, ccreadparameters
 ;
 ; @categories 
 ; Grid
 ;
-; @param startx {in}{required} scalar, x starting point
-; @param starty {in}{required} scalar, y starting point
-; @param stepxin {in}{required} scalar or vector: x direction step, must be > 0
-;             if vector nx is not used
-; @param stepyin {in}{required} scalar or vector: y direction step,
-;             could be > 0 (south to north) or < 0 (north to south)
-;             if vector ny is not used
-; @param nxin {in}{required} scalar, number of points in x direction
-; @param nyin {in}{required} scalar, number of points in y direction
-;
-; @keyword FULLCGRID activate to specify that you want to compute
-;       all the parameters 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
-;
-; @keyword GLAMBOUNDARY a 2 elements vector, [lon1,lon2], the longitude
-;       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
-;
-; @keyword MASK to specify the mask with a 2 or 3 dimension array
-;
-; @keyword ONEARTH {default= key_onearth=1}
-;       = 0 or 1 to force the manual definition of
-;       key_onearth (to specify if the data are on earth -> use longitude
-;       /latitude etc...).
-;       note that ONEARTH = 0 forces PERIODIC = 0, SHIFT = 0,
-;       and is cancelling GLAMBOUNDARY
-;
-; @keyword PERIODIC {default=key_periodic is automatically computed by using the first line of glamt}
-;       = 0 or 1 to force the manual definition of key_periodic.
-;
-; @keyword PLAIN force PERIODIC = 0, SHIFT = 0, STRIDE = [1, 1, 1] and
+; @param startx {in}{optional}{type=scalar}
+;       x starting point, optional if [XY]AXIS keyword is used
+;
+; @param starty {in}{optional}{type=scalar}
+;       y starting point, optional if [XY]AXIS keyword is used
+;
+; @param stepxin {in}{optional}{type=scalar or vector}
+;       x direction step, optional if [XY]AXIS keyword is used, must be > 0
+;       if stepxin is a vector nx is not used
+;
+; @param stepyin {in}{optional}{type=scalar or vector}
+;       y direction step, optional if [XY]AXIS keyword is used,
+;       could be > 0 (south to north) or < 0 (north to south)
+;       if stepyin is a vector ny is not used
+;
+; @param nxin {in}{optional}{type=scalar}
+;       number of points in x direction,
+;       optional if [XY]AXIS keyword is used or stepxin is a vector
+;
+; @param nyin {in}{optional}{type=scalar} 
+;       number of points in y direction,
+;       optional if [XY]AXIS keyword is used or stepyin is a vector
+;
+; @keyword FULLCGRID {default=0}{type=scalar: 0 or 1}
+;       Activate to specify that you want to compute all the C grid parameters:
+;       definition of glam[uv], gphi[uv], e1[uvf], e2[uvf], [uv]maskred and
+;       fmaskred[xy] will be add to the default computations
+;
+; @keyword GLAMBOUNDARY {default=those defined in the file}{type=2 elements vector}
+;       Longitude boundaries that should be used to visualize the data.
+;         lon2 > lon1
+;         lon2 - lon1 le 360
+;       By default, the common (cm_4mesh) variable key_shift will be automatically
+;       defined according to GLAMBOUNDARY.
+;
+; @keyword MASK {default=array of 1}{type=2D or 3D array}
+;       Specify the land(0)/sea(1) mask
+;
+; @keyword ONEARTH {default=1}{type=scalar: 0 or 1}
+;       Force the manual definition of data localization on the earth or not
+;          0) if the data are not on the earth
+;          1) if the data are on earth (in that case we can for example use
+;             the labels 'longitude', 'latitude' in plots).
+;       The resulting value will be stored in the common (cm_4mesh) variable key_onearth
+;       ONEARTH = 0 forces PERIODIC = 0, SHIFT = 0 and is cancelling GLAMBOUNDARY
+;
+; @keyword PERIODIC {default=computed by using the first line of glamt}{type=scalar: 0 or 1}
+;       Force the manual definition of the grid zonal periodicity.
+;       The resulting value will be stored in the common (cm_4mesh) variable key_periodic
+;       PERIODIC = 0 forces SHIFT = 0
+;
+; @keyword PLAIN {default=0}{type=scalar: 0 or 1}
+;       Force PERIODIC = 0, SHIFT = 0, STRIDE = [1, 1, 1] and
 ;       suppress the automatic redefinition of the domain in case of
 ;       x periodicity overlap, y periodicity overlap (ORCA type only)
 ;       and mask border to 0.
 ;
-; @keyword SHIFT = scalar to force the manual definition of key_shift. By
-;       default, key_shift is automatically computed according to
-;       glamboundary (when defined) by using the FIRST LINE of glamt. if
-;       key_periodic=0 then in any case key_shift = 0.
+; @keyword SHIFT {default=computed according to glamboundary}{type=scalar}
+;       Force the manual definition of the zonal shift that must be apply to the data. 
+;       The resulting value will be stored in the common (cm_4mesh) variable key_shift
+;       Note that if key_periodic=0 then in any case key_shift = 0.
 ;
 ; @keyword STRCALLING {type=string}
@@ -100,19 +119,18 @@
 ;       call computegrid (this is used by xxx.pro)
 ;
-; @keyword STRIDE {default=[1, 1, 1]} a 3 elements vector to specify the stride in x, y, z
-;       direction. The resulting value
-;       will be stored in the common (cm_4mesh) variable key_stride
-;
-; @keyword 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.
-;
-; @keyword 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).
+; @keyword STRIDE {default=[1, 1, 1]}{type=3 elements vector}
+;       Specify the stride in x, y and z direction. The resulting
+;       value will be stored in the common (cm_4mesh) variable key_stride
+;
+; @keyword XAXIS {type=1D or 2D array}
+;       Specify longitudes 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.
+;
+; @keyword YAXIS {type=1D or 2D array}
+;       Specify latitudes 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).
 ;
 ; @keyword XYINDEX activate to specify that the horizontal grid should
@@ -121,28 +139,57 @@
 ;       using this keyword forces key_onearth=0
 ;
-; @keyword XMINMESH {default=0L}
-; @keyword YMINMESH {default=0L}
-; @keyword ZMINMESH {default=0L}
-;       to define the common variables i[xyz]minmesh
-;       used to define the grid only in a zoomed part of the original
-;       grid. max value is [XYZ]MAXMESH
-;
-; @keyword XMAXMESH {default=jpiglo-1}
-; @keyword YMAXMESH {default=jpjglo-1}
-; @keyword ZMAXMESH {default=jpkglo-1}
-;       to define the common variables i[xyz]maxmesh
-;       used to define the grid only in a zoomed part of the original
-;       grid. 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 	
+; @keyword XMINMESH {default=0L}{type=scalar}
+;       Define common (cm_4mesh) variables ixminmesh used to define the localization
+;       of the first point of the grid along the x direction in a zoom of the original grid
+;
+; @keyword YMINMESH {default=0L}{type=scalar}
+;       Define common (cm_4mesh) variables iyminmesh used to define the localization
+;       of the first point of the grid along the y direction in a zoom of the original grid
+;
+; @keyword ZMINMESH {default=0L}{type=scalar}
+;       Define common (cm_4mesh) variables izminmesh used to define the localization
+;       of the first point of the grid along the z direction in a zoom of the original grid
+;
+; @keyword XMAXMESH {default=jpiglo-1}{type=scalar}
+;       Define common (cm_4mesh) variables ixmaxmesh used to define the localization
+;       of the last point of the grid along the x direction in a zoom of the original grid
+;
+; @keyword YMAXMESH {default=jpjglo-1}{type=scalar}
+;       Define common (cm_4mesh) variables iymaxmesh used to define the localization
+;       of the last point of the grid along the y direction in a zoom of the original grid
+;
+; @keyword ZMAXMESH {default=jpkglo-1}{type=scalar}
+;       Define common (cm_4mesh) variables izmaxmesh used to define the localization
+;       of the last point of the grid along the z direction in a zoom of the original grid
 ;
 ; @keyword FBASE2TBASE
-;
-; @keyword STRCALLING
-;
-; @keyword ZAXIS to specify the vertical axis with a 1 dimension
-;       array. Must be sorted in the increasing or deceasing order
-;
-; @keyword _EXTRA used to pass your keywords to the created function.
+;       Activate when the model is a C grid based on a F point
+;       (with a F point at the bottom-left corner and a T point at the
+;       upper-right corner). In this case, we ignore 
+;         - the first line of F and U points
+;         - the first column of F and V points
+;         - the last line of T and U points
+;         - the last column of T and V points.
+;      => we are back to a C grid based on T point as for OPA model.
+;
+; @keyword ROMSH {type=2D array}
+;       This array is the final bathymetry at RHO-points. It is stored in the common
+;       variable (cm_4mesh) romszinfos.h
+;       Used when the model is a ROMS C-grid with one more point
+;       in longitude for T and V grid and one more point in latitude
+;       for T and U grid. In this case, we ignore 
+;         - the last line of T and U points
+;         - the last column of T and V points.
+;      => we are back to a C grid based on T point as for OPA model.
+;       Note that activate ROMSH forces FULLCGRID=1
+;
+; @keyword STRCALLING {type=scalar string}
+;       Used by xxx...
+;
+; @keyword ZAXIS {type=1D}
+;       Specify the vertical axis. Must be sorted in the increasing or deceasing order
+;
+; @keyword _EXTRA 
+;    (not used in the present case...)
 ;
 ; @uses cm_4mesh cm_4data cm_4cal
@@ -156,5 +203,5 @@
 ;    to use cell_fill = 1).
 ;
-; @restrictions FUV points definition...
+; @restrictions FUV points definition is not exact if the grid is irregular
 ;
 ; @history Sebastien Masson (smasson\@lodyc.jussieu.fr)
@@ -179,5 +226,5 @@
                  , FULLCGRID = fullcgrid, XYINDEX = xyindex $
                  , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
-                 , _extra = ex
+                 , ROMSH = romsh, _extra = ex
 ;---------------------------------------------------------
 ;
@@ -276,4 +323,9 @@
   jpjglo = long(ny)
   jpkglo = long(nz)
+  IF keyword_set(romsh) THEN BEGIN
+    jpiglo = jpiglo - 1
+    jpjglo = jpjglo - 1
+    fullcgrid = 1
+  ENDIF
 ;
 ; impact of plain keyword:
@@ -484,5 +536,5 @@
                  , FULLCGRID = fullcgrid, XYINDEX = xyindex $
                  , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
-                 , _extra = ex
+                 , ROMSH = romsh, _extra = ex
     return
   ENDIF
@@ -910,4 +962,5 @@
       gdepw = gdept
     ENDELSE
+    IF keyword_set(romsh) THEN gdepw = gdept 
 ;
 ;====================================================
@@ -926,8 +979,13 @@
 ;
 ; default 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])
+  if NOT keyword_set(mask) then tmask = -1 ELSE tmask = mask
+  IF keyword_set(romsh) AND tmask[0] NE -1 THEN tmask = tmask[0:jpiglo-1, 0:jpjglo-1]
+  IF array_equal(size(tmask, /dimensions), [jpiglo, jpjglo]) EQ 1 AND jpkglo GT 1 THEN BEGIN
+    tmask = tmask[*]#replicate(1, jpkglo)
+    tmask = reform(tmask, jpiglo, jpjglo, jpkglo, /overwrite)
+  ENDIF
+;
+  if tmask[0] NE -1 then BEGIN
+    tmask = byte(tmask[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)
@@ -976,5 +1034,5 @@
                      , FULLCGRID = fullcgrid, XYINDEX = xyindex $
                      , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $
-                     , _extra = ex
+                     , ROMSH = romsh, _extra = ex
         return
   ENDIF
@@ -1131,4 +1189,8 @@
     gphiinfo = gphit
   ENDELSE
+  IF keyword_set(romsh) THEN $
+     romszinfos = {h:romsh[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh], zeta:replicate(0., jpi, jpj), theta_s:-1, theta_b:-1, hc:-1} $
+  ELSE romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1}
+
   ccmeshparameters = {filename:strcalling  $
           , glaminfo:float(string(glaminfo, format = '(E11.4)')) $
