Context Navigation

domdef.pro @ 493

Last change on this file since 493 was 493, checked in by pinsard, 10 years ago
fix some typos in comments
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 47.9 KB

Rev	Line
[2]	1	;+
	2	;
[226]	3	; @file_comments
	4	; Allows to extract a sub-domain of study by providing parameters
[142]	5	; needed for drawings (see outputs)
[2]	6	;
[142]	7	; @categories
[268]	8	; Grid
[2]	9	;
[142]	10	; @param Z1 {in}{optional}
	11	; For a 3d domain whose the horizontal part cover all glam
[226]	12	;
[142]	13	; @param Z2 {in}{optional}
	14	; For a 3d domain whose the horizontal part cover all gphi
[2]	15	;
[142]	16	; @param X1 {in}{optional}
	17	; Define the minimum longitude. (All levels are selected)
[2]	18	;
[142]	19	; @param X2 {in}{optional}
	20	; Define the maximum longitude. (All levels are selected)
[2]	21	;
[142]	22	; @param Y1 {in}{optional}
	23	; Define the minimum latitude. (All levels are selected)
[13]	24	;
[142]	25	; @param Y2 {in}{optional}
	26	; Define the maximum latitude. (All levels are selected)
[2]	27	;
[163]	28	; @keyword ENDPOINTS {type=vector}
[142]	29	; A four elements vector [x1,y1,x2,y2] used to specify
[268]	30	; that <pro>domdef</pro> must define the box used to make a plot
	31	; (<pro>pltz</pro>, <pro>pltt</pro>, <pro>plt1d</pro>) done strictly along the line (that can have any direction)
[142]	32	; starting at (x1, y1) ending at (x2, y2). When defining endpoints,
	33	; you must also define TYPE which define the type of plots
	34	; ('pltz', 'xt', 'yt', 'zt', 'x', 'y', 'z', 't') will used
	35	; ENDPOINTS keywords
[2]	36	;
[142]	37	; @keyword FINDALWAYS
[226]	38	; Force to redefine a box even when none point is find in the box.
[142]	39	; In this case, we select all the grid.
[2]	40	;
[163]	41	; @keyword GRIDTYPE {type=string or vector}
[226]	42	; It is a string or a vector of strings containing the grids's name
	43	; (Only determined by 'T','U','V','W','F') for which the calculation
	44	; must be done.
[163]	45	; For example: 'T' or ['T','U']
[2]	46	;
[142]	47	; @keyword MEMEINDICES
[268]	48	; It is possible that points T, U, V and F correspond to a same geographic
[142]	49	; box which do not concern the same array indexes. This is sometimes a
	50	; problem (or at least serious complications) in programs where several
[268]	51	; type of grid intervene (see <pro>norme</pro>, <pro>curl</pro>...).
	52	; Activate MEMEINDICES to force <pro>domdef</pro> to take same indexes - those
[254]	53	; of the grid T- for all other grids.
[2]	54	;
[226]	55	; @keyword INDEX
[268]	56	; We activate it if we want that all elements passed in input of
	57	; <pro>domdef</pro> refer to indexes of glam, gphi and gdep arrays rather
[254]	58	; than to values of these arrays.
[2]	59	;
[325]	60	; @keyword TYPE
	61	;
[226]	62	; @keyword XINDEX
[268]	63	; We activate it if we want that all elements passed in input of
	64	; <pro>domdef</pro>
[226]	65	; and concerning the X dimension refer to indexes of glam arrays rather
[142]	66	; than to values of these arrays.
[2]	67	;
[226]	68	; @keyword YINDEX
[268]	69	; We activate it if we want that all elements passed in input of
[254]	70	; <pro>domdef</pro>
[226]	71	; and concerning the X dimension refer to indexes of gphi arrays rather
[142]	72	; than to values of these arrays.
[226]	73	;
	74	; @keyword ZINDEX
[268]	75	; We activate it if we want that all elements passed in input of
[254]	76	; <pro>domdef</pro>
[226]	77	; and concerning the X dimension refer to indexes of gdep arrays rather
[142]	78	; than to values of these arrays.
	79	;
	80	; @uses
[370]	81	; <pro>common</pro>
[2]	82	;
[142]	83	; @history
[157]	84	; Sebastien Masson (smasson\@lodyc.jussieu.fr)
[142]	85	; 8/2/98
[325]	86	; rewrite everything, debug and speed-up Sebastien Masson April 2005
[2]	87	;
[142]	88	; @version
	89	; $Id$
[2]	90	;
[231]	91	; @todo
	92	; seb: output pas clair/ pas d'input required?
[2]	93	;-
[327]	94	PRO domdef, x1, x2, y1, y2, z1, z2, FINDALWAYS=findalways $
	95	, GRIDTYPE=gridtype, MEMEINDICES=memeindices $
	96	, XINDEX=xindex, YINDEX=yindex, ZINDEX=zindex $
	97	, ENDPOINTS=endpoints, TYPE=type $
	98	, INDEX=index, _EXTRA=ex
[114]	99	;
	100	compile_opt idl2, strictarrsubs
	101	;
[13]	102	@cm_4mesh
	103	IF NOT keyword_set(key_forgetold) THEN BEGIN
	104	@updatenew
	105	@updatekwd
	106	ENDIF
	107	;---------------------
[142]	108	tempsun = systime(1) ; For key_performance
[13]	109	;
	110	CASE N_PARAMS() OF
	111	0:
	112	1:
	113	2:
	114	4:
	115	6:
	116	ELSE:BEGIN
	117	ras = report('Bad number of parameter in the call of domdef.')
	118	RETURN
	119	END
	120	ENDCASE
	121	;
	122	IF keyword_set(endpoints) THEN BEGIN
[226]	123	IF NOT keyword_set(type) THEN BEGIN
[13]	124	dummy = report(['If domdef is used do find the box associated' $
	125	, 'to endpoints, you must also specify type keyword'])
	126	return
[226]	127	ENDIF
[13]	128	CASE N_PARAMS() OF
	129	0:
	130	1:boxzoom = [x1]
	131	2:boxzoom = [x1, x2]
	132	4:boxzoom = [x1, x2, y1, y2]
	133	6:boxzoom = [x1, x2, y1, y2, z1, z2]
	134	ENDCASE
	135	section, BOXZOOM = boxzoom, ENDPOINTS = endpoints, TYPE = type, /ONLYBOX
	136	return
[226]	137	ENDIF
[13]	138
[2]	139	;-------------------------------------------------------------------
[13]	140	; recall domdef when there is only one input parameter
[2]	141	;-------------------------------------------------------------------
[13]	142	;
	143	IF N_PARAMS() EQ 1 THEN BEGIN
	144	CASE n_elements(x1) OF
	145	2:domdef, x1[0], x1[1], FINDALWAYS = findalways, GRIDTYPE = gridtype, MEMEINDICES = memeindices, XINDEX = xindex, YINDEX = yindex, ZINDEX = zindex, INDEX = index, _extra = ex
	146	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
	147	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
	148	ELSE:BEGIN
	149	ras = report('Bad number of elements in x1')
	150	RETURN
	151	END
	152	ENDCASE
[226]	153	RETURN
[13]	154	ENDIF
	155	;-------------------------------------------------------------------
	156	; default definitions and checks
	157	;-------------------------------------------------------------------
[469]	158	IF NOT keyword_set(gridtype) THEN gdtype = ['T', 'U', 'V', 'W', 'F'] $
	159	ELSE gdtype = strupcase(gridtype)
	160	IF keyword_set(memeindices) THEN gdtype = ['T', gdtype]
	161	IF finite(glamu[0]) eq 0 THEN gdtype = gdtype[where(gdtype NE 'U')]
	162	IF finite(glamv[0]) eq 0 THEN gdtype = gdtype[where(gdtype NE 'V')]
[13]	163	; default definitions
	164	lon1t = 99999. & lon2t = -99999. & lat1t = 99999. & lat2t = -99999.
	165	lon1u = 99999. & lon2u = -99999. & lat1u = 99999. & lat2u = -99999.
	166	lon1v = 99999. & lon2v = -99999. & lat1v = 99999. & lat2v = -99999.
	167	lon1f = 99999. & lon2f = -99999. & lat1f = 99999. & lat2f = -99999.
	168	vert1t = 99999. & vert2t = -99999. & vert1w = 99999. & vert2w = -99999.
	169	;
[226]	170	IF jpj EQ 1 THEN BEGIN
[469]	171	IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1 THEN BEGIN
[74]	172	glamt = reform(glamt, jpi, jpj, /over)
	173	gphit = reform(gphit, jpi, jpj, /over)
[226]	174	ENDIF
[469]	175	IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
[74]	176	glamu = reform(glamu, jpi, jpj, /over)
	177	gphiu = reform(gphiu, jpi, jpj, /over)
[226]	178	ENDIF
[469]	179	IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
[74]	180	glamv = reform(glamv, jpi, jpj, /over)
	181	gphiv = reform(gphiv, jpi, jpj, /over)
[226]	182	ENDIF
[469]	183	IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
[74]	184	glamf = reform(glamf, jpi, jpj, /over)
	185	gphif = reform(gphif, jpi, jpj, /over)
[226]	186	ENDIF
[74]	187	ENDIF
	188	;
[13]	189	IF N_PARAMS() EQ 2 THEN GOTO, vertical
	190	;
	191	;-------------------------------------------------------------------
	192	;-------------------------------------------------------------------
	193	; define all horizontal parameters ...
[296]	194	; lon1 and lon2
	195	; lat1 and lat2
[13]	196	; firstx[tuvf], lastx[tuvf], nx[tuvf]
	197	;-------------------------------------------------------------------
	198	; check if the grid is defined for U and V points. If not, take care
[469]	199	; of the cases gdtype eq 'U' or 'V'
[13]	200	;-------------------------------------------------------------------
	201	errstatus = 0
[469]	202	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
[13]	203	firstxu = !values.f_nan
	204	lastxu = !values.f_nan
	205	nxu = !values.f_nan
[469]	206	okgrid = where(gdtype NE 'U', count)
	207	IF count NE 0 THEN gdtype = gdtype[okgrid] $
[13]	208	ELSE errstatus = report('U grid is undefined. Impossible to call domdef with vargid = ''U''')
	209	ENDIF
	210	;
[469]	211	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
[13]	212	firstxv = !values.f_nan
	213	lastxv = !values.f_nan
	214	nxv = !values.f_nan
[469]	215	okgrid = where(gdtype NE 'V', count)
	216	IF count NE 0 THEN gdtype = gdtype[okgrid] $
[13]	217	ELSE errstatus = report('V grid is undefined. Impossible to call domdef with vargid = ''V''')
	218	ENDIF
	219	IF errstatus EQ -1 THEN return
	220	;-------------------------------------------------------------------
	221	;-------------------------------------------------------------------
	222	; horizontal domain defined with lon1, lon2, lat1 and lat2
	223	;-------------------------------------------------------------------
	224	;-------------------------------------------------------------------
	225	IF N_PARAMS() EQ 0 $
	226	OR ( (N_PARAMS() EQ 4 OR N_PARAMS() EQ 6) $
	227	AND NOT keyword_set(xindex) AND NOT keyword_set(yindex) AND NOT keyword_set(index) ) THEN BEGIN
	228	IF N_PARAMS() EQ 0 THEN BEGIN
[226]	229	; find lon1 and lon2 the longitudinal boundaries of the full domain
[469]	230	IF (where(gdtype eq 'T'))[0] NE -1 THEN lon1t = min(glamt, max = lon2t)
	231	IF (where(gdtype eq 'W'))[0] NE -1 AND (where(gdtype eq 'T'))[0] EQ -1 THEN lon1t = min(glamt, max = lon2t)
	232	IF (where(gdtype eq 'U'))[0] NE -1 THEN lon1u = min(glamu, max = lon2u)
	233	IF (where(gdtype eq 'V'))[0] NE -1 THEN lon1v = min(glamv, max = lon2v)
	234	IF (where(gdtype eq 'F'))[0] NE -1 THEN lon1f = min(glamf, max = lon2f)
[13]	235	lon1 = min([lon1t, lon1u, lon1v, lon1f])
	236	lon2 = max([lon2t, lon2u, lon2v, lon2f])
[226]	237	; find lat1 and lat2 the latitudinal boundaries of the full domain
[469]	238	IF (where(gdtype eq 'T'))[0] NE -1 THEN lat1t = min(gphit, max = lat2t)
	239	IF (where(gdtype eq 'W'))[0] NE -1 AND (where(gdtype eq 'T'))[0] EQ -1 THEN lat1t = min(gphit, max = lat2t)
	240	IF (where(gdtype eq 'U'))[0] NE -1 THEN lat1u = min(gphiu, max = lat2u)
	241	IF (where(gdtype eq 'V'))[0] NE -1 THEN lat1v = min(gphiv, max = lat2v)
	242	IF (where(gdtype eq 'F'))[0] NE -1 THEN lat1f = min(gphif, max = lat2f)
[13]	243	lat1 = min([lat1t, lat1u, lat1v, lat1f])
	244	lat2 = max([lat2t, lat2u, lat2v, lat2f])
[226]	245	ENDIF ELSE BEGIN
[13]	246	lon1 = min([x1, x2], max = lon2)
	247	lat1 = min([y1, y2], max = lat2)
	248	ENDELSE
	249	; find firstxt, firstxt, nxt and nyt according to lon1, lon2, lat1 and lat2
[469]	250	IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1 THEN BEGIN
[13]	251	dom = where( (glamt GE lon1) AND (glamt LE lon2) $
	252	AND (gphit GE lat1) AND (gphit LE lat2) )
	253	IF (dom[0] EQ -1) THEN BEGIN
	254	IF keyword_set(findalways) THEN BEGIN
[280]	255	ras = report('WARNING, empty T points box... we get the neighbors to define a new box...')
[74]	256	neig1 = neighbor(lon1, lat1, glamt, gphit, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	257	neig2 = neighbor(lon2, lat2, glamt, gphit, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	258	CASE N_PARAMS() OF
[469]	259	4:domdef, min(glamt[neig1]), max(glamt[neig2]), min(gphit[neig1]), max(gphit[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
	260	6:domdef, min(glamt[neig1]), max(glamt[neig2]), min(gphit[neig1]), max(gphit[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
[74]	261	ENDCASE
[13]	262	RETURN
	263	ENDIF ELSE BEGIN
	264	ras = report('WARNING, The box does not contain any T points...')
	265	firstxt = -1 & lastxt = -1 & nxt = 0
	266	firstyt = -1 & lastyt = -1 & nyt = 0
	267	ENDELSE
	268	ENDIF ELSE BEGIN
	269	jyt = dom / jpi
	270	ixt = temporary(dom) MOD jpi
	271	firstxt = min(temporary(ixt), max = lastxt)
	272	firstyt = min(temporary(jyt), max = lastyt)
	273	nxt = lastxt - firstxt + 1
	274	nyt = lastyt - firstyt + 1
[2]	275	ENDELSE
[13]	276	ENDIF
	277	; find firstxu, firstxu, firstyu, firstyu, nxu and nyu
	278	; according to lon1, lon2, lat1 and lat2
[469]	279	IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
[13]	280	IF keyword_set(memeindices) THEN BEGIN
	281	firstxu = firstxt & lastxu = lastxt & nxu = nxt
	282	firstyu = firstyt & lastyu = lastyt & nyu = nyt
	283	ENDIF ELSE BEGIN
	284	dom = where( (glamu GE lon1) AND (glamu LE lon2) $
	285	AND (gphiu GE lat1) AND (gphiu LE lat2) )
	286	IF (dom[0] EQ -1) THEN BEGIN
	287	IF keyword_set(findalways) THEN BEGIN
[493]	288	; if t grid parameters already defined, we use them...
[74]	289	CASE 1 OF
[469]	290	(where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
[240]	291	ras = report('WARNING, empty U points box... we use the same index as T points...')
[74]	292	firstxu = firstxt & lastxu = lastxt & nxu = nxt
	293	firstyu = firstyt & lastyu = lastyt & nyu = nyt
	294	END
	295	ELSE:BEGIN
[280]	296	ras = report('WARNING, empty U points box... we get the neighbors to define a new box...')
[74]	297	neig1 = neighbor(lon1, lat1, glamu, gphiu, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	298	neig2 = neighbor(lon2, lat2, glamu, gphiu, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	299	CASE N_PARAMS() OF
[469]	300	4:domdef, min(glamu[neig1]), max(glamu[neig2]), min(gphiu[neig1]), max(gphiu[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
	301	6:domdef, min(glamu[neig1]), max(glamu[neig2]), min(gphiu[neig1]), max(gphiu[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
[74]	302	ENDCASE
	303	RETURN
	304	END
	305	ENDCASE
[13]	306	ENDIF ELSE BEGIN
	307	ras = report('WARNING, The box does not contain any U points...')
	308	firstxu = -1 & lastxu = -1 & nxu = 0
	309	firstyu = -1 & lastyu = -1 & nyu = 0
	310	ENDELSE
	311	ENDIF ELSE BEGIN
	312	jyu = dom / jpi
	313	ixu = temporary(dom) MOD jpi
	314	firstxu = min(temporary(ixu), max = lastxu)
	315	firstyu = min(temporary(jyu), max = lastyu)
	316	nxu = lastxu - firstxu + 1
	317	nyu = lastyu - firstyu + 1
	318	ENDELSE
[2]	319	ENDELSE
[13]	320	ENDIF
	321	; find firstxv, firstxv, firstyv, firstyv, nxv and nyv
	322	; according to lon1, lon2, lat1 and lat2
[469]	323	IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
[13]	324	IF keyword_set(memeindices) THEN BEGIN
	325	firstxv = firstxt & lastxv = lastxt & nxv = nxt
	326	firstyv = firstyt & lastyv = lastyt & nyv = nyt
	327	ENDIF ELSE BEGIN
	328	dom = where( (glamv GE lon1) AND (glamv LE lon2) $
	329	AND (gphiv GE lat1) AND (gphiv LE lat2) )
	330	IF (dom[0] EQ -1) THEN BEGIN
	331	IF keyword_set(findalways) THEN BEGIN
[74]	332	CASE 1 OF
[469]	333	(where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
[240]	334	ras = report('WARNING, empty V points box... we use the same index as T points...')
[74]	335	firstxv = firstxt & lastxv = lastxt & nxv = nxt
	336	firstyv = firstyt & lastyv = lastyt & nyv = nyt
[226]	337	END
[469]	338	(where(gdtype eq 'U'))[0] NE -1:BEGIN
[240]	339	ras = report('WARNING, empty V points box... we use the same index as U points...')
[74]	340	firstxv = firstxu & lastxv = lastxu & nxv = nxu
	341	firstyv = firstyu & lastyv = lastyu & nyv = nyu
[226]	342	END
[74]	343	ELSE:BEGIN
[280]	344	ras = report('WARNING, empty V points box... we get the neighbors to define a new box...')
[74]	345	neig1 = neighbor(lon1, lat1, glamv, gphiv, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	346	neig2 = neighbor(lon2, lat2, glamv, gphiv, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	347	CASE N_PARAMS() OF
[469]	348	4:domdef, min(glamv[neig1]), max(glamv[neig2]), min(gphiv[neig1]), max(gphiv[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
	349	6:domdef, min(glamv[neig1]), max(glamv[neig2]), min(gphiv[neig1]), max(gphiv[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
[74]	350	ENDCASE
	351	RETURN
	352	END
	353	ENDCASE
[13]	354	ENDIF ELSE BEGIN
	355	ras = report('WARNING, The box does not contain any V points...')
	356	firstxv = -1 & lastxv = -1 & nxv = 0
	357	firstyv = -1 & lastyv = -1 & nyv = 0
	358	ENDELSE
	359	ENDIF ELSE BEGIN
	360	jyv = dom / jpi
	361	ixv = temporary(dom) MOD jpi
	362	firstxv = min(temporary(ixv), max = lastxv)
	363	firstyv = min(temporary(jyv), max = lastyv)
	364	nxv = lastxv - firstxv + 1
	365	nyv = lastyv - firstyv + 1
	366	ENDELSE
[2]	367	ENDELSE
[13]	368	ENDIF
	369	; find firstxf, firstxf, firstyf, firstyf, nxf and nyf
	370	; according to lon1, lon2, lat1 and lat2
[469]	371	IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
[13]	372	IF keyword_set(memeindices) THEN BEGIN
	373	firstxf = firstxt & lastxf = lastxt & nxf = nxt
	374	firstyf = firstyt & lastyf = lastyt & nyf = nyt
	375	ENDIF ELSE BEGIN
	376	dom = where( (glamf GE lon1) AND (glamf LE lon2) $
	377	AND (gphif GE lat1) AND (gphif LE lat2) )
	378	IF (dom[0] EQ -1) THEN BEGIN
	379	IF keyword_set(findalways) THEN BEGIN
[74]	380	CASE 1 OF
[469]	381	(where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
[240]	382	ras = report('WARNING, empty F points box... we use the same index as T points...')
[74]	383	firstxf = firstxt & lastxf = lastxt & nxf = nxt
	384	firstyf = firstyt & lastyf = lastyt & nyf = nyt
[226]	385	END
[469]	386	(where(gdtype eq 'U'))[0] NE -1:BEGIN
[240]	387	ras = report('WARNING, empty F points box... we use the same index as U points...')
[74]	388	firstxf = firstxu & lastxf = lastxu & nxf = nxu
	389	firstyf = firstyu & lastyf = lastyu & nyf = nyu
[226]	390	END
[469]	391	(where(gdtype eq 'V'))[0] NE -1:BEGIN
[240]	392	ras = report('WARNING, empty F points box... we use the same index as V points...')
[74]	393	firstxf = firstxv & lastxf = lastxv & nxf = nxv
	394	firstyf = firstyv & lastyf = lastyv & nyf = nyv
[226]	395	END
[74]	396	ELSE:BEGIN
[280]	397	ras = report('WARNING, empty F points box... we get the neighbors to define a new box...')
[74]	398	neig1 = neighbor(lon1, lat1, glamf, gphif, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	399	neig2 = neighbor(lon2, lat2, glamf, gphif, sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	400	CASE N_PARAMS() OF
[469]	401	4:domdef, min(glamf[neig1]), max(glamf[neig2]), min(gphif[neig1]), max(gphif[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
	402	6:domdef, min(glamf[neig1]), max(glamf[neig2]), min(gphif[neig1]), max(gphif[neig2]), z1, z2, GDTYPE = gdtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
[74]	403	ENDCASE
	404	RETURN
	405	END
[226]	406	ENDCASE
[13]	407	ENDIF ELSE BEGIN
	408	ras = report('WARNING, The box does not contain any F points...')
	409	firstxf = -1 & lastxf = -1 & nxf = 0
	410	firstyf = -1 & lastyf = -1 & nyf = 0
	411	ENDELSE
	412	ENDIF ELSE BEGIN
	413	jyf = dom / jpi
	414	ixf = temporary(dom) MOD jpi
	415	firstxf = min(temporary(ixf), max = lastxf)
	416	firstyf = min(temporary(jyf), max = lastyf)
	417	nxf = lastxf - firstxf + 1
	418	nyf = lastyf - firstyf + 1
	419	ENDELSE
[2]	420	ENDELSE
[13]	421	ENDIF
	422	;
	423	ENDIF ELSE BEGIN
	424	CASE 1 OF
[2]	425	;-------------------------------------------------------------------
	426	;-------------------------------------------------------------------
[13]	427	; horizontal domain defined with the X and Y indexes
[2]	428	;-------------------------------------------------------------------
	429	;-------------------------------------------------------------------
[13]	430	(keyword_set(xindex) AND keyword_set(yindex)) OR keyword_set(index):BEGIN
[393]	431	fstx = min(long([x1, x2]), max = lstx)
	432	fsty = min(long([y1, y2]), max = lsty)
[13]	433	IF fstx LT 0 OR lstx GE jpi THEN BEGIN
	434	ras = report('Bad definition of X1 or X2')
	435	return
	436	ENDIF
	437	IF fsty LT 0 OR lsty GE jpj THEN BEGIN
	438	ras = report('Bad definition of Y1 or Y2')
	439	return
	440	ENDIF
	441	nx = lstx - fstx + 1
	442	ny = lsty - fsty + 1
	443	; find lon1t, lon2t, lat1t, lat2t, firstxt, firstxt, nxt and nyt
	444	; according to x1, x2, y1, y2
[469]	445	IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype eq 'W'))[0] NE -1 THEN BEGIN
[13]	446	lon1t = min(glamt[fstx:lstx, fsty:lsty], max = lon2t)
	447	lat1t = min(gphit[fstx:lstx, fsty:lsty], max = lat2t)
	448	firstxt = fstx & lastxt = lstx
	449	firstyt = fsty & lastyt = lsty
	450	nxt = nx & nyt = ny
	451	ENDIF
	452	; find lon1u, lon2u, lat1u, lat2u, firstxu, firstxu, nxu and nyu
	453	; according to x1, x2, y1, y2
[469]	454	IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
[13]	455	lon1u = min(glamu[fstx:lstx, fsty:lsty], max = lon2u)
	456	lat1u = min(gphiu[fstx:lstx, fsty:lsty], max = lat2u)
	457	firstxu = fstx & lastxu = lstx
	458	firstyu = fsty & lastyu = lsty
	459	nxu = nx & nyu = ny
	460	ENDIF
	461	; find lon1v, lon2v, lat1v, lat2v, firstxv, firstxv, nxv and nyv
	462	; according to x1, x2, y1, y2
[469]	463	IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
[13]	464	lon1v = min(glamv[fstx:lstx, fsty:lsty], max = lon2v)
	465	lat1v = min(gphiv[fstx:lstx, fsty:lsty], max = lat2v)
	466	firstxv = fstx & lastxv = lstx
	467	firstyv = fsty & lastyv = lsty
	468	nxv = nx & nyv = ny
[226]	469	ENDIF
[13]	470	; find lon1f, lon2f, lat1f, lat2f, firstxf, firstxf, nxf and nyf
	471	; according to x1, x2, y1, y2
[469]	472	IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
[13]	473	lon1f = min(glamf[fstx:lstx, fsty:lsty], max = lon2f)
	474	lat1f = min(gphif[fstx:lstx, fsty:lsty], max = lat2f)
	475	firstxf = fstx & lastxf = lstx
	476	firstyf = fsty & lastyf = lsty
	477	nxf = nx & nyf = ny
	478	ENDIF
	479	lon1 = min([lon1t, lon1u, lon1v, lon1f])
	480	lon2 = max([lon2t, lon2u, lon2v, lon2f])
	481	lat1 = min([lat1t, lat1u, lat1v, lat1f])
	482	lat2 = max([lat2t, lat2u, lat2v, lat2f])
	483	END
	484	;-------------------------------------------------------------------
	485	;-------------------------------------------------------------------
	486	; horizontal domain defined with the X index and lat1/lat2
	487	;-------------------------------------------------------------------
	488	;-------------------------------------------------------------------
	489	keyword_set(xindex):BEGIN
[393]	490	fstx = min(long([x1, x2]), max = lstx)
[13]	491	IF fstx LT 0 OR lstx GE jpi THEN BEGIN
	492	ras = report('Bad definition of X1 or X2')
	493	return
	494	ENDIF
	495	nx = lstx - fstx + 1
	496	lat1 = min([y1, y2], max = lat2)
	497	; find lon1t, lon2t, firstxt, firstxt, firstyt, firstyt, nxt
	498	; and nyt according to x1, x2, lat1 and lat2
[469]	499	IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1 THEN BEGIN
[13]	500	firstxt = fstx & lastxt = lstx & nxt = nx
	501	dom = where( (gphit[fstx:lstx, ] GE lat1) AND (gphit[fstx:lstx, ] LE lat2) )
	502	IF (dom[0] EQ -1) THEN BEGIN
	503	IF keyword_set(findalways) THEN BEGIN
[280]	504	ras = report('WARNING, empty T points box... we get the neighbors to define a new box...')
[74]	505	neig1 = neighbor(lon1, lat1, glamt[fstx:lstx, ], gphit[fstx:lstx, ], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	506	neig2 = neighbor(lon2, lat2, glamt[fstx:lstx, ], gphit[fstx:lstx, ], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	507	CASE N_PARAMS() OF
[469]	508	4:domdef, fstx, lstx, min((gphit[fstx:lstx, ])[neig1]), max((gphit[fstx:lstx, ])[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
	509	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
[74]	510	ENDCASE
[13]	511	RETURN
[2]	512	ENDIF ELSE BEGIN
[13]	513	ras = report('WARNING, The box does not contain any T points...')
	514	firstyt = -1 & lastyt = -1 & nyt = 0
[2]	515	ENDELSE
[13]	516	ENDIF ELSE BEGIN
	517	jyt = temporary(dom) / nx
	518	firstyt = min(temporary(jyt), max = lastyt)
	519	nyt = lastyt - firstyt + 1
	520	ENDELSE
	521	IF nyt NE 0 THEN lon1t = min(glamt[firstxt:lastxt, firstyt:lastyt], max = lon2t)
	522	ENDIF
	523	; find lon1u, lon2u, firstxu, firstxu, firstyu, firstyu, nxu
	524	; and nyu according to x1, x2, lat1 and lat2
[469]	525	IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
[13]	526	firstxu = fstx & lastxu = lstx & nxu = nx
	527	IF keyword_set(memeindices) THEN BEGIN
	528	firstyu = firstyt & lastyu = lastyt & nyu = nyt
	529	ENDIF ELSE BEGIN
	530	dom = where( (gphiu[fstx:lstx, ] GE lat1) AND (gphiu[fstx:lstx, ] LE lat2) )
	531	IF (dom[0] EQ -1) THEN BEGIN
	532	IF keyword_set(findalways) THEN BEGIN
[74]	533	CASE 1 OF
[469]	534	(where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
[240]	535	ras = report( 'WARNING, empty U points box... we use the same index as T points...')
[74]	536	firstyu = firstyt & lastyu = lastyt & nyu = nyt
[226]	537	END
[74]	538	ELSE:BEGIN
[280]	539	ras = report('WARNING, empty U points box... we get the neighbors to define a new box...')
[74]	540	neig1 = neighbor(lon1, lat1, glamu[fstx:lstx, ], gphiu[fstx:lstx, ], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	541	neig2 = neighbor(lon2, lat2, glamu[fstx:lstx, ], gphiu[fstx:lstx, ], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	542	CASE N_PARAMS() OF
[469]	543	4:domdef, fstx, lstx, min((gphiu[fstx:lstx, ])[neig1]), max((gphiu[fstx:lstx, ])[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
	544	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
[74]	545	ENDCASE
	546	RETURN
	547	END
[226]	548	ENDCASE
[13]	549	ENDIF ELSE BEGIN
	550	ras = report('WARNING, The box does not contain any U points...')
	551	firstyu = -1 & lastyu = -1 & nyu = 0
	552	ENDELSE
	553	ENDIF ELSE BEGIN
	554	jyu = temporary(dom) / nx
	555	firstyu = min(temporary(jyu), max = lastyu)
	556	nyu = lastyu - firstyu + 1
	557	ENDELSE
	558	ENDELSE
	559	IF nyu NE 0 THEN lon1u = min(glamu[firstxu:lastxu, firstyu:lastyu], max = lon2u)
	560	ENDIF
	561	; find lon1v, lon2v, firstxv, firstxv, firstyv, firstyv, nxv
	562	; and nyv according to x1, x2, lat1 and lat2
[469]	563	IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
[13]	564	firstxv = fstx & lastxv = lstx & nxv = nx
	565	IF keyword_set(memeindices) THEN BEGIN
	566	firstyv = firstyt & lastyv = lastyt & nyv = nyt
	567	ENDIF ELSE BEGIN
	568	dom = where( (gphiv[fstx:lstx, ] GE lat1) AND (gphiv[fstx:lstx, ] LE lat2) )
	569	IF (dom[0] EQ -1) THEN BEGIN
	570	IF keyword_set(findalways) THEN BEGIN
[74]	571	CASE 1 OF
[469]	572	(where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
[240]	573	ras = report('WARNING, empty V points box... we use the same index as T points...')
[74]	574	firstyv = firstyt & lastyv = lastyt & nyv = nyt
[226]	575	END
[469]	576	(where(gdtype eq 'U'))[0] NE -1:BEGIN
[240]	577	ras = report('WARNING, empty V points box... we use the same index as U points...')
[74]	578	firstyv = firstyu & lastyv = lastyu & nyv = nyu
[226]	579	END
[74]	580	ELSE:BEGIN
[280]	581	ras = report('WARNING, empty V points box... we get the neighbors to define a new box...')
[74]	582	neig1 = neighbor(lon1, lat1, glamv[fstx:lstx, ], gphiv[fstx:lstx, ], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	583	neig2 = neighbor(lon2, lat2, glamv[fstx:lstx, ], gphiv[fstx:lstx, ], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	584	CASE N_PARAMS() OF
[469]	585	4:domdef, fstx, lstx, min((gphiv[fstx:lstx, ])[neig1]), max((gphiv[fstx:lstx, ])[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
	586	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
[74]	587	ENDCASE
	588	RETURN
[226]	589	END
	590	ENDCASE
[13]	591	ENDIF ELSE BEGIN
	592	ras = report('WARNING, The box does not contain any V points...')
	593	firstyv = -1 & lastyv = -1 & nyv = 0
	594	ENDELSE
	595	ENDIF ELSE BEGIN
	596	jyv = temporary(dom) / nx
	597	firstyv = min(temporary(jyv), max = lastyv)
	598	nyv = lastyv - firstyv + 1
	599	ENDELSE
	600	ENDELSE
	601	IF nyv NE 0 THEN lon1v = min(glamv[firstxv:lastxv, firstyv:lastyv], max = lon2v)
	602	ENDIF
	603	; find lon1f, lon2f, firstxf, firstxf, firstyf, firstyf, nxf
	604	; and nyf according to x1, x2, lat1 and lat2
[469]	605	IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
[13]	606	firstxf = fstx & lastxf = lstx & nxf = nx
	607	IF keyword_set(memeindices) THEN BEGIN
	608	firstyf = firstyt & lastyf = lastyt & nyf = nyt
	609	ENDIF ELSE BEGIN
	610	dom = where( (gphif[fstx:lstx, ] GE lat1) AND (gphif[fstx:lstx, ] LE lat2) )
	611	IF (dom[0] EQ -1) THEN BEGIN
	612	IF keyword_set(findalways) THEN BEGIN
[74]	613	CASE 1 OF
[469]	614	(where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
[240]	615	ras = report('WARNING, empty F points box... we use the same index as T points...')
[74]	616	firstyf = firstyt & lastyf = lastyt & nyf = nyt
[226]	617	END
[469]	618	(where(gdtype eq 'U'))[0] NE -1:BEGIN
[240]	619	ras = report('WARNING, empty F points box... we use the same index as U points...')
[74]	620	firstyf = firstyu & lastyf = lastyu & nyf = nyu
[226]	621	END
[469]	622	(where(gdtype eq 'V'))[0] NE -1:BEGIN
[240]	623	ras = report('WARNING, empty F points box... we use the same index as V points...')
[74]	624	firstyf = firstyv & lastyf = lastyv & nyf = nyv
[226]	625	END
[74]	626	ELSE:BEGIN
[280]	627	ras = report('WARNING, empty F points box... we get the neighbors to define a new box...')
[74]	628	neig1 = neighbor(lon1, lat1, glamf[fstx:lstx, ], gphif[fstx:lstx, ], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	629	neig2 = neighbor(lon2, lat2, glamf[fstx:lstx, ], gphif[fstx:lstx, ], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	630	CASE N_PARAMS() OF
[469]	631	4:domdef, fstx, lstx, min((gphif[fstx:lstx, ])[neig1]), max((gphif[fstx:lstx, ])[neig2]), GDTYPE = gdtype, MEMEINDICES = memeindices, /XINDEX, ZINDEX = zindex, _extra = ex
	632	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
[74]	633	ENDCASE
	634	RETURN
	635	END
[226]	636	ENDCASE
[13]	637	ENDIF ELSE BEGIN
	638	ras = report('WARNING, The box does not contain any F points...')
	639	firstyf = -1 & lastyf = -1 & nyf = 0
	640	ENDELSE
	641	ENDIF ELSE BEGIN
	642	jyf = temporary(dom) / nx
	643	firstyf = min(temporary(jyf), max = lastyf)
	644	nyf = lastyf - firstyf + 1
	645	ENDELSE
	646	ENDELSE
	647	IF nyf NE 0 THEN lon1f = min(glamf[firstxf:lastxf, firstyf:lastyf], max = lon2f)
	648	ENDIF
	649	lon1 = min([lon1t, lon1u, lon1v, lon1f])
	650	lon2 = max([lon2t, lon2u, lon2v, lon2f])
	651	END
[2]	652	;-------------------------------------------------------------------
	653	;-------------------------------------------------------------------
[13]	654	; horizontal domain defined with the Y index and lon1/lon2
[2]	655	;-------------------------------------------------------------------
	656	;-------------------------------------------------------------------
[13]	657	keyword_set(yindex):BEGIN
[393]	658	fsty = min(long([y1, y2]), max = lsty)
[13]	659	IF fsty LT 0 OR lsty GE jpj THEN BEGIN
	660	ras = report('Bad definition of Y1 or Y2')
	661	return
	662	ENDIF
	663	ny = lsty - fsty + 1
	664	lon1 = min([x1, x2], max = lon2)
	665	; find lat1t, lat2t, firstxt, firstxt, firstyt, firstyt, nxt
	666	; and nyt according to x1, x2, lon1 and lon2
[469]	667	IF (where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1 THEN BEGIN
[13]	668	firstyt = fsty & lastyt = lsty & nyt = ny
	669	dom = where( (glamt[, fsty:lsty] GE lon1) AND (glamt[, fsty:lsty] LE lon2) )
	670	IF (dom[0] EQ -1) THEN BEGIN
	671	IF keyword_set(findalways) THEN BEGIN
[280]	672	ras = report('WARNING, empty T points box... we get the neighbors to define a new box...')
[74]	673	neig1 = neighbor(lon1, lat1, glamt[, fsty:lsty], gphit[, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	674	neig2 = neighbor(lon2, lat2, glamt[, fsty:lsty], gphit[, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	675	CASE N_PARAMS() OF
[469]	676	4:domdef, min((glamt[, fsty:lsty])[neig1]), max((glamt[, fsty:lsty])[neig2]), fsty, lsty, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
	677	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
[74]	678	ENDCASE
[13]	679	RETURN
	680	ENDIF ELSE BEGIN
	681	ras = report('WARNING, The box does not contain any T points...')
	682	firstxt = -1 & lastxt = -1 & nxt = 0
	683	ENDELSE
	684	ENDIF ELSE BEGIN
	685	jxt = temporary(dom) MOD jpi
	686	firstxt = min(temporary(jxt), max = lastxt)
	687	nxt = lastxt - firstxt + 1
	688	ENDELSE
	689	IF nxt NE 0 THEN lat1t = min(gphit[firstxt:lastxt, firstyt:lastyt], max = lat2t)
	690	ENDIF
	691	; find lat1u, lat2u, firstxu, firstxu, firstyu, firstyu, nxu
	692	; and nyu according to x1, x2, lon1 and lon2
[469]	693	IF (where(gdtype eq 'U'))[0] NE -1 THEN BEGIN
[13]	694	firstyu = fsty & lastyu = lsty & nyu = ny
	695	IF keyword_set(memeindices) THEN BEGIN
[492]	696	firstxu = firstxt & lastxu = lastxt & nxu = nxt
[13]	697	ENDIF ELSE BEGIN
	698	dom = where( (glamu[, fsty:lsty] GE lon1) AND (glamu[, fsty:lsty] LE lon2) )
	699	IF (dom[0] EQ -1) THEN BEGIN
	700	IF keyword_set(findalways) THEN BEGIN
[74]	701	CASE 1 OF
[469]	702	(where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
[240]	703	ras = report('WARNING, empty U points box... we use the same index as T points...')
[74]	704	firstxu = firstxt & lastxu = lastxt & nxu = nxt
[226]	705	END
[74]	706	ELSE:BEGIN
[280]	707	ras = report('WARNING, empty U points box... we get the neighbors to define a new box...')
[74]	708	neig1 = neighbor(lon1, lat1, glamu[, fsty:lsty], gphiu[, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	709	neig2 = neighbor(lon2, lat2, glamu[, fsty:lsty], gphiu[, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	710	CASE N_PARAMS() OF
[469]	711	4:domdef, min((glamu[, fsty:lsty])[neig1]), max((glamu[, fsty:lsty])[neig2]), fsty, lsty, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
	712	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
[74]	713	ENDCASE
	714	RETURN
	715	END
[226]	716	ENDCASE
[13]	717	ENDIF ELSE BEGIN
	718	ras = report('WARNING, The box does not contain any U points...')
	719	firstxu = -1 & lastxu = -1 & nxu = 0
	720	ENDELSE
	721	ENDIF ELSE BEGIN
	722	jxu = temporary(dom) MOD jpi
	723	firstxu = min(temporary(jxu), max = lastxu)
	724	nxu = lastxu - firstxu + 1
	725	ENDELSE
	726	ENDELSE
	727	IF nxu NE 0 THEN lat1u = min(gphiu[firstxu:lastxu, firstyu:lastyu], max = lat2u)
	728	ENDIF
	729	; find lat1v, lat2v, firstxv, firstxv, firstyv, firstyv, nxv
	730	; and nyv according to x1, x2, lon1 and lon2
[469]	731	IF (where(gdtype eq 'V'))[0] NE -1 THEN BEGIN
[13]	732	firstyv = fsty & lastyv = lsty & nyv = ny
	733	IF keyword_set(memeindices) THEN BEGIN
[492]	734	firstxv = firstxt & lastxv = lastxt & nxv = nxt
[13]	735	ENDIF ELSE BEGIN
	736	dom = where( (glamv[, fsty:lsty] GE lon1) AND (glamv[, fsty:lsty] LE lon2) )
	737	IF (dom[0] EQ -1) THEN BEGIN
	738	IF keyword_set(findalways) THEN BEGIN
[74]	739	CASE 1 OF
[469]	740	(where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
[240]	741	ras = report('WARNING, empty V points box... we use the same index as T points...')
[74]	742	firstxv = firstxt & lastxv = lastxt & nxv = nxt
[226]	743	END
[469]	744	(where(gdtype eq 'U'))[0] NE -1:BEGIN
[240]	745	ras = report('WARNING, empty V points box... we use the same index as U points...')
[74]	746	firstxv = firstxu & lastxv = lastxu & nxv = nxu
[226]	747	END
[74]	748	ELSE:BEGIN
[280]	749	ras = report('WARNING, empty V points box... we get the neighbors to define a new box...')
[74]	750	neig1 = neighbor(lon1, lat1, glamv[, fsty:lsty], gphiv[, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	751	neig2 = neighbor(lon2, lat2, glamv[, fsty:lsty], gphiv[, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	752	CASE N_PARAMS() OF
[469]	753	4:domdef, min((glamv[, fsty:lsty])[neig1]), max((glamv[, fsty:lsty])[neig2]), fsty, lsty, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
	754	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
[74]	755	ENDCASE
	756	RETURN
[226]	757	END
	758	ENDCASE
[13]	759	ENDIF ELSE BEGIN
	760	ras = report('WARNING, The box does not contain any V points...')
	761	firstxv = -1 & lastxv = -1 & nxv = 0
	762	ENDELSE
	763	ENDIF ELSE BEGIN
	764	jxv = temporary(dom) MOD jpi
	765	firstxv = min(temporary(jxv), max = lastxv)
	766	nxv = lastxv - firstxv + 1
	767	ENDELSE
	768	ENDELSE
	769	IF nxv NE 0 THEN lat1v = min(gphiv[firstxv:lastxv, firstyv:lastyv], max = lat2v)
	770	ENDIF
	771	; find lat1f, lat2f, firstxf, firstxf, firstyf, firstyf, nxf
	772	; and nyf according to x1, x2, lon1 and lon2
[469]	773	IF (where(gdtype eq 'F'))[0] NE -1 THEN BEGIN
[13]	774	firstyf = fsty & lastyf = lsty & nyf = ny
	775	IF keyword_set(memeindices) THEN BEGIN
[492]	776	firstxf = firstxt & lastxf = lastxt & nxf = nxt
[13]	777	ENDIF ELSE BEGIN
	778	dom = where( (glamf[, fsty:lsty] GE lon1) AND (glamf[, fsty:lsty] LE lon2) )
	779	IF (dom[0] EQ -1) THEN BEGIN
	780	IF keyword_set(findalways) THEN BEGIN
[74]	781	CASE 1 OF
[469]	782	(where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype EQ 'W'))[0] NE -1:BEGIN
[240]	783	ras = report('WARNING, empty F points box... we use the same index as T points...')
[74]	784	firstxf = firstxt & lastxf = lastxt & nxf = nxt
[226]	785	END
[469]	786	(where(gdtype eq 'U'))[0] NE -1:BEGIN
[240]	787	ras = report('WARNING, empty F points box... we use the same index as U points...')
[74]	788	firstxf = firstxu & lastxf = lastxu & nxf = nxu
[226]	789	END
[469]	790	(where(gdtype eq 'V'))[0] NE -1:BEGIN
[240]	791	ras = report('WARNING, empty F points box... we use the same index as V points...')
[74]	792	firstxf = firstxv & lastxf = lastxv & nxf = nxv
[226]	793	END
[74]	794	ELSE:BEGIN
[280]	795	ras = report('WARNING, empty F points box... we get the neighbors to define a new box...')
[74]	796	neig1 = neighbor(lon1, lat1, glamf[, fsty:lsty], gphif[, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	797	neig2 = neighbor(lon2, lat2, glamf[, fsty:lsty], gphif[, fsty:lsty], sphere = keyword_set(key_onearth) * keyword_set(key_irregular))
	798	CASE N_PARAMS() OF
[469]	799	4:domdef, min((glamf[, fsty:lsty])[neig1]), max((glamf[, fsty:lsty])[neig2]), fsty, lsty, GDTYPE = gdtype, MEMEINDICES = memeindices, /YINDEX, ZINDEX = zindex, _extra = ex
	800	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
[74]	801	ENDCASE
	802	RETURN
[226]	803	END
	804	ENDCASE
[13]	805	ENDIF ELSE BEGIN
	806	ras = report('WARNING, The box does not contain any F points...')
	807	firstxf = -1 & lastyf = -1 & nxf = 0
	808	ENDELSE
	809	ENDIF ELSE BEGIN
	810	jxf = temporary(dom) MOD jpi
	811	firstxf = min(temporary(jxf), max = lastxf)
	812	nxf = lastxf - firstxf + 1
	813	ENDELSE
	814	ENDELSE
	815	IF nxf NE 0 THEN lat1f = min(gphif[firstxf:lastxf, firstyf:lastyf], max = lat2f)
	816	ENDIF
	817	lat1 = min([lat1t, lat1u, lat1v, lat1f])
	818	lat2 = max([lat2t, lat2u, lat2v, lat2f])
	819	END
	820	ENDCASE
[226]	821	ENDELSE
[2]	822	;-------------------------------------------------------------------
[226]	823	; The extracted domain is it regular or not?
[2]	824	;-------------------------------------------------------------------
[13]	825	CASE 1 OF
[469]	826	((where(gdtype eq 'T'))[0] NE -1 OR (where(gdtype eq 'W'))[0] NE -1) AND nxt NE 0 AND nyt NE 0:BEGIN
[226]	827	; to get faster, we first test the most basic cases before testing the
[13]	828	; full array.
	829	CASE 0 OF
	830	array_equal(glamt[firstxt:lastxt, firstyt], glamt[firstxt:lastxt, lastyt]):key_irregular = 1b
	831	array_equal(gphit[firstxt, firstyt:lastyt], gphit[lastxt, firstyt:lastyt]):key_irregular = 1b
	832	array_equal(glamt[firstxt:lastxt, firstyt:lastyt] $
	833	, glamt[firstxt:lastxt, firstyt]#replicate(1, nyt)):key_irregular = 1b
	834	array_equal(gphit[firstxt:lastxt, firstyt:lastyt] $
	835	, replicate(1, nxt)#(gphit[firstxt, firstyt:lastyt])[*]):key_irregular = 1b
	836	ELSE:key_irregular = 0b
	837	ENDCASE
	838	END
[469]	839	(where(gdtype eq 'U'))[0] NE -1 AND nxu NE 0 AND nyu NE 0:BEGIN
[13]	840	CASE 0 OF
	841	array_equal(glamu[firstxu:lastxu, firstyu], glamu[firstxu:lastxu, lastyu]):key_irregular = 1b
	842	array_equal(gphiu[firstxu, firstyu:lastyu], gphiu[lastxu, firstyu:lastyu]):key_irregular = 1b
	843	array_equal(glamu[firstxu:lastxu, firstyu:lastyu] $
	844	, glamu[firstxu:lastxu, firstyu]#replicate(1, nyu)):key_irregular = 1b
	845	array_equal(gphiu[firstxu:lastxu, firstyu:lastyu] $
	846	, replicate(1, nxu)#(gphiu[firstxu, firstyu:lastyu])[*]):key_irregular = 1b
	847	ELSE:key_irregular = 0b
	848	ENDCASE
	849	END
[469]	850	(where(gdtype eq 'V'))[0] NE -1 AND nxv NE 0 AND nyv NE 0:BEGIN
[13]	851	CASE 0 OF
	852	array_equal(glamv[firstxv:lastxv, firstyv], glamv[firstxv:lastxv, lastyv]):key_irregular = 1b
	853	array_equal(gphiv[firstxv, firstyv:lastyv], gphiv[lastxv, firstyv:lastyv]):key_irregular = 1b
	854	array_equal(glamv[firstxv:lastxv, firstyv:lastyv] $
	855	, glamv[firstxv:lastxv, firstyv]#replicate(1, nyv)):key_irregular = 1b
	856	array_equal(gphiv[firstxv:lastxv, firstyv:lastyv] $
	857	, replicate(1, nxv)#(gphiv[firstxv, firstyv:lastyv])[*]):key_irregular = 1b
	858	ELSE:key_irregular = 0b
	859	ENDCASE
	860	END
[469]	861	(where(gdtype eq 'F'))[0] NE -1 AND nxf NE 0 AND nyf NE 0:BEGIN
[13]	862	CASE 0 OF
	863	array_equal(glamf[firstxf:lastxf, firstyf], glamf[firstxf:lastxf, lastyf]):key_irregular = 1b
	864	array_equal(gphif[firstxf, firstyf:lastyf], gphif[lastxf, firstyf:lastyf]):key_irregular = 1b
	865	array_equal(glamf[firstxf:lastxf, firstyf:lastyf] $
	866	, glamf[firstxf:lastxf, firstyf]#replicate(1, nyf)):key_irregular = 1b
	867	array_equal(gphif[firstxf:lastxf, firstyf:lastyf] $
	868	, replicate(1, nxf)#(gphif[firstxf, firstyf:lastyf])[*]):key_irregular = 1b
	869	ELSE:key_irregular = 0b
	870	ENDCASE
	871	END
	872	ELSE:
	873	ENDCASE
	874	;
[2]	875	;-------------------------------------------------------------------
	876	;-------------------------------------------------------------------
[13]	877	; define all vertical parameters ...
	878	; vert1, vert2
	879	; firstz[tw], lastz[tw], nz[tw]
[2]	880	;-------------------------------------------------------------------
	881	;-------------------------------------------------------------------
[13]	882	;
	883	vertical:
	884	;
[2]	885	;-------------------------------------------------------------------
[13]	886	; vertical domain defined with vert1, vert2
[2]	887	;-------------------------------------------------------------------
[13]	888	IF NOT (keyword_set(zindex) OR keyword_set(index)) THEN BEGIN
[296]	889	; define vert1 and vert2
[13]	890	CASE N_PARAMS() OF
	891	2:vert1 = min([x1, x2], max = vert2)
	892	6:vert1 = min([z1, z2], max = vert2)
	893	ELSE:BEGIN
[469]	894	IF (inter(byte(gdtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN $
[13]	895	vert1t = min(gdept, max = vert2t)
[469]	896	IF (where(gdtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN $
[13]	897	vert1w = min(gdepw, max = vert2w)
	898	vert1 = min([vert1t, vert1w])
	899	vert2 = max([vert2t, vert2w])
	900	END
	901	ENDCASE
	902	; define firstzt, firstzt, nzt
[469]	903	IF (inter(byte(gdtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN
[13]	904	domz = where(gdept ge vert1 and gdept le vert2, nzt)
	905	IF nzt NE 0 THEN BEGIN
	906	firstzt = domz[0]
	907	lastzt = domz[nzt-1]
[226]	908	ENDIF ELSE BEGIN
[13]	909	ras = report('WARNING, The box does not contain any T level...')
	910	firstzt = -1
	911	lastzt = -1
[226]	912	ENDELSE
[13]	913	ENDIF
	914	; define firstzw, firstzw, nzw
[469]	915	IF (where(gdtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
[13]	916	IF keyword_set(memeindices) THEN BEGIN
	917	firstzw = firstzt
	918	lastzw = lastzt
	919	nzw = nzt
[226]	920	ENDIF ELSE BEGIN
[13]	921	domz = where(gdepw ge vert1 and gdepw le vert2, nzw)
	922	IF nzw NE 0 THEN BEGIN
	923	firstzw = domz[0]
	924	lastzw = domz[nzw-1]
[226]	925	ENDIF ELSE BEGIN
[13]	926	ras = report('WARNING, The box does not contain any W level...')
	927	firstzw = -1
	928	lastzw = -1
[226]	929	ENDELSE
	930	ENDELSE
[13]	931	ENDIF
[2]	932	;-------------------------------------------------------------------
[13]	933	; vertical domain defined with the Z index
[2]	934	;-------------------------------------------------------------------
[226]	935	ENDIF ELSE BEGIN
[13]	936	CASE N_PARAMS() OF
[393]	937	2:fstz = min(long([x1, x2]), max = lstz)
[13]	938	4:return
[393]	939	6:fstz = min(long([z1, z2]), max = lstz)
[13]	940	ENDCASE
	941	IF fstz LT 0 OR lstz GE jpk THEN BEGIN
	942	ras = report('Bad definition of X1, X2, Z1 or Z2')
	943	return
	944	ENDIF
	945	nz = lstz - fstz + 1
	946	; find vert1t, vert2t, firstzt, firstzt, nzt
	947	; according to (x1, x2) or (z1, z2)
[469]	948	IF (where(gdtype eq 'T'))[0] NE -1 THEN BEGIN
[13]	949	vert1t = min(gdept[fstz:lstz], max = vert2t)
	950	firstzt = fstz & lastzt = lstz & nzt = nz
[226]	951	ENDIF
[13]	952	; find vert1w, vert2w, firstzw, firstzw, nzw
	953	; according to (x1, x2) or (z1, z2)
[469]	954	IF (where(gdtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
[13]	955	vert1w = min(gdepw[fstz:lstz], max = vert2w)
	956	firstzw = fstz & lastzw = lstz & nzw = nz
[226]	957	ENDIF
[13]	958	vert1 = min([vert1t, vert1w])
	959	vert2 = max([vert2t, vert2w])
	960	ENDELSE
	961	;
[2]	962	;-------------------------------------------------------------------
[13]	963	IF NOT keyword_set(key_forgetold) THEN BEGIN
	964	@updateold
[226]	965	ENDIF
[2]	966	;-------------------------------------------------------------------
[226]	967	if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun
[2]	968	;-------------------------------------------------------------------
	969
	970	;-------------------------------------------------------------------
[13]	971	return
[2]	972	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format