Context Navigation

domdef.pro @ 114

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

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

Download in other formats:

Original Format