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

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

english and nicer header (3a)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.4 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 on matrixes
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 matrixes can be 2 or 4d.
26; Beware, to discern differents 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 hase not the same
33; size and refered to different indexes. In this case, arrays are recut on
34; common indexesand the domain is redifined to match with these common
35; indexes. To avoid these recuts, 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 pour marcher avec un domaine reduit
45;
46;                      21/5/1999: missing values at !values.f_nan
47;periodicite
48;
49; @version
50; $Id$
51;
52;-
53;------------------------------------------------------------
54;------------------------------------------------------------
55;------------------------------------------------------------
56FUNCTION curl, uu, vv
57;
58  compile_opt idl2, strictarrsubs
59;
60@common
61   tempsun = systime(1)         ; To 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; We find common points between U and V
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 of U and V on the appropriated domain
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; calculation of the curl
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; Edging put at !values.f_nan
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; For the graphic drawing
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 of U and V on the appropriated domain
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; Calculation of the curl
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; Calculation of the curl
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; extraction of U and V on the appropriated domain
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;------------------------------------------------------------
255         case 1 of
256            (size(u))[0] NE 2 OR (size(v))[0] NE 2: return,  -1
257            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
258             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
259               if nxu NE nx then $
260                if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
261               IF nxv NE nx THEN $
262                if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
263               IF nyu NE ny THEN $
264                if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
265               IF nyv NE ny THEN $
266                if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
267            END
268            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
269             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
270               u = u[indice2d]
271               v = v[indice2d]
272            END
273            ELSE:return,  -1
274         endcase
275;------------------------------------------------------------
276; Calculation of the curl
277;------------------------------------------------------------
278         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
279         terreu = where(coefu EQ 0)
280         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
281         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
282         terrev = where(coefv EQ 0)
283         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
284         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
285;
286         zu = u*temporary(coefu)
287         zv = v*temporary(coefv)
288
289         psi =  (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1))
290         psi = tabf*psi
291
292;------------------------------------------------------------
293; Edging put at !values.f_nan
294;------------------------------------------------------------
295         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
296            psi[0, *] = !values.f_nan
297            psi[nx-1, *] = !values.f_nan
298         endif
299         psi[*, 0] = !values.f_nan
300         psi[*, ny-1] = !values.f_nan
301;
302         if n_elements(valmask) EQ 0 THEN valmask = 1e20
303         terref =  where(tabf EQ 0)
304         if terref[0] NE -1 then psi[temporary(terref)] = valmask
305;------------------------------------------------------------
306; for the graphic drawing
307;------------------------------------------------------------
308         domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
309         if keyword_set(direc) then psi = moyenne(psi,direc,/nan)
310
311      END
312;----------------------------------------------------------------------------
313   endcase
314;------------------------------------------------------------
315   if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun
316;------------------------------------------------------------
317   vargrid = 'F'
318   varname = 'vorticity'
319   return, psi
320end
Note: See TracBrowser for help on using the repository browser.