begin
;************************************************
; open netCDF file
;************************************************
;  a = addfile("ne30pg3_gmted2010_bedmachine_nc3000_NoAniso_Laplace0100_20220516_smooth_over_ocean.nc","r")
;  a = addfile("/glade/p/cgd/amp/pel/topo/cubedata/gmted2010_modis_bedmahcine-ncube3000-220518.nc","r")
  a = addfile("/glade/scratch/pel/topowopts_new/regression-test-data/gmted2010_bedmachine-ncube0540-220518.nc","r")
;  b = addfile("/glade/p/cgd/amp/pel/topo/cubedata/gmted2010_bedmachine-ncube3000.nc","r")
;  b = addfile("gmted2010_modis-ncube3000-stitch.nc","r")
;  b = addfile("gmted2010_bedmachine-ncube3000.nc","r")
  ;
  ; SCRIP file
  ;
;  sf =  addfile(scripFile,"r")
  cellfill = False
;************************************************
; Read in Regression Coef
;************************************************
;  PHIS  = a->terr
;  PHIS_ref=b->terr
  PHIS  = a->LANDFRAC
;  PHIS_ref=b->LANDFRAC
;  PHIS = PHIS-PHIS_ref
;  PHIS = a->GBXAR
;************************************************
; create plot
;************************************************

  wks = gsn_open_wks("pdf","raster")         ; send graphics to PNG file
;  gsn_define_colormap(wks,"BlAqGrYeOrReVi200")   ; choose colormap
  gsn_define_colormap(wks,"cmocean_balance")   ; choose colormap


  res                     = True         ; plot modifications desired
  res@gsnMaximize         = False         ; Maximize size of plot in frame
  res@cnFillOn            = True         ; Turn on contour fill
                                         ; use "CellFill" and "RasterFill"
  res@cnLinesOn           = False        ; Turn off contour lines
  res@cnLineLabelsOn      = False        ; Turn off contour line labels
  res@lbLabelAutoStride   = True         ; Clean up labelbar labels.

  res@tiMainString = "Latest bednmachine topo update (Rene; 5/18/2022)"
  res@gsnStringFontHeightF = 0.015
  if (cellfill) then
    lat1d   = sf->grid_center_lat
    lon1d   = sf->grid_center_lon
    latvert = sf->grid_corner_lat
    lonvert = sf->grid_corner_lon
	
    if (lonvert@units.eq."radians") then
      r2d  = get_r2d("double")
      latvert = r2d*latvert
      lonvert = r2d*lonvert
    end if

    res@cnFillMode = "CellFill"
    res@sfXCellBounds = lonvert
    res@sfYCellBounds = latvert

    delete([/sf,latvert,lonvert/])
  else
    lat1d = a->lat
    lon1d = a->lon
    res@cnFillMode          ="RasterFill"
  end if   

  res@sfXArray            = lon1d        ; Required to tell NCL where to
  res@sfYArray            = lat1d        ; overlay data on globe.

;===================================
; you can have the contour lines on in raster mode, but their thickness
; actually make the plot look like is was contoured normally.

;  res@cnFillMode       = "RasterFill"       ; Raster Mode
;================================
; these three resources make the continents look nice. The first two
; make them color, and the later adds continental outlines when
; raster mode is used.

;  res@cnLevelSelectionMode = "ManualLevels"	; manually set the contour levels with the following 3 resources
;  res@cnMinLevelValF  = -10			; set the minimum contour level
;  res@cnMaxLevelValF  = 10			; set the maximum contour level
;  res@cnLevelSpacingF = 1

;  res@cnMinLevelValF  =  0.0                     ; set the minimum contour level
;  res@cnMaxLevelValF  = 1.0                      ; set the maximum contour level
;  res@cnLevelSpacingF = 0.001


;
  res@mpCenterLonF     = 180                ; set map center at 180

  res@mpMinLatF      = 15
  res@mpMaxLatF      = 40
  res@mpMinLonF      = 60
  res@mpMaxLonF      = 100
                                            ; must be in colormap  
  plot = gsn_csm_contour_map_ce(wks,PHIS, res) ; create plot

end