source: trunk/ToBeReviewed/GRILLE/tracegrille.pro @ 13

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:tracegrille
;
; PURPOSE:dessine la grille
;
; CATEGORY:
;
; CALLING SEQUENCE:tracegrille
; 
; INPUTS:glam et gphi, les tableaux 1d ou 2d des position en
; longitude/latitude des points de la grille. Si glam et gphi ne sont
; pas specifies, trace la grille specifiee par vargrid, sur le domaine
; definit par le dernier domdef.
;
; KEYWORD PARAMETERS:
;
;        XSTRIDE un entier pour specifier qu''on ne veut tracer
;        qu''une ligne de i constant tout les xstride points
;
;        YSTRIDE un entier pour specifier qu''on ne veut tracer
;        qu''une ligne de j constant tout les ystride points
;
;        /OCEAN: pour ne tracer la grille que sur les points oceans
; 
;        /EARTH: pour ne tracer la grille que sur les points terre
;
;        /RMOUT:select to remove all cell having one corner out of the
;        plot boundaries (!x.range, !y.range)
;
;        + tous les mots clefs de la procedure PLOTS
;
; OUTPUTS:none
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; EXAMPLE:
;
;     IDL> plt,indgen(jpi,jpj),/nocontour,/nofill
;     IDL> vargrid='T'
;     IDL> tracegrille,/ocean,color=20
;     IDL> tracegrille,/earth,color=80
;
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
PRO tracegrille, glamin, gphiin, OCEAN = ocean, EARTH = earth $
                 , XSTRIDE = xstride, YSTRIDE = ystride, RMOUT = rmout $
                 , _extra = extra
;---------------------------------------------------------
@cm_4mesh
@cm_4data
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
  ENDIF
;---------------------------------------------------------
  tempsun = systime(1)          ; pour 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 key_gridtype 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 key_gridtype 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