Context Navigation

div.pro @ 134

Last change on this file since 134 was 134, checked in by navarro, 18 years ago
change *.pro file properties (del eof-style, del executable, set keywords Id
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 14.2 KB

Line
1	;------------------------------------------------------------
2	;------------------------------------------------------------
3	;------------------------------------------------------------
4	;+
5	; NAME:div
6	;
7	; PURPOSE:calcule la divergence d'un champ 2D
8	;
9	; CATEGORY:calcule sur les matrices
10	;
11	; CALLING SEQUENCE:res=div(u,v)
12	;
13	; INPUTS:
14	; u et v deux matrices representant les coordonnes d''un
15	; champ de vecteur
16	;
17	; KEYWORD PARAMETERS:
18	;
19	; OUTPUTS:res: une matrice 2d
20	;
21	; COMMON BLOCKS:
22	; common.pro
23	;
24	; SIDE EFFECTS:
25	;
26	; RESTRICTIONS:
27	; les matrices u et v peuvent de 2 a 4 dimensions.
28	; attention pour distinger les differents configurations de u et v
29	; (xy, xyz, xyt, xyzt), on regarde la variable du common
30	; -time qui contient le calendrier en jour julien d''IDL auquel
31	; se rapportent u et v ansi que la variable
32	; -jpt qui est le nombre de pas de temps a considerer ds time.
33	; les tableaux u et v sont decoupes sur le meme domaine
34	; geographique. A cause du decalage des grilles T, U, V et F il est
35	; possiible que ces 2 tableaux n''aient pas la meme taille et se
36	; repportent a des indices differents. Si tel est le cas les tableaux
37	; sont redecoupes sur les indices qu'ils ont en commun et le dommaine
38	; est redefinit pour qu'il colle a ces indices communs.
39	; pour eviter ces redecoupes utiliser le mot cles /memeindice ds
40	; domdef.pro
41	;
42	; les points sur le bord du dessin sont mis a !values.f_nan
43	;
44	; EXAMPLE:
45	;
46	; MODIFICATION HISTORY:Guillaume Roullet (grlod@ipsl.jussieu.fr)
47	; Creation : printemps 1998
48	; Sebastien Masson (smasson@lodyc.jussieu.fr)
49	; adaptation pour marcher avec un domaine reduit
50	; 12/1/2000
51	;-
52	;------------------------------------------------------------
53	;------------------------------------------------------------
54	;------------------------------------------------------------
55	FUNCTION div, uu, vv
56	;
57	compile_opt idl2, strictarrsubs
58	;
59	tempsun = systime(1) ; pour key_performance
60	@common
61	;
62	IF finite(glamu[0])finite(gphiu[0])finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
63	return, report(['This version of div is based on Arakawa C-grid.' $
64	, 'U and V grids must therefore be defined'])
65	;
66	u = litchamp(uu)
67	v = litchamp(vv)
68	;
69	date1 = time[0]
70	if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
71	if (size(u))[0] NE (size(v))[0] then return, -1
72
73	;------------------------------------------------------------
74	; on trouve les points que u et v ont en communs
75	;------------------------------------------------------------
76	indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
77	indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
78	indicex = inter(indicexu, indicexv)
79	indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
80	indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
81	indicey = inter(indiceyu, indiceyv)
82	nx = n_elements(indicex)
83	ny = n_elements(indicey)
84	indice2d = lindgen(jpi, jpj)
85	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
86	case 1 of
87	;----------------------------------------------------------------------------
88	;----------------------------------------------------------------------------
89	;xyz
90	;----------------------------------------------------------------------------
91	(size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN
92	;------------------------------------------------------------
93	; extraction de u et v sur le domaine qui convient
94	;------------------------------------------------------------
95	case 1 of
96	(size(v))[0] NE 3: return, -1
97	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
98	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
99	case 1 of
100	nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, , ] ELSE u = u[1: nx, , ]
101	nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, , ] ELSE v = v[1: nx, , ]
102	nyu NE ny:if indicey[0] EQ firstyu then u = u[, 0:ny-1, ] ELSE u = u[, 1: ny, ]
103	nyv NE ny:if indicey[0] EQ firstyv then v = v[, 0:ny-1, ] ELSE v = v[, 1: ny, ]
104	ELSE :
105	endcase
106	END
107	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
108	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
109	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
110	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
111	END
112	ELSE:BEGIN
113	zdiv = -1
114	GOTO, sortie
115	end
116
117	endcase
118	;------------------------------------------------------------
119	; calcul de la divergence
120	;------------------------------------------------------------
121	zu = (e2u[indice2d])[*]#replicate(1, nzt)
122	zu = reform(zu, nx, ny, nzt, /over)
123	zu = temporary(u)*temporary(zu) $
124	*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
125	terreu = where(zu EQ 0)
126	if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
127	;
128	zv = (e1v[indice2d])[*]#replicate(1, nzt)
129	zv = reform(zv, nx, ny, nzt, /over)
130	zv = temporary(v)*temporary(zv) $
131	*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
132	terrev = where(zv EQ 0)
133	if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
134	;
135	zdiv = 1e6/(e1t[indice2d]*e2t[indice2d])
136	zdiv = (zdiv)[*]#replicate(1, nzt)
137	zdiv = reform(zdiv, nx, ny, nzt, /over)
138	zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $
139	*tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
140	;------------------------------------------------------------
141	; mise a !values.f_nan de la bordure
142	;------------------------------------------------------------
143	if NOT keyword_set(key_periodic) OR nx NE jpi then begin
144	zdiv[0, , ] = !values.f_nan
145	zdiv[nx-1, , ] = !values.f_nan
146	endif
147	zdiv[, 0, ] = !values.f_nan
148	zdiv[, ny-1, ] = !values.f_nan
149	;
150	zdiv = temporary(zdiv)
151	;
152	if n_elements(valmask) EQ 0 THEN valmask = 1e20
153	terre = where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
154	if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
155	;------------------------------------------------------------
156	; pour le trace graphique
157	;------------------------------------------------------------
158	vargrid = 'T'
159	varname = 'div'
160	varunits = '1e6*s-1'
161	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
162	if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan)
163	END
164	;----------------------------------------------------------------------------
165	;----------------------------------------------------------------------------
166	;xyt
167	;----------------------------------------------------------------------------
168	;----------------------------------------------------------------------------
169	date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN
170	;------------------------------------------------------------
171	; extraction de u et v sur le domaine qui convient
172	;------------------------------------------------------------
173	case 1 of
174	(size(u))[0] NE 3 OR (size(v))[0] NE 3: return, -1
175	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
176	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
177	case 1 of
178	nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, , ] ELSE u = u[1: nx, , ]
179	nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, , ] ELSE v = v[1: nx, , ]
180	nyu NE ny:if indicey[0] EQ firstyu then u = u[, 0:ny-1, ] ELSE u = u[, 1: ny, ]
181	nyv NE ny:if indicey[0] EQ firstyv then v = v[, 0:ny-1, ] ELSE v = v[, 1: ny, ]
182	ELSE :
183	endcase
184	END
185	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
186	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
187	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
188	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
189	END
190	ELSE:return, -1
191	endcase
192	;------------------------------------------------------------
193	; calcul de la divergence
194	;------------------------------------------------------------
195	zu = e2u[indice2d](umask())[indice2d+jpijpj*firstzt]
196	terreu = where(zu EQ 0)
197	if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
198	zu = (zu)[*]#replicate(1, jpt)
199	zu = reform(zu, nx, ny, jpt, /over)
200	zu = temporary(u)*temporary(zu)
201	;
202	zv = e1v[indice2d](vmask())[indice2d+jpijpj*firstzt]
203	terrev = where(zv EQ 0)
204	if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
205	zv = (zv)[*]#replicate(1, jpt)
206	zv = reform(zv, nx, ny, jpt, /over)
207	zv = temporary(v)*temporary(zv)
208	;
209	zdiv = 1e6tmask[indice2d+jpijpjfirstzt]/(e1t[indice2d]e2t[indice2d])
210	zdiv = (zdiv)[*]#replicate(1, jpt)
211	zdiv = reform(zdiv, nx, ny, jpt, /over)
212	terre = where(zdiv EQ 0)
213	zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0))
214	;------------------------------------------------------------
215	; mise a !values.f_nan de la bordure
216	;------------------------------------------------------------
217	if NOT keyword_set(key_periodic) OR nx NE jpi then begin
218	zdiv[0, , ] = !values.f_nan
219	zdiv[nx-1, , ] = !values.f_nan
220	endif
221	zdiv[, 0, ] = !values.f_nan
222	zdiv[, ny-1, ] = !values.f_nan
223	;
224	if n_elements(valmask) EQ 0 THEN valmask = 1e20
225	if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
226	;------------------------------------------------------------
227	; pour le trace graphique
228	;------------------------------------------------------------
229	vargrid = 'T'
230	varname = 'div'
231	varunits = '1e6*s-1'
232	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
233	if keyword_set(direc) then zdiv = grossemoyenne(zdiv,direc,/nan)
234	END
235	;----------------------------------------------------------------------------
236	;----------------------------------------------------------------------------
237	;xyzt
238	;----------------------------------------------------------------------------
239	;----------------------------------------------------------------------------
240	date1 NE date2 AND (size(u))[0] EQ 4:BEGIN
241	return, report('non code!')
242	END
243	;----------------------------------------------------------------------------
244	;----------------------------------------------------------------------------
245	;xy
246	;----------------------------------------------------------------------------
247	;----------------------------------------------------------------------------
248	ELSE:BEGIN ;xy
249	indice3d = lindgen(jpi, jpj, jpk)
250	indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt]
251	;------------------------------------------------------------
252	; extraction de u et v sur le domaine qui convient
253	;------------------------------------------------------------
254	case 1 of
255	(size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN
256	zdiv = -1
257	GOTO, sortie
258	end
259	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
260	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
261	case 1 of
262	nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, ] ELSE u = u[1: nx, ]
263	nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, ] ELSE v = v[1: nx, ]
264	nyu NE ny:if indicey[0] EQ firstyu then u = u[, 0:ny-1] ELSE u = u[, 1: ny]
265	nyv NE ny:if indicey[0] EQ firstyv then v = v[, 0:ny-1] ELSE v = v[, 1: ny]
266	ELSE :
267	endcase
268	END
269	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
270	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
271	u = u[indice2d]
272	v = v[indice2d]
273	END
274	ELSE:return, -1
275	endcase
276	;------------------------------------------------------------
277	; calcul de la divergence
278	;------------------------------------------------------------
279	zu = temporary(u)e2u[indice2d](umask())[indice3d]
280	terreu = where(zu EQ 0)
281	if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
282	zv = temporary(v)e1v[indice2d](vmask())[indice3d]
283	terrev = where(zv EQ 0)
284	if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
285	zdiv = zu-shift(zu, 1, 0)+zv-shift(zv, 0, 1)
286	zdiv = temporary(zdiv)tmask[indice3d]/(e1t[indice2d]e2t[indice2d])
287	;------------------------------------------------------------
288	; mise a !values.f_nan de la bordure
289	;------------------------------------------------------------
290	if NOT keyword_set(key_periodic) OR nx NE jpi then begin
291	zdiv[0, *] = !values.f_nan
292	zdiv[nx-1, *] = !values.f_nan
293	endif
294	zdiv[*, 0] = !values.f_nan
295	zdiv[*, ny-1] = !values.f_nan
296	;
297	zdiv = temporary(zdiv)*1e6
298	;
299	if n_elements(valmask) EQ 0 THEN valmask = 1e20
300	terre = where(tmask[indice3d] EQ 0)
301	if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
302	;------------------------------------------------------------
303	; pour le trace graphique
304	;------------------------------------------------------------
305	vargrid = 'T'
306	varname = 'div'
307	varunits = '1e6*s-1'
308	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
309	if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan)
310
311	END
312	endcase
313
314	;------------------------------------------------------------
315	sortie:
316	if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun
317
318	return, zdiv
319	end

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

Download in other formats:

Original Format