source: trunk/SRC/ToBeReviewed/GRILLE/decoupeterre.pro @ 76

Last change on this file since 76 was 14, checked in by pinsard, 18 years ago

upgrade of COULEURS according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : directories
  • Property svn:executable set to *
File size: 10.1 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME:decoupeterre
6;
7; PURPOSE:tres semblable a grille. Ici qd vargrid ne 'T' ou 'W' alors
8; pour le trace il faut recuperer Tmask, glamt, gphit et le tableau de
9; triangulation sur le sous domaine considere. La specificite de
10; decoupeterre par rapport a grille, c''est que l''on prend ds la
11; mesure du possible un sous domaine juste un peu plus grand que celui
12; definit par domdef de facon a etre sur que le masque que l''on trace
13; recouvrira bien tout le dessin.
14;
15; CATEGORY:pour plt
16;
17; CALLING SEQUENCE:decoupeterre, mask, glam, gphi, z, nx, ny, nz, TRI = tri
18;
19; INPUTS:
20;
21; KEYWORD PARAMETERS:
22;         TRI si ce mot clef sert a obtenir grace a grille la
23;         triangulation qui se rapporte a la grille mais uniquement
24;         sur la partie du zoom. ce tableau de triangulation reduit
25;         est passe ds la variable que l''on a egalee a tri.par ex:
26;         grille,...,tri=triangulation_reduite. ne mot clef est
27;         utilise dans plt.pro
28;
29;         /WDEPTH: to specify that the field is at W depth instad of T
30;         depth (automatically activated if vargrid eq 'W')
31;
32;
33; OUTPUTS:le masque et ses coordonnees
34;
35; COMMON BLOCKS:
36;       common.pro
37;
38; SIDE EFFECTS:
39;
40; RESTRICTIONS:
41;
42; EXAMPLE:
43;
44; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
45;                       24/2/99
46;-
47;------------------------------------------------------------
48;------------------------------------------------------------
49;------------------------------------------------------------
50PRO decoupeterre, mask, glam, gphi, gdep, TYPE = type, TRI = tri, INDICEZOOM = indicezoom, COINMONTE = coinmonte, COINDESCEND = coindescend, WDEPTH = wdepth, REALSECTION = realsection, USETRI = usetri, _extra = ex
51;---------------------------------------------------------
52@cm_4mesh
53@cm_4data
54  IF NOT keyword_set(key_forgetold) THEN BEGIN
55@updatenew
56  ENDIF
57;---------------------------------------------------------
58  tempsun = systime(1)          ; pour key_performance
59;------------------------------------------------------------
60  if vargrid EQ 'W' then wdepth = 1
61;------------------------------------------------------------
62;------------------------------------------------------------
63; horizontal parameters
64;------------------------------------------------------------
65; if possible extent the domain according to the grid type
66; default case
67  case vargrid of
68    'U':BEGIN
69      firstx = 0 > (min([firstxt, firstxu])-1)
70      lastx = (max([lastxt, lastxu])+1) < (jpi-1)
71      firsty = 0 > (min([firstyt, firstyu])-1)
72      lasty = (max([lastyt, lastyu])+1) < (jpj-1)
73    end
74    'V':BEGIN
75      firstx = 0 > (min([firstxt, firstxv])-1)
76      lastx = (max([lastxt, lastxv])+1) < (jpi-1)
77      firsty = 0 > (min([firstyt, firstyv])-1)
78      lasty = (max([lastyt, lastyv])+1) < (jpj-1)
79    end
80    'F':BEGIN
81      firstx = 0 > (min([firstxt, firstxf])-1)
82      lastx = (max([lastxt, lastxf])+1) < (jpi-1)
83      firsty = 0 > (min([firstyt, firstyf])-1)
84      lasty = (max([lastyt, lastyf])+1) < (jpj-1)
85    END
86    ELSE:BEGIN
87      firstx = firstxt
88      lastx = lastxt
89      firsty = firstyt
90      lasty = lastyt
91    END
92  ENDCASE
93; however for the vertical section we don''t whant to extent the domain
94; in the direction perpendicular to the vertical section
95  if type EQ 'xz' then begin
96    case vargrid of
97      'U':BEGIN
98        firsty = firstyu
99        lasty = lastyu
100      END
101      'V':BEGIN
102        firsty = firstyv
103        lasty = lastyv
104      END
105      'F':BEGIN
106        firsty = firstyf
107        lasty = lastyf
108      END
109      ELSE:
110    ENDCASE
111  endif
112;
113  if type EQ 'yz' then begin
114    case vargrid of
115      'U':BEGIN
116        firstx = firstxu
117        lastx = lastxu
118      END
119      'V':BEGIN
120        firstx = firstxv
121        lastx = lastxv
122      END
123      'F':BEGIN
124        firstx = firstxf
125        lastx = lastxf
126      END
127      ELSE:
128    ENDCASE
129  endif
130;
131  nx = lastx-firstx+1
132  ny = lasty-firsty+1
133;------------------------------------------------------------
134; horizontal grid
135;------------------------------------------------------------
136; default case
137  glam = glamt[firstx:lastx, firsty:lasty]
138  gphi = gphit[firstx:lastx, firsty:lasty]
139; for the vertical section, 2 cases:
140; 1) keyword_set(realsection) eq 1: then we use drawsectionbottom.pro
141; that directly draw the bottom using polyfill and plot. In this case,
142; the position of the mask must given by the edge of the mask
143; cells. As we are considering the edges of the mask cells, glam or
144; gphi as one more point than the mask (if it is possible)
145; 2) keyword_set(realsection) eq 0: then we draw the mask by contouring
146; the mask with contour (contouring the level=.5). In this case, the
147; mask position must be given by the middle on the mask cells.
148  CASE  type OF
149    'xz':BEGIN
150      if keyword_set(realsection) EQ 0 then begin
151        if vargrid EQ 'V' OR vargrid EQ 'F' then $
152          glam = glamv[firstx:lastx, firsty:lasty]
153      ENDIF ELSE BEGIN          ; to drawsectionbottom
154        if vargrid EQ 'V' OR vargrid EQ 'F' OR finite(glamu[0]) EQ 0 then $
155          glam = glamf[0 > (firstx-1):lastx, firsty:lasty] $
156        ELSE glam = glamu[0 > (firstx-1):lastx, firsty:lasty]
157      ENDELSE
158    END
159    'yz':BEGIN
160      if keyword_set(realsection) EQ 0 then begin
161        if vargrid EQ 'U' OR vargrid EQ 'F' then $
162          gphi = gphiu[firstx:lastx, firsty:lasty]
163      ENDIF ELSE BEGIN          ; to drawsectionbottom 
164        if vargrid EQ 'U' OR vargrid EQ 'F' OR finite(gphiv[0]) EQ 0 then $
165          gphi = gphif[firstx:lastx, 0 > (firsty-1):lasty] $
166        ELSE gphi = gphiv[firstx:lastx, 0 > (firsty-1):lasty]
167      ENDELSE
168    END
169    ELSE:
170  ENDCASE
171;------------------------------------------------------------
172; vertical boundaries
173;------------------------------------------------------------
174  if keyword_set(wdepth)  then begin
175    firstz = 0 > (min([firstzt, firstzw])-1)
176    lastz = (max([lastzt, lastzw])+1) < (jpk-1)
177  ENDIF ELSE BEGIN
178    firstz = firstzt
179    lastz = lastzt
180  ENDELSE
181  nz = lastz-firstz+1
182;------------------------------------------------------------
183; mask
184;------------------------------------------------------------
185  case type of
186    'xy':BEGIN
187      mask = tmask[firstx:lastx, firsty:lasty, firstz]
188      profond = firstz NE 0
189    END
190; for the verical section, we have to choose the right mask according
191; to the grid point and to the direction of the section
192    'xz':BEGIN
193      if vargrid EQ 'V' OR vargrid EQ 'F' then begin
194        mask = (vmask())[firstx:lastx, firstyv:lastyv, firstz:lastz]
195      ENDIF ELSE mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
196    END
197    'yz':BEGIN
198      if vargrid EQ 'U' OR vargrid EQ 'F' then begin
199        mask = (umask())[firstxu:lastxu, firsty:lasty, firstz:lastz]
200      ENDIF ELSE mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
201    END
202    ELSE:mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
203  endcase
204;------------------------------------------------------------
205; vertical axis
206;------------------------------------------------------------
207; when we do a real section we directly plot the gdepw
208; (in drawsectionbottom.pro) instead of contouring the mask at 0.5 at
209; gdept
210  IF keyword_set(realsection) EQ 0 then gdep = gdept[firstz:lastz] $
211  ELSE BEGIN
212;     if lastz EQ jpk-1 then $
213; we add some fictive very deep level that will not be used but that is
214; necessary to avoid array size bugs in draw bottom section
215;      gdep = [gdepw[firstz+1:lastz], 2*gdept[jpk-1]] $
216;    ELSE gdep = gdepw[firstz+1:lastz+1]
217     gdep = gdepw[firstz:lastz]
218; special case when we are using the partial steps in the vertical
219; section that are only 1 point wide.
220; in that case, the z axis is a 2d array and we modify the depth of
221; the last level ocean with hdepw that is the real depth of the bottom.
222    CASE 1 OF
223      keyword_set(key_partialstep) and type EQ 'xz' $
224        AND ny EQ 1 AND keyword_set(realsection):BEGIN
225        bottom = total(mask, 3)
226        good = where(bottom NE 0 AND bottom NE nz+1)
227        bottom = lindgen(nx)+(bottom)*nx
228        IF good[0] NE -1 THEN BEGIN
229          bottom = bottom[good]
230          gdep = replicate(1, nx)#gdep
231          truegdep = hdepw[firstx:lastx, firsty:lasty]
232          gdep[bottom] = truegdep[good]
233        ENDIF
234      END
235      keyword_set(key_partialstep) and type EQ 'yz' $
236        AND nx EQ 1 AND keyword_set(realsection):BEGIN
237        bottom = total(mask, 3)
238        good = where(bottom NE 0 AND bottom NE nz+1)
239        bottom = lindgen(ny)+(bottom)*ny
240        IF good[0] NE -1 THEN BEGIN
241          bottom = bottom[good]
242          gdep = replicate(1, ny)#gdep
243          truegdep = hdepw[firstx:lastx, firsty:lasty]
244          gdep[bottom] = truegdep[good]
245        ENDIF
246      END
247      ELSE:
248    ENDCASE
249  ENDELSE
250;------------------------------------------------------------
251; vecteur triangulation Qd TRI est active
252;------------------------------------------------------------
253  IF arg_present(TRI) then $
254    if triangles_list[0] EQ -1 OR usetri LT 1 then tri = -1 ELSE BEGIN
255; si on est en train de tracer un niveau profond on refait la
256; triangulation
257    if keyword_set(profond) then begin
258      tri = triangule(mask, coinmonte = coinmonte, coindescend = coindescend, _extra = ex)
259      indicezoom = (lindgen(jpi, jpj))[firstx:lastx, firsty:lasty]
260    ENDIF ELSE BEGIN
261; sinon on recupere la partie de triangulation qui nous interesse et
262; on la numerote convenablement!
263      if nx EQ jpi AND ny EQ jpj then tri = triangles_list ELSE BEGIN
264        msk = bytarr(jpi, jpj)
265        msk[firstx:lastx, firsty:lasty] = 1
266        ind = where( msk[triangles_list[0, *]] EQ 1 $
267                     AND msk[triangles_list[1, *]] EQ 1 $
268                     AND msk[triangles_list[2, *]] EQ 1 )
269        tri = triangles_list[*, ind]-(firstx+firsty*jpi)
270        y = tri/jpi
271        x = tri-y*jpi
272        tri = x+y*nx
273      ENDELSE
274    ENDELSE
275  ENDELSE
276;-------------------------------------------------------------------
277  if keyword_set(key_performance) THEN print, 'temps decoupeterre', systime(1)-tempsun
278;------------------------------------------------------------
279  return
280end
281
Note: See TracBrowser for help on using the repository browser.