Index: /issm/trunk-jpl/src/jl/classes.jl
===================================================================
--- /issm/trunk-jpl/src/jl/classes.jl	(revision 26522)
+++ /issm/trunk-jpl/src/jl/classes.jl	(revision 26522)
@@ -0,0 +1,166 @@
+using Printf
+
+abstract type AbstractMesh end
+mutable struct Mesh2dTriangle <: AbstractMesh
+	numberofvertices::Int64
+	numberofelements::Int64
+	x::Vector{Float64}
+	y::Vector{Float64}
+	elements::Matrix{Int64}
+	segments::Matrix{Int64}
+	vertexonboundary::Vector{Bool}
+end
+function Mesh2dTriangle() #{{{
+	return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int64}(undef, 0, 0), Matrix{Int64}(undef, 0, 0), Vector{Bool}(undef,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) = ")
+		@printf "%19s: " name
+		if isa(a,String)
+			println(io, a)
+		elseif length(a)>1
+			if !isempty(a)
+				println(io, typeof(a), " of size ", size(a))
+			else
+				println(io,"empty")
+			end
+		else
+			println(io, a)
+		end
+	end
+end# }}}
+mutable struct Mesh3dPrism{T} <: AbstractMesh
+	numberofvertices::Int64
+	numberofelements::Int64
+	numberoflayers::Int64
+	x::Vector{Float64}
+	y::Vector{Float64}
+	z::Vector{Float64}
+	elements::Matrix{Int64}
+	segments::Matrix{Int64}
+	vertexonboundary::Vector{Bool}
+end
+function Mesh3dPrism() #{{{
+	return Mesh3dPrism( 0, 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Matrix{Int64}(undef, 0, 0), Matrix{Int64}(undef, 0, 0), Vector{Bool}(undef,0))
+end# }}}
+
+mutable struct Geometry
+	surface::Vector{Float64}
+	base::Vector{Float64}
+	thickness::Vector{Float64}
+	bed::Vector{Float64}
+end
+function Geometry() #{{{
+	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 Initialization
+	vx::Vector{Float64}
+	vy::Vector{Float64}
+end
+function Initialization() #{{{
+	return Initialization( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
+end# }}}
+
+mutable struct Stressbalance
+	spcvx::Vector{Float64}
+	spcvy::Vector{Float64}
+	restol::Float64
+	reltol::Float64
+	abstol::Float64
+	maxiter::Int64
+end
+function Stressbalance() #{{{
+	return Stressbalance( Vector{Float64}(undef,0), Vector{Float64}(undef,0), 1.e-4, 0.01, 10., 100)
+end# }}}
+function Base.show(io::IO, this::Stressbalance)# {{{
+
+	println(io,typeof(this),":")
+	for name in fieldnames(typeof(this))
+		a=getfield(this,name)
+		#print(io,"   $(name) = ")
+		@printf "%19s: " name
+		if isa(a,String)
+			println(io, a)
+		elseif length(a)>1
+			if !isempty(a)
+				println(io, typeof(a), " of size ", size(a))
+			else
+				println(io,"empty")
+			end
+		else
+			println(io, a)
+		end
+	end
+end# }}}
+
+mutable struct Constants
+	g::Float64
+	yts::Float64
+end
+function Constants() #{{{
+	return Constants( 9.81,  365*24*3600.)
+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::Int64
+	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
+	stressbalance::Stressbalance
+	constants::Constants
+end
+function model() #{{{
+	return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants())
+end#}}}
+function Base.show(io::IO, md::model)# {{{
+
+	compact = get(io, :compact, false)
+
+	println(io,"Model:")
+	@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"
+	@printf "%19s: %-22s -- %s\n" "stressbalance" typeof(md.stressbalance) "stress balance parameters"
+	@printf "%19s: %-22s -- %s\n" "constants" typeof(md.constants) "physical constants"
+
+end# }}}
+
+
Index: /issm/trunk-jpl/src/jl/issm.jl
===================================================================
--- /issm/trunk-jpl/src/jl/issm.jl	(revision 26521)
+++ /issm/trunk-jpl/src/jl/issm.jl	(revision 26522)
@@ -1,237 +1,15 @@
-using Printf
+# Julia version of the Ice-sheet and Sea-level System Model
+#
+# Author: Mathieu Morlighem
+# email:  mathieu.morlighem@dartmouth.edu
 
-#Define classes used by model struct
-abstract type AbstractMesh end
-mutable struct Mesh2dTriangle <: AbstractMesh
-	numberofvertices::Int64
-	numberofelements::Int64
-	x::Vector{Float64}
-	y::Vector{Float64}
-	elements::Matrix{Int64}
+module ISSM
+
+include("classes.jl")
+include("utils.jl")
+include("triangle.jl")
+include("parameterization.jl")
+include("solve.jl")
+#include("plotmodel.jl")
+
 end
-function Mesh2dTriangle() #{{{
-	return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int64}(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::Int64
-	numberofelements::Int64
-	numberoflayers::Int64
-	x::Vector{Float64}
-	y::Vector{Float64}
-	z::Vector{Float64}
-	elements::Matrix{Int64}
-end
-function Mesh3dPrism() #{{{
-	return Mesh3dPrism( 0, 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Matrix{Int64}(undef, 0, 0))
-end# }}}
-mutable struct Geometry
-	surface::Vector{Float64}
-	base::Vector{Float64}
-	thickness::Vector{Float64}
-	bed::Vector{Float64}
-end
-function Geometry() #{{{
-	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 Initialization
-	vx::Vector{Float64}
-	vy::Vector{Float64}
-end
-function Initialization() #{{{
-	return Initialization( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
-end# }}}
-mutable struct Stressbalance
-	spcvx::Vector{Float64}
-	spcvy::Vector{Float64}
-	restol::Float64
-	reltol::Float64
-	abstol::Float64
-	maxiter::Int64
-end
-function Stressbalance() #{{{
-	return Stressbalance( Vector{Float64}(undef,0), Vector{Float64}(undef,0), 0., 0., 0., 0)
-end# }}}
-mutable struct Constants
-	g::Float64
-	yts::Float64
-end
-function Constants() #{{{
-	return Constants( 9.81,  365*24*3600.)
-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::Int64
-	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
-	stressbalance::Stressbalance
-	constants::Constants
-end
-function model() #{{{
-	return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants())
-end#}}}
-function Base.show(io::IO, md::model)# {{{
-
-	compact = get(io, :compact, false)
-
-	println(io,"Model:")
-	@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"
-	@printf "%19s: %-22s -- %s\n" "stressbalance" typeof(md.stressbalance) "stress balance parameters"
-	@printf "%19s: %-22s -- %s\n" "constants" typeof(md.constants) "physical constants"
-
-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))
-				data = collect(data)
-				if cols == 1
-					data = vec(data)
-				end
-			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# }}}
-function InterpFromMeshToMesh2d(index::Array,x::Vector,y::Vector,data::Vector,xout::Vector,yout::Vector) #{{{
-	#prepare input arrays
-	nods = Cint(length(x))
-	nels = Cint(size(index,1))
-	nods_interp = Cint(length(xout))
-	Cindex=Array{Cint,1}(undef,length(index))
-	for i in 1:size(index,1)
-		for j in 1:3
-			Cindex[(i-1)*3+j] = Int32(index[i,j])
-		end
-	end
-	Cx    = Array{Cdouble,1}(undef,nods)
-	Cy    = Array{Cdouble,1}(undef,nods)
-	Cdata = Array{Cdouble,1}(undef,nods)
-	for i in 1:nods
-		Cx[i]    = x[i]
-		Cy[i]    = y[i]
-		Cdata[i] = data[i]
-	end
-	Cxout = Array{Cdouble,1}(undef,nods_interp)
-	Cyout = Array{Cdouble,1}(undef,nods_interp)
-	for i in 1:nods_interp
-		Cxout[i] = xout[i]
-		Cyout[i] = yout[i]
-	end
-
-	Cdataout = Vector{Float64}(undef,nods_interp)
-
-	#This is not working....
-	rc=ccall( (:InterpFromMeshToMesh2dx,"libISSMCore"),
-				Cint, (Ptr{Ptr{Cdouble}},Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint),
-				Ref(Ref(Cdataout)), Ref(Cindex), Ref(Cx), Ref(Cy), nods, nels,
-				Ref(Cdata), nods, 1, Ref(Cxout), Ref(Cyout), nods_interp)
-
-	#Process output
-	dataout = Vector{Float64}(undef,nods_interp)
-	for i in 1:nods_interp
-		dataout[i] = Cdataout[i]
-	end
-
-	return dataout
-end #}}}
-function solve(md::model,solution::String) #{{{
-
-	if solution=="sb" || solution=="Stressbalance"
-		solutionstring = "StressbalanceSolution"
-	else
-		error("solutionstring "*solution*" not supported!");
-	end
-
-	IssmCore(md)
-
-	return md
-end #}}}
Index: /issm/trunk-jpl/src/jl/parameterization.jl
===================================================================
--- /issm/trunk-jpl/src/jl/parameterization.jl	(revision 26521)
+++ /issm/trunk-jpl/src/jl/parameterization.jl	(revision 26522)
@@ -1,3 +1,3 @@
-include("issm.jl")
+include("classes.jl")
 include("exp.jl")
 
Index: /issm/trunk-jpl/src/jl/plotmodel.jl
===================================================================
--- /issm/trunk-jpl/src/jl/plotmodel.jl	(revision 26521)
+++ /issm/trunk-jpl/src/jl/plotmodel.jl	(revision 26522)
@@ -2,6 +2,8 @@
 
 using GLMakie
+include("classes.jl")
+using .ISSM
 
-function plotmodel(md::model,data)
+function plotmodel(md::ISSM.model,data::Vector)
 
 	vertexcolor  = :black
Index: /issm/trunk-jpl/src/jl/solve.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve.jl	(revision 26521)
+++ /issm/trunk-jpl/src/jl/solve.jl	(revision 26522)
@@ -1,3 +1,3 @@
-include("issm.jl")
+include("classes.jl")
 include("issmenums.jl")
 
@@ -209,5 +209,4 @@
 
 	#Solve
-	error("not there yet")
 	Stressbalance(femmodel)
 
@@ -317,6 +316,25 @@
 end#}}}
 function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{Constraint},md::model) #{{{
-	error("stop")
-end#}}}
+
+	#load constraints from model
+	spcvx = md.stressbalance.spcvx
+	spcvy = md.stressbalance.spcvy
+
+	count = 1
+	for i in 1:md.mesh.numberofvertices
+		if ~isnan(spcvx[i])
+			push!(constraints,Constraint(count,i,1,spcvx[i]))
+			count+=1
+		end
+		if ~isnan(spcvy[i])
+			push!(constraints,Constraint(count,i,2,spcvy[i]))
+			count+=1
+		end
+	end
+
+end#}}}
+function Stressbalance(femmodel::FemModel)
+	error("STOP")
+end
 
 #Element functions
Index: /issm/trunk-jpl/src/jl/test101.jl
===================================================================
--- /issm/trunk-jpl/src/jl/test101.jl	(revision 26521)
+++ /issm/trunk-jpl/src/jl/test101.jl	(revision 26522)
@@ -1,10 +1,9 @@
 #!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia
-include("triangle.jl")
-include("parameterization.jl")
-include("solve.jl")
+include("issm.jl")
+using .ISSM
 
-md = model()
-md = triangle(md,"../../test/Exp/Square.exp",180000.)
-md = setmask(md,"all","")
+md = ISSM.model()
+md = ISSM.triangle(md,"../../test/Exp/Square.exp",180000.)
+md = ISSM.setmask(md,"all","")
 
 #Geometry
@@ -21,9 +20,9 @@
 
 #Initial velocity
-x     = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","x")
-y     = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","y")
-vx    = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vx")
-vy    = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vy")
-index = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","index")
+x     = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","x")
+y     = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","y")
+vx    = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","vx")
+vy    = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","vy")
+index = ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","index")
 md.initialization.vx=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
 md.initialization.vy=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y)
@@ -36,3 +35,10 @@
 md.stressbalance.abstol=NaN
 
-md=solve(md,"Stressbalance")
+#Boundary conditions
+md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices)
+md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices)
+pos = findall(md.mesh.vertexonboundary)
+md.stressbalance.spcvx[pos] .= 0.0
+md.stressbalance.spcvy[pos] .= 0.0
+
+md=ISSM.solve(md,"Stressbalance")
Index: /issm/trunk-jpl/src/jl/triangle.jl
===================================================================
--- /issm/trunk-jpl/src/jl/triangle.jl	(revision 26521)
+++ /issm/trunk-jpl/src/jl/triangle.jl	(revision 26522)
@@ -51,5 +51,5 @@
 
 using Printf #needed for sprintf
-include("issm.jl")
+include("classes.jl")
 include("exp.jl")
 
@@ -162,4 +162,7 @@
 	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
+	segments  = convert(Array{Cint,2},    Base.unsafe_wrap(Array, ctio_out.segmentlist,  (2,Int(ctio_out.numberofsegments)), own=true))' .+1
+	
+	#assign output
 	md.mesh = Mesh2dTriangle()
 	md.mesh.numberofvertices = ctio_out.numberofpoints
@@ -168,4 +171,9 @@
 	md.mesh.y                = points[:,2]
 	md.mesh.elements         = triangles
+	md.mesh.segments         = segments
+
+	#post processing
+	md.mesh.vertexonboundary = zeros(Bool,md.mesh.numberofvertices)
+	md.mesh.vertexonboundary[md.mesh.segments] .= true
 
    return md
Index: /issm/trunk-jpl/src/jl/utils.jl
===================================================================
--- /issm/trunk-jpl/src/jl/utils.jl	(revision 26522)
+++ /issm/trunk-jpl/src/jl/utils.jl	(revision 26522)
@@ -0,0 +1,108 @@
+#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))
+				data = collect(data)
+				if cols == 1
+					data = vec(data)
+				end
+			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# }}}
+function InterpFromMeshToMesh2d(index::Array,x::Vector,y::Vector,data::Vector,xout::Vector,yout::Vector) #{{{
+	#prepare input arrays
+	nods = Cint(length(x))
+	nels = Cint(size(index,1))
+	nods_interp = Cint(length(xout))
+	Cindex=Array{Cint,1}(undef,length(index))
+	for i in 1:size(index,1)
+		for j in 1:3
+			Cindex[(i-1)*3+j] = Int32(index[i,j])
+		end
+	end
+	Cx    = Array{Cdouble,1}(undef,nods)
+	Cy    = Array{Cdouble,1}(undef,nods)
+	Cdata = Array{Cdouble,1}(undef,nods)
+	for i in 1:nods
+		Cx[i]    = x[i]
+		Cy[i]    = y[i]
+		Cdata[i] = data[i]
+	end
+	Cxout = Array{Cdouble,1}(undef,nods_interp)
+	Cyout = Array{Cdouble,1}(undef,nods_interp)
+	for i in 1:nods_interp
+		Cxout[i] = xout[i]
+		Cyout[i] = yout[i]
+	end
+
+	Cdataout = Vector{Float64}(undef,nods_interp)
+
+	#This is not working....
+	rc=ccall( (:InterpFromMeshToMesh2dx,"libISSMCore"),
+				Cint, (Ptr{Ptr{Cdouble}},Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint),
+				Ref(Ref(Cdataout)), Ref(Cindex), Ref(Cx), Ref(Cy), nods, nels,
+				Ref(Cdata), nods, 1, Ref(Cxout), Ref(Cyout), nods_interp)
+
+	#Process output
+	dataout = Vector{Float64}(undef,nods_interp)
+	for i in 1:nods_interp
+		dataout[i] = Cdataout[i]
+	end
+
+	return dataout
+end #}}}
+function solve(md::model,solution::String) #{{{
+
+	if solution=="sb" || solution=="Stressbalance"
+		solutionstring = "StressbalanceSolution"
+	else
+		error("solutionstring "*solution*" not supported!");
+	end
+
+	IssmCore(md)
+
+	return md
+end #}}}