source: trunk/SRC/Grid/computegrid.pro

;+
;
; @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
;
; 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}{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 YREVERSE=0, ZREVERSE=0, 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 {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}
;       a string containing the calling command used to
;       call computegrid (this is used by <pro>xxx</pro>)
;
; @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 decreasing order (along each column if 2d array).
;
; @keyword 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
;
; @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
;       Note that if XMAXMESH < 0 then ixmaxmesh is defined as ixmaxmesh = jpiglo -1 + xmaxmesh
;
; @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
;       Note that if YMAXMESH < 0 then iymaxmesh is defined as iymaxmesh = jpjglo -1 + ymaxmesh
;
; @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
;       Note that if ZMAXMESH < 0 then izmaxmesh is defined as izmaxmesh = jpkglo -1 + maxmesh
;
; @keyword FBASE2TBASE
;       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 V points
;         - the last  line of T and U points
;         - if the grid is not x-periodic, the first column of F and U points
;         - if the grid is not x-periodic, 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 in that case, key_gridtype = 'c_f' and not 'c' (-> used in read_ncdf)
;       Note that activate FBASE2TBASE forces FULLCGRID=1
;
; @keyword FFBASE2TBASE
;       Activate when the model is a C grid based on a F point
;       (with a F point at the bottom-left corner and again a F point at the
;       upper-right corner). In this case, we ignore
;         - the first line of F and V points
;         - if the grid is not x-periodic, the first column of F and U points
;       => we are back to a C grid based on T point as for OPA model.
;       Note that in that case, key_gridtype = 'c_ff' and not 'c' (-> used in read_ncdf)
;       Note that activate FFBASE2TBASE forces FULLCGRID=1
;
; @keyword UBASE2TBASE
;       Activate when the model is a C grid based on a U point
;       (with a U point at the bottom-left corner and a T point at the
;       upper-right corner). In this case, we ignore
;         - if the grid is not x-periodic, the first column of F and U points
;         - if the grid is not x-periodic, 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 in that case, key_gridtype = 'c_u' and not 'c' (-> used in read_ncdf)
;       Note that activate UBASE2TBASE forces FULLCGRID=1
;
; @keyword VBASE2TBASE
;       Activate when the model is a C grid based on a V point
;       (with a V 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 V points
;         - the last  line of T and U points
;       => we are back to a C grid based on T point as for OPA model.
;       Note that in that case, key_gridtype = 'c_v' and not 'c' (-> used in read_ncdf)
;       Note that activate VBASE2TBASE forces FULLCGRID=1
;
; @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 YREVERSE {default=computed according to gphit[0, 1:jpj-1] LT gphit[0, 0:jpj-2]}{type=scalar}
;       Force the manual definition of the y reverse that must be apply to the data.
;       The resulting value will be stored in the common (cm_4mesh) variable key_yreverse
;
; @keyword ZAXIS {type=1D}
;       Specify the vertical axis. Must be sorted in the increasing or decreasing order
;
; @keyword ZZAXIS {type=1D}
;       same a zaxis (overwrite zaxis definition). Needed when
;       computegrid is called by inincdf because inincdf as alreday a keyword with a
;       name sarting by zaxis: zaxisname
;
; @keyword ZREVERSE {default=computed according to gdept[0] GT gdept[1]}{type=scalar}
;       Force the manual definition of the z reverse that must be apply to the data.
;       The resulting value will be stored in the common (cm_4mesh) variable key_zreverse
;
; @keyword _EXTRA
; not used in the present case ...
;
; @uses
; <pro>cm_4mesh</pro>
; <pro>cm_4data</pro>
; <pro>cm_4cal</pro>
;

; @restrictions
; if the grid has x/y periodicity overlap and/or if
; the mask has 0 everywhere at the border (like a closed sea) and
; if (we did not activate /plain and xminmesh, xmaxmesh, yminmesh,
; ymaxmesh keywords are defined to their default values), we redefine
; xminmesh, xmaxmesh, yminmesh, ymaxmesh in order to remove the
; overlapping part and/or to open the domain (avoid it be forced
; to use cell_fill = 1).
;
; FUV points definition is not exact if the grid is irregular
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      2000-04-20
; Sept 2004, several bug fix to suit C grid type...
; Aug 2005, rewrite almost everything...
;
; @version
; $Id$
;
;-
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
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
  ENDIF
;---------------------------------------------------------
;------------------------------------------------------------
  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(zzaxis) NE 0 then zaxis = zzaxis
  if n_elements(zaxis) NE 0 then BEGIN
    CASE (size(zaxis))[0] OF
      0:nz = 1L
      1:nz = (size(zaxis))[1]
      ELSE:BEGIN
        ras = report( 'not coded')
        stop
      END
    ENDCASE
  ENDIF ELSE nz = 1L
;
;====================================================
; Others automatic definitions...
;====================================================
;
  IF keyword_set(romsh) THEN fullcgrid = 1
  CASE 1 OF
    keyword_set(ffbase2tbase):key_gridtype = 'c_ff'
    keyword_set(fbase2tbase):key_gridtype = 'c_f'
    keyword_set(ubase2tbase):key_gridtype = 'c_u'
    keyword_set(vbase2tbase):key_gridtype = 'c_v'
    else:key_gridtype = 'c'
  ENDCASE
  IF strlen(key_gridtype) NE 1 THEN fullcgrid = 1
;
  IF n_elements(xminmesh) NE 0 OR n_elements(xmaxmesh) NE 0 THEN BEGIN
    IF n_elements(xminmesh) EQ 0 THEN xminmesh = ixminmesh
    IF n_elements(xmaxmesh) EQ 0 THEN xmaxmesh = ixmaxmesh
    IF nx EQ jpi AND xminmesh EQ ixminmesh AND xmaxmesh EQ ixmaxmesh THEN xalreadycut = 1
  ENDIF
  IF keyword_set(xalreadycut) THEN BEGIN
    xmin = 0
    xmax = jpi - 1
    nxx = jpi
  ENDIF ELSE BEGIN
    jpiglo = long(nx)
    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 ixmaxmesh LT 0 THEN ixmaxmesh = jpiglo -1 + ixmaxmesh
    ixmaxmesh = 0 > ixmaxmesh < (jpiglo-1)
    ixminmesh = 0 > ixminmesh < ixmaxmesh
    jpi = ixmaxmesh-ixminmesh+1
    xmin = ixminmesh
    xmax = ixmaxmesh
    nxx = jpiglo
  ENDELSE

  IF n_elements(yminmesh) NE 0 OR n_elements(ymaxmesh) NE 0 THEN BEGIN
    IF n_elements(yminmesh) EQ 0 THEN yminmesh = iyminmesh
    IF n_elements(ymaxmesh) EQ 0 THEN ymaxmesh = iymaxmesh
    IF ny EQ jpj AND yminmesh EQ iyminmesh AND ymaxmesh EQ iymaxmesh THEN yalreadycut = 1
  ENDIF
  IF keyword_set(yalreadycut) THEN BEGIN
    ymin = 0
    ymax = jpj - 1
    nyy = jpj
  ENDIF ELSE BEGIN
    jpjglo = long(ny)
    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 key_gridtype EQ 'c_v' OR key_gridtype EQ 'c_f' THEN iymaxmesh = iymaxmesh-1
    IF iymaxmesh LT 0 THEN iymaxmesh = jpjglo -1 + iymaxmesh
    iymaxmesh = 0 > iymaxmesh < (jpjglo-1)
    iyminmesh = 0 > iyminmesh < iymaxmesh
    jpj = iymaxmesh-iyminmesh+1
    ymin = iyminmesh
    ymax = iymaxmesh
    nyy = jpjglo
  ENDELSE

  IF n_elements(zminmesh) NE 0 AND n_elements(zmaxmesh) NE 0 THEN BEGIN
    IF nz EQ jpk AND zminmesh EQ izminmesh AND zmaxmesh EQ izmaxmesh THEN zalreadycut = 1
  ENDIF
  IF keyword_set(zalreadycut) THEN BEGIN
    zmin = 0
    zmax = jpk - 1
    nzz = jpk
  ENDIF ELSE BEGIN
    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
    IF izmaxmesh LT 0 THEN izmaxmesh = jpkglo -1 + izmaxmesh
    izmaxmesh = 0 > izmaxmesh < (jpkglo-1)
    izminmesh = 0 > izminmesh < izmaxmesh
    jpk = izmaxmesh-izminmesh+1
    zmin = izminmesh
    zmax = izmaxmesh
    nzz = jpkglo
  ENDELSE
;
; impact of plain keyword:
;
  IF keyword_set(plain) THEN BEGIN
    yreverse = 0
    zreverse = 0
    periodic = 0
    shift = 0
    stride = [1, 1, 1]
  ENDIF
;
; avoid basics errors...
;
  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]
;
; 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
;
  IF n_elements(glamt) EQ nxx*nyy THEN glamt = reform(glamt, nxx, nyy, /over) $
  ELSE glamt = reform(glamt, nxx, /over)
  CASE size(glamt, /n_dimensions) OF
    1:BEGIN
      IF n_elements(glamt) EQ 1 THEN glamt = replicate(glamt[0], jpi, jpj) $
      ELSE glamt = glamt[xmin:xmax]#replicate(1, jpj)
    END
    2:glamt = glamt[xmin:xmax, ymin:ymax]
  ENDCASE
; keep 2d array even with degenerated 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
;
  IF n_elements(gphit) EQ nxx*nyy THEN gphit = reform(gphit, nxx, nyy, /over) $
  ELSE gphit = reform(gphit, nyy, /over)
  CASE size(gphit, /n_dimensions) OF
    1:BEGIN
      IF n_elements(gphit) EQ 1 THEN gphit = replicate(gphit[0], jpi, jpj) $
      ELSE gphit = replicate(1, jpi)#gphit[ymin:ymax]
    END
    2:gphit = gphit[xmin:xmax, ymin:ymax]
  ENDCASE
; keep 2d array even with degenerated dimension
  IF jpj EQ 1 THEN gphit = reform(gphit, jpi, jpj, /over)
;
;====================================================
; check y periodicity... Only according to ORCA grid
;====================================================
; 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

    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
; T pivot
        ymaxmesh = -1
        recall = 1
      END
      ixminmesh EQ 1l AND ixmaxmesh eq jpiglo-2 $
         AND array_equal(gphit[*, jpj-1], reverse(shift(gphit[*, jpj-3], -1))) EQ 1:BEGIN
; T pivot
        ymaxmesh = -1
        recall = 1
      END
      ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
         AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN
; F pivot
        ymaxmesh = -1
        recall = 1
      END
      ixminmesh EQ 1l AND ixmaxmesh eq jpiglo-2 $
         AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN
; F pivot
        ymaxmesh = -1
        recall = 1
      END
      ELSE:
    ENDCASE
  ENDIF
;
;====================================================
; check x periodicity...
;====================================================
  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 $
     AND ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 AND jpi GE 3 THEN BEGIN
    CASE 0 OF
      total((glamt[0, *] - glamt[jpi-2, *]) MOD 360) $
         + total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN
        xminmesh = 1
        xmaxmesh = -1
        recall = 1
      END
      total((glamt[0, *] - glamt[jpi-2, *]) MOD 360):BEGIN
        xminmesh = 1
        recall = 1
      END
      total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN
        xmaxmesh = -1
        recall = 1
      END
      total((glamt[0, *] - glamt[jpi-1, *]) MOD 360):BEGIN
        xmaxmesh = -1
        recall = 1
      END
      ELSE:
    ENDCASE
  ENDIF
;====================================================
; recall computegrid if needed...
;====================================================
  IF keyword_set(recall) THEN BEGIN
    computegrid, XAXIS = glamt, YAXIS = gphit, ZAXIS = zaxis $
                 , MASK = mask, GLAMBOUNDARY = glamboundary $
                 , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
                 , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
                 , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
                 , PERIODIC = periodic, SHIFT = shift, STRIDE = stride $
                 , FULLCGRID = fullcgrid, XYINDEX = xyindex $
                 , STRCALLING = strcalling $
                 , ROMSH = romsh, _extra = ex
    return
  ENDIF
;
;====================================================
; def key_yreverse
;====================================================
;
  IF n_elements(yreverse) EQ 0 THEN BEGIN
    IF jpj GT 1 THEN BEGIN
      IF total(gphit[0, 1:jpj-1] LT gphit[0, 0:jpj-2]) GT jpj/2 THEN key_yreverse = 1 ELSE key_yreverse = 0
    ENDIF ELSE key_yreverse = 0
  ENDIF ELSE key_yreverse = yreverse
  IF keyword_set(key_yreverse) THEN BEGIN
    gphit = reverse(gphit, 2)
    glamt = reverse(glamt, 2)
  ENDIF
;
;====================================================
; def of key_shift
;====================================================
;
; definition of key_shift by shifting 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
        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...'])
        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 (key_gridtype EQ 'c_u' OR key_gridtype EQ 'c_f') AND NOT keyword_set(key_periodic) 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
;
;====================================================
; 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
  ENDIF ELSE glamf = glamt + 0.5
;
  IF keyword_set(key_periodic) AND (max(glamf)-min(glamt)) GE 360 THEN BEGIN
    IF NOT keyword_set(glamboundary) THEN BEGIN
      bigger = where(glamf GE min(glamt)+360)
      glamf[bigger] = glamf[bigger]-360.
    ENDIF ELSE glamf = glamboundary[0] > temporary(glamf) < glamboundary[1]
  ENDIF
;
  IF jpj EQ 1 THEN glamf = reform(glamf, jpi, jpj, /over)
;
;====================================================
; 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.
;
  IF jpj EQ 1 THEN gphif = reform(gphif, jpi, jpj, /over)
;
;====================================================
; 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, *]
    ENDIF ELSE e1t = replicate(stepx, jpi, jpj)
  ENDIF ELSE e1t = replicate(1b, jpi, jpj)
;
  IF jpj EQ 1 THEN e1t = reform(e1t, jpi, jpj, /over)
;
;====================================================
; 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 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    
;
;====================================================
; 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
;
    IF jpj EQ 1 THEN e2u = reform(e2u, jpi, jpj, /over)
;
;====================================================
; 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
;
    IF jpj EQ 1 THEN e1v = reform(e1v, jpi, jpj, /over)
;
;====================================================
; 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
;
    IF jpj EQ 1 THEN e1f = reform(e1f, jpi, jpj, /over)
;
;====================================================
; 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
;
    IF jpj EQ 1 THEN e2f = reform(e2f, jpi, jpj, /over)
;
  ENDIF
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
;
;====================================================
; e1[tuvf] from degree to meters
;====================================================
;
  IF keyword_set(key_onearth)  THEN BEGIN
    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
;
  IF jpj EQ 1 THEN BEGIN
    e1t = reform(e1t, jpi, jpj, /over)
    IF keyword_set(fullcgrid) THEN BEGIN
      e1u = reform(e1u, jpi, jpj, /over)
      e1v = reform(e1v, jpi, jpj, /over)
      e1f = reform(e1f, jpi, jpj, /over)
    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[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
;
;====================================================
; 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
;====================================================
;
; default mask eq 1
  if NOT keyword_set(mask) then tmask = -1 ELSE tmask = mask
;
  if tmask[0] NE -1 then BEGIN
    tmask = byte(temporary(tmask))
    IF n_elements(mask) EQ nxx*nyy AND nzz GT 1 THEN BEGIN
      tmask = tmask[*]#replicate(1b, nzz)
      tmask = reform(tmask, nxx, nyy, nzz, /overwrite)
    ENDIF
    IF nxx EQ 1 OR nyy EQ 1 OR nzz EQ 1 THEN tmask = reform(tmask, nxx, nyy, nzz, /overwrite)
    tmask = tmask[xmin:xmax, ymin:ymax, zmin:zmax]
    IF jpi EQ 1 OR jpj EQ 1 OR jpk EQ 1 THEN 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, *]
    ENDIF
    fmaskredy = tmask[jpi-1, *, *]
    fmaskredx = tmask[*, jpj-1, *]
  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)
    ENDIF
    fmaskredy = replicate(1b, jpj, jpk)
    fmaskredx = replicate(1b, jpi, jpk)
  ENDELSE
;
  IF jpi GT 2 AND jpj GT 2 AND NOT keyword_set(plain) $
;;      AND ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $
;;      AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 $
     AND total(tmask[*, 0, *]) EQ 0 AND total(tmask[*, jpj-1, *]) EQ 0 $
     AND total(tmask[0, *, *]) EQ 0 AND total(tmask[jpi-1, *, *]) EQ 0 THEN BEGIN
    IF NOT keyword_set(key_periodic) THEN BEGIN
      xminmesh = ixminmesh + 1
      xmaxmesh = ixmaxmesh - 1
    ENDIF
    yminmesh = iyminmesh + 1
    ymaxmesh = iymaxmesh - 1
; come back to the original grid before calling computegrid with the
; new parameters...
    IF keyword_set(yreverse) THEN BEGIN
      gphit = reverse(gphit, 2)
      glamt = reverse(glamt, 2)
      tmask = reverse(tmask, 2)
      IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
    ENDIF
    IF keyword_set(zreverse) THEN BEGIN
      zaxis = reverse(zaxis)
      tmask = reverse(tmask, 3)
      IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
    ENDIF
    IF key_shift NE 0 THEN BEGIN
      glamt = shift(glamt, -key_shift, 0)
      gphit = shift(gphit, -key_shift, 0)
      tmask = shift(tmask, -key_shift, 0, 0)
      IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
    ENDIF
    IF ixminmesh NE 0 THEN BEGIN
      glamt = [fltarr(ixminmesh, jpj), glamt]
      gphit = [fltarr(ixminmesh, jpj), gphit]
      tmask = [fltarr(ixminmesh, jpj), tmask]
      jpi = jpi+ixminmesh
      IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
    ENDIF
    IF ixmaxmesh NE jpiglo-1 THEN BEGIN
      glamt = [glamt, fltarr(jpiglo-1-ixmaxmesh, jpj)]
      gphit = [gphit, fltarr(jpiglo-1-ixmaxmesh, jpj)]
      tmask = [tmask, fltarr(jpiglo-1-ixmaxmesh, jpj)]
      jpi = jpi+jpiglo-1-ixmaxmesh
      IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
    ENDIF
    IF iyminmesh NE 0 THEN BEGIN
      glamt = [[fltarr(jpi, iyminmesh)], [glamt]]
      gphit = [[fltarr(jpi, iyminmesh)], [gphit]]
      tmask = [[fltarr(jpi, iyminmesh)], [tmask]]
      jpj = jpj+iyminmesh
      IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
    ENDIF
    IF iymaxmesh NE jpjglo-1 THEN BEGIN
      glamt = [[glamt], [fltarr(jpi, jpjglo-1-iymaxmesh)]]
      gphit = [[gphit], [fltarr(jpi, jpjglo-1-iymaxmesh)]]
      tmask = [[tmask], [fltarr(jpi, jpjglo-1-iymaxmesh)]]
      jpj = jpj+jpjglo-1-iymaxmesh
      IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over)
    ENDIF
    computegrid, XAXIS = glamt, YAXIS = gphit, ZAXIS = zaxis $
                 , MASK = tmask, GLAMBOUNDARY = glamboundary $
                 , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
                 , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
                 , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
                 , ONEARTH = onearth, PERIODIC = key_periodic $
                 , PLAIN = plain, SHIFT = key_shift, STRIDE = key_stride $
                 , FULLCGRID = fullcgrid, XYINDEX = xyindex $
                 , UBASE2TBASE = ubase2tbase, VBASE2TBASE = vbase2tbase $
                 , FBASE2TBASE = fbase2tbase, FFBASE2TBASE = ffbase2tbase, STRCALLING = strcalling $
                 , ROMSH = romsh, _extra = ex
    return
  ENDIF
;
  IF NOT keyword_set(fullcgrid) THEN BEGIN
    umaskred = !values.f_nan
    vmaskred = !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
    ENDIF ELSE glamf = glamt + 0.5
    IF keyword_set(key_periodic) AND (max(glamf)-min(glamt)) GE 360 THEN BEGIN
      IF NOT keyword_set(glamboundary) THEN BEGIN
        bigger = where(glamf GE min(glamt)+360)
        glamf[bigger] = glamf[bigger]-360.
      ENDIF ELSE glamf = glamboundary[0] > temporary(glamf) < glamboundary[1]
    ENDIF
    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 jpj EQ 1 THEN BEGIN
      glamt = reform(glamt, jpi, jpj, /over)
      gphit = reform(gphit, jpi, jpj, /over)
      glamf = reform(glamf, jpi, jpj, /over)
      gphif = reform(gphif, jpi, jpj, /over)
      e1t = reform(e1t, jpi, jpj, /over)
      e2t = reform(e2t, jpi, jpj, /over)
    ENDIF
;
    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]]
      IF jpj EQ 1 THEN BEGIN
        glamu = reform(glamu, jpi, jpj, /over)
        gphiu = reform(gphiu, jpi, jpj, /over)
        e1u = reform(e1u, jpi, jpj, /over)
        e2u = reform(e2u, jpi, jpj, /over)
        glamv = reform(glamv, jpi, jpj, /over)
        gphiv = reform(gphiv, jpi, jpj, /over)
        e1v = reform(e1v, jpi, jpj, /over)
        e2v = reform(e2v, jpi, jpj, /over)
        e1f = reform(e1f, jpi, jpj, /over)
        e2f = reform(e2f, jpi, jpj, /over)
      ENDIF
    ENDIF
  ENDIF
;
;====================================================
; apply all the grid parameters
;====================================================
;
  @updateold
  domdef
;
;====================================================
; Triangulation
;====================================================
;
  IF total(tmask) EQ jpi*jpj*jpk $
     AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $
  ELSE BEGIN
; are we using ORCA2 ?
    IF jpiglo EQ 182 AND jpi EQ 180 AND jpjglo EQ 149 AND jpj EQ 148 THEN $
       triangles_list = triangule() ELSE triangles_list = triangule(/keep_cont)
  ENDELSE
;
;====================================================
; time axis (default definition)
;====================================================
;
  IF n_elements(time) EQ 0 OR n_elements(jpt) EQ 0 THEN BEGIN
    jpt = 1
    time = 0
  ENDIF
;
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updateold
  ENDIF
;====================================================
; grid parameters used by xxx
;====================================================
;
  IF NOT keyword_set(strcalling) THEN BEGIN
    IF n_elements(ccmeshparameters) EQ 0 THEN strcalling = 'computegrid' $
    ELSE strcalling = ccmeshparameters.filename
  ENDIF
  IF n_elements(glamt) GE 2 THEN BEGIN
    glaminfo = moment(glamt)
    IF finite(glaminfo[2]) EQ 0 THEN glaminfo = glaminfo[0:1]
    gphiinfo = moment(gphit)
    IF finite(gphiinfo[2]) EQ 0 THEN gphiinfo = gphiinfo[0:1]
  ENDIF ELSE BEGIN
    glaminfo = glamt
    gphiinfo = gphit
  ENDELSE
  IF keyword_set(romsh) THEN $
     romszinfos = {h:romsh[xmin:xmax, ymin:ymax], 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)')) $
                      , 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}
;------------------------------------------------------------
  IF keyword_set(key_performance) EQ 1 THEN $
     print, 'time computegrid', systime(1)-time1
;------------------------------------------------------------
  return
end