;----------------------------------------------------------------------
; This script reads corner edges off a SCRIP file and
; plots them on a map using three different methods, 
; simply to illustrate how different in speed the three
; methods are.
;
; Method 1 uses gsn_polyline in a loop, which can be a 
;   little faster than method 1 because the lines are
;   not taking up memory.
;
; Method 2 - uses gsn_add_polyline in a loop, which can
;   be the slowest method because each "lnid" line 
;   segment takes up memory.
;
; Method 3 - the fastest - uses the special "gsSegments"
;   resource to speed up line draws. No looping is required.
;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"

begin
  dir   = "/glade/p/cesm/cseg/mapping/grids/"
  fname = "ne16np4_110512_pentagons.nc"
  fgrid = addfile(dir+fname,"r")

  latctr = fgrid->grid_center_lat   ; center points, n
  lonctr = fgrid->grid_center_lon

  latcnr   = fgrid->grid_corner_lat
  loncnr   = fgrid->grid_corner_lon
  latcnr1d = ndtooned(latcnr)
  loncnr1d = ndtooned(loncnr)

  latcnr_dims = dimsizes(latcnr)
  ncnrs       = latcnr_dims(0)
  nedges      = latcnr_dims(1)

  wks = gsn_open_wks("png","seam_grid")

  res                       = True
  res@gsnDraw               = False
  res@gsnFrame              = False
  res@gsnMaximize           = True           ; Maximize size of plot in frame

  res@mpProjection          = "Orthographic"
  res@mpCenterLatF          =  40            ; Rotate map.
  res@mpCenterLonF          = 0
  res@mpGridAndLimbOn       = True           ; turn on lat/lon lines
  res@mpPerimOn             = False          ; turn off perimeter
  res@mpFillOn              = False
  res@mpOutlineOn           = True
; res@mpOutlineBoundarySets = "NoBoundaries"
  res@mpLimbLineColor       = "blue4"

  plot1 = gsn_csm_map(wks,res)    ; For method 1
  plot2 = gsn_csm_map(wks,res)    ; For method 2
  plot3 = gsn_csm_map(wks,res)    ; For method 3

;---Add the center locations to each plot
  mkres               = True
  mkres@gsMarkerIndex = 17              ; filled dot
  mkres@gsMarkerColor = "ForestGreen"
  mkid1 = gsn_add_polymarker(wks,plot1,lonctr,latctr,mkres)
  mkid2 = gsn_add_polymarker(wks,plot2,lonctr,latctr,mkres)
  mkid3 = gsn_add_polymarker(wks,plot3,lonctr,latctr,mkres)
  
;---Resources for grid lines
  lnres                     = True
  lnres@gsLineThicknessF    = 1.0             ; default is 1
  lnres@gsLineColor         = "black"           ; default is black.

;---Method 1
  start_time = get_cpu_time()
  draw(plot1)
  do n=0,ncnrs-1
    gsn_polyline(wks,plot1,loncnr(n,:),latcnr(n,:),lnres)
  end do
  frame(wks)
  end_time = get_cpu_time()

  print("Elapsed time for method 1 using gsn_polyline = " + \
        (end_time - start_time))

;---Method 2
  start_time = get_cpu_time()
  lnid1 = new(ncnrs,graphic)
  do n=0,ncnrs-1
    lnid1(n) = gsn_add_polyline(wks,plot2,loncnr(n,:),latcnr(n,:),lnres)
  end do
  draw(plot2)
  frame(wks)
  end_time = get_cpu_time()

  print("Elapsed time for method 2 using gsn_add_polyline = " + \
        (end_time - start_time))

;---Method 3
  segment_indices  = ispan(0,ncnrs*nedges-1,nedges) ; pen-up every 5th index
  lnres@gsSegments = segment_indices
  start_time = get_cpu_time()
  lnid3 = gsn_add_polyline(wks,plot3,loncnr1d,latcnr1d,lnres)
  draw(plot3)
  frame(wks)
  end_time = get_cpu_time()

  print("Elapsed time for method 3 using gsn_add_polyline/gsSegments = " + \
        (end_time - start_time))


end