Index: ../trunk-jpl/src/jl/plotmodel.jl
===================================================================
--- ../trunk-jpl/src/jl/plotmodel.jl	(revision 26491)
+++ ../trunk-jpl/src/jl/plotmodel.jl	(revision 26492)
@@ -27,7 +27,7 @@
 			fig, ax, h = Makie.mesh( [md.mesh.x md.mesh.y], md.mesh.elements, shading = false, color = data, colormap = jet)
 
 			#Add colorbar
-			Colorbar(fig[1, 2], h, width=25)
+			#Colorbar(fig[1, 2], h, width=25)
 		else
 			error("data of size "*string(length(data))*" not supported yet!")
 		end
Index: ../trunk-jpl/src/jl/issm.jl
===================================================================
--- ../trunk-jpl/src/jl/issm.jl	(revision 26491)
+++ ../trunk-jpl/src/jl/issm.jl	(revision 26492)
@@ -12,6 +12,19 @@
 function Mesh2dTriangle() #{{{
 	return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int32}(undef, 0, 0))
 end# }}}
+function Base.show(io::IO, this::Mesh2dTriangle)# {{{
+
+	println(io,typeof(this),":")
+	for name in fieldnames(typeof(this))
+		a=getfield(this,name)
+		print(io,"   $(name) = ")
+		if !isempty(a)
+			println(io, typeof(a), " of size ", size(a))
+		else
+			println(io,"empty")
+		end
+	end
+end# }}}
 mutable struct Mesh3dPrism{T} <: AbstractMesh
 	numberofvertices::Int32
 	numberofelements::Int32
@@ -40,14 +53,44 @@
 function Mask() #{{{
 	return Mask( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
 end# }}}
+mutable struct Initialization
+	vx::Vector{Float64}
+	vy::Vector{Float64}
+end
+function Initialization() #{{{
+	return Initialization( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
+end# }}}
+mutable struct Materials
+	rho_ice::Float64
+	rho_water::Float64
+	rho_freshwater::Float64
+	mu_water::Float64
+	heatcapacity::Float64
+	latentheat::Float64
+	thermalconductivity::Float64
+	temperateiceconductivity::Float64
+	effectiveconductivity_averaging::Int32
+	meltingpoint::Float64
+	beta::Float64
+	mixed_layer_capacity::Float64
+	thermal_exchange_velocity::Float64
+	rheology_B::Vector{Float64}
+	rheology_n::Vector{Float64}
+	rheology_law::String
+end
+function Materials() #{{{
+	return Materials(917., 1023., 1000., 0.001787, 2093., 3.34*10^5, 2.4, .24, 1, 273.15, 9.8*10^-8, 3974., 1.00*10^-4, Vector{Float64}(undef,0), Vector{Float64}(undef,0), "Cuffey")
+end# }}}
 
 mutable struct model
 	mesh::AbstractMesh
 	geometry::Geometry
 	mask::Mask
+	materials::Materials
+	initialization::Initialization
 end
 function model() #{{{
-	return model( Mesh2dTriangle(), Geometry(), Mask())
+	return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization())
 end#}}}
 function Base.show(io::IO, md::model)# {{{
 
@@ -57,12 +100,58 @@
 	@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"
+	@printf "%19s: %-22s -- %s\n" "materials" typeof(md.materials) "material properties"
+	@printf "%19s: %-22s -- %s\n" "initialization" typeof(md.initialization) "initial state"
 
 end# }}}
 
-# Probably actually want something more like:
-# function solve(mesh::AbstractMesh, data::AbstractDataset)
-#     do_something_with_mesh(mesh) # This will dispatch
-#     do_something_with_data(data) # This will dispatch
-#     # do something
-# end
+#utils
+function issmdir() #{{{
+	issmdir = ENV["ISSM_DIR"]
+
+	if isempty(issmdir)
+		error("Could not determine the location of ISSM")
+	else
+		return issmdir
+	end
+end#}}}
+function archread(filename::String,variablename::String) #{{{
+
+	#initialize variables
+	found = false
+
+	#open file
+	output = open(filename, "r") do f
+
+		while !eof(f)
+			reclen  = bswap(read(f, Int32))
+			rectype = bswap(read(f, Int32))
+			if rectype!=1
+				error("Expected variable of type string")
+			else
+				fieldname_length = bswap(read(f, Int32))
+				field_name = String(read(f, fieldname_length))
+			end
+			rec_length = bswap(read(f, Int32))
+			field_type = bswap(read(f, Int32))
+			if field_type==2
+				data = bswap(read(f, Float64))
+			elseif field_type==3
+				rows = bswap(read(f, Int32))
+				cols = bswap(read(f, Int32))
+				data = reinterpret(Float64, read(f, sizeof(Float64)*rows*cols))
+				data .= ntoh.(data)
+				data = reshape(data, (rows,cols))
+			else
+				error("Error: Encountered invalid field type when reading data.")
+			end
+
+			if field_name == variablename
+				found = true
+				return data
+			end
+		end
+	end
+
+	return output
+end# }}}