;+
;
; @file_comments
; Similar to grille.
; Here, when vargrid is not 'T' or 'W', we have to
; recuperate tmask, glamt, gphit and the array of triangulation on the
; considered sub-domain for the drawing.
; The specificity of decoupeterre, in comparaison with grille, is
; that we take, if possible, a sub-domain just a little bit bigger than the
; one defined by domdef in order to be
; sure that the mask we draw will cover over all the drawing.
;
; @categories
; Grid
;
; @param MASK
;
; @param GLAM
;
; @param GPHI
;
; @param GDEP
;
; @keyword TYPE
;
; @keyword INDICEZOOM
;
; @keyword COINMONTE
;
; @keyword COINDESCEND
;
; @keyword REALSECTION
;
; @keyword USETRI
;
; @keyword _EXTRA
; Used to pass keywords
;
; @keyword TRI
; This keyword serve to obtain, thanks to grille, the triangulation which
; refer to the grid but only on the part of the zoom. This array of triangulation
; is passed in the variable we have equate at TRI.
; For example: grille,...,tri=triangulation_reduite.
; This keyword is used in plt.
;
; @keyword WDEPTH
; To specify that the field is at W depth instead of T
; depth (automatically activated if vargrid eq 'W')
;
; @uses
; common.pro
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; 24/2/99
;
; @version
; $Id$
;
; @todo
; seb : manque tous les param et plein de keywords.
;
;-
PRO decoupeterre, mask, glam, gphi, gdep $
, TYPE=type, TRI=tri, INDICEZOOM=indicezoom $
, COINMONTE=coinmonte, COINDESCEND=coindescend $
, WDEPTH=wdepth, REALSECTION=realsection, USETRI=usetri $
, _EXTRA=ex
;
compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
ENDIF
;---------------------------------------------------------
tempsun = systime(1) ; For key_performance
;------------------------------------------------------------
if vargrid EQ 'W' then wdepth = 1
;------------------------------------------------------------
;------------------------------------------------------------
; horizontal parameters
;------------------------------------------------------------
; if possible extent the domain according to the grid type
; default case
case vargrid of
'U':BEGIN
firstx = 0 > (min([firstxt, firstxu])-1)
lastx = (max([lastxt, lastxu])+1) < (jpi-1)
firsty = 0 > (min([firstyt, firstyu])-1)
lasty = (max([lastyt, lastyu])+1) < (jpj-1)
end
'V':BEGIN
firstx = 0 > (min([firstxt, firstxv])-1)
lastx = (max([lastxt, lastxv])+1) < (jpi-1)
firsty = 0 > (min([firstyt, firstyv])-1)
lasty = (max([lastyt, lastyv])+1) < (jpj-1)
end
'F':BEGIN
firstx = 0 > (min([firstxt, firstxf])-1)
lastx = (max([lastxt, lastxf])+1) < (jpi-1)
firsty = 0 > (min([firstyt, firstyf])-1)
lasty = (max([lastyt, lastyf])+1) < (jpj-1)
END
ELSE:BEGIN
firstx = firstxt
lastx = lastxt
firsty = firstyt
lasty = lastyt
END
ENDCASE
; however for the vertical section we don''t want to extent the domain
; in the direction perpendicular to the vertical section
if type EQ 'xz' then begin
case vargrid of
'U':BEGIN
firsty = firstyu
lasty = lastyu
END
'V':BEGIN
firsty = firstyv
lasty = lastyv
END
'F':BEGIN
firsty = firstyf
lasty = lastyf
END
ELSE:
ENDCASE
endif
;
if type EQ 'yz' then begin
case vargrid of
'U':BEGIN
firstx = firstxu
lastx = lastxu
END
'V':BEGIN
firstx = firstxv
lastx = lastxv
END
'F':BEGIN
firstx = firstxf
lastx = lastxf
END
ELSE:
ENDCASE
endif
;
nx = lastx-firstx+1
ny = lasty-firsty+1
;------------------------------------------------------------
; horizontal grid
;------------------------------------------------------------
; default case
glam = glamt[firstx:lastx, firsty:lasty]
gphi = gphit[firstx:lastx, firsty:lasty]
; for the vertical section, 2 cases:
; 1) keyword_set(realsection) eq 1: then we use drawsectionbottom.pro
; that directly draw the bottom using polyfill and plot. In this case,
; the position of the mask must given by the edge of the mask
; cells. As we are considering the edges of the mask cells, glam or
; gphi as one more point than the mask (if it is possible)
; 2) keyword_set(realsection) eq 0: then we draw the mask by contouring
; the mask with contour (contouring the level=.5). In this case, the
; mask position must be given by the middle on the mask cells.
CASE type OF
'xz':BEGIN
if keyword_set(realsection) EQ 0 then begin
if vargrid EQ 'V' OR vargrid EQ 'F' then $
glam = glamv[firstx:lastx, firsty:lasty]
ENDIF ELSE BEGIN ; to drawsectionbottom
if vargrid EQ 'V' OR vargrid EQ 'F' OR finite(glamu[0]) EQ 0 then $
glam = glamf[0 > (firstx-1):lastx, firsty:lasty] $
ELSE glam = glamu[0 > (firstx-1):lastx, firsty:lasty]
ENDELSE
END
'yz':BEGIN
if keyword_set(realsection) EQ 0 then begin
if vargrid EQ 'U' OR vargrid EQ 'F' then $
gphi = gphiu[firstx:lastx, firsty:lasty]
ENDIF ELSE BEGIN ; to drawsectionbottom
if vargrid EQ 'U' OR vargrid EQ 'F' OR finite(gphiv[0]) EQ 0 then $
gphi = gphif[firstx:lastx, 0 > (firsty-1):lasty] $
ELSE gphi = gphiv[firstx:lastx, 0 > (firsty-1):lasty]
ENDELSE
END
ELSE:
ENDCASE
;------------------------------------------------------------
; vertical boundaries
;------------------------------------------------------------
if keyword_set(wdepth) then begin
firstz = 0 > (min([firstzt, firstzw])-1)
lastz = (max([lastzt, lastzw])+1) < (jpk-1)
ENDIF ELSE BEGIN
firstz = firstzt
lastz = lastzt
ENDELSE
nz = lastz-firstz+1
;------------------------------------------------------------
; mask
;------------------------------------------------------------
case type of
'xy':BEGIN
mask = tmask[firstx:lastx, firsty:lasty, firstz]
profond = firstz NE 0
END
; for the vertical section, we have to choose the right mask according
; to the grid point and to the direction of the section
'xz':BEGIN
if vargrid EQ 'V' OR vargrid EQ 'F' then begin
mask = (vmask())[firstx:lastx, firstyv:lastyv, firstz:lastz]
ENDIF ELSE mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
END
'yz':BEGIN
if vargrid EQ 'U' OR vargrid EQ 'F' then begin
mask = (umask())[firstxu:lastxu, firsty:lasty, firstz:lastz]
ENDIF ELSE mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
END
ELSE:mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
endcase
;------------------------------------------------------------
; vertical axis
;------------------------------------------------------------
; when we do a real section we directly plot the gdepw
; (in drawsectionbottom.pro) instead of contouring the mask at 0.5 at
; gdept
IF keyword_set(realsection) EQ 0 then gdep = gdept[firstz:lastz] $
ELSE BEGIN
; if lastz EQ jpk-1 then $
; we add some fictive very deep level that will not be used but that is
; necessary to avoid array size bugs in draw bottom section
; gdep = [gdepw[firstz+1:lastz], 2*gdept[jpk-1]] $
; ELSE gdep = gdepw[firstz+1:lastz+1]
gdep = gdepw[firstz:lastz]
; special case when we are using the partial steps in the vertical
; section that are only 1 point wide.
; in that case, the z axis is a 2d array and we modify the depth of
; the last level ocean with hdepw that is the real depth of the bottom.
CASE 1 OF
keyword_set(key_partialstep) and type EQ 'xz' $
AND ny EQ 1 AND keyword_set(realsection):BEGIN
bottom = total(mask, 3)
good = where(bottom NE 0 AND bottom NE nz+1)
bottom = lindgen(nx)+(bottom)*nx
IF good[0] NE -1 THEN BEGIN
bottom = bottom[good]
gdep = replicate(1, nx)#gdep
truegdep = hdepw[firstx:lastx, firsty:lasty]
gdep[bottom] = truegdep[good]
ENDIF
END
keyword_set(key_partialstep) and type EQ 'yz' $
AND nx EQ 1 AND keyword_set(realsection):BEGIN
bottom = total(mask, 3)
good = where(bottom NE 0 AND bottom NE nz+1)
bottom = lindgen(ny)+(bottom)*ny
IF good[0] NE -1 THEN BEGIN
bottom = bottom[good]
gdep = replicate(1, ny)#gdep
truegdep = hdepw[firstx:lastx, firsty:lasty]
gdep[bottom] = truegdep[good]
ENDIF
END
ELSE:
ENDCASE
ENDELSE
;------------------------------------------------------------
; Triangulation vector when TRI is activated.
;------------------------------------------------------------
IF arg_present(TRI) then $
if triangles_list[0] EQ -1 OR usetri LT 1 then tri = -1 ELSE BEGIN
; If we are tracing a deep level, we redo the triangulation
if keyword_set(profond) then begin
tri = triangule(mask, coinmonte = coinmonte, coindescend = coindescend, _extra = ex)
indicezoom = (lindgen(jpi, jpj))[firstx:lastx, firsty:lasty]
ENDIF ELSE BEGIN
; Otherwise, we recuperate the part of triangulation that interest us and we number them well!!
if nx EQ jpi AND ny EQ jpj then tri = triangles_list ELSE BEGIN
msk = bytarr(jpi, jpj)
msk[firstx:lastx, firsty:lasty] = 1
ind = where( msk[triangles_list[0, *]] EQ 1 $
AND msk[triangles_list[1, *]] EQ 1 $
AND msk[triangles_list[2, *]] EQ 1 )
tri = triangles_list[*, ind]-(firstx+firsty*jpi)
y = tri/jpi
x = tri-y*jpi
tri = x+y*nx
ENDELSE
ENDELSE
ENDELSE
;-------------------------------------------------------------------
if keyword_set(key_performance) THEN print, 'temps decoupeterre', systime(1)-tempsun
;------------------------------------------------------------
return
end