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

Last change on this file since 371 was 371, checked in by pinsard, 16 years ago

improvements of headers (alignments of IDL prompt in examples)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 17.0 KB
Line 
1;+
2;
3; @file_comments
4; Choose the grid which must be used to do the graph in function of
5; vargrid and send back corresponding parameters calculated in
6; <pro>domdef</pro> and reduced at the domain defined by
7; <pro>domdef</pro> (contrarily to
8; <pro>grandegrille</pro>)
9; BEWARE!! The choice of the grid is made from the value of the
10; global variable vargrid, which can be equal to 'T', 'U', 'V', 'W' ou 'F'.
11;
12; @categories
13; Grid
14;
15; @keyword TRI
16; This keyword serve to obtain the triangulation which
17; refer to the grid but only on the part of the zoom. This array of triangulation
18; is passed in the variable we have equate at TRI.
19; For example: grille,...,tri=triangulation_reduite.
20; This keyword is used in <pro>plt</pro>.
21;
22; @keyword WDEPTH
23; To specify that the field is at W depth instead of T
24; depth (automatically activated if vargrid eq 'W')
25;
26; @keyword FORPLT
27; In <pro>plt</pro>, we want that land points,
28; glam and gphi, be equal to glamt and
29; gphit regardless of the grid.
30;
31; @keyword NOTRI
32; Useful only when TRI is activated. In this case, grille send back -1 in the
33; variable tri even if the variable of the common triangles_list is defined
34; and different of -1
35;
36; @keyword _EXTRA
37; Used to pass keywords
38;
39; @keyword TOUT
40;
41; @param MASK {out}{optional}
42; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
43;
44; @param GLAM {out}{optional}
45; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
46;
47; @param GPHI {out}{optional}
48; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
49;
50; @param GDEP {out}{optional}
51; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
52;
53; @param NX {out}{optional}
54; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
55;
56; @param NY {out}{optional}
57; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
58;
59; @param NZ {out}{optional}
60; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
61;
62; @param FIRSTX {out}{optional}
63; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
64;
65; @param FIRSTY {out}{optional}
66; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
67;
68; @param FIRSTZ {out}{optional}
69; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
70;
71; @param LASTX {out}{optional}
72; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
73;
74; @param LASTY {out}{optional}
75; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
76;
77; @param LASTZ {out}{optional}
78; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
79;
80; @param E1 {out}{optional}
81; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
82;
83; @param E2 {out}{optional}
84; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
85;
86; @param E3 {out}{optional}
87; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
88;
89; @uses
90; <pro>cm_4mesh</pro>
91; <pro>cm_4data</pro>
92;
93; @restrictions
94; Use the variable vargrid
95;
96; @restrictions
97; Vargrid must be 'T', 'W', 'U', 'V' or 'F'
98;
99; @history
100;  Sebastien Masson (smasson\@lodyc.jussieu.fr)
101;                       12/2/1999
102;                       10/11/1999 /forplt
103;
104; @version
105; $Id$
106;
107; @todo
108; Comment ecrire la remarque sur les inputs?
109;
110;-
111PRO grille, mask, glam, gphi, gdep, nx, ny, nz $
112          , firstx, firsty, firstz, lastx, lasty, lastz $
113          , e1, e2, e3 $
114          , TRI=tri, NOTRI=notri, TOUT=tout, FORPLT=forplt $
115          , TYPE=type, WDEPTH=wdepth, _EXTRA=ex
116;
117  compile_opt idl2, strictarrsubs
118;
119@cm_4mesh
120@cm_4data
121@cm_4cal ; for jpt
122  IF NOT keyword_set(key_forgetold) THEN BEGIN
123@updatenew
124  ENDIF
125;---------------------
126  tempsun = systime(1)          ; For key_performance
127;------------------------------------------------------------
128  vargrid = strupcase(strmid(vargrid,0,/reverse_offset))
129;
130  if vargrid eq 'W' then wdepth = 1
131  if keyword_set(tout) then begin
132    savedbox = 1b
133    saveboxparam, 'boxparam4grille.dat'
134    domdef, gridtype = vargrid, _EXTRA = ex
135  endif
136  tempdeux = systime(1)         ; For key_performance =2
137;------------------------------------------------------------
138;------------------------------------------------------------
139  IF keyword_set(wdepth) THEN BEGIN
140    firstz = firstzw
141    lastz = lastzw
142    nz = nzw
143  ENDIF ELSE BEGIN
144    firstz = firstzt
145    lastz = lastzt
146    nz = nzt
147  ENDELSE
148;------------------------------------------------------------
149;------------------------------------------------------------
150  CASE 1 OF
151;------------------------------------------------------------
152; grid T and W
153;------------------------------------------------------------
154    vargrid eq 'T' OR vargrid eq 'W' : begin
155;scalars
156      nx = nxt
157      ny = nyt
158      firstx = firstxt
159      firsty = firstyt
160      lastx = lastxt
161      lasty = lastyt
162;2d vectors
163      IF arg_present(glam) THEN glam = glamt[firstx:lastx, firsty:lasty]
164      IF arg_present(gphi) THEN gphi = gphit[firstx:lastx, firsty:lasty]
165      IF arg_present(e1) THEN e1 = e1t[firstx:lastx, firsty:lasty]
166      IF arg_present(e2) THEN e2 = e2t[firstx:lastx, firsty:lasty]
167;3d vectors
168      IF keyword_set(forplt) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz] $
169      ELSE IF arg_present(mask) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
170    end
171;------------------------------------------------------------
172; grid U
173;------------------------------------------------------------
174    vargrid eq 'U': begin
175;scalars
176      nx = nxu
177      ny = nyu
178      firstx = firstxu
179      firsty = firstyu
180      lastx = lastxu
181      lasty = lastyu
182;2d vectors
183      IF arg_present(glam) THEN glam = glamu[firstx:lastx, firsty:lasty]
184      IF arg_present(gphi) THEN gphi = gphiu[firstx:lastx, firsty:lasty]
185      if keyword_set(forplt) then BEGIN
186        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
187        eastboarder = mask-shift(mask, 1, 0)*mask
188        westboarder = mask-shift(mask, -1, 0)*mask
189        if key_periodic NE 1 OR nx NE jpi then westboarder[nx-1, *] = 0b
190        tmp1 = shift(eastboarder, 0, 1)
191        tmp1[*, 0] = 0b
192        tmp2 = shift(eastboarder, 0, -1)
193        tmp2[*, ny-1] = 0b
194        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
195        eastboarder = temporary(eastboarder)+temporary(add)
196        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
197        tmp1[*, ny-1] = 1b
198        tmp1[*, 0] = 1b
199        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
200        if key_periodic NE 1 OR nx NE jpi then begin
201          tmp2[nx-1, *] = 1b
202          tmp2[0, *] = 0b
203        endif
204        no1 = temporary(tmp1)*temporary(tmp2)
205        tmp = temporary(eastboarder)*temporary(no1)*mask
206        mask[0:nx-2, *] = 0b
207        tmp = temporary(tmp)+temporary(mask)
208        tmp = where(tmp GE 1)
209        if tmp[0] NE -1 then begin
210          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
211          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
212        endif
213      ENDIF
214      IF arg_present(e1) THEN e1  = e1u[firstx:lastx, firsty:lasty]
215      IF arg_present(e2) THEN e2  = e2u[firstx:lastx, firsty:lasty]
216;3d vectors
217      IF keyword_set(forplt) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz] $
218      ELSE IF arg_present(mask) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
219    end
220;------------------------------------------------------------
221; grid V
222;------------------------------------------------------------
223    vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $
224      or vargrid eq 'V': begin
225;scalars
226      nx = nxv
227      ny = nyv
228      firstx = firstxv
229      firsty = firstyv
230      lastx = lastxv
231      lasty = lastyv
232;2d vectors
233      IF arg_present(glam) THEN glam = glamv[firstx:lastx, firsty:lasty]
234      IF arg_present(gphi) THEN gphi = gphiv[firstx:lastx, firsty:lasty]
235      if keyword_set(forplt) then BEGIN
236        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
237        northboarder = mask-shift(mask, 0, 1)*mask
238        southboarder = mask-shift(mask, 0, -1)*mask
239        southboarder[*, ny-1] = 0b
240        tmp1 = shift(northboarder, -1, 0)
241        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
242        tmp2 = shift(northboarder, 1, 0)
243        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
244        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
245        northboarder = temporary(northboarder)+temporary(add)
246        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
247        tmp1[*, ny-1] = 1b
248        tmp1[*, 0] = 0b
249        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
250        if key_periodic NE 1 OR nx NE jpi then begin
251          tmp2[nx-1, *] = 1b
252          tmp2[0, *] = 1b
253        endif
254        no1 = temporary(tmp1)*temporary(tmp2)
255        tmp = temporary(northboarder)*mask*temporary(no1)
256        mask[*, 0:ny-2] = 0b
257        tmp = temporary(tmp)+temporary(mask)
258        tmp = where(tmp GE 1)
259        if tmp[0] NE -1 then begin
260          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
261          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
262        endif
263      ENDIF
264      IF arg_present(e1) THEN e1  = e1v[firstx:lastx, firsty:lasty]
265      IF arg_present(e2) THEN e2  = e2v[firstx:lastx, firsty:lasty]
266;3d vecteurs
267      IF keyword_set(forplt) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz] $
268      ELSE IF arg_present(mask) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
269    end
270;------------------------------------------------------------
271; grid F
272;------------------------------------------------------------
273    vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $
274      or vargrid eq 'F': begin
275;scalars
276      nx = nxf
277      ny = nyf
278      firstx = firstxf
279      firsty = firstyf
280      lastx = lastxf
281      lasty = lastyf
282;2d vectors
283      IF arg_present(glam) THEN glam = glamf[firstx:lastx, firsty:lasty]
284      IF arg_present(gphi) THEN gphi = gphif[firstx:lastx, firsty:lasty]
285      if keyword_set(forplt) then BEGIN
286        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
287        eastboarder = mask-shift(mask, 1, 0)*mask
288        westboarder = mask-shift(mask, -1, 0)*mask
289        westboarder[nx-1, *] = 0b
290        northboarder = mask-shift(mask, 0, 1)*mask
291        southboarder = mask-shift(mask, 0, -1)*mask
292        southboarder[*, ny-1] = 0b
293        tmp1 = shift(northboarder, -1, 0)
294        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
295        tmp2 = shift(northboarder, 1, 0)
296        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
297        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
298        northboarder = temporary(northboarder)+temporary(add)
299        tmp1 = shift(eastboarder, 0, 1)
300        tmp1[*, 0] = 0b
301        tmp2 = shift(eastboarder, 0, -1)
302        tmp2[*, ny-1] = 0b
303        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
304        eastboarder = temporary(eastboarder)+temporary(add)
305        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
306        tmp1[*, ny-1] = 1b
307        tmp1[*, 0] = 1b
308        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
309        if key_periodic NE 1 OR nx NE jpi then begin
310          tmp2[nx-1, *] = 1b
311          tmp2[0, *] = 1b
312        endif
313        no1 = temporary(tmp1)*temporary(tmp2)
314        tmp = (temporary(northboarder)+temporary(eastboarder))*mask*temporary(no1)
315        mask[0:nx-2, *] = 0b
316        mask[*, 0:ny-2] = 0b
317        tmp = temporary(tmp)+temporary(mask)
318        tmp = where(tmp GE 1)
319        if tmp[0] NE -1 then begin
320          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
321          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
322        endif
323      ENDIF
324      IF arg_present(e1) THEN e1  = e1f[firstx:lastx, firsty:lasty]
325      IF arg_present(e2) THEN e2  = e2f[firstx:lastx, firsty:lasty]
326;3d vectors
327      IF keyword_set(forplt) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz] $
328      ELSE IF arg_present(mask) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
329    END
330;------------------------------------------------------------
331    ELSE:BEGIN
332      ras = report('Wrong definition of vargrid = '+vargrid+'. Only T, U, V, W or F are acceptable')
333      stop
334    END
335  ENDCASE
336  IF testvar(var = key_performance) EQ 2 THEN $
337    print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux
338;
339;------------------------------------------------------------
340;------------------------------------------------------------
341;------------------------------------------------------------
342; Variables refering to the vertical dimension
343;------------------------------------------------------------
344;------------------------------------------------------------
345;------------------------------------------------------------
346;
347;
348  tempdeux = systime(1)         ; For key_performance =2
349  if keyword_set(wdepth) then begin
350    gdep = gdepw[firstz:lastz]
351    e3 = e3w[firstz:lastz]
352  endif else begin
353    gdep = gdept[firstz:lastz]
354    e3 = e3t[firstz:lastz]
355  ENDELSE
356; for the vertical sections with partial steps
357  IF keyword_set(type) AND keyword_set(key_partialstep) THEN BEGIN
358    CASE 1 OF
359      type EQ 'xz' AND ny EQ 1:BEGIN
360        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
361        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
362        bottom = lindgen(nx)+(bottom-1l+keyword_set(wdepth))*nx
363        IF good[0] NE -1 THEN BEGIN
364          bottom = bottom[good]
365          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
366          gdep = replicate(1, nx)#gdep
367          if keyword_set(wdepth) THEN $
368            truegdep = hdepw[firstx:lastx, firsty:lasty] $
369          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
370          gdep[bottom] = truegdep[good]
371        ENDIF
372      END
373      type EQ 'yz' AND nx EQ 1:BEGIN
374        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
375        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
376        bottom = lindgen(ny)+(bottom-1l+keyword_set(wdepth))*ny
377        IF good[0] NE -1 THEN BEGIN
378          bottom = bottom[good]
379          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
380          gdep = replicate(1, ny)#gdep
381          if keyword_set(wdepth) THEN $
382            truegdep = hdepw[firstx:lastx, firsty:lasty] $
383          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
384          gdep[bottom] = truegdep[good]
385        ENDIF
386      END
387      ELSE:
388    ENDCASE
389  ENDIF
390; for the vertical sections with roms
391  IF keyword_set(type) AND n_elements(romszinfos) NE 0 THEN BEGIN
392    romsdp = romsdepth()
393    IF romsdp[0] NE -1 THEN BEGIN
394      IF jpt EQ 1 THEN BEGIN
395        CASE type OF
396          'z':gdep = moyenne(temporary(romsdp), 'xy')
397          'xz':gdep = moyenne(temporary(romsdp), 'y')
398          'yz':gdep = moyenne(temporary(romsdp), 'x')
399          'zt':gdep = moyenne(temporary(romsdp), 'xy')
400          ELSE:
401        ENDCASE
402      ENDIF ELSE BEGIN
403        CASE type OF
404          'z':gdep = moyenne(temporary(romsdp), 'xyt')
405          'xz':gdep = grossemoyenne(temporary(romsdp), 'yt')
406          'yz':gdep = grossemoyenne(temporary(romsdp), 'xt')
407          'zt':gdep = grossemoyenne(temporary(romsdp), 'xy')
408          ELSE:
409        ENDCASE
410      ENDELSE
411      IF size(gdep, /n_dimensions) EQ 2 THEN $
412         gdep = remplit(gdep, niter = 32000, mask = gdep LT 1.E19, /basique, /fillxdir)
413    ENDIF
414  ENDIF
415;
416  IF testvar(var = key_performance) EQ 2 THEN $
417    print, 'temps grille: Variables se rapportant a la dimension verticale ', systime(1)-tempdeux
418;------------------------------------------------------------
419; Triangulation vector when TRI is activated.
420;------------------------------------------------------------
421  if arg_present(TRI) then $
422    if triangles_list[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN
423    tempdeux = systime(1)       ; pour key_performance =2
424    msk = bytarr(jpi, jpj)
425    msk[firstx:lastx, firsty:lasty] = 1
426    ind = where( msk[triangles_list[0, *]]*msk[triangles_list[1, *]]*msk[triangles_list[2, *]] EQ 1 )
427    tri = triangles_list[*, ind]-(firstx+firsty*jpi)
428    y = tri/jpi
429    x = tri-y*jpi
430    tri = x+y*nx
431    IF testvar(var = key_performance) EQ 2 THEN $
432      print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux
433  ENDELSE
434;------------------------------------------------------------------
435; To make sure there is not any degenerated dimension (=1)
436;------------------------------------------------------------------
437;    mask=reform(mask, /over)
438;    glam=reform(glam, /over)
439;    gphi=reform(gphi, /over)
440;    gdep=reform(gdep, /over)
441;    e1=reform(e1, /over)
442;    e2=reform(e2, /over)
443;    e3=reform(e3, /over)
444
445  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grille.dat'
446  if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun
447
448;------------------------------------------------------------
449  IF NOT keyword_set(key_forgetold) THEN BEGIN
450@updateold
451  ENDIF
452;---------------------
453  return
454
455end
Note: See TracBrowser for help on using the repository browser.