source: trunk/divfred.pro @ 12

Last change on this file since 12 was 12, checked in by pinsard, 17 years ago

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