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

Last change on this file since 483 was 388, checked in by smasson, 16 years ago

introduce meridional and barotropic stream functions, see ticket:59

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