Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

div.pro @ 72

Last change on this file since 72 was 25, checked in by pinsard, 18 years ago
upgrade of CALCULS 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: 14.1 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	tempsun = systime(1) ; pour key_performance
57	@common
58	;
59	IF finite(glamu[0])finite(gphiu[0])finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
60	return, report(['This version of div is based on Arakawa C-grid.' $
61	, 'U and V grids must therefore be defined'])
62	;
63	u = litchamp(uu)
64	v = litchamp(vv)
65	;
66	date1 = time[0]
67	if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
68	if (size(u))[0] NE (size(v))[0] then return, -1
69
70	;------------------------------------------------------------
71	; on trouve les points que u et v ont en communs
72	;------------------------------------------------------------
73	indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
74	indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
75	indicex = inter(indicexu, indicexv)
76	indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
77	indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
78	indicey = inter(indiceyu, indiceyv)
79	nx = n_elements(indicex)
80	ny = n_elements(indicey)
81	indice2d = lindgen(jpi, jpj)
82	indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
83	case 1 of
84	;----------------------------------------------------------------------------
85	;----------------------------------------------------------------------------
86	;xyz
87	;----------------------------------------------------------------------------
88	(size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN
89	;------------------------------------------------------------
90	; extraction de u et v sur le domaine qui convient
91	;------------------------------------------------------------
92	case 1 of
93	(size(v))[0] NE 3: return, -1
94	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
95	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
96	case 1 of
97	nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, , ] ELSE u = u[1: nx, , ]
98	nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, , ] ELSE v = v[1: nx, , ]
99	nyu NE ny:if indicey[0] EQ firstyu then u = u[, 0:ny-1, ] ELSE u = u[, 1: ny, ]
100	nyv NE ny:if indicey[0] EQ firstyv then v = v[, 0:ny-1, ] ELSE v = v[, 1: ny, ]
101	ELSE :
102	endcase
103	END
104	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
105	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
106	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
107	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
108	END
109	ELSE:BEGIN
110	zdiv = -1
111	GOTO, sortie
112	end
113
114	endcase
115	;------------------------------------------------------------
116	; calcul de la divergence
117	;------------------------------------------------------------
118	zu = (e2u[indice2d])[*]#replicate(1, nzt)
119	zu = reform(zu, nx, ny, nzt, /over)
120	zu = temporary(u)*temporary(zu) $
121	*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
122	terreu = where(zu EQ 0)
123	if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
124	;
125	zv = (e1v[indice2d])[*]#replicate(1, nzt)
126	zv = reform(zv, nx, ny, nzt, /over)
127	zv = temporary(v)*temporary(zv) $
128	*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
129	terrev = where(zv EQ 0)
130	if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
131	;
132	zdiv = 1e6/(e1t[indice2d]*e2t[indice2d])
133	zdiv = (zdiv)[*]#replicate(1, nzt)
134	zdiv = reform(zdiv, nx, ny, nzt, /over)
135	zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $
136	*tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
137	;------------------------------------------------------------
138	; mise a !values.f_nan de la bordure
139	;------------------------------------------------------------
140	if NOT keyword_set(key_periodic) OR nx NE jpi then begin
141	zdiv(0, , ) = !values.f_nan
142	zdiv(nx-1, , ) = !values.f_nan
143	endif
144	zdiv(, 0, ) = !values.f_nan
145	zdiv(, ny-1, ) = !values.f_nan
146	;
147	zdiv = temporary(zdiv)
148	;
149	if n_elements(valmask) EQ 0 THEN valmask = 1e20
150	terre = where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
151	if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
152	;------------------------------------------------------------
153	; pour le trace graphique
154	;------------------------------------------------------------
155	vargrid = 'T'
156	varname = 'div'
157	varunits = '1e6*s-1'
158	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
159	if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan)
160	END
161	;----------------------------------------------------------------------------
162	;----------------------------------------------------------------------------
163	;xyt
164	;----------------------------------------------------------------------------
165	;----------------------------------------------------------------------------
166	date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN
167	;------------------------------------------------------------
168	; extraction de u et v sur le domaine qui convient
169	;------------------------------------------------------------
170	case 1 of
171	(size(u))[0] NE 3 OR (size(v))[0] NE 3: return, -1
172	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
173	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
174	case 1 of
175	nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, , ] ELSE u = u[1: nx, , ]
176	nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, , ] ELSE v = v[1: nx, , ]
177	nyu NE ny:if indicey[0] EQ firstyu then u = u[, 0:ny-1, ] ELSE u = u[, 1: ny, ]
178	nyv NE ny:if indicey[0] EQ firstyv then v = v[, 0:ny-1, ] ELSE v = v[, 1: ny, ]
179	ELSE :
180	endcase
181	END
182	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
183	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
184	u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
185	v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
186	END
187	ELSE:return, -1
188	endcase
189	;------------------------------------------------------------
190	; calcul de la divergence
191	;------------------------------------------------------------
192	zu = e2u[indice2d](umask())[indice2d+jpijpj*firstzt]
193	terreu = where(zu EQ 0)
194	if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
195	zu = (zu)[*]#replicate(1, jpt)
196	zu = reform(zu, nx, ny, jpt, /over)
197	zu = temporary(u)*temporary(zu)
198	;
199	zv = e1v[indice2d](vmask())[indice2d+jpijpj*firstzt]
200	terrev = where(zv EQ 0)
201	if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
202	zv = (zv)[*]#replicate(1, jpt)
203	zv = reform(zv, nx, ny, jpt, /over)
204	zv = temporary(v)*temporary(zv)
205	;
206	zdiv = 1e6tmask[indice2d+jpijpjfirstzt]/(e1t[indice2d]e2t[indice2d])
207	zdiv = (zdiv)[*]#replicate(1, jpt)
208	zdiv = reform(zdiv, nx, ny, jpt, /over)
209	terre = where(zdiv EQ 0)
210	zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0))
211	;------------------------------------------------------------
212	; mise a !values.f_nan de la bordure
213	;------------------------------------------------------------
214	if NOT keyword_set(key_periodic) OR nx NE jpi then begin
215	zdiv(0, , ) = !values.f_nan
216	zdiv(nx-1, , ) = !values.f_nan
217	endif
218	zdiv(, 0, ) = !values.f_nan
219	zdiv(, ny-1, ) = !values.f_nan
220	;
221	if n_elements(valmask) EQ 0 THEN valmask = 1e20
222	if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
223	;------------------------------------------------------------
224	; pour le trace graphique
225	;------------------------------------------------------------
226	vargrid = 'T'
227	varname = 'div'
228	varunits = '1e6*s-1'
229	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
230	if keyword_set(direc) then zdiv = grossemoyenne(zdiv,direc,/nan)
231	END
232	;----------------------------------------------------------------------------
233	;----------------------------------------------------------------------------
234	;xyzt
235	;----------------------------------------------------------------------------
236	;----------------------------------------------------------------------------
237	date1 NE date2 AND (size(u))[0] EQ 4:BEGIN
238	return, report('non code!')
239	END
240	;----------------------------------------------------------------------------
241	;----------------------------------------------------------------------------
242	;xy
243	;----------------------------------------------------------------------------
244	;----------------------------------------------------------------------------
245	ELSE:BEGIN ;xy
246	indice3d = lindgen(jpi, jpj, jpk)
247	indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt]
248	;------------------------------------------------------------
249	; extraction de u et v sur le domaine qui convient
250	;------------------------------------------------------------
251	case 1 of
252	(size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN
253	zdiv = -1
254	GOTO, sortie
255	end
256	(size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
257	(size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
258	case 1 of
259	nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, ] ELSE u = u[1: nx, ]
260	nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, ] ELSE v = v[1: nx, ]
261	nyu NE ny:if indicey[0] EQ firstyu then u = u[, 0:ny-1] ELSE u = u[, 1: ny]
262	nyv NE ny:if indicey[0] EQ firstyv then v = v[, 0:ny-1] ELSE v = v[, 1: ny]
263	ELSE :
264	endcase
265	END
266	(size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
267	(size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
268	u = u[indice2d]
269	v = v[indice2d]
270	END
271	ELSE:return, -1
272	endcase
273	;------------------------------------------------------------
274	; calcul de la divergence
275	;------------------------------------------------------------
276	zu = temporary(u)e2u[indice2d](umask())[indice3d]
277	terreu = where(zu EQ 0)
278	if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
279	zv = temporary(v)e1v[indice2d](vmask())[indice3d]
280	terrev = where(zv EQ 0)
281	if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
282	zdiv = zu-shift(zu, 1, 0)+zv-shift(zv, 0, 1)
283	zdiv = temporary(zdiv)tmask[indice3d]/(e1t[indice2d]e2t[indice2d])
284	;------------------------------------------------------------
285	; mise a !values.f_nan de la bordure
286	;------------------------------------------------------------
287	if NOT keyword_set(key_periodic) OR nx NE jpi then begin
288	zdiv(0, *) = !values.f_nan
289	zdiv(nx-1, *) = !values.f_nan
290	endif
291	zdiv(*, 0) = !values.f_nan
292	zdiv(*, ny-1) = !values.f_nan
293	;
294	zdiv = temporary(zdiv)*1e6
295	;
296	if n_elements(valmask) EQ 0 THEN valmask = 1e20
297	terre = where(tmask[indice3d] EQ 0)
298	if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
299	;------------------------------------------------------------
300	; pour le trace graphique
301	;------------------------------------------------------------
302	vargrid = 'T'
303	varname = 'div'
304	varunits = '1e6*s-1'
305	domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
306	if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan)
307
308	END
309	endcase
310
311	;------------------------------------------------------------
312	sortie:
313	if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun
314
315	return, zdiv
316	end

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

Download in other formats:

Original Format