Context Navigation

domdef.pro @ 157

Last change on this file since 157 was 157, checked in by navarro, 18 years ago
header improvements + xxx doc
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 49.1 KB

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

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

Download in other formats:

Original Format