source: trunk/CALCULS/curl.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: 14.5 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME:curl
6;
7; PURPOSE:calcule la composante verticale du rotationnel d''un champ
8; de vecteur horizontaux
9;
10; CATEGORY:calcule sur les matrices
11;
12; CALLING SEQUENCE:res=curl(u,v)
13;
14; INPUTS:
15;       u et v deux matrices representant les coordonnes d''un
16;       champ de vecteur
17;
18; KEYWORD PARAMETERS:
19;
20; OUTPUTS:res: une matrice 2d
21;
22; COMMON BLOCKS:
23;       common.pro
24;
25; SIDE EFFECTS:
26;
27; RESTRICTIONS:
28; les matrices u et v peuvent de 2 a 4 dimensions.
29; attention pour distinger les differents configurations de u et v
30; (xy, xyz, xyt, xyzt), on regarde la variable du common
31;        -time qui contient le calendrier en jour julien d''IDL auquel
32;        se rapportent u et v ansi que la variable
33;        -jpt qui est le nombre de pas de temps a considerer ds time.
34; les tableaux u et v sont decoupes sur le meme domaine
35; geographique. A cause du decalage des grilles T, U, V et F il est
36; possiible que ces 2 tableaux n''aient pas la meme taille et se
37; repportent a des indices differents. Si tel est le cas les tableaux
38; sont redecoupes sur les indices qu'ils ont en commun et le dommaine
39; est redefinit pour qu'il colle a ces indices communs.
40; pour eviter ces redecoupes utiliser le mot cles /memeindice ds
41; domdef.pro
42;
43;
44; les points sur le bord du dessin sont mis a !values.f_nan
45;
46; EXAMPLE:
47;
48; MODIFICATION HISTORY:Guillaume Roullet (grlod@ipsl.jussieu.fr)
49;
50;                      Sebastien Masson (smasson@lodyc.jussieu.fr)
51;                      adaptation pour marcher avec un domaine reduit
52;
53;                      21/5/1999: valeurs manquantes a !values.f_nan
54;periodicite
55;-
56;------------------------------------------------------------
57;------------------------------------------------------------
58;------------------------------------------------------------
59FUNCTION curl, uu, vv
60@common
61   tempsun = systime(1)         ; pour key_performance
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))[premierxu:premierxu+nxu-1]
74   indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1]
75   indicex = inter(indicexu, indicexv)
76   indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1]
77   indiceyv = (lindgen(jpj))[premieryv:premieryv+nyv-1]
78   indicey = inter(indiceyu, indiceyv)
79   nx = n_elements(indicex)
80   ny = n_elements(indicey)
81   case 1 of
82;----------------------------------------------------------------------------
83;----------------------------------------------------------------------------
84;xyz
85;----------------------------------------------------------------------------
86      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN
87         indice2d = lindgen(jpi, jpj)
88         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
89;------------------------------------------------------------
90; extraction de u et v sur le domaine qui convient
91;------------------------------------------------------------
92         case 1 of
93            (size(u))[0] NE 3 OR (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 premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
98               nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
99               nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
100               nyv NE ny:if indicey[0] EQ premieryv 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:return,  -1
110         endcase
111;------------------------------------------------------------
112; calcul du rotationnel
113;------------------------------------------------------------
114         coefu = (e1u[indice2d])[*]#replicate(1, nzt)
115         coefu = reform(coefu, nx, ny, nzt, /over)
116         coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]
117         terreu = where(coefu EQ 0)
118         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
119         
120         coefv = (e2v[indice2d])[*]#replicate(1, nzt)
121         coefv = reform(coefv, nx, ny, nzt, /over)
122         coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]
123         terrev = where(coefv EQ 0)
124         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
125
126         tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]
127         div = (e1f[indice2d]*e2f[indice2d])[*]#replicate(1, nzt)
128         div = reform(div, nx, ny, nzt, /over)
129         tabf = tabf/div
130;
131         zu = u*temporary(coefu)
132         zv = v*temporary(coefv)
133
134         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
135         psi = tabf*psi
136;------------------------------------------------------------
137; mise a !values.f_nan de la bordure
138;------------------------------------------------------------
139         if NOT keyword_set(key_periodique)  OR nx NE jpi then begin
140            psi(0, *, *) = !values.f_nan
141            psi(nx-1, *, *) = !values.f_nan
142         endif
143         psi(*, 0, *) = !values.f_nan
144         psi(*, ny-1, *) = !values.f_nan
145;
146         if n_elements(valmask) EQ 0 THEN valmask = 1e20
147         terref =  where(tabf EQ 0)
148         if terref[0] NE -1 then psi[temporary(terref)] = valmask
149;------------------------------------------------------------
150; pour le trace graphique
151;------------------------------------------------------------
152         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f']
153         if keyword_set(direc) then psi = moyenne(psi,direc,/nan, boite = boite)
154
155;----------------------------------------------------------------------------
156      END
157;----------------------------------------------------------------------------
158;----------------------------------------------------------------------------
159;xyt
160;----------------------------------------------------------------------------
161;----------------------------------------------------------------------------
162      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN
163         indice2d = lindgen(jpi, jpj)
164         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
165;------------------------------------------------------------
166; extraction de u et v sur le domaine qui convient
167;------------------------------------------------------------
168         case 1 of
169            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
170             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
171               if nxu NE nx then $
172                if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
173               IF nxv NE nx THEN $
174                if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
175               IF nyu NE ny THEN $
176                if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
177               IF nyv NE ny THEN $
178                if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
179            END
180            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
181             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
182               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
183               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
184            END
185            ELSE:BEGIN
186               print, 'problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs'
187               return, -1
188            end
189         endcase
190;----------------------------------------------------------------------------
191; calcul du rotationnel
192;----------------------------------------------------------------------------
193         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*premierzt]
194         terreu = where(coefu EQ 0)
195         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
196         coefu = temporary(coefu[*])#replicate(1, jpt)
197         coefu = reform(coefu, nx, ny, jpt, /over)
198;
199         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*premierzt]
200         terrev = where(coefv EQ 0)
201         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
202         coefv = temporary(coefv[*])#replicate(1, jpt)
203         coefv = reform(coefv, nx, ny, jpt, /over)
204;
205         tabf = (fmask())[indice2d+jpi*jpj*premierzt]/(e1f[indice2d]*e2f[indice2d])
206         tabf = temporary(tabf[*])#replicate(1, jpt)
207         tabf = reform(tabf, nx, ny, jpt, /over)
208;------------------------------------------------------------
209; calcul du rotationnel
210;------------------------------------------------------------
211         zu = u*temporary(coefu)
212         zv = v*temporary(coefv)
213;
214         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
215         psi = tabf*psi
216;------------------------------------------------------------
217; mise a !values.f_nan de la bordure
218;------------------------------------------------------------
219         if NOT keyword_set(key_periodique) OR nx NE jpi then begin
220            psi(0, *, *) = !values.f_nan
221            psi(nx-1, *, *) = !values.f_nan
222         endif
223         psi(*, 0, *) = !values.f_nan
224         psi(*, ny-1, *) = !values.f_nan
225         if n_elements(valmask) EQ 0 THEN valmask = 1e20
226         terref =  where(tabf EQ 0)
227         if terref[0] NE -1 then psi[temporary(terref)] = valmask
228;----------------------------------------------------------------------------
229         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f']
230         if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan, boite = boite)
231;----------------------------------------------------------------------------
232      END
233;----------------------------------------------------------------------------
234;----------------------------------------------------------------------------
235;xyzt
236;----------------------------------------------------------------------------
237;----------------------------------------------------------------------------
238      date1 NE date2 AND (size(u))[0] EQ 4:BEGIN
239         return, report('non code!')
240      END
241;----------------------------------------------------------------------------
242;----------------------------------------------------------------------------
243;xy
244;----------------------------------------------------------------------------
245;----------------------------------------------------------------------------
246      ELSE:BEGIN                ;xy
247         indice2d = lindgen(jpi, jpj)
248         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
249;------------------------------------------------------------
250; extraction de u et v sur le domaine qui convient
251;------------------------------------------------------------
252         case 1 of
253            (size(u))[0] NE 2 OR (size(v))[0] NE 2: return,  -1
254            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
255             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
256               if nxu NE nx then $
257                if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
258               IF nxv NE nx THEN $
259                if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
260               IF nyu NE ny THEN $
261                if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
262               IF nyv NE ny THEN $
263                if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
264            END
265            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
266             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
267               u = u[indice2d]
268               v = v[indice2d]
269            END
270            ELSE:return,  -1
271         endcase
272;------------------------------------------------------------
273; calcul du rotationnel
274;------------------------------------------------------------
275         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*premierzt]
276         terreu = where(coefu EQ 0)
277         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
278         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*premierzt]
279         terrev = where(coefv EQ 0)
280         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
281         tabf = (fmask())[indice2d+jpi*jpj*premierzt]/(e1f[indice2d]*e2f[indice2d])
282;
283         zu = u*temporary(coefu)
284         zv = v*temporary(coefv)
285
286         psi =  (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1))
287         psi = tabf*psi
288
289;------------------------------------------------------------
290; mise a !values.f_nan de la bordure
291;------------------------------------------------------------
292         if  NOT keyword_set(key_periodique) OR nx NE jpi then begin
293            psi(0, *) = !values.f_nan
294            psi(nx-1, *) = !values.f_nan
295         endif
296         psi(*, 0) = !values.f_nan
297         psi(*, ny-1) = !values.f_nan
298;
299         if n_elements(valmask) EQ 0 THEN valmask = 1e20
300         terref =  where(tabf EQ 0)
301         if terref[0] NE -1 then psi[temporary(terref)] = valmask
302;------------------------------------------------------------
303; pour le trace graphique
304;------------------------------------------------------------
305         domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f']
306         if keyword_set(direc) then psi = moyenne(psi,direc,/nan, boite = boite)
307
308
309      END
310;----------------------------------------------------------------------------
311   endcase
312;------------------------------------------------------------
313   if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun
314;------------------------------------------------------------
315   vargrid = 'F'
316   varname = 'vorticity'
317   return, psi
318end
Note: See TracBrowser for help on using the repository browser.