begin

;---Initialize the graphics

; wks_type = "x11"

; wks_type = "ncgm"

; wks_type = "ps"
; wks_type@wkOrientation = "portrait"
; wks_type@wkPaperWidthF  =  8.5  ; in inches
; wks_type@wkPaperHeightF = 11.0  ; in inches

; wks_type = "pdf"
; wks_type@wkPaperSize = "A4"

  wks_type = "png"
  wks_type@wkWidth  = 1000
  wks_type@wkHeight =  800

  plot_name = "mpas_grids"

  wks = gsn_open_wks(wks_type,plot_name)     ; send graphics to PNG file

; ----------- set map resources and prepare map -----------------------

  res                             = True
; res@gsnMaximize                 = True

  res@mpShapeMode                 = "FreeAspect"
  res@vpWidthF                    = 0.80
  res@vpHeightF                   = 0.75

  res@gsnFrame                    = False        ; don't advance the frame yet
  res@gsnDraw                     = False        ; don't draw the plot yet

  res@mpFillOn                    = False        ; turn off gray fill
  res@mpOutlineBoundarySets       = "National"   ; turn on country boundaries
  res@mpGeophysicalLineColor      = "Navy"       ; color of cont. outlines
  res@mpGeophysicalLineThicknessF = 2.0          ; lines separating land/ocean
  res@mpNationalLineThicknessF    = 2.0          ; interior boundaries

  res@pmTickMarkDisplayMode       = "Always"

  res@mpOutlineBoundarySets       = "GeophysicalAndUSStates" ; state boundaries
; res@mpOutlineBoundarySets       = "AllBoundaries"
  res@mpDataBaseVersion           = "MediumRes"
; res@mpDataBaseVersion           = "Ncarg4_1"
  res@mpDataSetName               = "Earth..4"
  res@mpFillBoundarySets          = "National"
  res@mpLandFillColor             = "transparent"
  res@mpOutlineDrawOrder          = "PostDraw"

  res@mpMaxLatF                   = 43           ; choose subregion           
  res@mpMinLatF                   = 31.5
  res@mpMaxLonF                   = -113
  res@mpMinLonF                   = -127.5

  map = gsn_csm_map(wks,res)

; draw(map)     ; just to test if map is plotted
; frame(wks)

; ------------ work on plotting MPAS grid cells -------------------------

;--Open MPAS file

  filename = "../x4x5/x4x5.794626.rotated.grid.static.nc"
  f        = addfile(filename,"r")

  RAD2DEG   = get_r2d("float") ; Radian to Degree

;---Read edge and lat/lon information

  verticesOnEdge = f->verticesOnEdge
  lonCell        = f->lonCell * RAD2DEG
  latCell        = f->latCell * RAD2DEG
  lonVertex      = f->lonVertex * RAD2DEG
  latVertex      = f->latVertex * RAD2DEG

  res@mpProjection         = "CylindricalEquidistant"
  res@mpMinLonF            = -130
  res@mpMaxLonF            = -125
  res@mpMinLatF            =  30
  res@mpMaxLatF            =  35
  res@pmTickMarkDisplayMode = "Always"   ; better map tickmarks

  map = gsn_csm_map(wks,res)  ; Create the map, don't draw it.

;---Code to attach MPAS edge lines to the existing map

  lnres                     = True
  lnres@gsLineThicknessF    = 0.10             ; default is 1
  lnres@gsLineColor         = "NavyBlue"       ; default is black.
  lnres@gsLineThicknessF    = 2

;---This is the code for the MPAS grid edges

  esizes = getfilevardimsizes(f,"latEdge")
  nedges = esizes(0)
  print("Number of edges = " + nedges)

  ecx = new((/nedges,2/),double)
  ecy = new((/nedges,2/),double)

  ecx(:,0) = lonVertex(verticesOnEdge(:,0)-1)
  ecx(:,1) = lonVertex(verticesOnEdge(:,1)-1)
  ecy(:,0) = latVertex(verticesOnEdge(:,0)-1)
  ecy(:,1) = latVertex(verticesOnEdge(:,1)-1)

  ii0 = ind((abs(ecx(:,0)-ecx(:,1)).gt.180.and.(ecx(:,0).gt.ecx(:,1))))
  ii1 = ind((abs(ecx(:,0)-ecx(:,1)).gt.180.and.(ecx(:,0).lt.ecx(:,1))))
  ecx(ii0,0) = ecx(ii0,0) - 360.0
  ecx(ii1,1) = ecx(ii1,1) - 360.0

; Attach the polylines using special "gsSegments" resource. This
; is *much* faster than attaching every line individually.

  start_cpu_time = get_cpu_time()

  print("Attaching the polylines...")
  lnres@gsSegments = ispan(0,nedges * 2,2)

  poly = gsn_add_polyline(wks,map,ndtooned(ecx),ndtooned(ecy),lnres)

  draw(map)
  frame(wks)

  end_cpu_time = get_cpu_time()
  print("CPU elapsed time = " + (end_cpu_time-start_cpu_time))

end