Context Navigation

domdef.pro @ 205

Last change on this file since 205 was 163, checked in by navarro, 18 years ago
header improvements : type of parameters and keywords, default values, spell checking + idldoc assistant (IDL online_help)
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 48.3 KB

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

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

Download in other formats:

Original Format