Context Navigation

norme.pro @ 2

Last change on this file since 2 was 2, checked in by opalod, 22 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 21.6 KB

Line
1	;------------------------------------------------------------
2	;------------------------------------------------------------
3	;------------------------------------------------------------
4	;+
5	; NAME:norme
6	;
7	; PURPOSE: calcule la norme d''un champ de vecteurs, puis fait une
8	; moyenne eventuelle.
9	; Rq1: le champ de vecteur peut etre, 2d:xy, 3d: xyz ou xyt,
10	; 4d: xyzt
11	; Rq2: le calcul de la norme est fait avant l''eventuelle
12	; moyenne spatiale ou temporelle car la moyenne de la norme
13	; n''est pas egale a la norme des moyennes.
14	;
15	; CATEGORY: calcul de post traitement
16	;
17	; CALLING SEQUENCE:res=norme(champ_de_vecteurs)
18	;
19	; INPUTS:un tableau 2d, 3d ou 4d
20	;
21	; KEYWORD PARAMETERS:
22	;
23	; BOITE: boite sur laquelle moyenner (par defaut le domaine
24	; selectionner par le dernier domdef effectue)
25	;
26	; DIREC:'t' 'x' 'y' 'z' 'xy' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt'
27	; 'xzt' 'yzt' 'xyzt' directions selon lesquelles effectuer les
28	; moyennes
29	;
30	; OUTPUTS:tableau a tracer avec plt, pltz ou pltt.
31	;
32	; COMMON BLOCKS:
33	; common.pro
34	;
35	; SIDE EFFECTS:
36	;
37	; La norme est calculee aux points T. Pour faire ce calcul, on
38	; moyenne les champs U et V aux points T avant de calculer la norme.
39	; Au bord des cotes et du domaine, on ne peut pas calculer les
40	; champs U et V aux points T, ces points sont donc a la valeur
41	; !values.f_nan.
42	;
43	; lorsqu''on fait le calcul sur un domaine geographique reduit,
44	; les champs U et V ne comprennent pas forcement le meme nombre
45	; de points. Dans ce cas on redecoupe U et V pour ne garder que
46	; les points en commun. Au passage on refait un domdef qui
47	; redefinit un domaine geographique sur lequel les champs U et V
48	; sont extraits sur les meme points.
49	;
50	; RESTRICTIONS:
51	;
52	; pour savoir a quel type de tableau on a a faire, on teste la
53	; taille de celui-ci et les dates donnees par time[0] et
54	; time[jpt-1] pour savoir si il y a une dimension
55	; temporelle. Avant de lancer norme s''assurer que time et jpt
56	; sont bien definis comme il faut!
57	;
58	; EXAMPLE:
59	;
60	; pour calculer la moyenne de la norme des courants sur tout le
61	; dommaine entre 0 et 50:
62	; IDL> res=norme(un,vn,boite=[0,50],dir='xyz')
63	;
64	; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
65	; 9/6/1999
66	;-
67	;------------------------------------------------------------
68	;------------------------------------------------------------
69	;------------------------------------------------------------
70	FUNCTION norme,composanteu, composantev, BOITE=boite,DIREC=direc
71	@common
72	tempsun = systime(1) ; pour key_performance
73	;------------------------------------------------------------
74	if NOT keyword_set(direc) then direc = 0
75	; construction de u et v aux pts T
76	u = litchamp(composanteu)
77	v = litchamp(composantev)
78	date1 = time[0]
79	if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
80
81	if (size(u))[0] NE (size(v))[0] then return, -1
82
83	vargrid='T'
84	varname = 'norme '
85	valmask = 1e20
86	;
87	grilleu = litchamp(composanteu, /grid)
88	if grilleu EQ '' then grilleu = 'U'
89	grillev = litchamp(composantev, /grid)
90	if grillev EQ '' then grillev = 'V'
91	IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1
92	IF grilleu EQ 'T' AND grillev EQ 'T' THEN BEGIN
93	interpolle = 0
94	return, report('cas non code mais facile a faire!')
95	ENDIF ELSE interpolle = 1
96	if keyword_set(inverse) then begin
97	rien = u
98	u = v
99	v = rien
100	endif
101
102
103	;------------------------------------------------------------
104	; on trouve les points que u et v ont en communs
105	;------------------------------------------------------------
106	indicexu = (lindgen(jpi))[premierxu:premierxu+nxu-1]
107	indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1]
108	indicex = inter(indicexu, indicexv)
109	indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1]
110	indiceyv = (lindgen(jpj))[premieryv:premieryv+nyv-1]
111	indicey = inter(indiceyu, indiceyv)
112	nx = n_elements(indicex)
113	ny = n_elements(indicey)
114	;----------------------------------------------------------------------------
115	case 1 of
116	;----------------------------------------------------------------------------
117	;----------------------------------------------------------------------------
118	;xyz
119	;----------------------------------------------------------------------------
120	;----------------------------------------------------------------------------
121	(size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN
122	;----------------------------------------------------------------------------
123	indice2d = lindgen(jpi, jpj)
124	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
125	indice3d = lindgen(jpi, jpj, jpk)
126	indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,premierzt:dernierzt]
127	;------------------------------------------------------------
128	; extraction de u et v sur le domaine qui convient
129	;------------------------------------------------------------
130	case 1 of
131	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
132	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
133	case (size(u))[3] OF
134	nzt:BEGIN
135	if nxu NE nx then $
136	if indicex[0] EQ premierxu then u = u[0:nx-1,,] ELSE u = u[1: nx, ,]
137	IF nxv NE nx THEN $
138	if indicex[0] EQ premierxv then v = v[0:nx-1,,] ELSE v = v[1: nx, ,]
139	IF nyu NE ny THEN $
140	if indicey[0] EQ premieryu then u = u[,0:ny-1,] ELSE u = u[, 1: ny,]
141	IF nyv NE ny THEN $
142	if indicey[0] EQ premieryv then v = v[,0:ny-1,] ELSE v = v[, 1: ny,]
143	end
144	jpk:BEGIN
145	if nxu NE nx then $
146	if indicex[0] EQ premierxu then u = u[0:nx-1, ,premierzt:dernierzt] ELSE u = u[1: nx, ,premierzt:dernierzt]
147	IF nxv NE nx THEN $
148	if indicex[0] EQ premierxv then v = v[0:nx-1, ,premierzt:dernierzt] ELSE v = v[1: nx, ,premierzt:dernierzt]
149	IF nyu NE ny THEN $
150	if indicey[0] EQ premieryu then u = u[, 0:ny-1,premierzt:dernierzt] ELSE u = u[, 1: ny,premierzt:dernierzt]
151	IF nyv NE ny THEN $
152	if indicey[0] EQ premieryv then v = v[, 0:ny-1,premierzt:dernierzt] ELSE v = v[, 1: ny,premierzt:dernierzt]
153	end
154	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
155	endcase
156	END
157	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
158	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN
159	u = u[indice3d]
160	v = v[indice3d]
161	END
162	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
163	endcase
164	;------------------------------------------------------------------
165	; on reform u et v pour s'assurer qu'aucune dimension n'a ete ecrasee
166	;------------------------------------------------------------------
167	if nzt EQ 1 then begin
168	u = reform(u, nx, ny, nzt, /over)
169	v = reform(v, nx, ny, nzt, /over)
170	endif
171	;------------------------------------------------------------------
172	; construction de u et v aux pts T
173	;-----------------------------------------------------------
174	a=u(0,,)
175	u=(u+shift(u,1,0,0))/2.
176	if NOT keyword_set(key_periodique) OR nx NE jpi then u(0,,)=a
177	a=v(,0,)
178	v=(v+shift(v,0,1,0))/2.
179	if NOT keyword_set(key_periodique) OR nx NE jpi then v(,0,)=a
180	;----------------------------------------------------------------------------
181	; attribution du mask et des tableau de longitude et latitude
182	;----------------------------------------------------------------------------
183	mask = tmask[indice3d]
184	if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
185	;-----------------------------------------------------------
186	if n_elements(valmask) EQ 0 THEN valmask = 1e20
187	landu = where(u GE valmask/10)
188	if landu[0] NE -1 then u[landu] = 0
189	landv = where(v GE valmask/10)
190	if landv[0] NE -1 then v[landv] = 0
191	res=sqrt(u^2+v^2)
192	if NOT keyword_set(key_periodique) OR nx NE jpi then res(0,, )=!values.f_nan
193	res(,0, )=!values.f_nan
194	mask = where(mask eq 0)
195	IF mask[0] NE -1 THEN res(mask) = valmask
196	; moyennes en tous genres
197	oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
198	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme
199	if keyword_set(direc) then res = moyenne(res,direc,/nan, boite = boite)
200	;-----------------------------------------------------------
201	END
202	;----------------------------------------------------------------------------
203	;----------------------------------------------------------------------------
204	;xyt
205	;----------------------------------------------------------------------------
206	;----------------------------------------------------------------------------
207	date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN
208	indice2d = lindgen(jpi, jpj)
209	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
210	;------------------------------------------------------------
211	; extraction de u et v sur le domaine qui convient
212	;------------------------------------------------------------
213	case 1 of
214	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
215	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
216	if nxu NE nx then $
217	if indicex[0] EQ premierxu then u = u[0:nx-1, , ] ELSE u = u[1: nx, , ]
218	IF nxv NE nx THEN $
219	if indicex[0] EQ premierxv then v = v[0:nx-1, , ] ELSE v = v[1: nx, , ]
220	IF nyu NE ny THEN $
221	if indicey[0] EQ premieryu then u = u[, 0:ny-1, ] ELSE u = u[, 1: ny, ]
222	IF nyv NE ny THEN $
223	if indicey[0] EQ premieryv then v = v[, 0:ny-1, ] ELSE v = v[, 1: ny, ]
224	END
225	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
226	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
227	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
228	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
229	END
230	ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
231	endcase
232	;------------------------------------------------------------------
233	; construction de u et v aux pts T
234	;-----------------------------------------------------------
235	a=u(0,,)
236	u=(u+shift(u,1,0,0))/2.
237	if NOT keyword_set(key_periodique) OR nx NE jpi then u(0,,)=a
238	a=v(,0,)
239	v=(v+shift(v,0,1,0))/2.
240	if NOT keyword_set(key_periodique) OR nx NE jpi then v(,0,)=a
241	;----------------------------------------------------------------------------
242	; attribution du mask et des tableau de longitude et latitude
243	; on recupere la grille complette pour etablir un grand mask etendu ds les 4
244	; directions pour couvrir les points pour lesquels un pt terre a ete pris en
245	; compte (faire un petit dessin...)
246	;----------------------------------------------------------------------------
247	mask = tmask[indice2d+jpijpj(niveau-1)]
248	if ny EQ 1 then mask = reform(mask, nx, ny, /over)
249	;-----------------------------------------------------------
250	; construction de terre qui contient tous les point a masquer
251	;-----------------------------------------------------------
252	if n_elements(valmask) EQ 0 THEN valmask = 1e20
253	landu = where(u GE valmask/10)
254	if landu[0] NE -1 then u[landu] = 0
255	landv = where(v GE valmask/10)
256	if landv[0] NE -1 then v[landv] = 0
257	res=sqrt(u^2+v^2)
258	if NOT keyword_set(key_periodique) OR nx NE jpi then res(0,, )=!values.f_nan
259	res(,0, )=!values.f_nan
260	mask = where(mask eq 0)
261	IF mask[0] NE -1 THEN BEGIN
262	coeftps = lindgen(jpt)nxny
263	coeftps = replicate(1, n_elements(mask))#coeftps
264	mask = (temporary(mask))[*]#replicate(1, jpt)
265	mask =temporary(mask[]) + temporary(coeftps[])
266	res(temporary(mask)) = valmask
267	ENDIF
268	; moyennes en tous genres
269	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme
270	if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boite = boite)
271	END
272	;----------------------------------------------------------------------------
273	;----------------------------------------------------------------------------
274	;xyzt
275	;----------------------------------------------------------------------------
276	;----------------------------------------------------------------------------
277	date1 NE date2 AND (size(u))[0] EQ 4:BEGIN
278	indice2d = lindgen(jpi, jpj)
279	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
280	indice3d = lindgen(jpi, jpj, jpk)
281	indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,premierzt:dernierzt]
282	;------------------------------------------------------------
283	; extraction de u et v sur le domaine qui convient
284	;------------------------------------------------------------
285	case 1 of
286	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
287	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
288	case (size(u))[3] OF
289	nzt:BEGIN
290	if nxu NE nx then $
291	if indicex[0] EQ premierxu then u = u[0:nx-1,,,] ELSE u = u[1: nx, ,,]
292	IF nxv NE nx THEN $
293	if indicex[0] EQ premierxv then v = v[0:nx-1,,,] ELSE v = v[1: nx, ,,]
294	IF nyu NE ny THEN $
295	if indicey[0] EQ premieryu then u = u[,0:ny-1,,] ELSE u = u[, 1: ny,,]
296	IF nyv NE ny THEN $
297	if indicey[0] EQ premieryv then v = v[,0:ny-1,,] ELSE v = v[, 1: ny,,]
298	end
299	jpk:BEGIN
300	if nxu NE nx then $
301	if indicex[0] EQ premierxu then u = u[0:nx-1, ,premierzt:dernierzt,] ELSE u = u[1: nx, ,premierzt:dernierzt,]
302	IF nxv NE nx THEN $
303	if indicex[0] EQ premierxv then v = v[0:nx-1, ,premierzt:dernierzt,] ELSE v = v[1: nx, ,premierzt:dernierzt,]
304	IF nyu NE ny THEN $
305	if indicey[0] EQ premieryu then u = u[, 0:ny-1,premierzt:dernierzt,] ELSE u = u[, 1: ny,premierzt:dernierzt,]
306	IF nyv NE ny THEN $
307	if indicey[0] EQ premieryv then v = v[, 0:ny-1,premierzt:dernierzt,] ELSE v = v[, 1: ny,premierzt:dernierzt,]
308	end
309	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
310	endcase
311	END
312	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
313	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN
314	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,premierzt:dernierzt, *]
315	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,premierzt:dernierzt, *]
316	END
317	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
318	endcase
319	;------------------------------------------------------------------
320	; construction de u et v aux pts T
321	;-----------------------------------------------------------
322	a=u(0,,,*)
323	u=(u+shift(u,1,0,0,0))/2.
324	if NOT keyword_set(key_periodique) OR nx NE jpi then u(0,,,*)=a
325	a=v(,0,,*)
326	v=(v+shift(v,0,1,0,0))/2.
327	if NOT keyword_set(key_periodique) OR nx NE jpi then v(,0,,*)=a
328	;----------------------------------------------------------------------------
329	; attribution du mask et des tableau de longitude et latitude
330	;----------------------------------------------------------------------------
331	mask = tmask[indice3d]
332	if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
333	;-----------------------------------------------------------
334	if n_elements(valmask) EQ 0 THEN valmask = 1e20
335	landu = where(u GE valmask/10)
336	if landu[0] NE -1 then u[landu] = 0
337	landv = where(v GE valmask/10)
338	if landv[0] NE -1 then v[landv] = 0
339	res=sqrt(u^2+v^2)
340	if NOT keyword_set(key_periodique) OR nx NE jpi then res(0,, , *)=!values.f_nan
341	res(,0, , *)=!values.f_nan
342	mask = where(mask eq 0)
343	IF mask[0] NE -1 THEN BEGIN
344	coeftps = lindgen(jpt)nxny*nzt
345	coeftps = replicate(1, n_elements(mask))#coeftps
346	mask = (temporary(mask))[*]#replicate(1, jpt)
347	mask =temporary(mask[]) + temporary(coeftps[])
348	res(temporary(mask)) = valmask
349	ENDIF
350	; moyennes en tous genres
351	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme
352	if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boite = boite)
353	END
354	;----------------------------------------------------------------------------
355	;----------------------------------------------------------------------------
356	;xy
357	;----------------------------------------------------------------------------
358	;----------------------------------------------------------------------------
359	ELSE:BEGIN ;xy
360	indice2d = lindgen(jpi, jpj)
361	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
362	;------------------------------------------------------------
363	; extraction de u et v sur le domaine qui convient
364	;------------------------------------------------------------
365	case 1 of
366	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
367	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
368	if nxu NE nx then $
369	if indicex[0] EQ premierxu then u = u[0:nx-1, ] ELSE u = u[1: nx, ]
370	IF nxv NE nx THEN $
371	if indicex[0] EQ premierxv then v = v[0:nx-1, ] ELSE v = v[1: nx, ]
372	IF nyu NE ny THEN $
373	if indicey[0] EQ premieryu then u = u[, 0:ny-1] ELSE u = u[, 1: ny]
374	IF nyv NE ny THEN $
375	if indicey[0] EQ premieryv then v = v[, 0:ny-1] ELSE v = v[, 1: ny]
376	END
377	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
378	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
379	u = u[indice2d]
380	v = v[indice2d]
381	END
382	ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
383	endcase
384	;------------------------------------------------------------------
385	; on reform u et v pour s'assurer qu'aucune dimension n'a ete ecrasee
386	;------------------------------------------------------------------
387	if ny EQ 1 then begin
388	u = reform(u, nx, ny, /over)
389	v = reform(v, nx, ny, /over)
390	endif
391	;------------------------------------------------------------------
392	; construction de u et v aux pts T
393	;-----------------------------------------------------------
394	a=u(0,*)
395	u=(u+shift(u,1,0))/2.
396	if NOT keyword_set(key_periodique) OR nx NE jpi then u(0,*)=a
397	a=v(*,0)
398	v=(v+shift(v,0,1))/2.
399	if NOT keyword_set(key_periodique) OR nx NE jpi then v(*,0)=a
400	;----------------------------------------------------------------------------
401	; attribution du mask et des tableau de longitude et latitude
402	; on recupere la grille complette pour etablir un grand mask etendu ds les 4
403	; directions pour couvrir les points pour lesquels un pt terre a ete pris en
404	; compte (faire un petit dessin...)
405	;----------------------------------------------------------------------------
406	mask = tmask[indice2d+jpijpj(niveau-1)]
407	if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over)
408	;-----------------------------------------------------------
409	; construction de terre qui contient tous les point a masquer
410	;-----------------------------------------------------------
411	if n_elements(valmask) EQ 0 THEN valmask = 1e20
412	landu = where(u GE valmask/10)
413	if landu[0] NE -1 then u[landu] = 0
414	landv = where(v GE valmask/10)
415	if landv[0] NE -1 then v[landv] = 0
416	res=sqrt(u^2+v^2)
417	if NOT keyword_set(key_periodique) OR nx NE jpi then res(0,*)=!values.f_nan
418	res(*,0)=!values.f_nan
419	mask = where(mask eq 0)
420	IF mask[0] NE -1 THEN res(mask) = valmask
421	; moyennes en tous genres
422	oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
423	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme
424	if keyword_set(direc) then res = moyenne(res,direc,/nan, boite = boite)
425	END
426	;----------------------------------------------------------------------------
427	endcase
428	;------------------------------------------------------------
429	if keyword_set(key_performance) THEN print, 'temps norme', systime(1)-tempsun
430	return, res
431	end

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

Download in other formats:

Original Format