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