Index: ../trunk-jpl/src/jl/md/parameterization.jl
===================================================================
--- ../trunk-jpl/src/jl/md/parameterization.jl	(revision 26712)
+++ ../trunk-jpl/src/jl/md/parameterization.jl	(revision 26713)
@@ -16,8 +16,8 @@
 """
 function setmask(md::model,floatingicename::String,groundedicename::String)
 
-	elementonfloatingice = FlagElements( md, floatingicename);
-	elementongroundedice = FlagElements( md, groundedicename);
+	elementonfloatingice = FlagElements( md, floatingicename)
+	elementongroundedice = FlagElements( md, groundedicename)
 
 	elementonfloatingice = convert( Array{Float64}, (elementonfloatingice.>0) .&  (elementongroundedice.==0.))
 	elementongroundedice = convert( Array{Float64}, elementonfloatingice.==0.)
@@ -44,7 +44,9 @@
 	elseif region == "all"
 		flags = ones(md.mesh.numberofelements)
 	else
-		flags = InContour(md.mesh.x,md.mesh.y,region);
+		xcenter = md.mesh.x[md.mesh.elements]*[1;1;1]/3
+		ycenter = md.mesh.y[md.mesh.elements]*[1;1;1]/3
+		flags = ContourToNodes(xcenter, ycenter, region, 2.)
 	end
 
 	return flags
Index: ../trunk-jpl/src/jl/test101.jl
===================================================================
--- ../trunk-jpl/src/jl/test101.jl	(revision 26712)
+++ ../trunk-jpl/src/jl/test101.jl	(revision 26713)
@@ -25,12 +25,13 @@
 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)
+index = Int.(ISSM.archread(ISSM.issmdir()*"/test/Data/SquareShelfConstrained.arch","index"))
+md.initialization.vx=ISSM.InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y,0.0)
+md.initialization.vy=ISSM.InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y,0.0)
 
 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
Index: ../trunk-jpl/src/jl/test501.jl
===================================================================
--- ../trunk-jpl/src/jl/test501.jl	(revision 26712)
+++ ../trunk-jpl/src/jl/test501.jl	(revision 26713)
@@ -5,35 +5,31 @@
 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 and geometry
+x         = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","x")
+y         = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","y")
+vx_obs    = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","vx_obs")
+vy_obs    = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","vy_obs")
+index     = Int.(ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","index"))
+surface   = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","surface")
+thickness = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","thickness")
+bed       = ISSM.archread(ISSM.issmdir()*"/test/Data/Pig.arch","bed")
+md.initialization.vx=ISSM.InterpFromMeshToMesh2d(index, x, y, vx_obs, md.mesh.x, md.mesh.y, 0.0)
+md.initialization.vy=ISSM.InterpFromMeshToMesh2d(index, x, y, vy_obs, md.mesh.x, md.mesh.y, 0.0)
+md.geometry.surface = ISSM.InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y, 0.0)
+md.geometry.thickness = ISSM.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
+pos = findall(md.mask.ocean_levelset.<0)
+md.geometry.bed[pos] = ISSM.InterpFromMeshToMesh2d(index, x, y, bed, md.mesh.x[pos], md.mesh.y[pos])
 
-#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.friction.coefficient=50*ones(md.mesh.numberofvertices)
 
 md.stressbalance.restol=0.05
-md.stressbalance.reltol=0.05
+md.stressbalance.reltol=1.0
 md.stressbalance.abstol=NaN
 
 #Boundary conditions