source: trunk/SRC/ToBeReviewed/CALCULS/curl.pro @ 163

Last change on this file since 163 was 163, checked in by navarro, 18 years ago

header improvements : type of parameters and keywords, default values, spell checking + idldoc assistant (IDL online_help)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.3 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5;
6; @file_comments
7; Calculate the vertical component of the curl of a field of horizontal vectors
8;
9; @categories
10; Calculation
11;
12; @param UU
13; Matrix representing coordinates of a field of vectors
14;
15; @param VV
16; Matrix representing coordinates of a field of vectors
17;
18; @returns RES
19; A 2d matrix
20;
21; @uses
22; common.pro
23;
24; @restrictions
25; U and V matrices can be 2 or 4d.
26; Beware, to discern different configuration of U and V (xy, xyz, xyt, xyzt),
27; we look at the variable of the common
28;        -time which contain the calendar in IDL julian days to which U and
29; V refered to, in the same way as the variable
30;        -jpt which is the number of time's step to consider in time.
31; U and V arrays ae cut in the same geographic domain. Because of the gap of
32; T, U, V and F grids, it is possible that these two arrays has not the same
33; size and refered to different indexes. In this case, arrays are re-cut on
34; common indexes and the domain is redefined to match with these common
35; indexes. To avoid these re-cuts, use the keyword /memeindice in domdef.pro
36;
37;
38; Points on the drawing edge are at !values.f_nan
39;
40; @history
41; Guillaume Roullet (grlod\@ipsl.jussieu.fr)
42;
43; Sebastien Masson (smasson\@lodyc.jussieu.fr)
44; adaptation to work with a reduce domain
45; 21/5/1999: missing values at !values.f_nan
46;
47; @version
48; $Id$
49;
50;-
51;------------------------------------------------------------
52;------------------------------------------------------------
53;------------------------------------------------------------
54FUNCTION curl, uu, vv
55;
56  compile_opt idl2, strictarrsubs
57;
58@common
59   tempsun = systime(1)         ; To key_performance
60;
61   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
62     return, report(['This version of curl is based on Arakawa C-grid.' $
63                     , 'U and V grids must therefore be defined'])
64;
65   u = litchamp(uu)
66   v = litchamp(vv)
67
68   date1 = time[0]
69   if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
70   if (size(u))[0] NE (size(v))[0] then return,  -1
71
72;------------------------------------------------------------
73; We find common points between U and V
74;------------------------------------------------------------
75   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
76   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
77   indicex = inter(indicexu, indicexv)
78   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
79   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
80   indicey = inter(indiceyu, indiceyv)
81   nx = n_elements(indicex)
82   ny = n_elements(indicey)
83   case 1 of
84;----------------------------------------------------------------------------
85;----------------------------------------------------------------------------
86;xyz
87;----------------------------------------------------------------------------
88      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN
89         indice2d = lindgen(jpi, jpj)
90         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
91;------------------------------------------------------------
92; extraction of U and V on the appropriated domain
93;------------------------------------------------------------
94         case 1 of
95            (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1
96            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
97             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
98            case 1 of
99               nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
100               nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
101               nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
102               nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
103               ELSE :
104            endcase
105            END
106            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
107             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
108               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
109               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
110            END
111            ELSE:return,  -1
112         endcase
113;------------------------------------------------------------
114; calculation of the curl
115;------------------------------------------------------------
116         coefu = (e1u[indice2d])[*]#replicate(1, nzt)
117         coefu = reform(coefu, nx, ny, nzt, /over)
118         coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
119         terreu = where(coefu EQ 0)
120         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
121         
122         coefv = (e2v[indice2d])[*]#replicate(1, nzt)
123         coefv = reform(coefv, nx, ny, nzt, /over)
124         coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
125         terrev = where(coefv EQ 0)
126         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
127
128         tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
129         div = (e1f[indice2d]*e2f[indice2d])[*]#replicate(1, nzt)
130         div = reform(div, nx, ny, nzt, /over)
131         tabf = tabf/div
132;
133         zu = u*temporary(coefu)
134         zv = v*temporary(coefv)
135
136         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
137         psi = tabf*psi
138;------------------------------------------------------------
139; Edging put at !values.f_nan
140;------------------------------------------------------------
141         if NOT keyword_set(key_periodic)  OR nx NE jpi then begin
142            psi[0, *, *] = !values.f_nan
143            psi[nx-1, *, *] = !values.f_nan
144         endif
145         psi[*, 0, *] = !values.f_nan
146         psi[*, ny-1, *] = !values.f_nan
147;
148         if n_elements(valmask) EQ 0 THEN valmask = 1e20
149         terref =  where(tabf EQ 0)
150         if terref[0] NE -1 then psi[temporary(terref)] = valmask
151;------------------------------------------------------------
152; For the graphic drawing
153;------------------------------------------------------------
154         domdef, (glagmt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
155         if keyword_set(direc) then psi = moyenne(psi,direc,/nan)
156
157;----------------------------------------------------------------------------
158      END
159;----------------------------------------------------------------------------
160;----------------------------------------------------------------------------
161;xyt
162;----------------------------------------------------------------------------
163;----------------------------------------------------------------------------
164      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN
165         indice2d = lindgen(jpi, jpj)
166         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
167;------------------------------------------------------------
168; extraction of U and V on the appropriated domain
169;------------------------------------------------------------
170         case 1 of
171            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
172             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
173               if nxu NE nx then $
174                if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
175               IF nxv NE nx THEN $
176                if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
177               IF nyu NE ny THEN $
178                if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
179               IF nyv NE ny THEN $
180                if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
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:BEGIN
188               print, 'problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs'
189               return, -1
190            end
191         endcase
192;----------------------------------------------------------------------------
193; Calculation of the curl
194;----------------------------------------------------------------------------
195         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
196         terreu = where(coefu EQ 0)
197         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
198         coefu = temporary(coefu[*])#replicate(1, jpt)
199         coefu = reform(coefu, nx, ny, jpt, /over)
200;
201         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
202         terrev = where(coefv EQ 0)
203         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
204         coefv = temporary(coefv[*])#replicate(1, jpt)
205         coefv = reform(coefv, nx, ny, jpt, /over)
206;
207         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
208         tabf = temporary(tabf[*])#replicate(1, jpt)
209         tabf = reform(tabf, nx, ny, jpt, /over)
210;------------------------------------------------------------
211; Calculation of the curl
212;------------------------------------------------------------
213         zu = u*temporary(coefu)
214         zv = v*temporary(coefv)
215;
216         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
217         psi = tabf*psi
218;------------------------------------------------------------
219; extraction of U and V on the appropriated domain
220;------------------------------------------------------------
221         if NOT keyword_set(key_periodic) OR nx NE jpi then begin
222            psi[0, *, *] = !values.f_nan
223            psi[nx-1, *, *] = !values.f_nan
224         endif
225         psi[*, 0, *] = !values.f_nan
226         psi[*, ny-1, *] = !values.f_nan
227         if n_elements(valmask) EQ 0 THEN valmask = 1e20
228         terref =  where(tabf EQ 0)
229         if terref[0] NE -1 then psi[temporary(terref)] = valmask
230;----------------------------------------------------------------------------
231         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
232         if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan)
233;----------------------------------------------------------------------------
234      END
235;----------------------------------------------------------------------------
236;----------------------------------------------------------------------------
237;xyzt
238;----------------------------------------------------------------------------
239;----------------------------------------------------------------------------
240      date1 NE date2 AND (size(u))[0] EQ 4:BEGIN
241         return, report('non code!')
242      END
243;----------------------------------------------------------------------------
244;----------------------------------------------------------------------------
245;xy
246;----------------------------------------------------------------------------
247;----------------------------------------------------------------------------
248      ELSE:BEGIN                ;xy
249         indice2d = lindgen(jpi, jpj)
250         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
251;------------------------------------------------------------
252;------------------------------------------------------------
253         case 1 of
254            (size(u))[0] NE 2 OR (size(v))[0] NE 2: return,  -1
255            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
256             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
257               if nxu NE nx then $
258                if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
259               IF nxv NE nx THEN $
260                if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
261               IF nyu NE ny THEN $
262                if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
263               IF nyv NE ny THEN $
264                if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
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; Calculation of the curl
275;------------------------------------------------------------
276         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
277         terreu = where(coefu EQ 0)
278         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
279         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
280         terrev = where(coefv EQ 0)
281         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
282         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
283;
284         zu = u*temporary(coefu)
285         zv = v*temporary(coefv)
286
287         psi =  (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1))
288         psi = tabf*psi
289
290;------------------------------------------------------------
291; Edging put at !values.f_nan
292;------------------------------------------------------------
293         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
294            psi[0, *] = !values.f_nan
295            psi[nx-1, *] = !values.f_nan
296         endif
297         psi[*, 0] = !values.f_nan
298         psi[*, ny-1] = !values.f_nan
299;
300         if n_elements(valmask) EQ 0 THEN valmask = 1e20
301         terref =  where(tabf EQ 0)
302         if terref[0] NE -1 then psi[temporary(terref)] = valmask
303;------------------------------------------------------------
304; for the graphic drawing
305;------------------------------------------------------------
306         domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
307         if keyword_set(direc) then psi = moyenne(psi,direc,/nan)
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.