Context Navigation

domdef.pro @ 268

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

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

Download in other formats:

Original Format