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

;+
;
; @file_comments
; Choose the grid which must be used to do the graph in function of
; vargrid and send back corresponding parameters calculated in
; <pro>domdef</pro> and reduced at the domain defined by
; <pro>domdef</pro> (contrarily to 
; <pro>grandegrille</pro>)
; BEWARE!! The choice of the grid is made from the value of the
; global variable vargrid, which can be equal to 'T', 'U', 'V', 'W' ou 'F'.
;
; @categories
; Grid
;
; @keyword TRI
; This keyword serve to obtain 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')
;
; @keyword FORPLT
; In <pro>plt</pro>, we want that land points,
; glam and gphi, be equal to glamt and
; gphit regardless of the grid.
;
; @keyword NOTRI
; Useful only when TRI is activated. In this case, grille send back -1 in the
; variable tri even if the variable of the common triangles_list is defined
; and different of -1
;
; @keyword _EXTRA
; Used to pass keywords
;
; @keyword TOUT
;
; @param MASK {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param GLAM {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param GPHI {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param GDEP {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param NX {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param NY {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param NZ {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param FIRSTX {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param FIRSTY {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param FIRSTZ {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param LASTX {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param LASTY {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param LASTZ {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param E1 {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param E2 {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @param E3 {out}{optional}
; For the definition, see <pro>domdef</pro> and the management of subdomains on the web.
;
; @uses
; <pro>cm_4mesh</pro>
; <pro>cm_4data</pro>
;
; @restrictions
; Use the variable vargrid
;
; @restrictions
; Vargrid must be 'T', 'W', 'U', 'V' or 'F'
;
; @history
;  Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                       12/2/1999
;                       10/11/1999 /forplt
;
; @version
; $Id$
;
; @todo
; Comment ecrire la remarque sur les inputs?
;
;-
PRO grille, mask, glam, gphi, gdep, nx, ny, nz $
          , firstx, firsty, firstz, lastx, lasty, lastz $
          , e1, e2, e3 $
          , TRI=tri, NOTRI=notri, TOUT=tout, FORPLT=forplt $
          , TYPE=type, WDEPTH=wdepth, _EXTRA=ex
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal ; for jpt
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
  ENDIF
;---------------------
  tempsun = systime(1)          ; For key_performance
;------------------------------------------------------------
  vargrid = strupcase(strmid(vargrid,0,/reverse_offset))
;
  if vargrid eq 'W' then wdepth = 1
  if keyword_set(tout) then begin
    savedbox = 1b
    saveboxparam, 'boxparam4grille.dat'
    domdef, gridtype = vargrid, _EXTRA = ex
  endif
  tempdeux = systime(1)         ; For key_performance =2
;------------------------------------------------------------
;------------------------------------------------------------
  IF keyword_set(wdepth) THEN BEGIN
    firstz = firstzw
    lastz = lastzw
    nz = nzw
  ENDIF ELSE BEGIN
    firstz = firstzt
    lastz = lastzt
    nz = nzt
  ENDELSE
;------------------------------------------------------------
;------------------------------------------------------------
  CASE 1 OF
;------------------------------------------------------------
; grid T and W
;------------------------------------------------------------
    vargrid eq 'T' OR vargrid eq 'W' : begin
;scalars
      nx = nxt
      ny = nyt
      firstx = firstxt
      firsty = firstyt
      lastx = lastxt
      lasty = lastyt
;2d vectors
      IF arg_present(glam) THEN glam = glamt[firstx:lastx, firsty:lasty]
      IF arg_present(gphi) THEN gphi = gphit[firstx:lastx, firsty:lasty]
      IF arg_present(e1) THEN e1 = e1t[firstx:lastx, firsty:lasty]
      IF arg_present(e2) THEN e2 = e2t[firstx:lastx, firsty:lasty]
;3d vectors
      IF keyword_set(forplt) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz] $
      ELSE IF arg_present(mask) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
    end
;------------------------------------------------------------
; grid U
;------------------------------------------------------------
    vargrid eq 'U': begin
;scalars
      nx = nxu
      ny = nyu
      firstx = firstxu
      firsty = firstyu
      lastx = lastxu
      lasty = lastyu
;2d vectors
      IF arg_present(glam) THEN glam = glamu[firstx:lastx, firsty:lasty]
      IF arg_present(gphi) THEN gphi = gphiu[firstx:lastx, firsty:lasty]
      if keyword_set(forplt) then BEGIN
        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
        eastboarder = mask-shift(mask, 1, 0)*mask
        westboarder = mask-shift(mask, -1, 0)*mask
        if key_periodic NE 1 OR nx NE jpi then westboarder[nx-1, *] = 0b
        tmp1 = shift(eastboarder, 0, 1)
        tmp1[*, 0] = 0b
        tmp2 = shift(eastboarder, 0, -1)
        tmp2[*, ny-1] = 0b
        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
        eastboarder = temporary(eastboarder)+temporary(add)
        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
        tmp1[*, ny-1] = 1b
        tmp1[*, 0] = 1b
        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
        if key_periodic NE 1 OR nx NE jpi then begin
          tmp2[nx-1, *] = 1b
          tmp2[0, *] = 0b
        endif
        no1 = temporary(tmp1)*temporary(tmp2)
        tmp = temporary(eastboarder)*temporary(no1)*mask
        mask[0:nx-2, *] = 0b
        tmp = temporary(tmp)+temporary(mask)
        tmp = where(tmp GE 1)
        if tmp[0] NE -1 then begin
          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
        endif
      ENDIF
      IF arg_present(e1) THEN e1  = e1u[firstx:lastx, firsty:lasty]
      IF arg_present(e2) THEN e2  = e2u[firstx:lastx, firsty:lasty]
;3d vectors
      IF keyword_set(forplt) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz] $
      ELSE IF arg_present(mask) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
    end
;------------------------------------------------------------
; grid V
;------------------------------------------------------------
    vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $
      or vargrid eq 'V': begin
;scalars
      nx = nxv
      ny = nyv
      firstx = firstxv
      firsty = firstyv
      lastx = lastxv
      lasty = lastyv
;2d vectors
      IF arg_present(glam) THEN glam = glamv[firstx:lastx, firsty:lasty]
      IF arg_present(gphi) THEN gphi = gphiv[firstx:lastx, firsty:lasty]
      if keyword_set(forplt) then BEGIN
        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
        northboarder = mask-shift(mask, 0, 1)*mask
        southboarder = mask-shift(mask, 0, -1)*mask
        southboarder[*, ny-1] = 0b
        tmp1 = shift(northboarder, -1, 0)
        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
        tmp2 = shift(northboarder, 1, 0)
        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
        northboarder = temporary(northboarder)+temporary(add)
        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
        tmp1[*, ny-1] = 1b
        tmp1[*, 0] = 0b
        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
        if key_periodic NE 1 OR nx NE jpi then begin
          tmp2[nx-1, *] = 1b
          tmp2[0, *] = 1b
        endif
        no1 = temporary(tmp1)*temporary(tmp2)
        tmp = temporary(northboarder)*mask*temporary(no1)
        mask[*, 0:ny-2] = 0b
        tmp = temporary(tmp)+temporary(mask)
        tmp = where(tmp GE 1)
        if tmp[0] NE -1 then begin
          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
        endif
      ENDIF
      IF arg_present(e1) THEN e1  = e1v[firstx:lastx, firsty:lasty]
      IF arg_present(e2) THEN e2  = e2v[firstx:lastx, firsty:lasty]
;3d vecteurs
      IF keyword_set(forplt) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz] $
      ELSE IF arg_present(mask) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
    end
;------------------------------------------------------------
; grid F
;------------------------------------------------------------
    vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $
      or vargrid eq 'F': begin
;scalars
      nx = nxf
      ny = nyf
      firstx = firstxf
      firsty = firstyf
      lastx = lastxf
      lasty = lastyf
;2d vectors
      IF arg_present(glam) THEN glam = glamf[firstx:lastx, firsty:lasty]
      IF arg_present(gphi) THEN gphi = gphif[firstx:lastx, firsty:lasty]
      if keyword_set(forplt) then BEGIN
        mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz]
        eastboarder = mask-shift(mask, 1, 0)*mask
        westboarder = mask-shift(mask, -1, 0)*mask
        westboarder[nx-1, *] = 0b
        northboarder = mask-shift(mask, 0, 1)*mask
        southboarder = mask-shift(mask, 0, -1)*mask
        southboarder[*, ny-1] = 0b
        tmp1 = shift(northboarder, -1, 0)
        if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b
        tmp2 = shift(northboarder, 1, 0)
        if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b
        add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder)
        northboarder = temporary(northboarder)+temporary(add)
        tmp1 = shift(eastboarder, 0, 1)
        tmp1[*, 0] = 0b
        tmp2 = shift(eastboarder, 0, -1)
        tmp2[*, ny-1] = 0b
        add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder))
        eastboarder = temporary(eastboarder)+temporary(add)
        tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b
        tmp1[*, ny-1] = 1b
        tmp1[*, 0] = 1b
        tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b
        if key_periodic NE 1 OR nx NE jpi then begin
          tmp2[nx-1, *] = 1b
          tmp2[0, *] = 1b
        endif
        no1 = temporary(tmp1)*temporary(tmp2)
        tmp = (temporary(northboarder)+temporary(eastboarder))*mask*temporary(no1)
        mask[0:nx-2, *] = 0b
        mask[*, 0:ny-2] = 0b
        tmp = temporary(tmp)+temporary(mask)
        tmp = where(tmp GE 1)
        if tmp[0] NE -1 then begin
          glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp]
          gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp]
        endif
      ENDIF
      IF arg_present(e1) THEN e1  = e1f[firstx:lastx, firsty:lasty]
      IF arg_present(e2) THEN e2  = e2f[firstx:lastx, firsty:lasty]
;3d vectors
      IF keyword_set(forplt) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz] $
      ELSE IF arg_present(mask) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
    END
;------------------------------------------------------------
    ELSE:BEGIN
      ras = report('Wrong definition of vargrid = '+vargrid+'. Only T, U, V, W or F are acceptable')
      stop
    END
  ENDCASE
  IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux
;
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
; Variables refering to the vertical dimension
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;
;
  tempdeux = systime(1)         ; For key_performance =2
  if keyword_set(wdepth) then begin
    gdep = gdepw[firstz:lastz]
    e3 = e3w[firstz:lastz]
  endif else begin
    gdep = gdept[firstz:lastz]
    e3 = e3t[firstz:lastz]
  ENDELSE
; for the vertical sections with partial steps
  IF keyword_set(type) AND keyword_set(key_partialstep) THEN BEGIN
    CASE 1 OF
      type EQ 'xz' AND ny EQ 1:BEGIN
        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
        bottom = lindgen(nx)+(bottom-1l+keyword_set(wdepth))*nx
        IF good[0] NE -1 THEN BEGIN
          bottom = bottom[good]
          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
          gdep = replicate(1, nx)#gdep
          if keyword_set(wdepth) THEN $
            truegdep = hdepw[firstx:lastx, firsty:lasty] $
          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
          gdep[bottom] = truegdep[good]
        ENDIF
      END
      type EQ 'yz' AND nx EQ 1:BEGIN
        bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3)
        good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth))
        bottom = lindgen(ny)+(bottom-1l+keyword_set(wdepth))*ny
        IF good[0] NE -1 THEN BEGIN
          bottom = bottom[good]
          IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw)
          gdep = replicate(1, ny)#gdep
          if keyword_set(wdepth) THEN $
            truegdep = hdepw[firstx:lastx, firsty:lasty] $
          ELSE truegdep = hdept[firstx:lastx, firsty:lasty]
          gdep[bottom] = truegdep[good]
        ENDIF
      END
      ELSE:
    ENDCASE
  ENDIF
; for the vertical sections with roms
  IF keyword_set(type) AND n_elements(romszinfos) NE 0 THEN BEGIN
    romsdp = romsdepth()
    IF romsdp[0] NE -1 THEN BEGIN
      IF jpt EQ 1 THEN BEGIN
        CASE type OF
          'z':gdep = moyenne(temporary(romsdp), 'xy')
          'xz':gdep = moyenne(temporary(romsdp), 'y')
          'yz':gdep = moyenne(temporary(romsdp), 'x')
          'zt':gdep = moyenne(temporary(romsdp), 'xy')
          ELSE:
        ENDCASE
      ENDIF ELSE BEGIN
        CASE type OF
          'z':gdep = moyenne(temporary(romsdp), 'xyt')
          'xz':gdep = grossemoyenne(temporary(romsdp), 'yt')
          'yz':gdep = grossemoyenne(temporary(romsdp), 'xt')
          'zt':gdep = grossemoyenne(temporary(romsdp), 'xy')
          ELSE:
        ENDCASE
      ENDELSE
      IF size(gdep, /n_dimensions) EQ 2 THEN $
         gdep = remplit(gdep, niter = 32000, mask = gdep LT 1.E19, /basique, /fillxdir)
    ENDIF
  ENDIF
;
  IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps grille: Variables se rapportant a la dimension verticale ', systime(1)-tempdeux
;------------------------------------------------------------
; Triangulation vector when TRI is activated.
;------------------------------------------------------------
  if arg_present(TRI) then $
    if triangles_list[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN
    tempdeux = systime(1)       ; pour key_performance =2
    msk = bytarr(jpi, jpj)
    msk[firstx:lastx, firsty:lasty] = 1
    ind = where( msk[triangles_list[0, *]]*msk[triangles_list[1, *]]*msk[triangles_list[2, *]] EQ 1 )
    tri = triangles_list[*, ind]-(firstx+firsty*jpi)
    y = tri/jpi
    x = tri-y*jpi
    tri = x+y*nx
    IF testvar(var = key_performance) EQ 2 THEN $
      print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux
  ENDELSE
;------------------------------------------------------------------
; To make sure there is not any degenerated dimension (=1)
;------------------------------------------------------------------
;    mask=reform(mask, /over)
;    glam=reform(glam, /over)
;    gphi=reform(gphi, /over)
;    gdep=reform(gdep, /over)
;    e1=reform(e1, /over)
;    e2=reform(e2, /over)
;    e3=reform(e3, /over)

  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grille.dat'
  if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun

;------------------------------------------------------------
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updateold
  ENDIF
;---------------------
  return

end