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

begin
  f = addfile("ctcbay.nc","r")

  wks = gsn_open_wks("png","chspkbay-cell")
  gsn_define_colormap(wks,"rainbow+gray")

  ele = f->ele
  lon = f->lon
  lat = f->lat
  depth = f->depth
  ;printVarSummary(ele)
  ;printVarSummary(lon)

  ;
  ; create arrays with the lat, lon, and data values  of the points defining the triangles

  xb = new(dimsizes(ele),typeof(lon))
  yb = new(dimsizes(ele),typeof(lat))
  data = new(dimsizes(ele),typeof(depth))

  do i = 0,2
    xb(:,i) = lon(ele(:,i)-1)
    yb(:,i) = lat(ele(:,i)-1)
    data(:,i) = depth(ele(:,i)-1)
  end do 

  ; the original grid defines the data on the points defining the corners of the triangles (the nodes)
  ; but for cell fill we need the values at the center of the triangles.
  ; set the data value to the average of the values at the nodes : roughly approximating the cell center

  data_cell = dim_avg(data)

  res                   = True
  res@gsnMaximize       = True
  res@gsnSpreadColors   = True
  res@gsnSpreadColorEnd = -2         ; make sure not to include gray

  ; currently it is necessary to set sfXArray and sfYArray even though they are not used.
  ; But I am setting them to the aveage of the node values anyway -- roughly approximating the center of the triangle where the data values are set
  res@sfXArray          = dim_avg(xb)
  res@sfYArray          = dim_avg(yb)
  res@sfXCellBounds = xb
  res@sfYCellBounds = yb

  res@cnFillOn          = True 
  res@cnFillMode        = "cellfill"
  res@cnCellFillEdgeColor = "black"
;  res@cnLinesOn         = False
;  res@cnLineLabelsOn         = False

  res@lbOrientation     = "Vertical"
  if(pre_ncl_610()) then
    res@lbLabelAutoStride = True
  end if
  contour = gsn_csm_contour(wks,data_cell,res)

  getvalues contour@data
    "sfXCActualStartF" : xs
    "sfXCActualEndF"   : xe
    "sfYCActualStartF" : ys
    "sfYCActualEndF"   : ye
  end getvalues

  res@mpDataBaseVersion     = "MediumRes"     ; "MediumRes" will run faster.
  res@mpLimitMode           = "LatLon"
  res@mpMinLonF             = xs -0.25
  res@mpMaxLonF             = xe +0.25
  res@mpMinLatF             = ys -0.25
  res@mpMaxLatF             = ye +0.25

  res@pmTickMarkDisplayMode = "Always"

  ;map = gsn_csm_contour_map(wks,data_cell,res)
  ; draw the grid only ; contour does not like it if there are no data values or they are constant
  res@cnFillOpacityF = 0.0
  res@lbLabelBarOn = False
  map = gsn_csm_contour_map(wks,data_cell,res)

end