source: trunk/SRC/ToBeReviewed/GRILLE/tracegrille.pro @ 325

;+
;
; @file_comments
; Draw the grid
;
; @categories
; Grid
;
; @param GLAMIN {in}{optional}{type=1d or 2d array}{default=glam specified by vargrid, on the domain defined by <pro>domdef</pro>}
; Longitude of points of the grid.
;
; @param GPHIIN {in}{optional}{type=1d or 2d array}{default=gphi specified by vargrid, on the domain defined by <pro>domdef</pro>}
; Latitude of points of the grid.
;
; @keyword XSTRIDE {type=integer}{default=1}
; It specify that we want to trace only one line of
; constant i every xstride points
;
; @keyword YSTRIDE {type=integer}{default=1}
; It specify that we want to trace only one line of
; constant j every xstride points
;
; @keyword OCEAN
; To trace the grid only on ocean points.
;
; @keyword EARTH
; To trace the grid only on land points.
;
; @keyword RMOUT
; Select to remove all cell having one corner out of the
; plot boundaries (!x.range, !y.range)
;
; @keyword _EXTRA
; Used to pass keywords
;
; @uses
; common.pro
;
; @examples
;
;     IDL> plt,indgen(jpi,jpj),/nocontour,/nofill
;     IDL> vargrid='T'
;     IDL> tracegrille,/ocean,color=20
;     IDL> tracegrille,/earth,color=80
;
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
PRO tracegrille, glamin, gphiin, OCEAN = ocean, EARTH = earth $
                 , XSTRIDE = xstride, YSTRIDE = ystride, RMOUT = rmout $
                 , _EXTRA = extra
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
  ENDIF
;---------------------------------------------------------
  tempsun = systime(1)          ; For key_performance
; to avoid warning message
  oldexcept = !except
  !except = 0
  if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c'
;
  if n_elements(glamin) * n_elements(gphiin) EQ 0 then BEGIN
    grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz
    IF keyword_set(ocean) AND strmid(key_gridtype, 0, 1) EQ 'c' THEN BEGIN
; we reduce the mask to take into account the point located ON the coastline.
      CASE vargrid OF
        'U':BEGIN
          mask = tmask[firstx:lastx, firsty:lasty]
          IF NOT keyword_set(key_periodic) OR nx NE jpi $
            THEN tmpx = mask[nx-1, *]
          mask = (mask+shift(mask, -1, 0)) < 1
          IF NOT keyword_set(key_periodic) OR nx NE jpi $
            THEN mask[nx-1, *] = temporary(tmpx)
        END
        'V':BEGIN
          mask = tmask[firstx:lastx, firsty:lasty]
          tmpy = mask[*, ny-1]
          mask = (mask+shift(mask, 0, -1)) < 1
          mask[*, ny-1] = temporary(tmpy)
        END
        'F':BEGIN
          mask = tmask[firstx:lastx, firsty:lasty]
          IF NOT keyword_set(key_periodic) OR nx NE jpi $
            THEN tmpx = mask[nx-1, *]
          tmpy = mask[*, ny-1]
          mask = (mask+shift(mask, -1, 0)+shift(mask, 0, -1)+shift(mask, -1, -1)) < 1
          mask[*, ny-1] = temporary(tmpy)
          IF NOT keyword_set(key_periodic) OR nx NE jpi $
            THEN mask[nx-1, *] = temporary(tmpx)
        END
        ELSE:
      ENDCASE
    ENDIF
  ENDIF ELSE BEGIN
    glam = glamin
    gphi = gphiin
    IF (size(glam))[0] EQ 1 AND (size(gphi))[0] EQ 1 THEN BEGIN
      nx = n_elements(glam)
      ny = n_elements(gphi)
      glam = glam#replicate(1, ny)
      gphi = replicate(1, nx)#gphi
    ENDIF ELSE BEGIN
      nx = (size(glam))[1]
      ny = (size(glam))[2]
    ENDELSE
  ENDELSE
  if n_elements(mask) EQ 0 then mask = replicate(1b, nx, ny)
  if (size(mask))[0] EQ 3 then mask = mask[*, *, 0]
;
  IF keyword_set(RMOUT) THEN BEGIN
    out = where(glam GT max(!x.range) OR glam LT min(!x.range) $
                OR gphi GT max(!y.range) OR gphi LT min(!y.range))
    IF out[0] NE -1 THEN BEGIN
      glam[out] = !values.f_nan
      gphi[out] = !values.f_nan
    ENDIF
  ENDIF
;
  IF keyword_set(ocean) then BEGIN
    earth = where(mask EQ 0)
    if earth[0] NE -1 then begin
      glam[earth] = !values.f_nan
      gphi[earth] = !values.f_nan
    ENDIF
    earth = 0
  ENDIF
;
  IF keyword_set(earth) THEN BEGIN
    ocean = where(mask EQ 1)
    if ocean[0] NE -1 then begin
      glam[ocean] = !values.f_nan
      gphi[ocean] = !values.f_nan
    ENDIF
    ocean = 0
  ENDIF
;
  if NOT keyword_set(xstride) then xstride = 1
  if NOT keyword_set(ystride) then ystride = 1
  case strmid(key_gridtype, 0, 1) of
    'c':BEGIN
      for i = 0, ny-1, ystride do begin
        plots,  glam[*, i], gphi[*, i], _extra = extra
      endfor
      for i = 0, nx-1, xstride do begin
        plots,  glam[i, *], gphi[i, *], _extra = extra
      endfor
    END
    'e':BEGIN
      shifted = glam[0, 0] LT glam[0, 1]
      glam2 = glam+(glam[1]-glam[0])/2.
      if shifted then begin
        for i = 0, ny-2 do BEGIN
          xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*]
          yy = (transpose([[gphi[*, i]], [gphi[*, i+1]]]))[*]
          plots, xx[0:2*nx-2], yy[0:2*nx-2], _extra = extra
        ENDFOR
      ENDIF ELSE BEGIN
        for i = 1, ny-1 do BEGIN
          xx = (transpose([[glam[*, i]], [glam2[*, i]]]))[*]
          yy = (transpose([[gphi[*, i]], [gphi[*, i-1]]]))[*]
          plots, xx[0:2*nx-2], yy[0:2*nx-2], _extra = extra
        ENDFOR
      ENDELSE
      for i = 1, (ny-1)/2 do $
        plots, [glam[0, 2*i-1], glam[0, 2*i]] $
        , [gphi[0, 2*i-1], gphi[0, 2*i]], _extra = extra
      for i = 0, (ny-2)/2 do $
        plots, [glam[nx-1, 2*i], glam[nx-1, 2*i+1]] $
        , [gphi[nx-1, 2*i], gphi[nx-1, 2*i+1]], _extra = extra
    END
  endcase

  if keyword_set(key_performance) THEN print, 'temps trace grille', systime(1)-tempsun
  !except = oldexcept

  return
end