source: trunk/SRC/ToBeReviewed/GRILLE/grille.pro @ 114

Last change on this file since 114 was 114, checked in by smasson, 18 years ago

new compilation options (compile_opt idl2, strictarrsubs) in each routine
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 15.0 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME:grille
6;
7; PURPOSE: choisit la grille qui doit etre utilisee pour faire le graphe en
8; fonction de vargrid et renvoie les parametres correspondants calcules ds
9; domdef.pro et reduit au domaine definit par domdef (contrairement a
10; grandegrille.pro)
11;
12; CATEGORY:
13;
14; CALLING SEQUENCE:
15;  grille,mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz,e1,e2,e3
16;
17; INPUTS:rien. ATTENTION les choix de la grille se fait a partir de la
18; valeur de la variable globale vargrid, qui peut etre egale a 'T',
19; 'U', 'V', 'W' ou 'F'.
20;
21; KEYWORD PARAMETERS:
22;         TRI si ce mot clef sert a obtenir grace a grille la
23;         triangulation qui se rapporte a la grille mais uniquement
24;         sur la partie du zoom. ce tableau de triangulation reduit
25;         est passe ds la variable que l''on a egalee a tri.par ex:
26;         grille,...,tri=triangulation_reduite. ne mot clef est
27;         utilise dans plt.pro
28;
29;         /FORPLT: ds plt on veut que sur les points terres, glam et
30;         gphi soit egale a glamt et gphit quelle que soit la grille.
31;
32;         /NOTRI: utile seulement qd TRI est active. dans ce cas
33;         grille retourne -1 ds la variable tri meme si la variable du
34;         common triangles_list est definie et differente de -1
35;
36;         /WDEPTH: to specify that the field is at W depth instad of T
37;         depth (automatically activated if vargrid eq 'W')
38;
39; OUTPUTS:mask,glam,gphi,gdep,nx,ny,nz,firstx,firsty,firstz,
40;         lastx,lasty,lastz,e1,e2,e3
41;
42;         pour leur definition cf domdef et la gestion des sous
43;         domaines sur le web
44;
45;         Rq: ces outputs sont optionnels, si je veux recuperer que
46;         mask, glam et gphi il suffit de taper grille, mask, glam, gphi
47;
48; COMMON BLOCKS: cm_4mesh and cm_4data
49;
50; SIDE EFFECTS: utilise la variable globale vargird
51;
52; RESTRICTIONS: vargrid doit etre 'T', 'W', 'U', 'V' ou 'F'
53;
54; EXAMPLE:
55;
56; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
57;                       12/2/1999
58;                       10/11/1999 /forplt
59;-
60;------------------------------------------------------------
61;------------------------------------------------------------
62;------------------------------------------------------------
63pro grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, TRI = tri, NOTRI = notri, TOUT = tout, FORPLT = forplt, IFPLTZ = ifpltz, WDEPTH = wdepth, _EXTRA = ex
64;------------------------------------------------------------
65; include commons
66;
67  compile_opt idl2, strictarrsubs
68;
69@cm_4mesh
70@cm_4data
71  IF NOT keyword_set(key_forgetold) THEN BEGIN
72@updatenew
73  ENDIF
74;---------------------
75  tempsun = systime(1)          ; pour key_performance
76;------------------------------------------------------------
77  vargrid = strupcase(strmid(vargrid,0,/reverse_offset))
78;
79  if vargrid eq 'W' then wdepth = 1
80  if keyword_set(tout) then begin
81    savedbox = 1b
82    saveboxparam, 'boxparam4grille.dat'
83    domdef, gridtype = vargrid, _EXTRA = ex
84  endif
85  tempdeux = systime(1)         ; pour key_performance =2
86;------------------------------------------------------------
87;------------------------------------------------------------
88  IF keyword_set(wdepth) THEN BEGIN
89    firstz = firstzw
90    lastz = lastzw
91    nz = nzw
92  ENDIF ELSE BEGIN
93    firstz = firstzt
94    lastz = lastzt
95    nz = nzt
96  ENDELSE
97;------------------------------------------------------------
98;------------------------------------------------------------
99  CASE 1 OF
100;------------------------------------------------------------
101; grille T and W
102;------------------------------------------------------------
103    vargrid eq 'T' OR vargrid eq 'W' : begin
104;scalaires
105      nx = nxt
106      ny = nyt
107      firstx = firstxt
108      firsty = firstyt
109      lastx = lastxt
110      lasty = lastyt
111;vecteurs 2d
112      IF arg_present(glam) THEN glam = glamt[firstx:lastx, firsty:lasty]
113      IF arg_present(gphi) THEN gphi = gphit[firstx:lastx, firsty:lasty]
114      IF arg_present(e1) THEN e1 = e1t[firstx:lastx, firsty:lasty]
115      IF arg_present(e2) THEN e2 = e2t[firstx:lastx, firsty:lasty]
116;vecteurs 3d
117      IF keyword_set(forplt) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz] $
118      ELSE IF arg_present(mask) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
119    end
120;------------------------------------------------------------
121; grille U
122;------------------------------------------------------------
123    vargrid eq 'U': begin
124;scalaires
125      nx = nxu
126      ny = nyu
127      firstx = firstxu
128      firsty = firstyu
129      lastx = lastxu
130      lasty = lastyu
131;vecteurs 2d
132      IF arg_present(glam) THEN glam = glamu[firstx:lastx, firsty:lasty]
133      IF arg_present(gphi) THEN gphi = gphiu[firstx:lastx, firsty:lasty]
134      if keyword_set(forplt) then BEGIN
135        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
136        eastboarder = mask-shift(mask, 1, 0)*mask
137        westboarder = mask-shift(mask, -1, 0)*mask
138        if key_periodic NE 1 OR nx NE jpi then westboarder[nx-1, *] = 0b
139        tmp1 = shift(eastboarder, 0, 1)
140        tmp1[*, 0] = 0b
141        tmp2 = shift(eastboarder, 0, -1)
142        tmp2[*, ny-1] = 0b
143        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
144        eastboarder = temporary(eastboarder)+temporary(add)
145        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
146        tmp1[*, ny-1] = 1b
147        tmp1[*, 0] = 1b
148        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
149        if key_periodic NE 1 OR nx NE jpi then begin
150          tmp2[nx-1, *] = 1b
151          tmp2[0, *] = 0b
152        endif
153        no1 = temporary(tmp1)*temporary(tmp2)
154        tmp = temporary(eastboarder)*temporary(no1)*mask
155        mask[0:nx-2, *] = 0b
156        tmp = temporary(tmp)+temporary(mask)
157        tmp = where(tmp GE 1)
158        if tmp[0] NE -1 then begin
159          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
160          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
161        endif
162      ENDIF
163      IF arg_present(e1) THEN e1  = e1u[firstx:lastx, firsty:lasty]
164      IF arg_present(e2) THEN e2  = e2u[firstx:lastx, firsty:lasty]
165;vecteurs 3d
166      IF keyword_set(forplt) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz] $
167      ELSE IF arg_present(mask) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
168    end
169;------------------------------------------------------------
170; grille V
171;------------------------------------------------------------
172    vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $
173      or vargrid eq 'V': begin
174;scalaires
175      nx = nxv
176      ny = nyv
177      firstx = firstxv
178      firsty = firstyv
179      lastx = lastxv
180      lasty = lastyv
181;vecteurs 2d
182      IF arg_present(glam) THEN glam = glamv[firstx:lastx, firsty:lasty]
183      IF arg_present(gphi) THEN gphi = gphiv[firstx:lastx, firsty:lasty]
184      if keyword_set(forplt) then BEGIN
185        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
186        northboarder = mask-shift(mask, 0, 1)*mask
187        southboarder = mask-shift(mask, 0, -1)*mask
188        southboarder[*, ny-1] = 0b
189        tmp1 = shift(northboarder, -1, 0)
190        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
191        tmp2 = shift(northboarder, 1, 0)
192        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
193        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
194        northboarder = temporary(northboarder)+temporary(add)
195        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
196        tmp1[*, ny-1] = 1b
197        tmp1[*, 0] = 0b
198        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
199        if key_periodic NE 1 OR nx NE jpi then begin
200          tmp2[nx-1, *] = 1b
201          tmp2[0, *] = 1b
202        endif
203        no1 = temporary(tmp1)*temporary(tmp2)
204        tmp = temporary(northboarder)*mask*temporary(no1)
205        mask[*, 0:ny-2] = 0b
206        tmp = temporary(tmp)+temporary(mask)
207        tmp = where(tmp GE 1)
208        if tmp[0] NE -1 then begin
209          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
210          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
211        endif
212      ENDIF
213      IF arg_present(e1) THEN e1  = e1v[firstx:lastx, firsty:lasty]
214      IF arg_present(e2) THEN e2  = e2v[firstx:lastx, firsty:lasty]
215;vecteurs 3d
216      IF keyword_set(forplt) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz] $
217      ELSE IF arg_present(mask) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
218    end
219;------------------------------------------------------------
220; grille F
221;------------------------------------------------------------
222    vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $
223      or vargrid eq 'F': begin
224;scalaires
225      nx = nxf
226      ny = nyf
227      firstx = firstxf
228      firsty = firstyf
229      lastx = lastxf
230      lasty = lastyf
231;vecteurs 2d
232      IF arg_present(glam) THEN glam = glamf[firstx:lastx, firsty:lasty]
233      IF arg_present(gphi) THEN gphi = gphif[firstx:lastx, firsty:lasty]
234      if keyword_set(forplt) then BEGIN
235        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
236        eastboarder = mask-shift(mask, 1, 0)*mask
237        westboarder = mask-shift(mask, -1, 0)*mask
238        westboarder[nx-1, *] = 0b
239        northboarder = mask-shift(mask, 0, 1)*mask
240        southboarder = mask-shift(mask, 0, -1)*mask
241        southboarder[*, ny-1] = 0b
242        tmp1 = shift(northboarder, -1, 0)
243        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
244        tmp2 = shift(northboarder, 1, 0)
245        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
246        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
247        northboarder = temporary(northboarder)+temporary(add)
248        tmp1 = shift(eastboarder, 0, 1)
249        tmp1[*, 0] = 0b
250        tmp2 = shift(eastboarder, 0, -1)
251        tmp2[*, ny-1] = 0b
252        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
253        eastboarder = temporary(eastboarder)+temporary(add)
254        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
255        tmp1[*, ny-1] = 1b
256        tmp1[*, 0] = 1b
257        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
258        if key_periodic NE 1 OR nx NE jpi then begin
259          tmp2[nx-1, *] = 1b
260          tmp2[0, *] = 1b
261        endif
262        no1 = temporary(tmp1)*temporary(tmp2)
263        tmp = (temporary(northboarder)+temporary(eastboarder))*mask*temporary(no1)
264        mask[0:nx-2, *] = 0b
265        mask[*, 0:ny-2] = 0b
266        tmp = temporary(tmp)+temporary(mask)
267        tmp = where(tmp GE 1)
268        if tmp[0] NE -1 then begin
269          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
270          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
271        endif
272      ENDIF
273      IF arg_present(e1) THEN e1  = e1f[firstx:lastx, firsty:lasty]
274      IF arg_present(e2) THEN e2  = e2f[firstx:lastx, firsty:lasty]
275;vecteurs 3d
276      IF keyword_set(forplt) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz] $
277      ELSE IF arg_present(mask) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
278    END
279;------------------------------------------------------------
280    ELSE:BEGIN
281      ras = report('Wrong definition of vargrid = '+vargrid+'. Only T, U, V, W or F are acceptable')
282      stop
283    END
284  ENDCASE
285  IF testvar(var = key_performance) EQ 2 THEN $
286    print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux
287;
288;------------------------------------------------------------
289;------------------------------------------------------------
290;------------------------------------------------------------
291; Variables se rapportant a la dimension verticale
292;------------------------------------------------------------
293;------------------------------------------------------------
294;------------------------------------------------------------
295;
296;
297  tempdeux = systime(1)         ; pour key_performance =2
298  if keyword_set(wdepth) then begin
299    gdep = gdepw[firstz:lastz]
300    e3 = e3w[firstz:lastz]
301  endif else begin
302    gdep = gdept[firstz:lastz]
303    e3 = e3t[firstz:lastz]
304  ENDELSE
305; for the vertical sections with partial steps
306  IF keyword_set(ifpltz) AND keyword_set(key_partialstep) THEN BEGIN
307    CASE 1 OF
308      ifpltz EQ 'xz' AND ny EQ 1:BEGIN
309        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
310        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
311        bottom = lindgen(nx)+(bottom-1l+keyword_set(wdepth))*nx
312        IF good[0] NE -1 THEN BEGIN
313          bottom = bottom[good]
314          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
315          gdep = replicate(1, nx)#gdep
316          if keyword_set(wdepth) THEN $
317            truegdep = hdepw[firstx:lastx, firsty:lasty] $
318          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
319          gdep[bottom] = truegdep[good]
320        ENDIF
321      END
322      ifpltz EQ 'yz' AND nx EQ 1:BEGIN
323        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
324        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
325        bottom = lindgen(ny)+(bottom-1l+keyword_set(wdepth))*ny
326        IF good[0] NE -1 THEN BEGIN
327          bottom = bottom[good]
328          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
329          gdep = replicate(1, ny)#gdep
330          if keyword_set(wdepth) THEN $
331            truegdep = hdepw[firstx:lastx, firsty:lasty] $
332          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
333          gdep[bottom] = truegdep[good]
334        ENDIF
335      END
336      ELSE:
337    ENDCASE
338  ENDIF
339  IF testvar(var = key_performance) EQ 2 THEN $
340    print, 'temps grille: Variables se rapportant a la dimension verticale ', systime(1)-tempdeux
341;------------------------------------------------------------
342; vecteur triangulation Qd TRI est active
343;------------------------------------------------------------
344  if arg_present(TRI) then $
345    if triangles_list[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN
346    tempdeux = systime(1)       ; pour key_performance =2
347    msk = bytarr(jpi, jpj)
348    msk[firstx:lastx, firsty:lasty] = 1
349    ind = where( msk[triangles_list[0, *]]*msk[triangles_list[1, *]]*msk[triangles_list[2, *]] EQ 1 )
350    tri = triangles_list[*, ind]-(firstx+firsty*jpi)
351    y = tri/jpi
352    x = tri-y*jpi
353    tri = x+y*nx
354    IF testvar(var = key_performance) EQ 2 THEN $
355      print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux
356  ENDELSE
357;------------------------------------------------------------------
358; pour s'assurer qu'il n'y a pas de dimension degenerees (=1)
359;-------------------------------------------------------------------
360;    mask=reform(mask, /over)
361;    glam=reform(glam, /over)
362;    gphi=reform(gphi, /over)
363;    gdep=reform(gdep, /over)
364;    e1=reform(e1, /over)
365;    e2=reform(e2, /over)
366;    e3=reform(e3, /over)
367
368  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grille.dat'
369  if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun
370
371;------------------------------------------------------------
372  IF NOT keyword_set(key_forgetold) THEN BEGIN
373@updateold
374  ENDIF
375;---------------------
376  return
377
378end
379
380
381
382
383
384
385
Note: See TracBrowser for help on using the repository browser.