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

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

typo and links in header in some pro files

  • 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
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; cm_4mesh
91; cm_4data
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;-
111;
112PRO 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
113;
114  compile_opt idl2, strictarrsubs
115;
116@cm_4mesh
117@cm_4data
118@cm_4cal ; for jpt
119  IF NOT keyword_set(key_forgetold) THEN BEGIN
120@updatenew
121  ENDIF
122;---------------------
123  tempsun = systime(1)          ; For key_performance
124;------------------------------------------------------------
125  vargrid = strupcase(strmid(vargrid,0,/reverse_offset))
126;
127  if vargrid eq 'W' then wdepth = 1
128  if keyword_set(tout) then begin
129    savedbox = 1b
130    saveboxparam, 'boxparam4grille.dat'
131    domdef, gridtype = vargrid, _EXTRA = ex
132  endif
133  tempdeux = systime(1)         ; For key_performance =2
134;------------------------------------------------------------
135;------------------------------------------------------------
136  IF keyword_set(wdepth) THEN BEGIN
137    firstz = firstzw
138    lastz = lastzw
139    nz = nzw
140  ENDIF ELSE BEGIN
141    firstz = firstzt
142    lastz = lastzt
143    nz = nzt
144  ENDELSE
145;------------------------------------------------------------
146;------------------------------------------------------------
147  CASE 1 OF
148;------------------------------------------------------------
149; grid T and W
150;------------------------------------------------------------
151    vargrid eq 'T' OR vargrid eq 'W' : begin
152;scalars
153      nx = nxt
154      ny = nyt
155      firstx = firstxt
156      firsty = firstyt
157      lastx = lastxt
158      lasty = lastyt
159;2d vectors
160      IF arg_present(glam) THEN glam = glamt[firstx:lastx, firsty:lasty]
161      IF arg_present(gphi) THEN gphi = gphit[firstx:lastx, firsty:lasty]
162      IF arg_present(e1) THEN e1 = e1t[firstx:lastx, firsty:lasty]
163      IF arg_present(e2) THEN e2 = e2t[firstx:lastx, firsty:lasty]
164;3d vectors
165      IF keyword_set(forplt) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz] $
166      ELSE IF arg_present(mask) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
167    end
168;------------------------------------------------------------
169; grid U
170;------------------------------------------------------------
171    vargrid eq 'U': begin
172;scalars
173      nx = nxu
174      ny = nyu
175      firstx = firstxu
176      firsty = firstyu
177      lastx = lastxu
178      lasty = lastyu
179;2d vectors
180      IF arg_present(glam) THEN glam = glamu[firstx:lastx, firsty:lasty]
181      IF arg_present(gphi) THEN gphi = gphiu[firstx:lastx, firsty:lasty]
182      if keyword_set(forplt) then BEGIN
183        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
184        eastboarder = mask-shift(mask, 1, 0)*mask
185        westboarder = mask-shift(mask, -1, 0)*mask
186        if key_periodic NE 1 OR nx NE jpi then westboarder[nx-1, *] = 0b
187        tmp1 = shift(eastboarder, 0, 1)
188        tmp1[*, 0] = 0b
189        tmp2 = shift(eastboarder, 0, -1)
190        tmp2[*, ny-1] = 0b
191        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
192        eastboarder = temporary(eastboarder)+temporary(add)
193        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
194        tmp1[*, ny-1] = 1b
195        tmp1[*, 0] = 1b
196        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
197        if key_periodic NE 1 OR nx NE jpi then begin
198          tmp2[nx-1, *] = 1b
199          tmp2[0, *] = 0b
200        endif
201        no1 = temporary(tmp1)*temporary(tmp2)
202        tmp = temporary(eastboarder)*temporary(no1)*mask
203        mask[0:nx-2, *] = 0b
204        tmp = temporary(tmp)+temporary(mask)
205        tmp = where(tmp GE 1)
206        if tmp[0] NE -1 then begin
207          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
208          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
209        endif
210      ENDIF
211      IF arg_present(e1) THEN e1  = e1u[firstx:lastx, firsty:lasty]
212      IF arg_present(e2) THEN e2  = e2u[firstx:lastx, firsty:lasty]
213;3d vectors
214      IF keyword_set(forplt) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz] $
215      ELSE IF arg_present(mask) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
216    end
217;------------------------------------------------------------
218; grid V
219;------------------------------------------------------------
220    vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $
221      or vargrid eq 'V': begin
222;scalars
223      nx = nxv
224      ny = nyv
225      firstx = firstxv
226      firsty = firstyv
227      lastx = lastxv
228      lasty = lastyv
229;2d vectors
230      IF arg_present(glam) THEN glam = glamv[firstx:lastx, firsty:lasty]
231      IF arg_present(gphi) THEN gphi = gphiv[firstx:lastx, firsty:lasty]
232      if keyword_set(forplt) then BEGIN
233        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
234        northboarder = mask-shift(mask, 0, 1)*mask
235        southboarder = mask-shift(mask, 0, -1)*mask
236        southboarder[*, ny-1] = 0b
237        tmp1 = shift(northboarder, -1, 0)
238        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
239        tmp2 = shift(northboarder, 1, 0)
240        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
241        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
242        northboarder = temporary(northboarder)+temporary(add)
243        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
244        tmp1[*, ny-1] = 1b
245        tmp1[*, 0] = 0b
246        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
247        if key_periodic NE 1 OR nx NE jpi then begin
248          tmp2[nx-1, *] = 1b
249          tmp2[0, *] = 1b
250        endif
251        no1 = temporary(tmp1)*temporary(tmp2)
252        tmp = temporary(northboarder)*mask*temporary(no1)
253        mask[*, 0:ny-2] = 0b
254        tmp = temporary(tmp)+temporary(mask)
255        tmp = where(tmp GE 1)
256        if tmp[0] NE -1 then begin
257          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
258          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
259        endif
260      ENDIF
261      IF arg_present(e1) THEN e1  = e1v[firstx:lastx, firsty:lasty]
262      IF arg_present(e2) THEN e2  = e2v[firstx:lastx, firsty:lasty]
263;3d vecteurs
264      IF keyword_set(forplt) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz] $
265      ELSE IF arg_present(mask) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
266    end
267;------------------------------------------------------------
268; grid F
269;------------------------------------------------------------
270    vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $
271      or vargrid eq 'F': begin
272;scalars
273      nx = nxf
274      ny = nyf
275      firstx = firstxf
276      firsty = firstyf
277      lastx = lastxf
278      lasty = lastyf
279;2d vectors
280      IF arg_present(glam) THEN glam = glamf[firstx:lastx, firsty:lasty]
281      IF arg_present(gphi) THEN gphi = gphif[firstx:lastx, firsty:lasty]
282      if keyword_set(forplt) then BEGIN
283        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
284        eastboarder = mask-shift(mask, 1, 0)*mask
285        westboarder = mask-shift(mask, -1, 0)*mask
286        westboarder[nx-1, *] = 0b
287        northboarder = mask-shift(mask, 0, 1)*mask
288        southboarder = mask-shift(mask, 0, -1)*mask
289        southboarder[*, ny-1] = 0b
290        tmp1 = shift(northboarder, -1, 0)
291        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
292        tmp2 = shift(northboarder, 1, 0)
293        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
294        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
295        northboarder = temporary(northboarder)+temporary(add)
296        tmp1 = shift(eastboarder, 0, 1)
297        tmp1[*, 0] = 0b
298        tmp2 = shift(eastboarder, 0, -1)
299        tmp2[*, ny-1] = 0b
300        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
301        eastboarder = temporary(eastboarder)+temporary(add)
302        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
303        tmp1[*, ny-1] = 1b
304        tmp1[*, 0] = 1b
305        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
306        if key_periodic NE 1 OR nx NE jpi then begin
307          tmp2[nx-1, *] = 1b
308          tmp2[0, *] = 1b
309        endif
310        no1 = temporary(tmp1)*temporary(tmp2)
311        tmp = (temporary(northboarder)+temporary(eastboarder))*mask*temporary(no1)
312        mask[0:nx-2, *] = 0b
313        mask[*, 0:ny-2] = 0b
314        tmp = temporary(tmp)+temporary(mask)
315        tmp = where(tmp GE 1)
316        if tmp[0] NE -1 then begin
317          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
318          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
319        endif
320      ENDIF
321      IF arg_present(e1) THEN e1  = e1f[firstx:lastx, firsty:lasty]
322      IF arg_present(e2) THEN e2  = e2f[firstx:lastx, firsty:lasty]
323;3d vectors
324      IF keyword_set(forplt) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz] $
325      ELSE IF arg_present(mask) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
326    END
327;------------------------------------------------------------
328    ELSE:BEGIN
329      ras = report('Wrong definition of vargrid = '+vargrid+'. Only T, U, V, W or F are acceptable')
330      stop
331    END
332  ENDCASE
333  IF testvar(var = key_performance) EQ 2 THEN $
334    print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux
335;
336;------------------------------------------------------------
337;------------------------------------------------------------
338;------------------------------------------------------------
339; Variables refering to the vertical dimension
340;------------------------------------------------------------
341;------------------------------------------------------------
342;------------------------------------------------------------
343;
344;
345  tempdeux = systime(1)         ; For key_performance =2
346  if keyword_set(wdepth) then begin
347    gdep = gdepw[firstz:lastz]
348    e3 = e3w[firstz:lastz]
349  endif else begin
350    gdep = gdept[firstz:lastz]
351    e3 = e3t[firstz:lastz]
352  ENDELSE
353; for the vertical sections with partial steps
354  IF keyword_set(type) AND keyword_set(key_partialstep) THEN BEGIN
355    CASE 1 OF
356      type EQ 'xz' AND ny EQ 1:BEGIN
357        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
358        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
359        bottom = lindgen(nx)+(bottom-1l+keyword_set(wdepth))*nx
360        IF good[0] NE -1 THEN BEGIN
361          bottom = bottom[good]
362          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
363          gdep = replicate(1, nx)#gdep
364          if keyword_set(wdepth) THEN $
365            truegdep = hdepw[firstx:lastx, firsty:lasty] $
366          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
367          gdep[bottom] = truegdep[good]
368        ENDIF
369      END
370      type EQ 'yz' AND nx EQ 1:BEGIN
371        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
372        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
373        bottom = lindgen(ny)+(bottom-1l+keyword_set(wdepth))*ny
374        IF good[0] NE -1 THEN BEGIN
375          bottom = bottom[good]
376          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
377          gdep = replicate(1, ny)#gdep
378          if keyword_set(wdepth) THEN $
379            truegdep = hdepw[firstx:lastx, firsty:lasty] $
380          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
381          gdep[bottom] = truegdep[good]
382        ENDIF
383      END
384      ELSE:
385    ENDCASE
386  ENDIF
387; for the vertical sections with roms
388  IF keyword_set(type) AND n_elements(romszinfos) NE 0 THEN BEGIN
389    romsdp = romsdepth()
390    IF romsdp[0] NE -1 THEN BEGIN
391      IF jpt EQ 1 THEN BEGIN
392        CASE type OF
393          'z':gdep = moyenne(temporary(romsdp), 'xy')
394          'xz':gdep = moyenne(temporary(romsdp), 'y')
395          'yz':gdep = moyenne(temporary(romsdp), 'x')
396          'zt':gdep = moyenne(temporary(romsdp), 'xy')
397          ELSE:
398        ENDCASE
399      ENDIF ELSE BEGIN
400        CASE type OF
401          'z':gdep = moyenne(temporary(romsdp), 'xyt')
402          'xz':gdep = grossemoyenne(temporary(romsdp), 'yt')
403          'yz':gdep = grossemoyenne(temporary(romsdp), 'xt')
404          'zt':gdep = grossemoyenne(temporary(romsdp), 'xy')
405          ELSE:
406        ENDCASE
407      ENDELSE
408      IF size(gdep, /n_dimensions) EQ 2 THEN $
409         gdep = remplit(gdep, niter = 32000, mask = gdep LT 1.E19, /basique, /fillxdir)
410    ENDIF
411  ENDIF
412;
413  IF testvar(var = key_performance) EQ 2 THEN $
414    print, 'temps grille: Variables se rapportant a la dimension verticale ', systime(1)-tempdeux
415;------------------------------------------------------------
416; Triangulation vector when TRI is activated.
417;------------------------------------------------------------
418  if arg_present(TRI) then $
419    if triangles_list[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN
420    tempdeux = systime(1)       ; pour key_performance =2
421    msk = bytarr(jpi, jpj)
422    msk[firstx:lastx, firsty:lasty] = 1
423    ind = where( msk[triangles_list[0, *]]*msk[triangles_list[1, *]]*msk[triangles_list[2, *]] EQ 1 )
424    tri = triangles_list[*, ind]-(firstx+firsty*jpi)
425    y = tri/jpi
426    x = tri-y*jpi
427    tri = x+y*nx
428    IF testvar(var = key_performance) EQ 2 THEN $
429      print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux
430  ENDELSE
431;------------------------------------------------------------------
432; To make sure there is not any degenerated dimension (=1)
433;------------------------------------------------------------------
434;    mask=reform(mask, /over)
435;    glam=reform(glam, /over)
436;    gphi=reform(gphi, /over)
437;    gdep=reform(gdep, /over)
438;    e1=reform(e1, /over)
439;    e2=reform(e2, /over)
440;    e3=reform(e3, /over)
441
442  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grille.dat'
443  if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun
444
445;------------------------------------------------------------
446  IF NOT keyword_set(key_forgetold) THEN BEGIN
447@updateold
448  ENDIF
449;---------------------
450  return
451
452end
Note: See TracBrowser for help on using the repository browser.