;+
;
; @file_comments
; Allows to extract a sub-domain of study by providing parameters
; needed for drawings (see outputs)
;
; @categories
; Grid
;
; @param Z1 {in}{optional}
; For a 3d domain whose the horizontal part cover all glam
;
; @param Z2 {in}{optional}
; For a 3d domain whose the horizontal part cover all gphi
;
; @param X1 {in}{optional}
; Define the minimum longitude. (All levels are selected)
;
; @param X2 {in}{optional}
; Define the maximum longitude. (All levels are selected)
;
; @param Y1 {in}{optional}
; Define the minimum latitude. (All levels are selected)
;
; @param Y2 {in}{optional}
; Define the maximum latitude. (All levels are selected)
;
; @keyword ENDPOINTS {type=vector}
; A four elements vector [x1,y1,x2,y2] used to specify
; that <pro>domdef</pro> must define the box used to make a plot
; (<pro>pltz</pro>, <pro>pltt</pro>, <pro>plt1d</pro>) done strictly along the line (that can have any direction)
; starting at (x1, y1) ending at (x2, y2). When defining endpoints,
; you must also define TYPE which define the type of plots
; ('pltz', 'xt', 'yt', 'zt', 'x', 'y', 'z', 't') will used
; ENDPOINTS keywords
;
; @keyword FINDALWAYS
; Force to redefine a box even when none point is find in the box.
; In this case, we select all the grid.
;
; @keyword GRIDTYPE {type=string or vector}
; It is a string or a vector of strings containing the grids's name
; (Only determined by 'T','U','V','W','F') for which the calculation
; must be done.
;  For example: 'T' or ['T','U']
;
; @keyword MEMEINDICES
; It is possible that points T, U, V and F correspond to a same geographic
; box which do not concern the same array indexes. This is sometimes a
; problem (or at least serious complications) in programs where several
; type of grid intervene (see <pro>norme</pro>, <pro>curl</pro>...).
; Activate MEMEINDICES to force <pro>domdef</pro> to take same indexes - those
; of the grid T- for all other grids.
;
; @keyword INDEX
; We activate it if we want that all elements passed in input of
; <pro>domdef</pro> refer to indexes of glam, gphi and gdep arrays rather
; than to values of these arrays.
;
; @keyword TYPE 
;
; @keyword XINDEX
; We activate it if we want that all elements passed in input of
; <pro>domdef</pro>
; and concerning the X dimension refer to indexes of glam arrays rather
; than to values of these arrays.
;
; @keyword YINDEX
; We activate it if we want that all elements passed in input of
; <pro>domdef</pro>
; and concerning the X dimension refer to indexes of gphi arrays rather
; than to values of these arrays.
;
; @keyword ZINDEX
; We activate it if we want that all elements passed in input of
; <pro>domdef</pro>
; and concerning the X dimension refer to indexes of gdep arrays rather
; than to values of these arrays.
;
; @uses
; <pro>common</pro>
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;			8/2/98
; rewrite everything, debug and speed-up Sebastien Masson April 2005
;
; @version
; $Id$
;
; @todo
; seb: output pas clair/ pas d'input required?
;-
PRO domdef, x1, x2, y1, y2, z1, z2, FINDALWAYS=findalways $
            , GRIDTYPE=gridtype, MEMEINDICES=memeindices $
            , XINDEX=xindex, YINDEX=yindex, ZINDEX=zindex $
            , ENDPOINTS=endpoints, TYPE=type $
            , INDEX=index, _EXTRA=ex
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
  ENDIF
;---------------------
  tempsun = systime(1)          ; For key_performance
;
  CASE N_PARAMS() OF
    0:
    1:
    2:
    4:
    6:
    ELSE:BEGIN
      ras = report('Bad number of parameter in the call of domdef.')
      RETURN
    END
  ENDCASE
;
  IF keyword_set(endpoints) THEN BEGIN
    IF NOT keyword_set(type) THEN BEGIN
      dummy = report(['If domdef is used do find the box associated' $
                      , 'to endpoints, you must also specify type keyword'])
      return
    ENDIF
    CASE N_PARAMS() OF
      0:
      1:boxzoom = [x1]
      2:boxzoom = [x1, x2]
      4:boxzoom = [x1, x2, y1, y2]
      6:boxzoom = [x1, x2, y1, y2, z1, z2]
    ENDCASE
    section, BOXZOOM = boxzoom, ENDPOINTS = endpoints, TYPE = type, /ONLYBOX
    return
  ENDIF

;-------------------------------------------------------------------
; recall domdef when there is only one input parameter
;-------------------------------------------------------------------
;
  IF N_PARAMS() EQ 1 THEN BEGIN
    CASE n_elements(x1) OF
      2:domdef, x1[0], x1[1], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
      4:domdef, x1[0], x1[1], x1[2], x1[3], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
      6:domdef, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
      ELSE:BEGIN
        ras = report('Bad number of elements in x1')
        RETURN
      END
    ENDCASE
    RETURN
  ENDIF
;-------------------------------------------------------------------
; default definitions and checks
;-------------------------------------------------------------------
  IF NOT keyword_set(gridtype) THEN gdtype = ['T', 'U', 'V', 'W', 'F'] $
  ELSE gdtype = strupcase(gridtype)
  IF keyword_set(memeindices) THEN gdtype = ['T', gdtype]
  IF finite(glamu[0]) eq 0 THEN gdtype = gdtype[where(gdtype NE 'U')]
  IF finite(glamv[0]) eq 0 THEN gdtype = gdtype[where(gdtype NE 'V')]
; default definitions
  lon1t = 99999. & lon2t = -99999. & lat1t = 99999. & lat2t = -99999.
  lon1u = 99999. & lon2u = -99999. & lat1u = 99999. & lat2u = -99999.
  lon1v = 99999. & lon2v = -99999. & lat1v = 99999. & lat2v = -99999.
  lon1f = 99999. & lon2f = -99999. & lat1f = 99999. & lat2f = -99999.
  vert1t = 99999. & vert2t = -99999. & vert1w = 99999. & vert2w = -99999.
;
  IF jpj EQ 1 THEN BEGIN
    IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1 THEN BEGIN
      glamt = reform(glamt, jpi, jpj, /over)
      gphit = reform(gphit, jpi, jpj, /over)
    ENDIF
    IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
      glamu = reform(glamu, jpi, jpj, /over)
      gphiu = reform(gphiu, jpi, jpj, /over)
    ENDIF
    IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
      glamv = reform(glamv, jpi, jpj, /over)
      gphiv = reform(gphiv, jpi, jpj, /over)
    ENDIF
    IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
      glamf = reform(glamf, jpi, jpj, /over)
      gphif = reform(gphif, jpi, jpj, /over)
    ENDIF
  ENDIF
;
  IF N_PARAMS() EQ 2 THEN GOTO, vertical
;
;-------------------------------------------------------------------
;-------------------------------------------------------------------
; define all horizontal parameters ...
; lon1 and lon2
; lat1 and lat2
; firstx[tuvf], lastx[tuvf], nx[tuvf]
;-------------------------------------------------------------------
; check if the grid is defined for U and V points. If not, take care
; of the cases gdtype eq 'U' or 'V'
;-------------------------------------------------------------------
  errstatus = 0
  IF (finite(glamu[0]*gphiu[0]) EQ 0 OR n_elements(glamu) EQ 0 OR n_elements(gphiu) EQ 0) AND (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
    firstxu = !values.f_nan
    lastxu = !values.f_nan
    nxu = !values.f_nan
    okgrid = where(gdtype NE 'U', count)
    IF count NE 0 THEN gdtype = gdtype[okgrid] $
    ELSE errstatus = report('U grid is undefined. Impossible to call domdef with vargid = ''U''')
  ENDIF
;
  IF (finite(glamv[0]*gphiv[0]) EQ 0 OR n_elements(glamv) EQ 0 OR n_elements(gphiv) EQ 0) AND (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
    firstxv = !values.f_nan
    lastxv = !values.f_nan
    nxv = !values.f_nan
    okgrid = where(gdtype NE 'V', count)
    IF count NE 0 THEN gdtype = gdtype[okgrid] $
    ELSE errstatus = report('V grid is undefined. Impossible to call domdef with vargid = ''V''')
  ENDIF
  IF errstatus EQ -1 THEN return
;-------------------------------------------------------------------
;-------------------------------------------------------------------
; horizontal domain defined with lon1, lon2, lat1 and lat2
;-------------------------------------------------------------------
;-------------------------------------------------------------------
  IF N_PARAMS() EQ 0 $
    OR ( (N_PARAMS() EQ 4 OR N_PARAMS() EQ 6) $
         AND NOT keyword_set(xindex) AND NOT keyword_set(yindex) AND NOT keyword_set(index) ) THEN BEGIN
    IF N_PARAMS() EQ 0 THEN BEGIN
; find lon1 and lon2 the longitudinal boundaries of the full domain
      IF (where(gdtype eq 'T'))[0] NE -1 THEN lon1t = min(glamt, max = lon2t)
      IF (where(gdtype eq 'W'))[0] NE -1 AND (where(gdtype eq 'T'))[0] EQ -1 THEN lon1t = min(glamt, max = lon2t)
      IF (where(gdtype eq 'U'))[0] NE -1 THEN lon1u = min(glamu, max = lon2u)
      IF (where(gdtype eq 'V'))[0] NE -1 THEN lon1v = min(glamv, max = lon2v)
      IF (where(gdtype eq 'F'))[0] NE -1 THEN lon1f = min(glamf, max = lon2f)
      lon1 = min([lon1t, lon1u, lon1v, lon1f])
      lon2 = max([lon2t, lon2u, lon2v, lon2f])
; find lat1 and lat2 the latitudinal boundaries of the full domain
      IF (where(gdtype eq 'T'))[0] NE -1 THEN lat1t = min(gphit, max = lat2t)
      IF (where(gdtype eq 'W'))[0] NE -1 AND (where(gdtype eq 'T'))[0] EQ -1 THEN lat1t = min(gphit, max = lat2t)
      IF (where(gdtype eq 'U'))[0] NE -1 THEN lat1u = min(gphiu, max = lat2u)
      IF (where(gdtype eq 'V'))[0] NE -1 THEN lat1v = min(gphiv, max = lat2v)
      IF (where(gdtype eq 'F'))[0] NE -1 THEN lat1f = min(gphif, max = lat2f)
      lat1 = min([lat1t, lat1u, lat1v, lat1f])
      lat2 = max([lat2t, lat2u, lat2v, lat2f])
    ENDIF ELSE BEGIN
      lon1 = min([x1, x2], max = lon2)
      lat1 = min([y1, y2], max = lat2)
    ENDELSE
; find firstxt, firstxt, nxt and nyt according to lon1, lon2, lat1 and lat2
    IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1 THEN BEGIN
      dom = where( (glamt GE lon1) AND (glamt LE lon2) $
                   AND (gphit GE lat1) AND (gphit LE lat2) )
      IF (dom[0] EQ -1) THEN BEGIN
        IF keyword_set(findalways) THEN BEGIN
          ras = report('WARNING, empty T points box... we get the neighbors to define a new box...')
          neig1 = neighbor(lon1, lat1, glamt, gphit, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
          neig2 = neighbor(lon2, lat2, glamt, gphit, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
          CASE N_PARAMS() OF
            4:domdef, min(glamt[neig1]), max(glamt[neig2]), min(gphit[neig1]), max(gphit[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
            6:domdef, min(glamt[neig1]), max(glamt[neig2]), min(gphit[neig1]), max(gphit[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
          ENDCASE
          RETURN
        ENDIF ELSE BEGIN
          ras = report('WARNING, The box does not contain any T points...')
          firstxt = -1 & lastxt = -1 & nxt = 0
          firstyt = -1 & lastyt = -1 & nyt = 0
        ENDELSE
      ENDIF ELSE BEGIN
        jyt = dom / jpi
        ixt = temporary(dom) MOD jpi
        firstxt = min(temporary(ixt), max = lastxt)
        firstyt = min(temporary(jyt), max = lastyt)
        nxt = lastxt - firstxt + 1
        nyt = lastyt - firstyt + 1
      ENDELSE
    ENDIF
; find firstxu, firstxu, firstyu, firstyu, nxu and nyu
; according to lon1, lon2, lat1 and lat2
    IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
      IF keyword_set(memeindices) THEN BEGIN
        firstxu = firstxt & lastxu = lastxt & nxu = nxt
        firstyu = firstyt & lastyu = lastyt & nyu = nyt
      ENDIF ELSE BEGIN
        dom = where( (glamu GE lon1) AND (glamu LE lon2) $
                     AND (gphiu GE lat1) AND (gphiu LE lat2) )
        IF (dom[0] EQ -1) THEN BEGIN
          IF keyword_set(findalways) THEN BEGIN
; if t grid parameters already defined, we use them...
            CASE 1 OF
              (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
                ras = report('WARNING, empty U points box... we use the same index as T points...')
                firstxu = firstxt & lastxu = lastxt & nxu = nxt
                firstyu = firstyt & lastyu = lastyt & nyu = nyt
              END
              ELSE:BEGIN
                ras = report('WARNING, empty U points box... we get the neighbors to define a new box...')
                neig1 = neighbor(lon1, lat1, glamu, gphiu, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                neig2 = neighbor(lon2, lat2, glamu, gphiu, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                CASE N_PARAMS() OF
                  4:domdef, min(glamu[neig1]), max(glamu[neig2]), min(gphiu[neig1]), max(gphiu[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
                  6:domdef, min(glamu[neig1]), max(glamu[neig2]), min(gphiu[neig1]), max(gphiu[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
                ENDCASE
                RETURN
              END
            ENDCASE
          ENDIF ELSE BEGIN
            ras = report('WARNING, The box does not contain any U points...')
            firstxu = -1 & lastxu = -1 & nxu = 0
            firstyu = -1 & lastyu = -1 & nyu = 0
          ENDELSE
        ENDIF ELSE BEGIN
          jyu = dom / jpi
          ixu = temporary(dom) MOD jpi
          firstxu = min(temporary(ixu), max = lastxu)
          firstyu = min(temporary(jyu), max = lastyu)
          nxu = lastxu - firstxu + 1
          nyu = lastyu - firstyu + 1
        ENDELSE
      ENDELSE
    ENDIF
; find firstxv, firstxv, firstyv, firstyv, nxv and nyv
;  according to lon1, lon2, lat1 and lat2
    IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
      IF keyword_set(memeindices) THEN BEGIN
        firstxv = firstxt & lastxv = lastxt & nxv = nxt
        firstyv = firstyt & lastyv = lastyt & nyv = nyt
      ENDIF ELSE BEGIN
        dom = where( (glamv GE lon1) AND (glamv LE lon2) $
                     AND (gphiv GE lat1) AND (gphiv LE lat2) )
        IF (dom[0] EQ -1) THEN BEGIN
          IF keyword_set(findalways) THEN BEGIN
            CASE 1 OF
              (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
                ras = report('WARNING, empty V points box... we use the same index as T points...')
                firstxv = firstxt & lastxv = lastxt & nxv = nxt
                firstyv = firstyt & lastyv = lastyt & nyv = nyt
              END
              (where(gdtype eq 'U'))[0] NE -1:BEGIN
                ras = report('WARNING, empty V points box... we use the same index as U points...')
                firstxv = firstxu & lastxv = lastxu & nxv = nxu
                firstyv = firstyu & lastyv = lastyu & nyv = nyu
              END
              ELSE:BEGIN
                ras = report('WARNING, empty V points box... we get the neighbors to define a new box...')
                neig1 = neighbor(lon1, lat1, glamv, gphiv, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                neig2 = neighbor(lon2, lat2, glamv, gphiv, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                CASE N_PARAMS() OF
                  4:domdef, min(glamv[neig1]), max(glamv[neig2]), min(gphiv[neig1]), max(gphiv[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
                  6:domdef, min(glamv[neig1]), max(glamv[neig2]), min(gphiv[neig1]), max(gphiv[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
                ENDCASE
                RETURN
              END
            ENDCASE
          ENDIF ELSE BEGIN
            ras = report('WARNING, The box does not contain any V points...')
            firstxv = -1 & lastxv = -1 & nxv = 0
            firstyv = -1 & lastyv = -1 & nyv = 0
          ENDELSE
        ENDIF ELSE BEGIN
          jyv = dom / jpi
          ixv = temporary(dom) MOD jpi
          firstxv = min(temporary(ixv), max = lastxv)
          firstyv = min(temporary(jyv), max = lastyv)
          nxv = lastxv - firstxv + 1
          nyv = lastyv - firstyv + 1
        ENDELSE
      ENDELSE
    ENDIF
; find firstxf, firstxf, firstyf, firstyf, nxf and nyf
; according to lon1, lon2, lat1 and lat2
    IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
      IF keyword_set(memeindices) THEN BEGIN
        firstxf = firstxt & lastxf = lastxt & nxf = nxt
        firstyf = firstyt & lastyf = lastyt & nyf = nyt
      ENDIF ELSE BEGIN
        dom = where( (glamf GE lon1) AND (glamf LE lon2) $
                     AND (gphif GE lat1) AND (gphif LE lat2) )
        IF (dom[0] EQ -1) THEN BEGIN
          IF keyword_set(findalways) THEN BEGIN
            CASE 1 OF
              (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
                ras = report('WARNING, empty F points box... we use the same index as T points...')
                firstxf = firstxt & lastxf = lastxt & nxf = nxt
                firstyf = firstyt & lastyf = lastyt & nyf = nyt
              END
              (where(gdtype eq 'U'))[0] NE -1:BEGIN
                ras = report('WARNING, empty F points box... we use the same index as U points...')
                firstxf = firstxu & lastxf = lastxu & nxf = nxu
                firstyf = firstyu & lastyf = lastyu & nyf = nyu
              END
              (where(gdtype eq 'V'))[0] NE -1:BEGIN
                ras = report('WARNING, empty F points box... we use the same index as V points...')
                firstxf = firstxv & lastxf = lastxv & nxf = nxv
                firstyf = firstyv & lastyf = lastyv & nyf = nyv
              END
              ELSE:BEGIN
                ras = report('WARNING, empty F points box... we get the neighbors to define a new box...')
                neig1 = neighbor(lon1, lat1, glamf, gphif, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                neig2 = neighbor(lon2, lat2, glamf, gphif, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                CASE N_PARAMS() OF
                  4:domdef, min(glamf[neig1]), max(glamf[neig2]), min(gphif[neig1]), max(gphif[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
                  6:domdef, min(glamf[neig1]), max(glamf[neig2]), min(gphif[neig1]), max(gphif[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
                ENDCASE
                RETURN
              END
            ENDCASE
          ENDIF ELSE BEGIN
            ras = report('WARNING, The box does not contain any F points...')
            firstxf = -1 & lastxf = -1 & nxf = 0
            firstyf = -1 & lastyf = -1 & nyf = 0
          ENDELSE
        ENDIF ELSE BEGIN
          jyf = dom / jpi
          ixf = temporary(dom) MOD jpi
          firstxf = min(temporary(ixf), max = lastxf)
          firstyf = min(temporary(jyf), max = lastyf)
          nxf = lastxf - firstxf + 1
          nyf = lastyf - firstyf + 1
        ENDELSE
      ENDELSE
    ENDIF
;
  ENDIF ELSE BEGIN
    CASE 1 OF
;-------------------------------------------------------------------
;-------------------------------------------------------------------
; horizontal domain defined with the X and Y indexes
;-------------------------------------------------------------------
;-------------------------------------------------------------------
      (keyword_set(xindex) AND keyword_set(yindex)) OR keyword_set(index):BEGIN
        fstx = min(long([x1, x2]), max = lstx)
        fsty = min(long([y1, y2]), max = lsty)
        IF fstx LT 0 OR lstx GE jpi THEN BEGIN
          ras = report('Bad definition of X1 or X2')
          return
        ENDIF
        IF fsty LT 0 OR lsty GE jpj THEN BEGIN
          ras = report('Bad definition of Y1 or Y2')
          return
        ENDIF
        nx = lstx - fstx + 1
        ny = lsty - fsty + 1
; find lon1t, lon2t, lat1t, lat2t, firstxt, firstxt, nxt and nyt
; according to x1, x2, y1, y2
        IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype eq 'W'))[0] NE -1 THEN BEGIN
          lon1t = min(glamt[fstx:lstx, fsty:lsty], max = lon2t)
          lat1t = min(gphit[fstx:lstx, fsty:lsty], max = lat2t)
          firstxt = fstx & lastxt = lstx
          firstyt = fsty & lastyt = lsty
          nxt = nx & nyt = ny
        ENDIF
; find lon1u, lon2u, lat1u, lat2u, firstxu, firstxu, nxu and nyu
; according to x1, x2, y1, y2
        IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
          lon1u = min(glamu[fstx:lstx, fsty:lsty], max = lon2u)
          lat1u = min(gphiu[fstx:lstx, fsty:lsty], max = lat2u)
          firstxu = fstx & lastxu = lstx
          firstyu = fsty & lastyu = lsty
          nxu = nx & nyu = ny
        ENDIF
; find lon1v, lon2v, lat1v, lat2v, firstxv, firstxv, nxv and nyv
; according to x1, x2, y1, y2
        IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
          lon1v = min(glamv[fstx:lstx, fsty:lsty], max = lon2v)
          lat1v = min(gphiv[fstx:lstx, fsty:lsty], max = lat2v)
          firstxv = fstx & lastxv = lstx
          firstyv = fsty & lastyv = lsty
          nxv = nx & nyv = ny
        ENDIF
; find lon1f, lon2f, lat1f, lat2f, firstxf, firstxf, nxf and nyf
; according to x1, x2, y1, y2
        IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
          lon1f = min(glamf[fstx:lstx, fsty:lsty], max = lon2f)
          lat1f = min(gphif[fstx:lstx, fsty:lsty], max = lat2f)
          firstxf = fstx & lastxf = lstx
          firstyf = fsty & lastyf = lsty
          nxf = nx & nyf = ny
        ENDIF
        lon1 = min([lon1t, lon1u, lon1v, lon1f])
        lon2 = max([lon2t, lon2u, lon2v, lon2f])
        lat1 = min([lat1t, lat1u, lat1v, lat1f])
        lat2 = max([lat2t, lat2u, lat2v, lat2f])
      END
;-------------------------------------------------------------------
;-------------------------------------------------------------------
; horizontal domain defined with the X index and lat1/lat2
;-------------------------------------------------------------------
;-------------------------------------------------------------------
      keyword_set(xindex):BEGIN
        fstx = min(long([x1, x2]), max = lstx)
        IF fstx LT 0 OR lstx GE jpi THEN BEGIN
          ras = report('Bad definition of X1 or X2')
          return
        ENDIF
        nx = lstx - fstx + 1
        lat1 = min([y1, y2], max = lat2)
; find lon1t, lon2t, firstxt, firstxt, firstyt, firstyt, nxt
; and nyt according to x1, x2, lat1 and lat2
        IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1 THEN BEGIN
          firstxt = fstx & lastxt = lstx & nxt = nx
          dom = where( (gphit[fstx:lstx, *] GE lat1) AND (gphit[fstx:lstx, *] LE lat2) )
          IF (dom[0] EQ -1) THEN BEGIN
            IF keyword_set(findalways) THEN BEGIN
              ras = report('WARNING, empty T points box... we get the neighbors to define a new box...')
              neig1 = neighbor(lon1, lat1, glamt[fstx:lstx, *], gphit[fstx:lstx, *], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
              neig2 = neighbor(lon2, lat2, glamt[fstx:lstx, *], gphit[fstx:lstx, *], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
              CASE N_PARAMS() OF
                4:domdef, fstx, lstx, min((gphit[fstx:lstx, *])[neig1]), max((gphit[fstx:lstx, *])[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
                6:domdef, fstx, lstx, min((gphit[fstx:lstx, *])[neig1]), max((gphit[fstx:lstx, *])[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
              ENDCASE
              RETURN
            ENDIF ELSE BEGIN
              ras = report('WARNING, The box does not contain any T points...')
              firstyt = -1 & lastyt = -1 & nyt = 0
            ENDELSE
          ENDIF ELSE BEGIN
            jyt = temporary(dom) / nx
            firstyt = min(temporary(jyt), max = lastyt)
            nyt = lastyt - firstyt + 1
          ENDELSE
          IF nyt NE 0 THEN lon1t = min(glamt[firstxt:lastxt, firstyt:lastyt], max = lon2t)
        ENDIF
; find lon1u, lon2u, firstxu, firstxu, firstyu, firstyu, nxu
; and nyu according to x1, x2, lat1 and lat2
        IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
          firstxu = fstx & lastxu = lstx & nxu = nx
          IF keyword_set(memeindices) THEN BEGIN
            firstyu = firstyt & lastyu = lastyt & nyu = nyt
          ENDIF ELSE BEGIN
            dom = where( (gphiu[fstx:lstx, *] GE lat1) AND (gphiu[fstx:lstx, *] LE lat2) )
            IF (dom[0] EQ -1) THEN BEGIN
              IF keyword_set(findalways) THEN BEGIN
                CASE 1 OF
                  (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
                    ras = report( 'WARNING, empty U points box... we use the same index as T points...')
                    firstyu = firstyt & lastyu = lastyt & nyu = nyt
                  END
                  ELSE:BEGIN
                    ras = report('WARNING, empty U points box... we get the neighbors to define a new box...')
                    neig1 = neighbor(lon1, lat1, glamu[fstx:lstx, *], gphiu[fstx:lstx, *], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    neig2 = neighbor(lon2, lat2, glamu[fstx:lstx, *], gphiu[fstx:lstx, *], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    CASE N_PARAMS() OF
                      4:domdef, fstx, lstx, min((gphiu[fstx:lstx, *])[neig1]), max((gphiu[fstx:lstx, *])[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
                      6:domdef, fstx, lstx, min((gphiu[fstx:lstx, *])[neig1]), max((gphiu[fstx:lstx, *])[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
                    ENDCASE
                    RETURN
                  END
                ENDCASE
              ENDIF ELSE BEGIN
                ras = report('WARNING, The box does not contain any U points...')
                firstyu = -1 & lastyu = -1 & nyu = 0
              ENDELSE
            ENDIF ELSE BEGIN
              jyu = temporary(dom) / nx
              firstyu = min(temporary(jyu), max = lastyu)
              nyu = lastyu - firstyu + 1
            ENDELSE
          ENDELSE
          IF nyu NE 0 THEN lon1u = min(glamu[firstxu:lastxu, firstyu:lastyu], max = lon2u)
        ENDIF
; find lon1v, lon2v, firstxv, firstxv, firstyv, firstyv, nxv
; and nyv according to x1, x2, lat1 and lat2
        IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
          firstxv = fstx & lastxv = lstx & nxv = nx
          IF keyword_set(memeindices) THEN BEGIN
            firstyv = firstyt & lastyv = lastyt & nyv = nyt
          ENDIF ELSE BEGIN
            dom = where( (gphiv[fstx:lstx, *] GE lat1) AND (gphiv[fstx:lstx, *] LE lat2) )
            IF (dom[0] EQ -1) THEN BEGIN
              IF keyword_set(findalways) THEN BEGIN
                CASE 1 OF
                  (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty V points box... we use the same index as T points...')
                    firstyv = firstyt & lastyv = lastyt & nyv = nyt
                  END
                  (where(gdtype eq 'U'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty V points box... we use the same index as U points...')
                    firstyv = firstyu & lastyv = lastyu & nyv = nyu
                  END
                  ELSE:BEGIN
                    ras = report('WARNING, empty V points box... we get the neighbors to define a new box...')
                    neig1 = neighbor(lon1, lat1, glamv[fstx:lstx, *], gphiv[fstx:lstx, *], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    neig2 = neighbor(lon2, lat2, glamv[fstx:lstx, *], gphiv[fstx:lstx, *], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    CASE N_PARAMS() OF
                      4:domdef, fstx, lstx, min((gphiv[fstx:lstx, *])[neig1]), max((gphiv[fstx:lstx, *])[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
                      6:domdef, fstx, lstx, min((gphiv[fstx:lstx, *])[neig1]), max((gphiv[fstx:lstx, *])[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
                    ENDCASE
                    RETURN
                  END
                ENDCASE
              ENDIF ELSE BEGIN
                ras = report('WARNING, The box does not contain any V points...')
                firstyv = -1 & lastyv = -1 & nyv = 0
              ENDELSE
            ENDIF ELSE BEGIN
              jyv = temporary(dom) / nx
              firstyv = min(temporary(jyv), max = lastyv)
              nyv = lastyv - firstyv + 1
            ENDELSE
          ENDELSE
          IF nyv NE 0 THEN lon1v = min(glamv[firstxv:lastxv, firstyv:lastyv], max = lon2v)
        ENDIF
; find lon1f, lon2f, firstxf, firstxf, firstyf, firstyf, nxf
; and nyf according to x1, x2, lat1 and lat2
        IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
          firstxf = fstx & lastxf = lstx & nxf = nx
          IF keyword_set(memeindices) THEN BEGIN
            firstyf = firstyt & lastyf = lastyt & nyf = nyt
          ENDIF ELSE BEGIN
            dom = where( (gphif[fstx:lstx, *] GE lat1) AND (gphif[fstx:lstx, *] LE lat2) )
            IF (dom[0] EQ -1) THEN BEGIN
              IF keyword_set(findalways) THEN BEGIN
                CASE 1 OF
                  (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty F points box... we use the same index as T points...')
                    firstyf = firstyt & lastyf = lastyt & nyf = nyt
                  END
                  (where(gdtype eq 'U'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty F points box... we use the same index as U points...')
                    firstyf = firstyu & lastyf = lastyu & nyf = nyu
                  END
                  (where(gdtype eq 'V'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty F points box... we use the same index as V points...')
                    firstyf = firstyv & lastyf = lastyv & nyf = nyv
                  END
                  ELSE:BEGIN
                    ras = report('WARNING, empty F points box... we get the neighbors to define a new box...')
                    neig1 = neighbor(lon1, lat1, glamf[fstx:lstx, *], gphif[fstx:lstx, *], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    neig2 = neighbor(lon2, lat2, glamf[fstx:lstx, *], gphif[fstx:lstx, *], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    CASE N_PARAMS() OF
                      4:domdef, fstx, lstx, min((gphif[fstx:lstx, *])[neig1]), max((gphif[fstx:lstx, *])[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
                      6:domdef, fstx, lstx, min((gphif[fstx:lstx, *])[neig1]), max((gphif[fstx:lstx, *])[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
                    ENDCASE
                    RETURN
                  END
                ENDCASE
              ENDIF ELSE BEGIN
                ras = report('WARNING, The box does not contain any F points...')
                firstyf = -1 & lastyf = -1 & nyf = 0
              ENDELSE
            ENDIF ELSE BEGIN
              jyf = temporary(dom) / nx
              firstyf = min(temporary(jyf), max = lastyf)
              nyf = lastyf - firstyf + 1
            ENDELSE
          ENDELSE
          IF nyf NE 0 THEN lon1f = min(glamf[firstxf:lastxf, firstyf:lastyf], max = lon2f)
        ENDIF
        lon1 = min([lon1t, lon1u, lon1v, lon1f])
        lon2 = max([lon2t, lon2u, lon2v, lon2f])
      END
;-------------------------------------------------------------------
;-------------------------------------------------------------------
; horizontal domain defined with the Y index and lon1/lon2
;-------------------------------------------------------------------
;-------------------------------------------------------------------
      keyword_set(yindex):BEGIN
        fsty = min(long([y1, y2]), max = lsty)
        IF fsty LT 0 OR lsty GE jpj THEN BEGIN
          ras = report('Bad definition of Y1 or Y2')
          return
        ENDIF
        ny = lsty - fsty + 1
        lon1 = min([x1, x2], max = lon2)
; find lat1t, lat2t, firstxt, firstxt, firstyt, firstyt, nxt
; and nyt according to x1, x2, lon1 and lon2
        IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1 THEN BEGIN
          firstyt = fsty & lastyt = lsty & nyt = ny
          dom = where( (glamt[*, fsty:lsty] GE lon1) AND (glamt[*, fsty:lsty] LE lon2) )
          IF (dom[0] EQ -1) THEN BEGIN
            IF keyword_set(findalways) THEN BEGIN
              ras = report('WARNING, empty T points box... we get the neighbors to define a new box...')
              neig1 = neighbor(lon1, lat1, glamt[*, fsty:lsty], gphit[*, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
              neig2 = neighbor(lon2, lat2, glamt[*, fsty:lsty], gphit[*, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
              CASE N_PARAMS() OF
                4:domdef, min((glamt[*, fsty:lsty])[neig1]), max((glamt[*, fsty:lsty])[neig2]), fsty, lsty, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
                6:domdef, min((glamt[*, fsty:lsty])[neig1]), max((glamt[*, fsty:lsty])[neig2]), fsty, lsty, z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
              ENDCASE
              RETURN
            ENDIF ELSE BEGIN
              ras = report('WARNING, The box does not contain any T points...')
              firstxt = -1 & lastxt = -1 & nxt = 0
            ENDELSE
          ENDIF ELSE BEGIN
            jxt = temporary(dom) MOD jpi
            firstxt = min(temporary(jxt), max = lastxt)
            nxt = lastxt - firstxt + 1
          ENDELSE
          IF nxt NE 0 THEN lat1t = min(gphit[firstxt:lastxt, firstyt:lastyt], max = lat2t)
        ENDIF
; find lat1u, lat2u, firstxu, firstxu, firstyu, firstyu, nxu
; and nyu according to x1, x2, lon1 and lon2
        IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
          firstyu = fsty & lastyu = lsty & nyu = ny
          IF keyword_set(memeindices) THEN BEGIN
            firstxu = firstxt & lastxu = lastxt & nxu = nxt
          ENDIF ELSE BEGIN
            dom = where( (glamu[*, fsty:lsty] GE lon1) AND (glamu[*, fsty:lsty] LE lon2) )
            IF (dom[0] EQ -1) THEN BEGIN
              IF keyword_set(findalways) THEN BEGIN
                CASE 1 OF
                  (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty U points box... we use the same index as T points...')
                    firstxu = firstxt & lastxu = lastxt & nxu = nxt
                  END
                  ELSE:BEGIN
                    ras = report('WARNING, empty U points box... we get the neighbors to define a new box...')
                    neig1 = neighbor(lon1, lat1, glamu[*, fsty:lsty], gphiu[*, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    neig2 = neighbor(lon2, lat2, glamu[*, fsty:lsty], gphiu[*, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    CASE N_PARAMS() OF
                      4:domdef, min((glamu[*, fsty:lsty])[neig1]), max((glamu[*, fsty:lsty])[neig2]), fsty, lsty, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
                      6:domdef, min((glamu[*, fsty:lsty])[neig1]), max((glamu[*, fsty:lsty])[neig2]), fsty, lsty, z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
                    ENDCASE
                    RETURN
                  END
                ENDCASE
              ENDIF ELSE BEGIN
                ras = report('WARNING, The box does not contain any U points...')
                firstxu = -1 & lastxu = -1 & nxu = 0
              ENDELSE
            ENDIF ELSE BEGIN
              jxu = temporary(dom) MOD jpi
              firstxu = min(temporary(jxu), max = lastxu)
              nxu = lastxu - firstxu + 1
            ENDELSE
          ENDELSE
          IF nxu NE 0 THEN lat1u = min(gphiu[firstxu:lastxu, firstyu:lastyu], max = lat2u)
        ENDIF
; find lat1v, lat2v, firstxv, firstxv, firstyv, firstyv, nxv
; and nyv according to x1, x2, lon1 and lon2
        IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
          firstyv = fsty & lastyv = lsty & nyv = ny
          IF keyword_set(memeindices) THEN BEGIN
            firstxv = firstxt & lastxv = lastxt & nxv = nxt
          ENDIF ELSE BEGIN
            dom = where( (glamv[*, fsty:lsty] GE lon1) AND (glamv[*, fsty:lsty] LE lon2) )
            IF (dom[0] EQ -1) THEN BEGIN
              IF keyword_set(findalways) THEN BEGIN
                CASE 1 OF
                  (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty V points box... we use the same index as T points...')
                    firstxv = firstxt & lastxv = lastxt & nxv = nxt
                  END
                  (where(gdtype eq 'U'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty V points box... we use the same index as U points...')
                    firstxv = firstxu & lastxv = lastxu & nxv = nxu
                  END
                  ELSE:BEGIN
                    ras = report('WARNING, empty V points box... we get the neighbors to define a new box...')
                    neig1 = neighbor(lon1, lat1, glamv[*, fsty:lsty], gphiv[*, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    neig2 = neighbor(lon2, lat2, glamv[*, fsty:lsty], gphiv[*, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    CASE N_PARAMS() OF
                      4:domdef, min((glamv[*, fsty:lsty])[neig1]), max((glamv[*, fsty:lsty])[neig2]), fsty, lsty, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
                      6:domdef, min((glamv[*, fsty:lsty])[neig1]), max((glamv[*, fsty:lsty])[neig2]), fsty, lsty, z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
                    ENDCASE
                    RETURN
                  END
                ENDCASE
              ENDIF ELSE BEGIN
                ras = report('WARNING, The box does not contain any V points...')
                firstxv = -1 & lastxv = -1 & nxv = 0
              ENDELSE
            ENDIF ELSE BEGIN
              jxv = temporary(dom) MOD jpi
              firstxv = min(temporary(jxv), max = lastxv)
              nxv = lastxv - firstxv + 1
            ENDELSE
          ENDELSE
          IF nxv NE 0 THEN lat1v = min(gphiv[firstxv:lastxv, firstyv:lastyv], max = lat2v)
        ENDIF
; find lat1f, lat2f, firstxf, firstxf, firstyf, firstyf, nxf
; and nyf according to x1, x2, lon1 and lon2
        IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
          firstyf = fsty & lastyf = lsty & nyf = ny
          IF keyword_set(memeindices) THEN BEGIN
            firstxf = firstxt & lastxf = lastxt & nxf = nxt
          ENDIF ELSE BEGIN
            dom = where( (glamf[*, fsty:lsty] GE lon1) AND (glamf[*, fsty:lsty] LE lon2) )
            IF (dom[0] EQ -1) THEN BEGIN
              IF keyword_set(findalways) THEN BEGIN
                CASE 1 OF
                  (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty F points box... we use the same index as T points...')
                    firstxf = firstxt & lastxf = lastxt & nxf = nxt
                  END
                  (where(gdtype eq 'U'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty F points box... we use the same index as U points...')
                    firstxf = firstxu & lastxf = lastxu & nxf = nxu
                  END
                  (where(gdtype eq 'V'))[0] NE -1:BEGIN
                    ras = report('WARNING, empty F points box... we use the same index as V points...')
                    firstxf = firstxv & lastxf = lastxv & nxf = nxv
                  END
                  ELSE:BEGIN
                    ras = report('WARNING, empty F points box... we get the neighbors to define a new box...')
                    neig1 = neighbor(lon1, lat1, glamf[*, fsty:lsty], gphif[*, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    neig2 = neighbor(lon2, lat2, glamf[*, fsty:lsty], gphif[*, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
                    CASE N_PARAMS() OF
                      4:domdef, min((glamf[*, fsty:lsty])[neig1]), max((glamf[*, fsty:lsty])[neig2]), fsty, lsty, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
                      6:domdef, min((glamf[*, fsty:lsty])[neig1]), max((glamf[*, fsty:lsty])[neig2]), fsty, lsty, z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
                    ENDCASE
                    RETURN
                  END
                ENDCASE
              ENDIF ELSE BEGIN
                ras = report('WARNING, The box does not contain any F points...')
                firstxf = -1 & lastyf = -1 & nxf = 0
              ENDELSE
            ENDIF ELSE BEGIN
              jxf = temporary(dom) MOD jpi
              firstxf = min(temporary(jxf), max = lastxf)
              nxf = lastxf - firstxf + 1
            ENDELSE
          ENDELSE
          IF nxf NE 0 THEN lat1f = min(gphif[firstxf:lastxf, firstyf:lastyf], max = lat2f)
        ENDIF
        lat1 = min([lat1t, lat1u, lat1v, lat1f])
        lat2 = max([lat2t, lat2u, lat2v, lat2f])
      END
    ENDCASE
  ENDELSE
;-------------------------------------------------------------------
; The extracted domain is it regular or not?
;-------------------------------------------------------------------
  CASE 1 OF
    ((where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype eq 'W'))[0] NE -1) AND nxt NE 0 AND nyt NE 0:BEGIN
; to get faster, we first test the most basic cases before testing the
; full array.
      CASE 0 OF
        array_equal(glamt[firstxt:lastxt, firstyt], glamt[firstxt:lastxt, lastyt]):key_irregular = 1b
        array_equal(gphit[firstxt, firstyt:lastyt], gphit[lastxt, firstyt:lastyt]):key_irregular = 1b
        array_equal(glamt[firstxt:lastxt, firstyt:lastyt] $
                    , glamt[firstxt:lastxt, firstyt]#replicate(1, nyt)):key_irregular = 1b
        array_equal(gphit[firstxt:lastxt, firstyt:lastyt] $
                    , replicate(1, nxt)#(gphit[firstxt, firstyt:lastyt])[*]):key_irregular = 1b
        ELSE:key_irregular = 0b
      ENDCASE
    END
    (where(gdtype eq 'U'))[0] NE -1 AND nxu NE 0 AND nyu NE 0:BEGIN
      CASE 0 OF
        array_equal(glamu[firstxu:lastxu, firstyu], glamu[firstxu:lastxu, lastyu]):key_irregular = 1b
        array_equal(gphiu[firstxu, firstyu:lastyu], gphiu[lastxu, firstyu:lastyu]):key_irregular = 1b
        array_equal(glamu[firstxu:lastxu, firstyu:lastyu] $
                    , glamu[firstxu:lastxu, firstyu]#replicate(1, nyu)):key_irregular = 1b
        array_equal(gphiu[firstxu:lastxu, firstyu:lastyu] $
                    , replicate(1, nxu)#(gphiu[firstxu, firstyu:lastyu])[*]):key_irregular = 1b
        ELSE:key_irregular = 0b
      ENDCASE
    END
    (where(gdtype eq 'V'))[0] NE -1 AND nxv NE 0 AND nyv NE 0:BEGIN
      CASE 0 OF
        array_equal(glamv[firstxv:lastxv, firstyv], glamv[firstxv:lastxv, lastyv]):key_irregular = 1b
        array_equal(gphiv[firstxv, firstyv:lastyv], gphiv[lastxv, firstyv:lastyv]):key_irregular = 1b
        array_equal(glamv[firstxv:lastxv, firstyv:lastyv] $
                    , glamv[firstxv:lastxv, firstyv]#replicate(1, nyv)):key_irregular = 1b
        array_equal(gphiv[firstxv:lastxv, firstyv:lastyv] $
                    , replicate(1, nxv)#(gphiv[firstxv, firstyv:lastyv])[*]):key_irregular = 1b
        ELSE:key_irregular = 0b
      ENDCASE
    END
    (where(gdtype eq 'F'))[0] NE -1 AND nxf NE 0 AND nyf NE 0:BEGIN
      CASE 0 OF
        array_equal(glamf[firstxf:lastxf, firstyf], glamf[firstxf:lastxf, lastyf]):key_irregular = 1b
        array_equal(gphif[firstxf, firstyf:lastyf], gphif[lastxf, firstyf:lastyf]):key_irregular = 1b
        array_equal(glamf[firstxf:lastxf, firstyf:lastyf] $
                    , glamf[firstxf:lastxf, firstyf]#replicate(1, nyf)):key_irregular = 1b
        array_equal(gphif[firstxf:lastxf, firstyf:lastyf] $
                    , replicate(1, nxf)#(gphif[firstxf, firstyf:lastyf])[*]):key_irregular = 1b
        ELSE:key_irregular = 0b
      ENDCASE
    END
    ELSE:
  ENDCASE
;
;-------------------------------------------------------------------
;-------------------------------------------------------------------
; define all vertical parameters ...
; vert1, vert2
; firstz[tw], lastz[tw], nz[tw]
;-------------------------------------------------------------------
;-------------------------------------------------------------------
;
  vertical:
;
;-------------------------------------------------------------------
; vertical domain defined with vert1, vert2
;-------------------------------------------------------------------
  IF NOT (keyword_set(zindex) OR keyword_set(index)) THEN BEGIN
; define vert1 and vert2
    CASE N_PARAMS() OF
      2:vert1 = min([x1, x2], max = vert2)
      6:vert1 = min([z1, z2], max = vert2)
      ELSE:BEGIN
        IF (inter(byte(gdtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN $
          vert1t = min(gdept, max = vert2t)
        IF (where(gdtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN $
          vert1w = min(gdepw, max = vert2w)
        vert1 = min([vert1t, vert1w])
        vert2 = max([vert2t, vert2w])
      END
    ENDCASE
; define firstzt, firstzt, nzt
    IF (inter(byte(gdtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN
      domz = where(gdept ge vert1 and gdept le vert2, nzt)
      IF nzt NE 0 THEN BEGIN
        firstzt = domz[0]
        lastzt = domz[nzt-1]
      ENDIF ELSE BEGIN
        ras = report('WARNING, The box does not contain any T level...')
        firstzt = -1
        lastzt = -1
      ENDELSE
    ENDIF
; define firstzw, firstzw, nzw
    IF (where(gdtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
      IF keyword_set(memeindices) THEN BEGIN
        firstzw = firstzt
        lastzw = lastzt
        nzw = nzt
      ENDIF ELSE BEGIN
        domz = where(gdepw ge vert1 and gdepw le vert2, nzw)
        IF nzw NE 0 THEN BEGIN
          firstzw = domz[0]
          lastzw = domz[nzw-1]
        ENDIF ELSE BEGIN
          ras = report('WARNING, The box does not contain any W level...')
          firstzw = -1
          lastzw = -1
        ENDELSE
      ENDELSE
    ENDIF
;-------------------------------------------------------------------
; vertical domain defined with the Z index
;-------------------------------------------------------------------
  ENDIF ELSE BEGIN
    CASE N_PARAMS() OF
      2:fstz = min(long([x1, x2]), max = lstz)
      4:return
      6:fstz = min(long([z1, z2]), max = lstz)
    ENDCASE
    IF fstz LT 0 OR lstz GE jpk THEN BEGIN
      ras = report('Bad definition of X1, X2, Z1 or Z2')
      return
    ENDIF
    nz = lstz - fstz + 1
; find vert1t, vert2t, firstzt, firstzt, nzt
; according to (x1, x2) or (z1, z2)
    IF (where(gdtype eq 'T'))[0] NE -1 THEN BEGIN
      vert1t = min(gdept[fstz:lstz], max = vert2t)
      firstzt = fstz & lastzt = lstz & nzt = nz
    ENDIF
; find vert1w, vert2w, firstzw, firstzw, nzw
; according to (x1, x2) or (z1, z2)
    IF (where(gdtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
      vert1w = min(gdepw[fstz:lstz], max = vert2w)
      firstzw = fstz & lastzw = lstz & nzw = nz
    ENDIF
    vert1 = min([vert1t, vert1w])
    vert2 = max([vert2t, vert2w])
  ENDELSE
;
;-------------------------------------------------------------------
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updateold
  ENDIF
;-------------------------------------------------------------------
  if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun
;-------------------------------------------------------------------

;-------------------------------------------------------------------
  return
end