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

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:decoupeterre
;
; PURPOSE:tres semblable a grille. Ici qd vargrid ne 'T' ou 'W' alors
; pour le trace il faut recuperer Tmask, glamt, gphit et le tableau de
; triangulation sur le sous domaine considere. La specificite de
; decoupeterre par rapport a grille, c''est que l''on prend ds la
; mesure du possible un sous domaine juste un peu plus grand que celui
; definit par domdef de facon a etre sur que le masque que l''on trace
; recouvrira bien tout le dessin.
;
; CATEGORY:pour plt
;
; CALLING SEQUENCE:decoupeterre, mask, glam, gphi, z, nx, ny, nz, TRI = tri
;
; INPUTS:
;
; KEYWORD PARAMETERS:
;         TRI si ce mot clef sert a obtenir grace a grille la
;         triangulation qui se rapporte a la grille mais uniquement
;         sur la partie du zoom. ce tableau de triangulation reduit
;         est passe ds la variable que l''on a egalee a tri.par ex:
;         grille,...,tri=triangulation_reduite. ne mot clef est
;         utilise dans plt.pro
;
;         /WDEPTH: to specify that the field is at W depth instad of T 
;         depth (automatically activated if vargrid eq 'W')
;
;
; OUTPUTS:le masque et ses coordonnees
;
; COMMON BLOCKS:
;       common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
;                       24/2/99
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
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)          ; pour 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 whant 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 verical 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
;------------------------------------------------------------
; vecteur triangulation Qd TRI est active
;------------------------------------------------------------
  IF arg_present(TRI) then $
    if triangles_list[0] EQ -1 OR usetri LT 1 then tri = -1 ELSE BEGIN
; si on est en train de tracer un niveau profond on refait la
; 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 
; sinon on recupere la partie de triangulation qui nous interesse et
; on la numerote convenablement!
      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