begin 

in = "./data/"
out = "./results/"
H_filename = "ice_thickness.nc"
H_file = addfile(in+H_filename,"r")
lonlat_filename = "greenland_4km_epsg3413_c170429.nc"
lonlat_file = addfile(in+lonlat_filename,"r")

H0            = H_file->thk
time_sim     = ispan(0,16100,100)
ny           =  dimsizes(H0&y1)
nx           =  dimsizes(H0&x1)
ntime        = dimsizes(time_sim)
lat        = lonlat_file->lat
lon       = lonlat_file->lon

printVarSummary(H0)

H = new((/ntime,ny,nx/),float)
H = H0*2000.

printMinMax(H, 0)

H = mask(H,H.lt.10.,False)

  res                       = True              ; plot mods desired
  res@gsnAddCyclic          = False
  res@mpLimitMode       = "Corners"
  res@mpLeftCornerLatF      = lat(0,0)
  res@mpLeftCornerLonF      = lon(0,0)
  res@mpRightCornerLatF     = lat(ny-1,nx-1)
  res@mpRightCornerLonF     = lon(ny-1,nx-1)
  res@trGridType   = "TriangularMesh"  
  res@mpProjection = "Stereographic"  ; Change the map projection.
  res@mpDataBaseVersion     = "mediumres"  
  res@mpCenterLonF = -45.0
  res@mpCenterLatF = 70.0
  res@tfDoNDCOverlay      =   True
  res@mpOutlineOn            = True
  res@cnFillOn               = True
  res@cnFillMode           = "RasterFill"
  ;res@gsnDraw                = False                  ; do not draw picture
  ;res@gsnFrame               = False
  res@cnInfoLabelOn = False        ; turn off contour info label
  res@lbAutoManage          = False             ; we control label bar
  res@pmLabelBarDisplayMode = "Always"          ; turns on label bar
  res@lbOrientation         = "Horizontal"      ; ncl default is vertical
  res@pmLabelBarSide        = "Bottom"          ; default is right
  res@lbLabelStride         = 1                 ; skip every other label
  res@pmLabelBarWidthF      = 0.3              ; default is shorter
  res@pmLabelBarHeightF     = 0.05               ; default is taller
  res@lbLabelFontHeightF    = .01              ; default is HUGE
  res@lbPerimOn             = False             ; default has box
  res@pmLabelBarOrthogonalPosF = 0.04
  
  res@lbTitleOn        = True                  ; turn on title
  res@lbTitleString    = "Ice thickness (m)"
  res@lbTitlePosition  = "Bottom"  
  ; title position  
  res@lbTitleFontHeightF= .01
  
  res@cnLinesOn              = False                   ;  contour lines
  res@cnLineLabelsOn         = False                   ; no contour labels
  res@gsnSpreadColors        = False                    ; use total colormap
  res@cnInfoLabelOn          = False              ; no contour info (right left)
  res@pmTickMarkDisplayMode  = "Never" 
  res@cnLevelSelectionMode = "ManualLevels"	; manurefy set the contour levels with the following 3 resources
  res@cnMinLevelValF  = 100.			; set the minimum contour level
  res@cnMaxLevelValF  = 3000.			; set the maximum contour level
  res@cnLevelSpacingF = 500.
  cmap = read_colormap_file("GMT_copper")
  res@cnFillPalette = cmap(5:45:5,:)   ; subset color map
  res@gsnMaximize = True

do i=0,ntime-1
	wks = gsn_open_wks("pdf","./results/image-"+sprinti("%03i",i))
        print("simulation time(" + i + ") = " + time_sim(i))
        res@gsnRightString = "Simulation time: " + time_sim(i) + " (yr)"
	ice = new((/ny,nx/),double)
	ice(:,:) = H(i,:,:)
	plot = new(1,graphic)
	plot(0) = gsn_csm_contour_map(wks,ice,res)
	delete(wks)
end do

cmd = "cd ./results/ ; pdflatex movie.tex"
system(cmd)

end