diff --git a/src/usr/utils.jl b/src/usr/utils.jl
index ce86e8c..2c5955c 100644
--- a/src/usr/utils.jl
+++ b/src/usr/utils.jl
@@ -266,6 +266,31 @@ function comparefields(fieldname::String, f1::String, f2::String)# {{{
 		@printf "%s differs\n" fieldname
 	end
 end #}}}
+function comparefields(fieldname::String, f1::Union{Bool, Float64, Int64}, f2::Union{Bool, Float64, Int64})# {{{
+	if f1!=f2
+		if ((!isnan(f1)) | (!isnan(f2)))
+			@printf "%s differs by %s \n" fieldname (f1-f2)
+		end
+	end
+end #}}}
+function comparefields(fieldname::String, f1::Union{Vector,Matrix}, f2::Union{Vector,Matrix})# {{{
+	if size(f1) != size(f2)
+		@printf "%s do not have the same size\n" fieldname
+	else
+		if any(f1 .!= f2)
+			# deal with NaN
+			pos1 = findall(isnan.(f1) .== false)
+			pos2 = findall(isnan.(f2) .== false)
+			if pos1 == pos2
+				if f1[pos1] != f2[pos2]
+					@printf "%s differs by norm2=%s\n" fieldname sum((f1[pos1]-f2[pos2]).^2)
+				end
+			else
+				@printf "%s differs, not the same size of NaN\n" fieldname
+			end
+		end
+	end
+end #}}}
 function comparefields(fieldname::String, f1::Union{Bool,Float64,Int64,Vector,Matrix}, f2::Union{Bool,Float64,Int64,Vector,Matrix})# {{{
 	if size(f1) != size(f2)
 		@printf "%s do not have the same size\n" fieldname
diff --git a/test/test501.jl b/test/test501.jl
index cdfb4bb..b81d589 100755
--- a/test/test501.jl
+++ b/test/test501.jl
@@ -19,15 +19,17 @@ md.initialization.vy=InterpFromMeshToMesh2d(index, x, y, vy_obs, md.mesh.x, md.m
 md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y, 0.0)
 md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y, 0.0)
 md.geometry.base=md.geometry.surface .- md.geometry.thickness
-md.geometry.bed =md.geometry.base
+md.geometry.bed = copy(md.geometry.base)
 pos = findall(md.mask.ocean_levelset.<0)
 md.geometry.bed[pos] = InterpFromMeshToMesh2d(index, x, y, bed, md.mesh.x[pos], md.mesh.y[pos])
+md.masstransport.min_thickness = 1.0
 
 md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices)
 md.materials.rheology_n=3*ones(md.mesh.numberofelements)
 md.friction.coefficient=50*ones(md.mesh.numberofvertices)
-md.friction.p=ones(md.mesh.numberofvertices)
-md.friction.q=ones(md.mesh.numberofvertices)
+md.friction.coefficient[findall(md.mask.ocean_levelset.<0.)] .= 0.
+md.friction.p=ones(md.mesh.numberofelements)
+md.friction.q=ones(md.mesh.numberofelements)
 
 md.stressbalance.restol=0.05
 md.stressbalance.reltol=1.0
@@ -45,13 +47,13 @@ md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices)
 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.stressbalance.spcvx[pos] .= md.initialization.vx[pos]
+md.stressbalance.spcvy[pos] .= md.initialization.vy[pos]
 
 md=solve(md, :Stressbalance)
 
 field_names =["Vx","Vy","Vel"]
-field_tolerances=[NaN,NaN,NaN]
+field_tolerances=[2.0e-1,2.0e-1,2.0e-1]
 field_values= [(md.results["StressbalanceSolution"]["Vx"]),
                (md.results["StressbalanceSolution"]["Vy"]),
                (md.results["StressbalanceSolution"]["Vel"]) ]