Context Navigation

domdef.pro @ 240

Last change on this file since 240 was 240, checked in by pinsard, 17 years ago
replace some print by some report in some .pro (continuation)
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 48.0 KB

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

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

Download in other formats:

Original Format