Context Navigation

norme.pro @ 232

Last change on this file since 232 was 232, checked in by pinsard, 17 years ago
improvements/corrections of some *.pro headers
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 21.6 KB

Line
1	;+
2	;
3	; @file_comments
4	; calculate the norm of a field of vectors, then make a possible average.
5	; Comment 1: The field of vector can be, 2d:xy, 3d: xyz or xyt,
6	; 4d: xyzt
7	; Comment 2:
8	; The calculation of the norm is made before the possible spatial or
9	; temporal average because the average of the norm is not equal to the
10	; norm of averages
11	;
12	; @categories
13	; Calculation
14	;
15	; @param COMPOSANTEU {in}{required}
16	; an 2d, 3d or 4d array
17	;
18	; @param COMPOSANTEV {in}{required}
19	; an 2d, 3d or 4d array
20	;
21	; @keyword BOXZOOM
22	; boxzoom on which do the average (by default the domain selected
23	; by the last domdef done)
24	;
25	; @keyword DIREC
26	; 't' 'x' 'y' 'z' 'xys' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt'
27	; 'xzt' 'yzt' 'xyzt' Direction on which do averages
28	;
29	; @returns
30	; Array to trace with plt, pltz or pltt.
31	;
32	; @uses
33	; common.pro
34	;
35	; @restrictions
36	; The norm is calculated on points TTo do this calculation, we average
37	; field U and Von points T before calculate the norme. At the edge of
38	; coast and of domain, we can not calculate fields U and V at points T,
39	; that is why these points are at value !values.f_nan.
40	;
41	; When we calculate on a reduce geographic domain, field U and V have not
42	; necessarily the same number of point. In this case, we recut U and V to
43	; keep only common points. We profit of this to redo a domdef which redefine
44	; a geographic domain on which fields U and V are extracted on same points
45	;
46	; To know what type of array we work with, we test its size and dates
47	; gave by time[0] and time[jpt-1] to know if thee is a temporal dimension.
48	; Before to start norme, make sure that time and jpt are defined how
49	; they have to!
50	;
51	; @examples
52	; To calculate the average of the norme of streams on all the domain
53	; between 0 et 50:
54	; IDL> res=norme(un,vn,boxzoom=[0,50],dir='xyz')
55	;
56	; @history
57	; Sebastien Masson (smasson\@lodyc.jussieu.fr)
58	; 9/6/1999
59	;
60	; @version
61	; $Id$
62	;
63	;-
64	;
65	FUNCTION norme, composanteu, composantev, BOXZOOM = boxzoom, DIREC = direc, _EXTRA = ex
66	;
67	compile_opt idl2, strictarrsubs
68	;
69	@cm_4mesh
70	@cm_4data
71	@cm_4cal
72	IF NOT keyword_set(key_forgetold) THEN BEGIN
73	@updatenew
74	@updatekwd
75	ENDIF
76	;---------------------------------------------------------
77	tempsun = systime(1) ; To key_performance
78	;
79	IF finite(glamu[0])finite(gphiu[0])finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
80	return, report(['This version of norme is based on Arakawa C-grid.' $
81	, 'U and V grids must therefore be defined'])
82	;
83	;------------------------------------------------------------
84	if keyword_set(boxzoom) then BEGIN
85	Case 1 Of
86	N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
87	N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
88	N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2]
89	N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]]
90	N_Elements(Boxzoom) Eq 6:bte = Boxzoom
91	Else: return, report('Mauvaise Definition de Boxzoom')
92	ENDCASE
93	domdef, boxzoom
94	ENDIF
95	;------------------------------------------------------------
96	if NOT keyword_set(direc) then direc = 0
97	; construction of u and v at points T
98	u = litchamp(composanteu)
99	v = litchamp(composantev)
100	date1 = time[0]
101	if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
102
103	if (size(u))[0] NE (size(v))[0] then return, -1
104
105	vargrid='T'
106	varname = 'norme '
107	valmask = 1e20
108	;
109	grilleu = litchamp(composanteu, /grid)
110	if grilleu EQ '' then grilleu = 'U'
111	grillev = litchamp(composantev, /grid)
112	if grillev EQ '' then grillev = 'V'
113	IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1
114	IF grilleu EQ 'T' AND grillev EQ 'T' THEN BEGIN
115	interpolle = 0
116	return, report('cas non code mais facile a faire!')
117	ENDIF ELSE interpolle = 1
118	if keyword_set(inverse) then begin
119	rien = u
120	u = v
121	v = rien
122	endif
123
124
125	;------------------------------------------------------------
126	; We find common points between u and v
127	;------------------------------------------------------------
128	indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
129	indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
130	indicex = inter(indicexu, indicexv)
131	indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
132	indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
133	indicey = inter(indiceyu, indiceyv)
134	nx = n_elements(indicex)
135	ny = n_elements(indicey)
136	;----------------------------------------------------------------------------
137	case 1 of
138	;----------------------------------------------------------------------------
139	;----------------------------------------------------------------------------
140	;xyz
141	;----------------------------------------------------------------------------
142	;----------------------------------------------------------------------------
143	(size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN
144	;----------------------------------------------------------------------------
145	indice2d = lindgen(jpi, jpj)
146	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
147	indice3d = lindgen(jpi, jpj, jpk)
148	indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt]
149	;------------------------------------------------------------
150	; extraction of u and v on the appropriated domain
151	;------------------------------------------------------------
152	case 1 of
153	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
154	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
155	case (size(u))[3] OF
156	nzt:BEGIN
157	if nxu NE nx then $
158	if indicex[0] EQ firstxu then u = u[0:nx-1,,] ELSE u = u[1: nx, ,]
159	IF nxv NE nx THEN $
160	if indicex[0] EQ firstxv then v = v[0:nx-1,,] ELSE v = v[1: nx, ,]
161	IF nyu NE ny THEN $
162	if indicey[0] EQ firstyu then u = u[,0:ny-1,] ELSE u = u[, 1: ny,]
163	IF nyv NE ny THEN $
164	if indicey[0] EQ firstyv then v = v[,0:ny-1,] ELSE v = v[, 1: ny,]
165	end
166	jpk:BEGIN
167	if nxu NE nx then $
168	if indicex[0] EQ firstxu then u = u[0:nx-1, ,firstzt:lastzt] ELSE u = u[1: nx, ,firstzt:lastzt]
169	IF nxv NE nx THEN $
170	if indicex[0] EQ firstxv then v = v[0:nx-1, ,firstzt:lastzt] ELSE v = v[1: nx, ,firstzt:lastzt]
171	IF nyu NE ny THEN $
172	if indicey[0] EQ firstyu then u = u[, 0:ny-1,firstzt:lastzt] ELSE u = u[, 1: ny,firstzt:lastzt]
173	IF nyv NE ny THEN $
174	if indicey[0] EQ firstyv then v = v[, 0:ny-1,firstzt:lastzt] ELSE v = v[, 1: ny,firstzt:lastzt]
175	end
176	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
177	endcase
178	END
179	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
180	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN
181	u = u[indice3d]
182	v = v[indice3d]
183	END
184	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
185	endcase
186	;------------------------------------------------------------------
187	; We reshape u and v to make sure that no dimension has been erased
188	;------------------------------------------------------------------
189	if nzt EQ 1 then begin
190	u = reform(u, nx, ny, nzt, /over)
191	v = reform(v, nx, ny, nzt, /over)
192	endif
193	;------------------------------------------------------------------
194	; construction of u and v at points T
195	;-----------------------------------------------------------
196	a=u[0,,]
197	u=(u+shift(u,1,0,0))/2.
198	if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,,]=a
199	a=v[,0,]
200	v=(v+shift(v,0,1,0))/2.
201	if NOT keyword_set(key_periodic) OR nx NE jpi then v[,0,]=a
202	;----------------------------------------------------------------------------
203	; attribution of the mask and of longitude and latitude arrays
204	;----------------------------------------------------------------------------
205	mask = tmask[indice3d]
206	if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
207	;-----------------------------------------------------------
208	if n_elements(valmask) EQ 0 THEN valmask = 1e20
209	landu = where(u GE valmask/10)
210	if landu[0] NE -1 then u[landu] = 0
211	landv = where(v GE valmask/10)
212	if landv[0] NE -1 then v[landv] = 0
213	res=sqrt(u^2+v^2)
214	if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,, ]=!values.f_nan
215	res[,0, ]=!values.f_nan
216	mask = where(mask eq 0)
217	IF mask[0] NE -1 THEN res[mask] = valmask
218	; All kind of average
219	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
220	if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
221	;-----------------------------------------------------------
222	END
223	;----------------------------------------------------------------------------
224	;----------------------------------------------------------------------------
225	;xyt
226	;----------------------------------------------------------------------------
227	;----------------------------------------------------------------------------
228	date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN
229	indice2d = lindgen(jpi, jpj)
230	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
231	;------------------------------------------------------------
232	; extraction of u and v on the appropriated domain
233	;------------------------------------------------------------
234	case 1 of
235	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
236	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
237	if nxu NE nx then $
238	if indicex[0] EQ firstxu then u = u[0:nx-1, , ] ELSE u = u[1: nx, , ]
239	IF nxv NE nx THEN $
240	if indicex[0] EQ firstxv then v = v[0:nx-1, , ] ELSE v = v[1: nx, , ]
241	IF nyu NE ny THEN $
242	if indicey[0] EQ firstyu then u = u[, 0:ny-1, ] ELSE u = u[, 1: ny, ]
243	IF nyv NE ny THEN $
244	if indicey[0] EQ firstyv then v = v[, 0:ny-1, ] ELSE v = v[, 1: ny, ]
245	END
246	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
247	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
248	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
249	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
250	END
251	ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
252	endcase
253	;------------------------------------------------------------------
254	; construction of u and v at points T
255	;-----------------------------------------------------------
256	a=u[0,,]
257	u=(u+shift(u,1,0,0))/2.
258	if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,,]=a
259	a=v[,0,]
260	v=(v+shift(v,0,1,0))/2.
261	if NOT keyword_set(key_periodic) OR nx NE jpi then v[,0,]=a
262	;----------------------------------------------------------------------------
263	; attribution of the mask and of longitude and latitude arrays.
264	; We recover the complete grid to establish a big mask extent in the four
265	; direction to cover pointsfor which a land point has been
266	; considerated (make a small drawing)
267	;----------------------------------------------------------------------------
268	mask = tmask[indice2d+jpijpjfirstzt]
269	if ny EQ 1 then mask = reform(mask, nx, ny, /over)
270	;-----------------------------------------------------------
271	; construction of land containing all points to mask
272	;-----------------------------------------------------------
273	if n_elements(valmask) EQ 0 THEN valmask = 1e20
274	landu = where(u GE valmask/10)
275	if landu[0] NE -1 then u[landu] = 0
276	landv = where(v GE valmask/10)
277	if landv[0] NE -1 then v[landv] = 0
278	res=sqrt(u^2+v^2)
279	if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,, ]=!values.f_nan
280	res[,0, ]=!values.f_nan
281	mask = where(mask eq 0)
282	IF mask[0] NE -1 THEN BEGIN
283	coeftps = lindgen(jpt)nxny
284	coeftps = replicate(1, n_elements(mask))#coeftps
285	mask = (temporary(mask))[*]#replicate(1, jpt)
286	mask =temporary(mask[]) + temporary(coeftps[])
287	res[temporary(mask)] = valmask
288	ENDIF
289	; moyennes en tous genres
290	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
291	if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
292	END
293	;----------------------------------------------------------------------------
294	;----------------------------------------------------------------------------
295	;xyzt
296	;----------------------------------------------------------------------------
297	;----------------------------------------------------------------------------
298	date1 NE date2 AND (size(u))[0] EQ 4:BEGIN
299	indice2d = lindgen(jpi, jpj)
300	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
301	indice3d = lindgen(jpi, jpj, jpk)
302	indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt]
303	;------------------------------------------------------------
304	; extraction of u and v on the appropriated domain
305	;------------------------------------------------------------
306	case 1 of
307	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
308	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
309	case (size(u))[3] OF
310	nzt:BEGIN
311	if nxu NE nx then $
312	if indicex[0] EQ firstxu then u = u[0:nx-1,,,] ELSE u = u[1: nx, ,,]
313	IF nxv NE nx THEN $
314	if indicex[0] EQ firstxv then v = v[0:nx-1,,,] ELSE v = v[1: nx, ,,]
315	IF nyu NE ny THEN $
316	if indicey[0] EQ firstyu then u = u[,0:ny-1,,] ELSE u = u[, 1: ny,,]
317	IF nyv NE ny THEN $
318	if indicey[0] EQ firstyv then v = v[,0:ny-1,,] ELSE v = v[, 1: ny,,]
319	end
320	jpk:BEGIN
321	if nxu NE nx then $
322	if indicex[0] EQ firstxu then u = u[0:nx-1, ,firstzt:lastzt,] ELSE u = u[1: nx, ,firstzt:lastzt,]
323	IF nxv NE nx THEN $
324	if indicex[0] EQ firstxv then v = v[0:nx-1, ,firstzt:lastzt,] ELSE v = v[1: nx, ,firstzt:lastzt,]
325	IF nyu NE ny THEN $
326	if indicey[0] EQ firstyu then u = u[, 0:ny-1,firstzt:lastzt,] ELSE u = u[, 1: ny,firstzt:lastzt,]
327	IF nyv NE ny THEN $
328	if indicey[0] EQ firstyv then v = v[, 0:ny-1,firstzt:lastzt,] ELSE v = v[, 1: ny,firstzt:lastzt,]
329	end
330	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
331	endcase
332	END
333	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
334	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN
335	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *]
336	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *]
337	END
338	ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
339	endcase
340	;------------------------------------------------------------------
341	; construction of u and v at points T
342	;-----------------------------------------------------------
343	a=u[0,,,*]
344	u=(u+shift(u,1,0,0,0))/2.
345	if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,,,*]=a
346	a=v[,0,,*]
347	v=(v+shift(v,0,1,0,0))/2.
348	if NOT keyword_set(key_periodic) OR nx NE jpi then v[,0,,*]=a
349	;----------------------------------------------------------------------------
350	; attribution of the mask and of longitude and latitude arrays
351	;----------------------------------------------------------------------------
352	mask = tmask[indice3d]
353	if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
354	;-----------------------------------------------------------
355	if n_elements(valmask) EQ 0 THEN valmask = 1e20
356	landu = where(u GE valmask/10)
357	if landu[0] NE -1 then u[landu] = 0
358	landv = where(v GE valmask/10)
359	if landv[0] NE -1 then v[landv] = 0
360	res=sqrt(u^2+v^2)
361	if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,, , *]=!values.f_nan
362	res[,0, , *]=!values.f_nan
363	mask = where(mask eq 0)
364	IF mask[0] NE -1 THEN BEGIN
365	coeftps = lindgen(jpt)nxny*nzt
366	coeftps = replicate(1, n_elements(mask))#coeftps
367	mask = (temporary(mask))[*]#replicate(1, jpt)
368	mask =temporary(mask[]) + temporary(coeftps[])
369	res[temporary(mask)] = valmask
370	ENDIF
371	; All kind of average
372	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
373	if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
374	END
375	;----------------------------------------------------------------------------
376	;----------------------------------------------------------------------------
377	;xy
378	;----------------------------------------------------------------------------
379	;----------------------------------------------------------------------------
380	ELSE:BEGIN ;xy
381	indice2d = lindgen(jpi, jpj)
382	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
383	;------------------------------------------------------------
384	; extraction of u and v on the appropriated domain
385	;------------------------------------------------------------
386	case 1 of
387	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
388	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
389	if nxu NE nx then $
390	if indicex[0] EQ firstxu then u = u[0:nx-1, ] ELSE u = u[1: nx, ]
391	IF nxv NE nx THEN $
392	if indicex[0] EQ firstxv then v = v[0:nx-1, ] ELSE v = v[1: nx, ]
393	IF nyu NE ny THEN $
394	if indicey[0] EQ firstyu then u = u[, 0:ny-1] ELSE u = u[, 1: ny]
395	IF nyv NE ny THEN $
396	if indicey[0] EQ firstyv then v = v[, 0:ny-1] ELSE v = v[, 1: ny]
397	END
398	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
399	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
400	u = u[indice2d]
401	v = v[indice2d]
402	END
403	ELSE:return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs')
404	endcase
405	;------------------------------------------------------------------
406	; We reshape u and v to make sure that no dimension has been erased
407	;------------------------------------------------------------------
408	if ny EQ 1 then begin
409	u = reform(u, nx, ny, /over)
410	v = reform(v, nx, ny, /over)
411	endif
412	;------------------------------------------------------------------
413	; construction of u and v at points T
414	;-----------------------------------------------------------
415	a=u[0,*]
416	u=(u+shift(u,1,0))/2.
417	if NOT keyword_set(key_periodic) OR nx NE jpi then u[0,*]=a
418	a=v[*,0]
419	v=(v+shift(v,0,1))/2.
420	if NOT keyword_set(key_periodic) OR nx NE jpi then v[*,0]=a
421	;----------------------------------------------------------------------------
422	; attribution of the mask and of longitude and latitude arrays.
423	; We recover the complete grid to establish a big mask extent in the four
424	; direction to cover pointsfor which a land point has been
425	; considerated (make a small drawing)
426	;----------------------------------------------------------------------------
427	mask = tmask[indice2d+jpijpjfirstzt]
428	if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over)
429	;-----------------------------------------------------------
430	; construction of land containing all points to mask
431	;-----------------------------------------------------------
432	if n_elements(valmask) EQ 0 THEN valmask = 1e20
433	landu = where(u GE valmask/10)
434	if landu[0] NE -1 then u[landu] = 0
435	landv = where(v GE valmask/10)
436	if landv[0] NE -1 then v[landv] = 0
437	res=sqrt(u^2+v^2)
438	if NOT keyword_set(key_periodic) OR nx NE jpi then res[0,*]=!values.f_nan
439	res[*,0]=!values.f_nan
440	mask = where(mask eq 0)
441	IF mask[0] NE -1 THEN res[mask] = valmask
442	; All kind of average
443	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
444	if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef)
445	END
446	;----------------------------------------------------------------------------
447	endcase
448	;------------------------------------------------------------
449	if keyword_set(key_performance) THEN print, 'temps norme', systime(1)-tempsun
450	return, res
451	end

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

Download in other formats:

Original Format