source: trunk/SRC/ToBeReviewed/TRIANGULATION/tracecote.pro @ 231

;+
;
; @file_comments
; Draw coasts in plt.
;
; @categories
; Graphics
;
; @keyword SURFACE_COASTLINE
; To draw the surface coast line instead of
; the coast line at level firstz[tw]. Useful only for deep
; plots!
;
; @keyword _EXTRA
; Used to pass keywords
;
; @uses
; common.pro
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                      30/9/1999
;
; @version
; $Id$
;
;-
;
PRO tracecote, SURFACE_COASTLINE = surface_coastline, _EXTRA = ex
;
  compile_opt idl2, strictarrsubs
;
@cm_4data
@cm_4mesh
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
  ENDIF
;--------------------------------------------------------------
   tempsun = systime(1)         ; For key_performance
   if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c'
;
; We enlarge a bit the frame defined by firsts..., lasts... in order to
; recuperate edges of the coast which are in the edging of the domain.
;
   tempdeux = systime(1)        ; For key_performance =2
   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)
   nx = lastx-firstx+1
   ny = lasty-firsty+1
; Which vertical level choose?
   IF keyword_set(surface_coastline) THEN firstz = 0 ELSE $
     IF strupcase(vargrid) eq 'W' THEN firstz = firstzw ELSE firstz = firstzt
; Attribution of the mask and of coordinates  delimiting limits of the land (coordinates f).
   mask = tmask[firstx:lastx, firsty:lasty, firstz]
   xf = glamf[firstx:lastx, firsty:lasty]
   yf = gphif[firstx:lastx, firsty:lasty] ;
   IF testvar(var = key_performance) EQ 2 THEN $
    print, 'temps tracecote: determiner mask xf yf', systime(1)-tempdeux
;
   if key_gridtype EQ 'e' then onemore = xf[0, 0] gT xf[0, 1]
; We pass in normalized coordinates to be able to become independant from the projection's
; type choosen and from the support on which we do the drawing (screen or postscript)
   z = convert_coord(xf[*],yf[*],/data,/to_normal)
   xf = reform(z[0, *], nx, ny)
   yf = reform(z[1, *], nx, ny)
   tempvar = SIZE(TEMPORARY(z))
;
; Beware, following the projection, some points x or y can become NaN (see point
; behind the earth in an orthographic projection).
;
; We put points to be eliminated at a very big value so that they will not pass the
; test with distanceseuil (see further).
;
   if (!map.projection LE 7 AND !map.projection NE 0) $
    OR !map.projection EQ 14 OR !map.projection EQ 15 OR !map.projection EQ 18 then begin
      ind = where(finite(xf*yf) EQ 0)
      IF ind[0] NE -1 THEN BEGIN
         xf[ind] = 1e5
         yf[ind] = 1e5
      ENDIF
   ENDIF
   ind = where(xf LT !p.position[0] OR xf GT !p.position[2])
   IF ind[0] NE -1 THEN xf[ind] = 1e5
   ind = where(yf LT !p.position[1] OR yf GT !p.position[3])
   IF ind[0] NE -1 THEN yf[ind] = 1e5
   tempvar = SIZE(TEMPORARY(ind)) ; we delete ind
;
   case strmid(key_gridtype, 0, 1) of
      'c':drawcoast_c, mask, xf, yf, nx, ny, _extra = ex
      'e':drawcoast_e, mask, xf, yf, nx, ny, onemore = onemore, _extra = ex
   endcase
   if keyword_set(key_performance) THEN print, 'temps tracecote', systime(1)-tempsun
   return
end