;+
;
; @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 domdef must define the box used to make a plot
; (pltz, pltt, plt1d) done strictly along the line (that can have any direction)
; starting at (x1, y1) ending at (x2, y2). When defining endpoints,
; you must also define TYPE which define the type of plots
; ('pltz', 'xt', 'yt', 'zt', 'x', 'y', 'z', 't') will used
; ENDPOINTS keywords
;
; @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 norme, curl...).
; Activate MEMEINDICES to force domdef 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
; domdef refer to indexes of glam, gphi and gdep arrays rather
; than to values of these arrays.
;
; @keyword XINDEX
; We activate it if we want that all elements passed in input of
; domdef
; 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
; domdef
; 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
; domdef
; and concerning the X dimension refer to indexes of gdep arrays rather
; than to values of these arrays.
;
; @uses
; common.pro
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; 8/2/98
; rewrite everything, debug and spee-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 gridtype = ['T', 'U', 'V', 'W', 'F'] $
ELSE gridtype = strupcase(gridtype)
IF keyword_set(memeindices) THEN gridtype = ['T', gridtype]
IF finite(glamu[0]) eq 0 THEN gridtype = gridtype[where(gridtype NE 'U')]
IF finite(glamv[0]) eq 0 THEN gridtype = gridtype[where(gridtype 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(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
glamt = reform(glamt, jpi, jpj, /over)
gphit = reform(gphit, jpi, jpj, /over)
ENDIF
IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
glamu = reform(glamu, jpi, jpj, /over)
gphiu = reform(gphiu, jpi, jpj, /over)
ENDIF
IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
glamv = reform(glamv, jpi, jpj, /over)
gphiv = reform(gphiv, jpi, jpj, /over)
ENDIF
IF (where(gridtype 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 et lon2
; lat1 et 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 gridtype eq 'U' or 'V'
;-------------------------------------------------------------------
errstatus = 0
IF (finite(glamu[0]*gphiu[0]) EQ 0 OR n_elements(glamu) EQ 0 OR n_elements(gphiu) EQ 0) AND (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
firstxu = !values.f_nan
lastxu = !values.f_nan
nxu = !values.f_nan
okgrid = where(gridtype NE 'U', count)
IF count NE 0 THEN gridtype = gridtype[okgrid] $
ELSE errstatus = report('U grid is undefined. Impossible to call domdef with vargid = ''U''')
ENDIF
;
IF (finite(glamv[0]*gphiv[0]) EQ 0 OR n_elements(glamv) EQ 0 OR n_elements(gphiv) EQ 0) AND (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
firstxv = !values.f_nan
lastxv = !values.f_nan
nxv = !values.f_nan
okgrid = where(gridtype NE 'V', count)
IF count NE 0 THEN gridtype = gridtype[okgrid] $
ELSE errstatus = report('V grid is undefined. Impossible to call domdef with vargid = ''V''')
ENDIF
IF errstatus EQ -1 THEN return
;-------------------------------------------------------------------
;-------------------------------------------------------------------
; horizontal domain defined with lon1, lon2, lat1 and lat2
;-------------------------------------------------------------------
;-------------------------------------------------------------------
IF N_PARAMS() EQ 0 $
OR ( (N_PARAMS() EQ 4 OR N_PARAMS() EQ 6) $
AND NOT keyword_set(xindex) AND NOT keyword_set(yindex) AND NOT keyword_set(index) ) THEN BEGIN
IF N_PARAMS() EQ 0 THEN BEGIN
; find lon1 and lon2 the longitudinal boundaries of the full domain
IF (where(gridtype eq 'T'))[0] NE -1 THEN lon1t = min(glamt, max = lon2t)
IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lon1t = min(glamt, max = lon2t)
IF (where(gridtype eq 'U'))[0] NE -1 THEN lon1u = min(glamu, max = lon2u)
IF (where(gridtype eq 'V'))[0] NE -1 THEN lon1v = min(glamv, max = lon2v)
IF (where(gridtype eq 'F'))[0] NE -1 THEN lon1f = min(glamf, max = lon2f)
lon1 = min([lon1t, lon1u, lon1v, lon1f])
lon2 = max([lon2t, lon2u, lon2v, lon2f])
; find lat1 and lat2 the latitudinal boundaries of the full domain
IF (where(gridtype eq 'T'))[0] NE -1 THEN lat1t = min(gphit, max = lat2t)
IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lat1t = min(gphit, max = lat2t)
IF (where(gridtype eq 'U'))[0] NE -1 THEN lat1u = min(gphiu, max = lat2u)
IF (where(gridtype eq 'V'))[0] NE -1 THEN lat1v = min(gphiv, max = lat2v)
IF (where(gridtype eq 'F'))[0] NE -1 THEN lat1f = min(gphif, max = lat2f)
lat1 = min([lat1t, lat1u, lat1v, lat1f])
lat2 = max([lat2t, lat2u, lat2v, lat2f])
ENDIF ELSE BEGIN
lon1 = min([x1, x2], max = lon2)
lat1 = min([y1, y2], max = lat2)
ENDELSE
; find firstxt, firstxt, nxt and nyt according to lon1, lon2, lat1 and lat2
IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
dom = where( (glamt GE lon1) AND (glamt LE lon2) $
AND (gphit GE lat1) AND (gphit LE lat2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
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]), GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
6:domdef, min(glamt[neig1]), max(glamt[neig2]), min(gphit[neig1]), max(gphit[neig2]), z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'U'))[0] NE -1 THEN BEGIN
IF keyword_set(memeindices) THEN BEGIN
firstxu = firstxt & lastxu = lastxt & nxu = nxt
firstyu = firstyt & lastyu = lastyt & nyu = nyt
ENDIF ELSE BEGIN
dom = where( (glamu GE lon1) AND (glamu LE lon2) $
AND (gphiu GE lat1) AND (gphiu LE lat2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
; if t grid parameters alreday defined, we use them...
CASE 1 OF
(where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype 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]), GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
6:domdef, min(glamu[neig1]), max(glamu[neig2]), min(gphiu[neig1]), max(gphiu[neig2]), z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'V'))[0] NE -1 THEN BEGIN
IF keyword_set(memeindices) THEN BEGIN
firstxv = firstxt & lastxv = lastxt & nxv = nxt
firstyv = firstyt & lastyv = lastyt & nyv = nyt
ENDIF ELSE BEGIN
dom = where( (glamv GE lon1) AND (glamv LE lon2) $
AND (gphiv GE lat1) AND (gphiv LE lat2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
CASE 1 OF
(where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype 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(gridtype 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]), GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
6:domdef, min(glamv[neig1]), max(glamv[neig2]), min(gphiv[neig1]), max(gphiv[neig2]), z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'F'))[0] NE -1 THEN BEGIN
IF keyword_set(memeindices) THEN BEGIN
firstxf = firstxt & lastxf = lastxt & nxf = nxt
firstyf = firstyt & lastyf = lastyt & nyf = nyt
ENDIF ELSE BEGIN
dom = where( (glamf GE lon1) AND (glamf LE lon2) $
AND (gphif GE lat1) AND (gphif LE lat2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
CASE 1 OF
(where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype 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(gridtype 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(gridtype 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]), GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
6:domdef, min(glamf[neig1]), max(glamf[neig2]), min(gphif[neig1]), max(gphif[neig2]), z1, z2, GRIDTYPE = gridtype, 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([x1, x2], max = lstx)
fsty = min([y1, y2], max = lsty)
IF fstx LT 0 OR lstx GE jpi THEN BEGIN
ras = report('Bad definition of X1 or X2')
return
ENDIF
IF fsty LT 0 OR lsty GE jpj THEN BEGIN
ras = report('Bad definition of Y1 or Y2')
return
ENDIF
nx = lstx - fstx + 1
ny = lsty - fsty + 1
; find lon1t, lon2t, lat1t, lat2t, firstxt, firstxt, nxt and nyt
; according to x1, x2, y1, y2
IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1 THEN BEGIN
lon1t = min(glamt[fstx:lstx, fsty:lsty], max = lon2t)
lat1t = min(gphit[fstx:lstx, fsty:lsty], max = lat2t)
firstxt = fstx & lastxt = lstx
firstyt = fsty & lastyt = lsty
nxt = nx & nyt = ny
ENDIF
; find lon1u, lon2u, lat1u, lat2u, firstxu, firstxu, nxu and nyu
; according to x1, x2, y1, y2
IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
lon1u = min(glamu[fstx:lstx, fsty:lsty], max = lon2u)
lat1u = min(gphiu[fstx:lstx, fsty:lsty], max = lat2u)
firstxu = fstx & lastxu = lstx
firstyu = fsty & lastyu = lsty
nxu = nx & nyu = ny
ENDIF
; find lon1v, lon2v, lat1v, lat2v, firstxv, firstxv, nxv and nyv
; according to x1, x2, y1, y2
IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
lon1v = min(glamv[fstx:lstx, fsty:lsty], max = lon2v)
lat1v = min(gphiv[fstx:lstx, fsty:lsty], max = lat2v)
firstxv = fstx & lastxv = lstx
firstyv = fsty & lastyv = lsty
nxv = nx & nyv = ny
ENDIF
; find lon1f, lon2f, lat1f, lat2f, firstxf, firstxf, nxf and nyf
; according to x1, x2, y1, y2
IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN
lon1f = min(glamf[fstx:lstx, fsty:lsty], max = lon2f)
lat1f = min(gphif[fstx:lstx, fsty:lsty], max = lat2f)
firstxf = fstx & lastxf = lstx
firstyf = fsty & lastyf = lsty
nxf = nx & nyf = ny
ENDIF
lon1 = min([lon1t, lon1u, lon1v, lon1f])
lon2 = max([lon2t, lon2u, lon2v, lon2f])
lat1 = min([lat1t, lat1u, lat1v, lat1f])
lat2 = max([lat2t, lat2u, lat2v, lat2f])
END
;-------------------------------------------------------------------
;-------------------------------------------------------------------
; horizontal domain defined with the X index and lat1/lat2
;-------------------------------------------------------------------
;-------------------------------------------------------------------
keyword_set(xindex):BEGIN
fstx = min([x1, x2], max = lstx)
IF fstx LT 0 OR lstx GE jpi THEN BEGIN
ras = report('Bad definition of X1 or X2')
return
ENDIF
nx = lstx - fstx + 1
lat1 = min([y1, y2], max = lat2)
; find lon1t, lon2t, firstxt, firstxt, firstyt, firstyt, nxt
; and nyt according to x1, x2, lat1 and lat2
IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
firstxt = fstx & lastxt = lstx & nxt = nx
dom = where( (gphit[fstx:lstx, *] GE lat1) AND (gphit[fstx:lstx, *] LE lat2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
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]), GRIDTYPE = gridtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
6:domdef, fstx, lstx, min((gphit[fstx:lstx, *])[neig1]), max((gphit[fstx:lstx, *])[neig2]), z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'U'))[0] NE -1 THEN BEGIN
firstxu = fstx & lastxu = lstx & nxu = nx
IF keyword_set(memeindices) THEN BEGIN
firstyu = firstyt & lastyu = lastyt & nyu = nyt
ENDIF ELSE BEGIN
dom = where( (gphiu[fstx:lstx, *] GE lat1) AND (gphiu[fstx:lstx, *] LE lat2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
CASE 1 OF
(where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype 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]), GRIDTYPE = gridtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
6:domdef, fstx, lstx, min((gphiu[fstx:lstx, *])[neig1]), max((gphiu[fstx:lstx, *])[neig2]), z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'V'))[0] NE -1 THEN BEGIN
firstxv = fstx & lastxv = lstx & nxv = nx
IF keyword_set(memeindices) THEN BEGIN
firstyv = firstyt & lastyv = lastyt & nyv = nyt
ENDIF ELSE BEGIN
dom = where( (gphiv[fstx:lstx, *] GE lat1) AND (gphiv[fstx:lstx, *] LE lat2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
CASE 1 OF
(where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype 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(gridtype 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]), GRIDTYPE = gridtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
6:domdef, fstx, lstx, min((gphiv[fstx:lstx, *])[neig1]), max((gphiv[fstx:lstx, *])[neig2]), z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'F'))[0] NE -1 THEN BEGIN
firstxf = fstx & lastxf = lstx & nxf = nx
IF keyword_set(memeindices) THEN BEGIN
firstyf = firstyt & lastyf = lastyt & nyf = nyt
ENDIF ELSE BEGIN
dom = where( (gphif[fstx:lstx, *] GE lat1) AND (gphif[fstx:lstx, *] LE lat2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
CASE 1 OF
(where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype 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(gridtype 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(gridtype 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]), GRIDTYPE = gridtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
6:domdef, fstx, lstx, min((gphif[fstx:lstx, *])[neig1]), max((gphif[fstx:lstx, *])[neig2]), z1, z2, GRIDTYPE = gridtype, 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([y1, y2], max = lsty)
IF fsty LT 0 OR lsty GE jpj THEN BEGIN
ras = report('Bad definition of Y1 or Y2')
return
ENDIF
ny = lsty - fsty + 1
lon1 = min([x1, x2], max = lon2)
; find lat1t, lat2t, firstxt, firstxt, firstyt, firstyt, nxt
; and nyt according to x1, x2, lon1 and lon2
IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
firstyt = fsty & lastyt = lsty & nyt = ny
dom = where( (glamt[*, fsty:lsty] GE lon1) AND (glamt[*, fsty:lsty] LE lon2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
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, GRIDTYPE = gridtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
6:domdef, min((glamt[*, fsty:lsty])[neig1]), max((glamt[*, fsty:lsty])[neig2]), fsty, lsty, z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'U'))[0] NE -1 THEN BEGIN
firstyu = fsty & lastyu = lsty & nyu = ny
IF keyword_set(memeindices) THEN BEGIN
firstxu = firstyt & lastxu = lastyt & nxu = nxt
ENDIF ELSE BEGIN
dom = where( (glamu[*, fsty:lsty] GE lon1) AND (glamu[*, fsty:lsty] LE lon2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
CASE 1 OF
(where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype 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, GRIDTYPE = gridtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
6:domdef, min((glamu[*, fsty:lsty])[neig1]), max((glamu[*, fsty:lsty])[neig2]), fsty, lsty, z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'V'))[0] NE -1 THEN BEGIN
firstyv = fsty & lastyv = lsty & nyv = ny
IF keyword_set(memeindices) THEN BEGIN
firstxv = firstyt & lastxv = lastyt & nxv = nxt
ENDIF ELSE BEGIN
dom = where( (glamv[*, fsty:lsty] GE lon1) AND (glamv[*, fsty:lsty] LE lon2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
CASE 1 OF
(where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype 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(gridtype 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, GRIDTYPE = gridtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
6:domdef, min((glamv[*, fsty:lsty])[neig1]), max((glamv[*, fsty:lsty])[neig2]), fsty, lsty, z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'F'))[0] NE -1 THEN BEGIN
firstyf = fsty & lastyf = lsty & nyf = ny
IF keyword_set(memeindices) THEN BEGIN
firstxf = firstyt & lastxf = lastyt & nxf = nxt
ENDIF ELSE BEGIN
dom = where( (glamf[*, fsty:lsty] GE lon1) AND (glamf[*, fsty:lsty] LE lon2) )
IF (dom[0] EQ -1) THEN BEGIN
IF keyword_set(findalways) THEN BEGIN
CASE 1 OF
(where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype 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(gridtype 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(gridtype 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, GRIDTYPE = gridtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
6:domdef, min((glamf[*, fsty:lsty])[neig1]), max((glamf[*, fsty:lsty])[neig2]), fsty, lsty, z1, z2, GRIDTYPE = gridtype, 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(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1) AND nxt NE 0 AND nyt NE 0:BEGIN
; to get faster, we first test the most basic cases 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(gridtype eq 'U'))[0] NE -1 AND nxu NE 0 AND nyu NE 0:BEGIN
CASE 0 OF
array_equal(glamu[firstxu:lastxu, firstyu], glamu[firstxu:lastxu, lastyu]):key_irregular = 1b
array_equal(gphiu[firstxu, firstyu:lastyu], gphiu[lastxu, firstyu:lastyu]):key_irregular = 1b
array_equal(glamu[firstxu:lastxu, firstyu:lastyu] $
, glamu[firstxu:lastxu, firstyu]#replicate(1, nyu)):key_irregular = 1b
array_equal(gphiu[firstxu:lastxu, firstyu:lastyu] $
, replicate(1, nxu)#(gphiu[firstxu, firstyu:lastyu])[*]):key_irregular = 1b
ELSE:key_irregular = 0b
ENDCASE
END
(where(gridtype eq 'V'))[0] NE -1 AND nxv NE 0 AND nyv NE 0:BEGIN
CASE 0 OF
array_equal(glamv[firstxv:lastxv, firstyv], glamv[firstxv:lastxv, lastyv]):key_irregular = 1b
array_equal(gphiv[firstxv, firstyv:lastyv], gphiv[lastxv, firstyv:lastyv]):key_irregular = 1b
array_equal(glamv[firstxv:lastxv, firstyv:lastyv] $
, glamv[firstxv:lastxv, firstyv]#replicate(1, nyv)):key_irregular = 1b
array_equal(gphiv[firstxv:lastxv, firstyv:lastyv] $
, replicate(1, nxv)#(gphiv[firstxv, firstyv:lastyv])[*]):key_irregular = 1b
ELSE:key_irregular = 0b
ENDCASE
END
(where(gridtype eq 'F'))[0] NE -1 AND nxf NE 0 AND nyf NE 0:BEGIN
CASE 0 OF
array_equal(glamf[firstxf:lastxf, firstyf], glamf[firstxf:lastxf, lastyf]):key_irregular = 1b
array_equal(gphif[firstxf, firstyf:lastyf], gphif[lastxf, firstyf:lastyf]):key_irregular = 1b
array_equal(glamf[firstxf:lastxf, firstyf:lastyf] $
, glamf[firstxf:lastxf, firstyf]#replicate(1, nyf)):key_irregular = 1b
array_equal(gphif[firstxf:lastxf, firstyf:lastyf] $
, replicate(1, nxf)#(gphif[firstxf, firstyf:lastyf])[*]):key_irregular = 1b
ELSE:key_irregular = 0b
ENDCASE
END
ELSE:
ENDCASE
;
;-------------------------------------------------------------------
;-------------------------------------------------------------------
; define all vertical parameters ...
; vert1, vert2
; firstz[tw], lastz[tw], nz[tw]
;-------------------------------------------------------------------
;-------------------------------------------------------------------
;
vertical:
;
;-------------------------------------------------------------------
; vertical domain defined with vert1, vert2
;-------------------------------------------------------------------
IF NOT (keyword_set(zindex) OR keyword_set(index)) THEN BEGIN
; define vert1 et vert2
CASE N_PARAMS() OF
2:vert1 = min([x1, x2], max = vert2)
6:vert1 = min([z1, z2], max = vert2)
ELSE:BEGIN
IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN $
vert1t = min(gdept, max = vert2t)
IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN $
vert1w = min(gdepw, max = vert2w)
vert1 = min([vert1t, vert1w])
vert2 = max([vert2t, vert2w])
END
ENDCASE
; define firstzt, firstzt, nzt
IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN
domz = where(gdept ge vert1 and gdept le vert2, nzt)
IF nzt NE 0 THEN BEGIN
firstzt = domz[0]
lastzt = domz[nzt-1]
ENDIF ELSE BEGIN
ras = report('WARNING, The box does not contain any T level...')
firstzt = -1
lastzt = -1
ENDELSE
ENDIF
; define firstzw, firstzw, nzw
IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
IF keyword_set(memeindices) THEN BEGIN
firstzw = firstzt
lastzw = lastzt
nzw = nzt
ENDIF ELSE BEGIN
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([x1, x2], max = lstz)
4:return
6:fstz = min([z1, z2], max = lstz)
ENDCASE
IF fstz LT 0 OR lstz GE jpk THEN BEGIN
ras = report('Bad definition of X1, X2, Z1 or Z2')
return
ENDIF
nz = lstz - fstz + 1
; find vert1t, vert2t, firstzt, firstzt, nzt
; according to (x1, x2) or (z1, z2)
IF (where(gridtype eq 'T'))[0] NE -1 THEN BEGIN
vert1t = min(gdept[fstz:lstz], max = vert2t)
firstzt = fstz & lastzt = lstz & nzt = nz
ENDIF
; find vert1w, vert2w, firstzw, firstzw, nzw
; according to (x1, x2) or (z1, z2)
IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
vert1w = min(gdepw[fstz:lstz], max = vert2w)
firstzw = fstz & lastzw = lstz & nzw = nz
ENDIF
vert1 = min([vert1t, vert1w])
vert2 = max([vert2t, vert2w])
ENDELSE
;
;-------------------------------------------------------------------
IF NOT keyword_set(key_forgetold) THEN BEGIN
@updateold
ENDIF
;-------------------------------------------------------------------
if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun
;-------------------------------------------------------------------
;-------------------------------------------------------------------
return
end