source: trunk/SRC/ToBeReviewed/CALCULS/div.pro @ 159

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

header improvements + xxx doc
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.0 KB
RevLine 
[2]1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5;
[142]6; @file_comments
7; calculation of the divergence of a 2d field
[2]8;
[142]9; @categories
[157]10; Calculation
[2]11;
[142]12; @param UU
13; Matrix representing coordinates of a field of vectors
[2]14;
[142]15; @param VV
16; Matrix representing coordinates of a field of vectors
[2]17;
[142]18; @returns RES
19; A 2d matrix
[2]20;
[142]21; @uses
22; common.pro
[2]23;
[142]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
[2]36;
37;
[142]38; Points on the drawing edge are at !values.f_nan
[2]39;
40;
[142]41; @history
[157]42; Guillaume Roullet (grlod\@ipsl.jussieu.fr)
[2]43;                      Creation : printemps 1998
[157]44;                      Sebastien Masson (smasson\@lodyc.jussieu.fr)
[2]45;                      adaptation pour marcher avec un domaine reduit
[142]46;                      12/1/2000;
47;
48; @version
49; $Id$
50;
[2]51;-
52;------------------------------------------------------------
53;------------------------------------------------------------
54;------------------------------------------------------------
55FUNCTION div, uu, vv
[114]56;
57  compile_opt idl2, strictarrsubs
58;
[142]59   tempsun = systime(1)         ; For key_performance
[2]60@common
[25]61;
62   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
63     return, report(['This version of div is based on Arakawa C-grid.' $
64                     , 'U and V grids must therefore be defined'])
65;
66   u = litchamp(uu)
67   v = litchamp(vv)
68;
[2]69   date1 = time[0]
70   if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
71   if (size(u))[0] NE (size(v))[0] then return,  -1
72
73;------------------------------------------------------------
[142]74; We find common points between U and V
[2]75;------------------------------------------------------------
[25]76   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
77   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
[2]78   indicex = inter(indicexu, indicexv)
[25]79   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
80   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
[2]81   indicey = inter(indiceyu, indiceyv)
82   nx = n_elements(indicex)
83   ny = n_elements(indicey)
84   indice2d = lindgen(jpi, jpj)
85   indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
86   case 1 of
87;----------------------------------------------------------------------------
88;----------------------------------------------------------------------------
89;xyz
90;----------------------------------------------------------------------------
91      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN
92;------------------------------------------------------------
[142]93; extraction of U and V on the appropriated domain
[2]94;------------------------------------------------------------
95         case 1 of
[25]96            (size(v))[0] NE 3: return,  -1
[2]97            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
98             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
99               case 1 of
[25]100                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
101                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
102                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
103                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
[2]104                  ELSE :
105               endcase
106            END
107            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
108             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
109               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
110               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
111            END
112            ELSE:BEGIN
113               zdiv = -1
114               GOTO, sortie
115            end
116           
117         endcase
118;------------------------------------------------------------
119; calcul de la divergence
120;------------------------------------------------------------
121         zu = (e2u[indice2d])[*]#replicate(1, nzt)
122         zu = reform(zu, nx, ny, nzt, /over)
123         zu = temporary(u)*temporary(zu) $
[25]124          *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
[2]125         terreu = where(zu EQ 0)
126         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
127;
128         zv = (e1v[indice2d])[*]#replicate(1, nzt)
129         zv = reform(zv, nx, ny, nzt, /over)
130         zv = temporary(v)*temporary(zv) $
[25]131          *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
[2]132         terrev = where(zv EQ 0)
133         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
134;
135         zdiv = 1e6/(e1t[indice2d]*e2t[indice2d])
136         zdiv = (zdiv)[*]#replicate(1, nzt)
137         zdiv = reform(zdiv, nx, ny, nzt, /over)
138         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $
[25]139          *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
[2]140;------------------------------------------------------------
[142]141; Edging put at !values.f_nan
[2]142;------------------------------------------------------------
[25]143         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
[114]144            zdiv[0, *, *] = !values.f_nan
145            zdiv[nx-1, *, *] = !values.f_nan
[2]146         endif
[114]147         zdiv[*, 0, *] = !values.f_nan
148         zdiv[*, ny-1, *] = !values.f_nan
[2]149;
150         zdiv = temporary(zdiv)
151;
152         if n_elements(valmask) EQ 0 THEN valmask = 1e20
[25]153         terre =  where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0)
[2]154         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
155;------------------------------------------------------------
[142]156; For the graphic drawing
[2]157;------------------------------------------------------------
158         vargrid = 'T'
159         varname = 'div'
160         varunits = '1e6*s-1'
[25]161         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
162         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan)
[2]163      END
164;----------------------------------------------------------------------------
165;----------------------------------------------------------------------------
166;xyt
167;----------------------------------------------------------------------------
168;----------------------------------------------------------------------------
169      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN
170;------------------------------------------------------------
[142]171; extraction of U and V on the appropriated domain
[2]172;------------------------------------------------------------
173         case 1 of
174            (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1
175            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
176             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
177               case 1 of
[25]178                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
179                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
180                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
181                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
[2]182                  ELSE :
183               endcase
184            END
185            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
186             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN
187               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
188               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
189            END
190            ELSE:return,  -1
191         endcase
192;------------------------------------------------------------
[142]193; Calculation of the divergence
[2]194;------------------------------------------------------------
[25]195         zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
[2]196         terreu = where(zu EQ 0)
197         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
198         zu = (zu)[*]#replicate(1, jpt)
199         zu = reform(zu, nx, ny, jpt, /over)
200         zu = temporary(u)*temporary(zu)
201;
[25]202         zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
[2]203         terrev = where(zv EQ 0)
204         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
205         zv = (zv)[*]#replicate(1, jpt)
206         zv = reform(zv, nx, ny, jpt, /over)
207         zv = temporary(v)*temporary(zv)
208;
[25]209         zdiv = 1e6*tmask[indice2d+jpi*jpj*firstzt]/(e1t[indice2d]*e2t[indice2d])
[2]210         zdiv = (zdiv)[*]#replicate(1, jpt)
211         zdiv = reform(zdiv, nx, ny, jpt, /over)
212         terre =  where(zdiv EQ 0)
213         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0))
214;------------------------------------------------------------
[142]215; Edging put at !values.f_nan
[2]216;------------------------------------------------------------
[25]217         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
[114]218            zdiv[0, *, *] = !values.f_nan
219            zdiv[nx-1, *, *] = !values.f_nan
[2]220         endif
[114]221         zdiv[*, 0, *] = !values.f_nan
222         zdiv[*, ny-1, *] = !values.f_nan
[2]223;
224         if n_elements(valmask) EQ 0 THEN valmask = 1e20
225         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
226;------------------------------------------------------------
[142]227; for the graphic drawing
[2]228;------------------------------------------------------------
229         vargrid = 'T'
230         varname = 'div'
231         varunits = '1e6*s-1'
[25]232         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
233         if keyword_set(direc) then  zdiv = grossemoyenne(zdiv,direc,/nan)
[2]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         indice3d = lindgen(jpi, jpj, jpk)
[25]250         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt]
[2]251;------------------------------------------------------------
[142]252; extraction of U and V on the appropriated domain
[2]253;------------------------------------------------------------
254         case 1 of
255            (size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN
256               zdiv = -1
257               GOTO, sortie
258            end
259            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
260             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN
261               case 1 of
[25]262                  nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
263                  nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
264                  nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
265                  nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
[2]266                  ELSE :
267               endcase
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;------------------------------------------------------------
[142]277; Calculation of the divergence
[2]278;------------------------------------------------------------
279         zu = temporary(u)*e2u[indice2d]*(umask())[indice3d]
280         terreu = where(zu EQ 0)
281         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
282         zv = temporary(v)*e1v[indice2d]*(vmask())[indice3d]
283         terrev = where(zv EQ 0)
284         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
285         zdiv = zu-shift(zu, 1, 0)+zv-shift(zv, 0, 1)
286         zdiv = temporary(zdiv)*tmask[indice3d]/(e1t[indice2d]*e2t[indice2d])
287;------------------------------------------------------------
[142]288; Edging put at !values.f_nan
[2]289;------------------------------------------------------------
[25]290         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
[114]291            zdiv[0, *] = !values.f_nan
292            zdiv[nx-1, *] = !values.f_nan
[2]293         endif
[114]294         zdiv[*, 0] = !values.f_nan
295         zdiv[*, ny-1] = !values.f_nan
[2]296;
297         zdiv = temporary(zdiv)*1e6
298;
299         if n_elements(valmask) EQ 0 THEN valmask = 1e20
300         terre =  where(tmask[indice3d] EQ 0)
301         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
302;------------------------------------------------------------
[142]303; for the graphic drawing
[2]304;------------------------------------------------------------
305         vargrid = 'T'
306         varname = 'div'
307         varunits = '1e6*s-1'
[25]308         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t']
309         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan)
[2]310
311      END
312   endcase
313   
314;------------------------------------------------------------
315sortie:
316   if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun
317   
318   return, zdiv
319end
Note: See TracBrowser for help on using the repository browser.