Context Navigation

domdef.pro @ 261

Last change on this file since 261 was 254, checked in by pinsard, 17 years ago
improvements/corrections of some *.pro headers
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 48.1 KB

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

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

Download in other formats:

Original Format