load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin
;
; This grid came to us via Dave Randall, Todd Ringler, and Ross Heikes
; of CSU. The data for this mesh were originally downloaded from:
; http://kiwi.atmos.colostate.edu/BUGS/geodesic/interpolate.html
;

;******************************************
; read in data
;******************************************
  r2d = 57.2957795             ; radians to degrees
  f = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/hswm_d000000p000.g2.nc","r")

  x  = f->grid_center_lon * r2d
  y  = f->grid_center_lat * r2d
  cx = f->grid_corner_lon * r2d
  cy = f->grid_corner_lat * r2d
  ke = f->kinetic_energy(2,:)

  wks = gsn_open_wks("ps","geo")
  gsn_define_colormap(wks,"gui_default")

  res                      = True            ; plot mods desired
;  res@gsnMaximize          = True           ; largest plot possible

;---This resource not needed in V6.1.0
  res@gsnSpreadColors      = True            ; use full range of colormap

  res@cnFillOn             = True            ; turn on color
  res@cnLinesOn            = False           ; turn off contour lines
  res@cnLineLabelsOn       = False           ; turn off line labels

  res@lbLabelAutoStride    = True            ; nice label bar stride

  res@tiMainString         = "Geodesic grid using triangular mesh conversion"

  res@mpProjection      = "Orthographic"     ; choose projection
  res@mpDataBaseVersion = "MediumRes"        ; change outline database
  res@mpCenterLatF      =  40                ; rotate globe
  res@mpCenterLonF      = -100
  res@mpGridAndLimbOn   = True               ; turn on lat/lon lines
  res@mpPerimOn         = False              ; turn off map perimeter

; these resources are what set the triangular mesh.
  res@sfXArray             = x
  res@sfYArray             = y
  res@sfXCellBounds        = cx
  res@sfYCellBounds        = cy

  map = gsn_csm_contour_map(wks,ke,res)
end