source: trunk/ToBeReviewed/CALCULS/curl.pro @ 25

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