Index: /issm/trunk-jpl/src/jl/exp.jl
===================================================================
--- /issm/trunk-jpl/src/jl/exp.jl	(revision 26471)
+++ /issm/trunk-jpl/src/jl/exp.jl	(revision 26472)
@@ -94,5 +94,5 @@
 			#read Info
 			A = readline(f); line += 1
-			if A[1:20]!="# Points Count Value"
+			if A[1:14]!="# Points Count"
 				println("line $(line): $(A)")
 				error("Unexpected exp file formatting") 
Index: /issm/trunk-jpl/src/jl/issm.jl
===================================================================
--- /issm/trunk-jpl/src/jl/issm.jl	(revision 26471)
+++ /issm/trunk-jpl/src/jl/issm.jl	(revision 26472)
@@ -34,11 +34,19 @@
 	return Geometry( Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0))
 end# }}}
+mutable struct Mask
+	ocean_levelset::Vector{Float64}
+	ice_levelset::Vector{Float64}
+end
+function Mask() #{{{
+	return Mask( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
+end# }}}
 
 mutable struct model
 	mesh::AbstractMesh
 	geometry::Geometry
+	mask::Mask
 end
 function model() #{{{
-	return model(Mesh2dTriangle(),Geometry())
+	return model( Mesh2dTriangle(), Geometry(), Mask())
 end#}}}
 function Base.show(io::IO, md::model)# {{{
@@ -49,4 +57,5 @@
 	@printf "%19s: %-22s -- %s\n" "mesh" typeof(md.mesh) "mesh properties"
 	@printf "%19s: %-22s -- %s\n" "geometry" typeof(md.geometry) "surface elevation, bedrock topography, ice thickness,..."
+	@printf "%19s: %-22s -- %s\n" "mask" typeof(md.mask) "defines grounded and floating regions"
 
 end# }}}
Index: /issm/trunk-jpl/src/jl/plotmodel.jl
===================================================================
--- /issm/trunk-jpl/src/jl/plotmodel.jl	(revision 26472)
+++ /issm/trunk-jpl/src/jl/plotmodel.jl	(revision 26472)
@@ -0,0 +1,23 @@
+
+#import ColorSchemes.leonardo
+#include("issm.jl")
+
+using GLMakie
+
+function plotmodel(md::model)
+
+	elementcolor = :yellow
+	vertexcolor  = :black
+	facetcolor   = :blue
+	showvertices = true
+	showfacets   = true
+
+	scene = Makie.mesh( [md.mesh.x md.mesh.y], md.mesh.elements, shading = false, color = elementcolor)
+
+	if showvertices
+		Makie.scatter!( [md.mesh.x md.mesh.y], markersize = 4, color = vertexcolor)
+	end
+
+
+	return scene
+end
Index: /issm/trunk-jpl/src/jl/triangle.jl
===================================================================
--- /issm/trunk-jpl/src/jl/triangle.jl	(revision 26471)
+++ /issm/trunk-jpl/src/jl/triangle.jl	(revision 26472)
@@ -1,4 +1,2 @@
-#see https://github.com/JuliaGeometry/Triangulate.jl/blob/master/src/plot.jl
-#see https://github.com/JuliaGeometry/MeshViz.jl/blob/master/src/simplemesh.jl
 
 #Class Triangle's triangulateio
@@ -52,6 +50,8 @@
 end# }}}
 
-using Printf, PyPlot
-using GLMakie
+using Printf #needed for sprintf
+include("issm.jl")
+include("exp.jl")
+
 """
 TRIANGLE - create model mesh using the triangle package
@@ -70,5 +70,5 @@
  - md=triangle(md,'DomainOutline.exp','Rifts.exp',1500);
 """
-function triangle(domainname::String,resolution::Float64) #{{{
+function triangle(md::model,domainname::String,resolution::Float64) #{{{
 
 	#read input file
@@ -161,21 +161,13 @@
 
 	#post process output
-	numberofvertices = ctio_out.numberofpoints
-	numberofelements = ctio_out.numberoftriangles
-	points           = convert(Array{Cdouble,2}, Base.unsafe_wrap(Array, ctio_out.pointlist,    (2,Int(ctio_out.numberofpoints)), own=true))'
-	elements         = convert(Array{Cint,2},    Base.unsafe_wrap(Array, ctio_out.trianglelist, (3,Int(ctio_out.numberoftriangles)), own=true))' .+1
+	points    = convert(Array{Cdouble,2}, Base.unsafe_wrap(Array, ctio_out.pointlist,    (2,Int(ctio_out.numberofpoints)), own=true))'
+	triangles = convert(Array{Cint,2},    Base.unsafe_wrap(Array, ctio_out.trianglelist, (3,Int(ctio_out.numberoftriangles)), own=true))' .+1
+	md.mesh = Mesh2dTriangle()
+	md.mesh.numberofvertices = ctio_out.numberofpoints
+	md.mesh.numberofelements = ctio_out.numberoftriangles
+	md.mesh.x                = points[:,1]
+	md.mesh.y                = points[:,2]
+	md.mesh.elements         = triangles
 
-	#plot
-	clf()
-	if false
-		poly(points, elements, strokewidth=1, shading=false,transparency=true)
-	else
-		fig = matplotlib[:pyplot][:figure]("2D Mesh Plot", figsize = (10,10))
-		ax = matplotlib[:pyplot][:axes]()
-		ax[:set_aspect]("equal")
-		tri = ax[:triplot](points[:,1], points[:,2], elements.-1)
-		setp(tri, linestyle = "-", linewidth = 1.5, marker = "None", markersize = 5., color = "blue")
-		fig[:canvas][:draw]()
-	end
-	
+   return md
 end#}}}