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

;+
;
; @file_comments
; Similar to <pro>grille</pro>. 
; 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 <pro>grille</pro>, is 
; that we take, if possible, a sub-domain just a little bit bigger than the 
; one defined by <pro>domdef</pro> 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 MASKDTA
; 2d array defining the mask that must be applied to the data instead
; of the default mask which is defined according to the grid (t/u/v/f mask)
;
; @keyword _EXTRA
; Used to pass keywords
;
; @keyword TRI
; This keyword serve to obtain, thanks to <pro>grille</pro>, 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 <pro>plt</pro>.
;
; @keyword WDEPTH
; To specify that the field is at W depth instead of T
; depth (automatically activated if vargrid eq 'W')
;
; @uses
; <pro>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, MASKDTA = maskdta $
                , 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
    firstzw = 0 > (min([firstzt, firstzw])-1)
    lastzw = (max([lastzt, lastzw])+1) < (jpk-1)
    nzw = lastzw - firstzw + 1
    firstz = firstzw
    lastz = lastzw
    nz = nzw
  ENDIF ELSE BEGIN
    firstz = firstzt
    lastz = lastzt
    nz = nzt
  ENDELSE
;------------------------------------------------------------
; mask
;------------------------------------------------------------
  IF keyword_set(maskdta) THEN BEGIN
    mask = byte(fitintobox(maskdta))
    IF keyword_set(wdepth) THEN BEGIN
      mask = shift(temporary(mask), 0, 1)
      mask[*, 0] = 0b
    ENDIF
  ENDIF ELSE BEGIN 
  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
ENDELSE
;------------------------------------------------------------
; 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.
; in domzgr:  hdepw(ji,jj) = gdepw(ji,jj,ik+1)
    CASE 1 OF
      keyword_set(key_partialstep) AND type EQ 'xz' $
        AND ny EQ 1 AND keyword_set(realsection):BEGIN
        bottom = total(mask, 3 - keyword_set(maskdta))
        good = where(bottom NE 0 AND bottom NE nz)
        IF good[0] NE -1 THEN BEGIN
          bottom = lindgen(nx)+(bottom)*nx
          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 - keyword_set(maskdta))
        good = where(bottom NE 0 AND bottom NE nz)
        IF good[0] NE -1 THEN BEGIN
          bottom = lindgen(ny)+(bottom)*ny
          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