Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26710)
+++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26711)
@@ -209,6 +209,4 @@
 			end
 		end
-		#error("S")
-		display(A)
 	end
 
Index: /issm/trunk-jpl/src/jl/test201.jl
===================================================================
--- /issm/trunk-jpl/src/jl/test201.jl	(revision 26711)
+++ /issm/trunk-jpl/src/jl/test201.jl	(revision 26711)
@@ -0,0 +1,51 @@
+#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia
+include("issm.jl")
+using .ISSM
+
+md = ISSM.model()
+md = ISSM.triangle(md,"../../test/Exp/Square.exp",50000.)
+md = ISSM.setmask(md,"all","")
+
+#Geometry
+hmin=300.
+hmax=1000.
+ymin=minimum(md.mesh.y)
+ymax=maximum(md.mesh.y)
+xmin=minimum(md.mesh.x)
+xmax=maximum(md.mesh.x)
+md.geometry.thickness = hmax .+ (hmin-hmax)*(md.mesh.y .- ymin)./(ymax-ymin) .+ 0.1*(hmin-hmax)*(md.mesh.x .- xmin)./(xmax-xmin)
+md.geometry.base      = -md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
+md.geometry.surface   = md.geometry.base+md.geometry.thickness
+md.geometry.bed       = md.geometry.base .-10
+
+#Initial velocity
+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)
+
+md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices)
+md.materials.rheology_n=3*ones(md.mesh.numberofelements)
+md.friction.coefficient=20*ones(md.mesh.numberofvertices)
+
+md.stressbalance.restol=0.1
+md.stressbalance.reltol=0.02
+md.stressbalance.abstol=NaN
+
+#Boundary conditions
+nodefront=ISSM.ContourToNodes(md.mesh.x,md.mesh.y,"../../test/Exp/SquareFront.exp",2.0) .& md.mesh.vertexonboundary
+md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices)
+md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices)
+pos = findall(md.mesh.vertexonboundary .& .~nodefront)
+md.mask.ice_levelset[findall(nodefront)] .= 0
+
+segmentsfront=md.mask.ice_levelset[md.mesh.segments[:,1:2]]==0
+segments = findall(vec(sum(Int64.(md.mask.ice_levelset[md.mesh.segments[:,1:2]].==0), dims=2)) .!=2)
+pos=md.mesh.segments[segments,1:2]
+md.stressbalance.spcvx[pos] .= 0.0
+md.stressbalance.spcvy[pos] .= 0.0
+
+md=ISSM.solve(md,"Stressbalance")
Index: /issm/trunk-jpl/src/jl/test501.jl
===================================================================
--- /issm/trunk-jpl/src/jl/test501.jl	(revision 26711)
+++ /issm/trunk-jpl/src/jl/test501.jl	(revision 26711)
@@ -0,0 +1,46 @@
+#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia
+include("issm.jl")
+using .ISSM
+
+md = ISSM.model()
+md = ISSM.triangle(md,"../../test/Exp/Pig.exp",10000.)
+md = ISSM.setmask(md,"../../test/Exp/PigShelves.exp","../../test/Exp/PigIslands.exp")
+error("not implemented yet")
+
+#Geometry
+hmin=300.
+hmax=1000.
+ymin=minimum(md.mesh.y)
+ymax=maximum(md.mesh.y)
+xmin=minimum(md.mesh.x)
+xmax=maximum(md.mesh.x)
+md.geometry.thickness = hmax .+ (hmin-hmax)*(md.mesh.y .- ymin)./(ymax-ymin) .+ 0.1*(hmin-hmax)*(md.mesh.x .- xmin)./(xmax-xmin)
+md.geometry.base      = -md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness .+20
+md.geometry.surface   = md.geometry.base+md.geometry.thickness
+md.geometry.bed       = md.geometry.base
+
+#Initial velocity
+#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)
+
+md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices)
+md.materials.rheology_n=3*ones(md.mesh.numberofelements)
+md.friction.coefficient=20*ones(md.mesh.numberofvertices)
+
+md.stressbalance.restol=0.05
+md.stressbalance.reltol=0.05
+md.stressbalance.abstol=NaN
+
+#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")