Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

domdef.pro @ 25

Last change on this file since 25 was 13, checked in by pinsard, 18 years ago
upgrade of GRILLE/Utilities according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 35.2 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 N_PARAMS() EQ 2 THEN GOTO, vertical
174	;
175	;-------------------------------------------------------------------
176	;-------------------------------------------------------------------
177	; define all horizontal parameters ...
178	; lon1 et lon2
179	; lat1 et lat2
180	; firstx[tuvf], lastx[tuvf], nx[tuvf]
181	;-------------------------------------------------------------------
182	; check if the grid is defined for U and V points. If not, take care
183	; of the cases gridtype eq 'U' or 'V'
184	;-------------------------------------------------------------------
185	errstatus = 0
186	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
187	firstxu = !values.f_nan
188	lastxu = !values.f_nan
189	nxu = !values.f_nan
190	okgrid = where(gridtype NE 'U', count)
191	IF count NE 0 THEN gridtype = gridtype[okgrid] $
192	ELSE errstatus = report('U grid is undefined. Impossible to call domdef with vargid = ''U''')
193	ENDIF
194	;
195	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
196	firstxv = !values.f_nan
197	lastxv = !values.f_nan
198	nxv = !values.f_nan
199	okgrid = where(gridtype NE 'V', count)
200	IF count NE 0 THEN gridtype = gridtype[okgrid] $
201	ELSE errstatus = report('V grid is undefined. Impossible to call domdef with vargid = ''V''')
202	ENDIF
203	IF errstatus EQ -1 THEN return
204	;-------------------------------------------------------------------
205	;-------------------------------------------------------------------
206	; horizontal domain defined with lon1, lon2, lat1 and lat2
207	;-------------------------------------------------------------------
208	;-------------------------------------------------------------------
209	IF N_PARAMS() EQ 0 $
210	OR ( (N_PARAMS() EQ 4 OR N_PARAMS() EQ 6) $
211	AND NOT keyword_set(xindex) AND NOT keyword_set(yindex) AND NOT keyword_set(index) ) THEN BEGIN
212	IF N_PARAMS() EQ 0 THEN BEGIN
213	; find lon1 and lon2 the longitudinal boudaries of the full domain
214	IF (where(gridtype eq 'T'))[0] NE -1 THEN lon1t = min(glamt, max = lon2t)
215	IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lon1t = min(glamt, max = lon2t)
216	IF (where(gridtype eq 'U'))[0] NE -1 THEN lon1u = min(glamu, max = lon2u)
217	IF (where(gridtype eq 'V'))[0] NE -1 THEN lon1v = min(glamv, max = lon2v)
218	IF (where(gridtype eq 'F'))[0] NE -1 THEN lon1f = min(glamf, max = lon2f)
219	lon1 = min([lon1t, lon1u, lon1v, lon1f])
220	lon2 = max([lon2t, lon2u, lon2v, lon2f])
221	; find lat1 and lat2 the latitudinal boudaries of the full domain
222	IF (where(gridtype eq 'T'))[0] NE -1 THEN lat1t = min(gphit, max = lat2t)
223	IF (where(gridtype eq 'W'))[0] NE -1 AND (where(gridtype eq 'T'))[0] EQ -1 THEN lat1t = min(gphit, max = lat2t)
224	IF (where(gridtype eq 'U'))[0] NE -1 THEN lat1u = min(gphiu, max = lat2u)
225	IF (where(gridtype eq 'V'))[0] NE -1 THEN lat1v = min(gphiv, max = lat2v)
226	IF (where(gridtype eq 'F'))[0] NE -1 THEN lat1f = min(gphif, max = lat2f)
227	lat1 = min([lat1t, lat1u, lat1v, lat1f])
228	lat2 = max([lat2t, lat2u, lat2v, lat2f])
229	ENDIF ELSE BEGIN
230	lon1 = min([x1, x2], max = lon2)
231	lat1 = min([y1, y2], max = lat2)
232	ENDELSE
233	; find firstxt, firstxt, nxt and nyt according to lon1, lon2, lat1 and lat2
234	IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
235	dom = where( (glamt GE lon1) AND (glamt LE lon2) $
236	AND (gphit GE lat1) AND (gphit LE lat2) )
237	IF (dom[0] EQ -1) THEN BEGIN
238	IF keyword_set(findalways) THEN BEGIN
239	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
240	IF N_PARAMS() EQ 6 THEN $
241	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
242	RETURN
243	ENDIF ELSE BEGIN
244	ras = report('WARNING, The box does not contain any T points...')
245	firstxt = -1 & lastxt = -1 & nxt = 0
246	firstyt = -1 & lastyt = -1 & nyt = 0
247	ENDELSE
248	ENDIF ELSE BEGIN
249	jyt = dom / jpi
250	ixt = temporary(dom) MOD jpi
251	firstxt = min(temporary(ixt), max = lastxt)
252	firstyt = min(temporary(jyt), max = lastyt)
253	nxt = lastxt - firstxt + 1
254	nyt = lastyt - firstyt + 1
255	ENDELSE
256	ENDIF
257	; find firstxu, firstxu, firstyu, firstyu, nxu and nyu
258	; according to lon1, lon2, lat1 and lat2
259	IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
260	IF keyword_set(memeindices) THEN BEGIN
261	firstxu = firstxt & lastxu = lastxt & nxu = nxt
262	firstyu = firstyt & lastyu = lastyt & nyu = nyt
263	ENDIF ELSE BEGIN
264	dom = where( (glamu GE lon1) AND (glamu LE lon2) $
265	AND (gphiu GE lat1) AND (gphiu LE lat2) )
266	IF (dom[0] EQ -1) THEN BEGIN
267	IF keyword_set(findalways) THEN BEGIN
268	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
269	IF N_PARAMS() EQ 6 THEN $
270	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
271	RETURN
272	ENDIF ELSE BEGIN
273	ras = report('WARNING, The box does not contain any U points...')
274	firstxu = -1 & lastxu = -1 & nxu = 0
275	firstyu = -1 & lastyu = -1 & nyu = 0
276	ENDELSE
277	ENDIF ELSE BEGIN
278	jyu = dom / jpi
279	ixu = temporary(dom) MOD jpi
280	firstxu = min(temporary(ixu), max = lastxu)
281	firstyu = min(temporary(jyu), max = lastyu)
282	nxu = lastxu - firstxu + 1
283	nyu = lastyu - firstyu + 1
284	ENDELSE
285	ENDELSE
286	ENDIF
287	; find firstxv, firstxv, firstyv, firstyv, nxv and nyv
288	; according to lon1, lon2, lat1 and lat2
289	IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
290	IF keyword_set(memeindices) THEN BEGIN
291	firstxv = firstxt & lastxv = lastxt & nxv = nxt
292	firstyv = firstyt & lastyv = lastyt & nyv = nyt
293	ENDIF ELSE BEGIN
294	dom = where( (glamv GE lon1) AND (glamv LE lon2) $
295	AND (gphiv GE lat1) AND (gphiv LE lat2) )
296	IF (dom[0] EQ -1) THEN BEGIN
297	IF keyword_set(findalways) THEN BEGIN
298	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
299	IF N_PARAMS() EQ 6 THEN $
300	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
301	RETURN
302	ENDIF ELSE BEGIN
303	ras = report('WARNING, The box does not contain any V points...')
304	firstxv = -1 & lastxv = -1 & nxv = 0
305	firstyv = -1 & lastyv = -1 & nyv = 0
306	ENDELSE
307	ENDIF ELSE BEGIN
308	jyv = dom / jpi
309	ixv = temporary(dom) MOD jpi
310	firstxv = min(temporary(ixv), max = lastxv)
311	firstyv = min(temporary(jyv), max = lastyv)
312	nxv = lastxv - firstxv + 1
313	nyv = lastyv - firstyv + 1
314	ENDELSE
315	ENDELSE
316	ENDIF
317	; find firstxf, firstxf, firstyf, firstyf, nxf and nyf
318	; according to lon1, lon2, lat1 and lat2
319	IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN
320	IF keyword_set(memeindices) THEN BEGIN
321	firstxf = firstxt & lastxf = lastxt & nxf = nxt
322	firstyf = firstyt & lastyf = lastyt & nyf = nyt
323	ENDIF ELSE BEGIN
324	dom = where( (glamf GE lon1) AND (glamf LE lon2) $
325	AND (gphif GE lat1) AND (gphif LE lat2) )
326	IF (dom[0] EQ -1) THEN BEGIN
327	IF keyword_set(findalways) THEN BEGIN
328	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
329	IF N_PARAMS() EQ 6 THEN $
330	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
331	RETURN
332	ENDIF ELSE BEGIN
333	ras = report('WARNING, The box does not contain any F points...')
334	firstxf = -1 & lastxf = -1 & nxf = 0
335	firstyf = -1 & lastyf = -1 & nyf = 0
336	ENDELSE
337	ENDIF ELSE BEGIN
338	jyf = dom / jpi
339	ixf = temporary(dom) MOD jpi
340	firstxf = min(temporary(ixf), max = lastxf)
341	firstyf = min(temporary(jyf), max = lastyf)
342	nxf = lastxf - firstxf + 1
343	nyf = lastyf - firstyf + 1
344	ENDELSE
345	ENDELSE
346	ENDIF
347	;
348	ENDIF ELSE BEGIN
349	CASE 1 OF
350	;-------------------------------------------------------------------
351	;-------------------------------------------------------------------
352	; horizontal domain defined with the X and Y indexes
353	;-------------------------------------------------------------------
354	;-------------------------------------------------------------------
355	(keyword_set(xindex) AND keyword_set(yindex)) OR keyword_set(index):BEGIN
356	fstx = min([x1, x2], max = lstx)
357	fsty = min([y1, y2], max = lsty)
358	IF fstx LT 0 OR lstx GE jpi THEN BEGIN
359	ras = report('Bad definition of X1 or X2')
360	return
361	ENDIF
362	IF fsty LT 0 OR lsty GE jpj THEN BEGIN
363	ras = report('Bad definition of Y1 or Y2')
364	return
365	ENDIF
366	nx = lstx - fstx + 1
367	ny = lsty - fsty + 1
368	; find lon1t, lon2t, lat1t, lat2t, firstxt, firstxt, nxt and nyt
369	; according to x1, x2, y1, y2
370	IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1 THEN BEGIN
371	lon1t = min(glamt[fstx:lstx, fsty:lsty], max = lon2t)
372	lat1t = min(gphit[fstx:lstx, fsty:lsty], max = lat2t)
373	firstxt = fstx & lastxt = lstx
374	firstyt = fsty & lastyt = lsty
375	nxt = nx & nyt = ny
376	ENDIF
377	; find lon1u, lon2u, lat1u, lat2u, firstxu, firstxu, nxu and nyu
378	; according to x1, x2, y1, y2
379	IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
380	lon1u = min(glamu[fstx:lstx, fsty:lsty], max = lon2u)
381	lat1u = min(gphiu[fstx:lstx, fsty:lsty], max = lat2u)
382	firstxu = fstx & lastxu = lstx
383	firstyu = fsty & lastyu = lsty
384	nxu = nx & nyu = ny
385	ENDIF
386	; find lon1v, lon2v, lat1v, lat2v, firstxv, firstxv, nxv and nyv
387	; according to x1, x2, y1, y2
388	IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
389	lon1v = min(glamv[fstx:lstx, fsty:lsty], max = lon2v)
390	lat1v = min(gphiv[fstx:lstx, fsty:lsty], max = lat2v)
391	firstxv = fstx & lastxv = lstx
392	firstyv = fsty & lastyv = lsty
393	nxv = nx & nyv = ny
394	ENDIF
395	; find lon1f, lon2f, lat1f, lat2f, firstxf, firstxf, nxf and nyf
396	; according to x1, x2, y1, y2
397	IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN
398	lon1f = min(glamf[fstx:lstx, fsty:lsty], max = lon2f)
399	lat1f = min(gphif[fstx:lstx, fsty:lsty], max = lat2f)
400	firstxf = fstx & lastxf = lstx
401	firstyf = fsty & lastyf = lsty
402	nxf = nx & nyf = ny
403	ENDIF
404	lon1 = min([lon1t, lon1u, lon1v, lon1f])
405	lon2 = max([lon2t, lon2u, lon2v, lon2f])
406	lat1 = min([lat1t, lat1u, lat1v, lat1f])
407	lat2 = max([lat2t, lat2u, lat2v, lat2f])
408	END
409	;-------------------------------------------------------------------
410	;-------------------------------------------------------------------
411	; horizontal domain defined with the X index and lat1/lat2
412	;-------------------------------------------------------------------
413	;-------------------------------------------------------------------
414	keyword_set(xindex):BEGIN
415	fstx = min([x1, x2], max = lstx)
416	IF fstx LT 0 OR lstx GE jpi THEN BEGIN
417	ras = report('Bad definition of X1 or X2')
418	return
419	ENDIF
420	nx = lstx - fstx + 1
421	lat1 = min([y1, y2], max = lat2)
422	; find lon1t, lon2t, firstxt, firstxt, firstyt, firstyt, nxt
423	; and nyt according to x1, x2, lat1 and lat2
424	IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
425	firstxt = fstx & lastxt = lstx & nxt = nx
426	dom = where( (gphit[fstx:lstx, ] GE lat1) AND (gphit[fstx:lstx, ] LE lat2) )
427	IF (dom[0] EQ -1) THEN BEGIN
428	IF keyword_set(findalways) THEN BEGIN
429	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
430	IF N_PARAMS() EQ 6 THEN $
431	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
432	RETURN
433	ENDIF ELSE BEGIN
434	ras = report('WARNING, The box does not contain any T points...')
435	firstyt = -1 & lastyt = -1 & nyt = 0
436	ENDELSE
437	ENDIF ELSE BEGIN
438	jyt = temporary(dom) / nx
439	firstyt = min(temporary(jyt), max = lastyt)
440	nyt = lastyt - firstyt + 1
441	ENDELSE
442	IF nyt NE 0 THEN lon1t = min(glamt[firstxt:lastxt, firstyt:lastyt], max = lon2t)
443	ENDIF
444	; find lon1u, lon2u, firstxu, firstxu, firstyu, firstyu, nxu
445	; and nyu according to x1, x2, lat1 and lat2
446	IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
447	firstxu = fstx & lastxu = lstx & nxu = nx
448	IF keyword_set(memeindices) THEN BEGIN
449	firstyu = firstyt & lastyu = lastyt & nyu = nyt
450	ENDIF ELSE BEGIN
451	dom = where( (gphiu[fstx:lstx, ] GE lat1) AND (gphiu[fstx:lstx, ] LE lat2) )
452	IF (dom[0] EQ -1) THEN BEGIN
453	IF keyword_set(findalways) THEN BEGIN
454	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
455	IF N_PARAMS() EQ 6 THEN $
456	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
457	RETURN
458	ENDIF ELSE BEGIN
459	ras = report('WARNING, The box does not contain any U points...')
460	firstyu = -1 & lastyu = -1 & nyu = 0
461	ENDELSE
462	ENDIF ELSE BEGIN
463	jyu = temporary(dom) / nx
464	firstyu = min(temporary(jyu), max = lastyu)
465	nyu = lastyu - firstyu + 1
466	ENDELSE
467	ENDELSE
468	IF nyu NE 0 THEN lon1u = min(glamu[firstxu:lastxu, firstyu:lastyu], max = lon2u)
469	ENDIF
470	; find lon1v, lon2v, firstxv, firstxv, firstyv, firstyv, nxv
471	; and nyv according to x1, x2, lat1 and lat2
472	IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
473	firstxv = fstx & lastxv = lstx & nxv = nx
474	IF keyword_set(memeindices) THEN BEGIN
475	firstyv = firstyt & lastyv = lastyt & nyv = nyt
476	ENDIF ELSE BEGIN
477	dom = where( (gphiv[fstx:lstx, ] GE lat1) AND (gphiv[fstx:lstx, ] LE lat2) )
478	IF (dom[0] EQ -1) THEN BEGIN
479	IF keyword_set(findalways) THEN BEGIN
480	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
481	IF N_PARAMS() EQ 6 THEN $
482	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
483	RETURN
484	ENDIF ELSE BEGIN
485	ras = report('WARNING, The box does not contain any V points...')
486	firstyv = -1 & lastyv = -1 & nyv = 0
487	ENDELSE
488	ENDIF ELSE BEGIN
489	jyv = temporary(dom) / nx
490	firstyv = min(temporary(jyv), max = lastyv)
491	nyv = lastyv - firstyv + 1
492	ENDELSE
493	ENDELSE
494	IF nyv NE 0 THEN lon1v = min(glamv[firstxv:lastxv, firstyv:lastyv], max = lon2v)
495	ENDIF
496	; find lon1f, lon2f, firstxf, firstxf, firstyf, firstyf, nxf
497	; and nyf according to x1, x2, lat1 and lat2
498	IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN
499	firstxf = fstx & lastxf = lstx & nxf = nx
500	IF keyword_set(memeindices) THEN BEGIN
501	firstyf = firstyt & lastyf = lastyt & nyf = nyt
502	ENDIF ELSE BEGIN
503	dom = where( (gphif[fstx:lstx, ] GE lat1) AND (gphif[fstx:lstx, ] LE lat2) )
504	IF (dom[0] EQ -1) THEN BEGIN
505	IF keyword_set(findalways) THEN BEGIN
506	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
507	IF N_PARAMS() EQ 6 THEN $
508	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
509	RETURN
510	ENDIF ELSE BEGIN
511	ras = report('WARNING, The box does not contain any F points...')
512	firstyf = -1 & lastyf = -1 & nyf = 0
513	ENDELSE
514	ENDIF ELSE BEGIN
515	jyf = temporary(dom) / nx
516	firstyf = min(temporary(jyf), max = lastyf)
517	nyf = lastyf - firstyf + 1
518	ENDELSE
519	ENDELSE
520	IF nyf NE 0 THEN lon1f = min(glamf[firstxf:lastxf, firstyf:lastyf], max = lon2f)
521	ENDIF
522	lon1 = min([lon1t, lon1u, lon1v, lon1f])
523	lon2 = max([lon2t, lon2u, lon2v, lon2f])
524	END
525	;-------------------------------------------------------------------
526	;-------------------------------------------------------------------
527	; horizontal domain defined with the Y index and lon1/lon2
528	;-------------------------------------------------------------------
529	;-------------------------------------------------------------------
530	keyword_set(yindex):BEGIN
531	fsty = min([y1, y2], max = lsty)
532	IF fsty LT 0 OR lsty GE jpj THEN BEGIN
533	ras = report('Bad definition of Y1 or Y2')
534	return
535	ENDIF
536	ny = lsty - fsty + 1
537	lon1 = min([x1, x2], max = lon2)
538	; find lat1t, lat2t, firstxt, firstxt, firstyt, firstyt, nxt
539	; and nyt according to x1, x2, lon1 and lon2
540	IF (where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype EQ 'W'))[0] NE -1 THEN BEGIN
541	firstyt = fsty & lastyt = lsty & nyt = ny
542	dom = where( (glamt[, fsty:lsty] GE lon1) AND (glamt[, fsty:lsty] LE lon2) )
543	IF (dom[0] EQ -1) THEN BEGIN
544	IF keyword_set(findalways) THEN BEGIN
545	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
546	IF N_PARAMS() EQ 6 THEN $
547	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
548	RETURN
549	ENDIF ELSE BEGIN
550	ras = report('WARNING, The box does not contain any T points...')
551	firstxt = -1 & lastxt = -1 & nxt = 0
552	ENDELSE
553	ENDIF ELSE BEGIN
554	jxt = temporary(dom) MOD jpi
555	firstxt = min(temporary(jxt), max = lastxt)
556	nxt = lastxt - firstxt + 1
557	ENDELSE
558	IF nxt NE 0 THEN lat1t = min(gphit[firstxt:lastxt, firstyt:lastyt], max = lat2t)
559	ENDIF
560	; find lat1u, lat2u, firstxu, firstxu, firstyu, firstyu, nxu
561	; and nyu according to x1, x2, lon1 and lon2
562	IF (where(gridtype eq 'U'))[0] NE -1 THEN BEGIN
563	firstyu = fsty & lastyu = lsty & nyu = ny
564	IF keyword_set(memeindices) THEN BEGIN
565	firstxu = firstyt & lastxu = lastyt & nxu = nxt
566	ENDIF ELSE BEGIN
567	dom = where( (glamu[, fsty:lsty] GE lon1) AND (glamu[, fsty:lsty] LE lon2) )
568	IF (dom[0] EQ -1) THEN BEGIN
569	IF keyword_set(findalways) THEN BEGIN
570	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
571	IF N_PARAMS() EQ 6 THEN $
572	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
573	RETURN
574	ENDIF ELSE BEGIN
575	ras = report('WARNING, The box does not contain any U points...')
576	firstxu = -1 & lastxu = -1 & nxu = 0
577	ENDELSE
578	ENDIF ELSE BEGIN
579	jxu = temporary(dom) MOD jpi
580	firstxu = min(temporary(jxu), max = lastxu)
581	nxu = lastxu - firstxu + 1
582	ENDELSE
583	ENDELSE
584	IF nxu NE 0 THEN lat1u = min(gphiu[firstxu:lastxu, firstyu:lastyu], max = lat2u)
585	ENDIF
586	; find lat1v, lat2v, firstxv, firstxv, firstyv, firstyv, nxv
587	; and nyv according to x1, x2, lon1 and lon2
588	IF (where(gridtype eq 'V'))[0] NE -1 THEN BEGIN
589	firstyv = fsty & lastyv = lsty & nyv = ny
590	IF keyword_set(memeindices) THEN BEGIN
591	firstxv = firstyt & lastxv = lastyt & nxv = nxt
592	ENDIF ELSE BEGIN
593	dom = where( (glamv[, fsty:lsty] GE lon1) AND (glamv[, fsty:lsty] LE lon2) )
594	IF (dom[0] EQ -1) THEN BEGIN
595	IF keyword_set(findalways) THEN BEGIN
596	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
597	IF N_PARAMS() EQ 6 THEN $
598	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
599	RETURN
600	ENDIF ELSE BEGIN
601	ras = report('WARNING, The box does not contain any V points...')
602	firstxv = -1 & lastxv = -1 & nxv = 0
603	ENDELSE
604	ENDIF ELSE BEGIN
605	jxv = temporary(dom) MOD jpi
606	firstxv = min(temporary(jxv), max = lastxv)
607	nxv = lastxv - firstxv + 1
608	ENDELSE
609	ENDELSE
610	IF nxv NE 0 THEN lat1v = min(gphiv[firstxv:lastxv, firstyv:lastyv], max = lat2v)
611	ENDIF
612	; find lat1f, lat2f, firstxf, firstxf, firstyf, firstyf, nxf
613	; and nyf according to x1, x2, lon1 and lon2
614	IF (where(gridtype eq 'F'))[0] NE -1 THEN BEGIN
615	firstyf = fsty & lastyf = lsty & nyf = ny
616	IF keyword_set(memeindices) THEN BEGIN
617	firstxf = firstyt & lastxf = lastyt & nxf = nxt
618	ENDIF ELSE BEGIN
619	dom = where( (glamf[, fsty:lsty] GE lon1) AND (glamf[, fsty:lsty] LE lon2) )
620	IF (dom[0] EQ -1) THEN BEGIN
621	IF keyword_set(findalways) THEN BEGIN
622	domdef, GRIDTYPE = gridtype, MEMEINDICES = memeindices, _extra = ex
623	IF N_PARAMS() EQ 6 THEN $
624	domdef, z1, z2, GRIDTYPE = gridtype, MEMEINDICES = memeindices, ZINDEX = zindex, _extra = ex
625	RETURN
626	ENDIF ELSE BEGIN
627	ras = report('WARNING, The box does not contain any F points...')
628	firstxf = -1 & lastyf = -1 & nxf = 0
629	ENDELSE
630	ENDIF ELSE BEGIN
631	jxf = temporary(dom) MOD jpi
632	firstxf = min(temporary(jxf), max = lastxf)
633	nxf = lastxf - firstxf + 1
634	ENDELSE
635	ENDELSE
636	IF nxf NE 0 THEN lat1f = min(gphif[firstxf:lastxf, firstyf:lastyf], max = lat2f)
637	ENDIF
638	lat1 = min([lat1t, lat1u, lat1v, lat1f])
639	lat2 = max([lat2t, lat2u, lat2v, lat2f])
640	END
641	ENDCASE
642	ENDELSE
643	;-------------------------------------------------------------------
644	; The extracted domain is it regular or not?
645	;-------------------------------------------------------------------
646	CASE 1 OF
647	((where(gridtype eq 'T'))[0] NE -1 OR (where(gridtype eq 'W'))[0] NE -1) AND nxt NE 0 AND nyt NE 0:BEGIN
648	; to get faster, we first test the most basic cases befor testing the
649	; full array.
650	CASE 0 OF
651	array_equal(glamt[firstxt:lastxt, firstyt], glamt[firstxt:lastxt, lastyt]):key_irregular = 1b
652	array_equal(gphit[firstxt, firstyt:lastyt], gphit[lastxt, firstyt:lastyt]):key_irregular = 1b
653	array_equal(glamt[firstxt:lastxt, firstyt:lastyt] $
654	, glamt[firstxt:lastxt, firstyt]#replicate(1, nyt)):key_irregular = 1b
655	array_equal(gphit[firstxt:lastxt, firstyt:lastyt] $
656	, replicate(1, nxt)#(gphit[firstxt, firstyt:lastyt])[*]):key_irregular = 1b
657	ELSE:key_irregular = 0b
658	ENDCASE
659	END
660	(where(gridtype eq 'U'))[0] NE -1 AND nxu NE 0 AND nyu NE 0:BEGIN
661	CASE 0 OF
662	array_equal(glamu[firstxu:lastxu, firstyu], glamu[firstxu:lastxu, lastyu]):key_irregular = 1b
663	array_equal(gphiu[firstxu, firstyu:lastyu], gphiu[lastxu, firstyu:lastyu]):key_irregular = 1b
664	array_equal(glamu[firstxu:lastxu, firstyu:lastyu] $
665	, glamu[firstxu:lastxu, firstyu]#replicate(1, nyu)):key_irregular = 1b
666	array_equal(gphiu[firstxu:lastxu, firstyu:lastyu] $
667	, replicate(1, nxu)#(gphiu[firstxu, firstyu:lastyu])[*]):key_irregular = 1b
668	ELSE:key_irregular = 0b
669	ENDCASE
670	END
671	(where(gridtype eq 'V'))[0] NE -1 AND nxv NE 0 AND nyv NE 0:BEGIN
672	CASE 0 OF
673	array_equal(glamv[firstxv:lastxv, firstyv], glamv[firstxv:lastxv, lastyv]):key_irregular = 1b
674	array_equal(gphiv[firstxv, firstyv:lastyv], gphiv[lastxv, firstyv:lastyv]):key_irregular = 1b
675	array_equal(glamv[firstxv:lastxv, firstyv:lastyv] $
676	, glamv[firstxv:lastxv, firstyv]#replicate(1, nyv)):key_irregular = 1b
677	array_equal(gphiv[firstxv:lastxv, firstyv:lastyv] $
678	, replicate(1, nxv)#(gphiv[firstxv, firstyv:lastyv])[*]):key_irregular = 1b
679	ELSE:key_irregular = 0b
680	ENDCASE
681	END
682	(where(gridtype eq 'F'))[0] NE -1 AND nxf NE 0 AND nyf NE 0:BEGIN
683	CASE 0 OF
684	array_equal(glamf[firstxf:lastxf, firstyf], glamf[firstxf:lastxf, lastyf]):key_irregular = 1b
685	array_equal(gphif[firstxf, firstyf:lastyf], gphif[lastxf, firstyf:lastyf]):key_irregular = 1b
686	array_equal(glamf[firstxf:lastxf, firstyf:lastyf] $
687	, glamf[firstxf:lastxf, firstyf]#replicate(1, nyf)):key_irregular = 1b
688	array_equal(gphif[firstxf:lastxf, firstyf:lastyf] $
689	, replicate(1, nxf)#(gphif[firstxf, firstyf:lastyf])[*]):key_irregular = 1b
690	ELSE:key_irregular = 0b
691	ENDCASE
692	END
693	ELSE:
694	ENDCASE
695	;
696	;-------------------------------------------------------------------
697	;-------------------------------------------------------------------
698	; define all vertical parameters ...
699	; vert1, vert2
700	; firstz[tw], lastz[tw], nz[tw]
701	;-------------------------------------------------------------------
702	;-------------------------------------------------------------------
703	;
704	vertical:
705	;
706	;-------------------------------------------------------------------
707	; vertical domain defined with vert1, vert2
708	;-------------------------------------------------------------------
709	IF NOT (keyword_set(zindex) OR keyword_set(index)) THEN BEGIN
710	; define vert1 et vert2
711	CASE N_PARAMS() OF
712	2:vert1 = min([x1, x2], max = vert2)
713	6:vert1 = min([z1, z2], max = vert2)
714	ELSE:BEGIN
715	IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN $
716	vert1t = min(gdept, max = vert2t)
717	IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN $
718	vert1w = min(gdepw, max = vert2w)
719	vert1 = min([vert1t, vert1w])
720	vert2 = max([vert2t, vert2w])
721	END
722	ENDCASE
723	; define firstzt, firstzt, nzt
724	IF (inter(byte(gridtype), byte(['T', 'U', 'V', 'F'])))[0] NE -1 THEN BEGIN
725	domz = where(gdept ge vert1 and gdept le vert2, nzt)
726	IF nzt NE 0 THEN BEGIN
727	firstzt = domz[0]
728	lastzt = domz[nzt-1]
729	ENDIF ELSE BEGIN
730	ras = report('WARNING, The box does not contain any T level...')
731	firstzt = -1
732	lastzt = -1
733	ENDELSE
734	ENDIF
735	; define firstzw, firstzw, nzw
736	IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
737	IF keyword_set(memeindices) THEN BEGIN
738	firstzw = firstzt
739	lastzw = lastzt
740	nzw = nzt
741	ENDIF ELSE BEGIN
742	domz = where(gdepw ge vert1 and gdepw le vert2, nzw)
743	IF nzw NE 0 THEN BEGIN
744	firstzw = domz[0]
745	lastzw = domz[nzw-1]
746	ENDIF ELSE BEGIN
747	ras = report('WARNING, The box does not contain any W level...')
748	firstzw = -1
749	lastzw = -1
750	ENDELSE
751	ENDELSE
752	ENDIF
753	;-------------------------------------------------------------------
754	; vertical domain defined with the Z index
755	;-------------------------------------------------------------------
756	ENDIF ELSE BEGIN
757	CASE N_PARAMS() OF
758	2:fstz = min([x1, x2], max = lstz)
759	4:return
760	6:fstz = min([z1, z2], max = lstz)
761	ENDCASE
762	IF fstz LT 0 OR lstz GE jpk THEN BEGIN
763	ras = report('Bad definition of X1, X2, Z1 or Z2')
764	return
765	ENDIF
766	nz = lstz - fstz + 1
767	; find vert1t, vert2t, firstzt, firstzt, nzt
768	; according to (x1, x2) or (z1, z2)
769	IF (where(gridtype eq 'T'))[0] NE -1 THEN BEGIN
770	vert1t = min(gdept[fstz:lstz], max = vert2t)
771	firstzt = fstz & lastzt = lstz & nzt = nz
772	ENDIF
773	; find vert1w, vert2w, firstzw, firstzw, nzw
774	; according to (x1, x2) or (z1, z2)
775	IF (where(gridtype eq 'W'))[0] NE -1 AND n_elements(gdepw) NE 0 THEN BEGIN
776	vert1w = min(gdepw[fstz:lstz], max = vert2w)
777	firstzw = fstz & lastzw = lstz & nzw = nz
778	ENDIF
779	vert1 = min([vert1t, vert1w])
780	vert2 = max([vert2t, vert2w])
781	ENDELSE
782	;
783	;-------------------------------------------------------------------
784	IF NOT keyword_set(key_forgetold) THEN BEGIN
785	@updateold
786	ENDIF
787	;-------------------------------------------------------------------
788	if keyword_set(key_performance) THEN print, 'temps domdef', systime(1)-tempsun
789	;-------------------------------------------------------------------
790
791	;-------------------------------------------------------------------
792	return
793	end

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

Download in other formats:

Original Format