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

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

modification of some headers (+some corrections) to prepare usage of the new idldoc

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