Context Navigation

domdef.pro @ 493

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

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

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

Download in other formats:

Original Format