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

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

improvements of headers (typo, links, paragraphes, etc)

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