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

Last change on this file since 231 was 231, checked in by pinsard, 17 years ago

improvements/corrections of some *.pro headers

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