Context Navigation

domdef.pro @ 76

Last change on this file since 76 was 74, checked in by smasson, 18 years ago
debug xxx and cie + clean data file + rm perldoc_idl
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 48.8 KB

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

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

Download in other formats:

Original Format