Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26704)
+++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26705)
@@ -44,6 +44,13 @@
 	FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum)
 	FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum)
-
-	FetchDataToInput(md,inputs,elements,md.friction.coefficient,FrictionCoefficientEnum)
+	FetchDataToInput(md,inputs,elements,md.mask.ice_levelset, MaskIceLevelsetEnum)
+	FetchDataToInput(md,inputs,elements,md.mask.ocean_levelset, MaskOceanLevelsetEnum)
+
+	#Deal with friction
+	if typeof(md.friction)==BuddFriction
+		FetchDataToInput(md,inputs,elements,md.friction.coefficient,FrictionCoefficientEnum)
+	else
+		error("Friction ", typeof(md.friction), " not supported yet")
+	end
 
 end#}}}
@@ -53,4 +60,11 @@
 	AddParam(parameters,md.stressbalance.abstol,StressbalanceAbstolEnum)
 	AddParam(parameters,md.stressbalance.maxiter,StressbalanceMaxiterEnum)
+
+	#Deal with friction
+	if typeof(md.friction)==BuddFriction
+		AddParam(parameters, 1, FrictionLawEnum)
+	else
+		error("Friction ", typeof(md.friction), " not supported yet")
+	end
 end#}}}
 
@@ -104,24 +118,27 @@
 
 	#Add basal friction
-	basis = Vector{Float64}(undef,numnodes)
-	friction = Friction(element)
-
-	#Start integrating
-	gauss = GaussTria(2)
-	for ig in 1:gauss.numgauss
-
-		Jdet = JacobianDeterminant(xyz_list, gauss)
-		NodalFunctions(element, basis, gauss, ig, P1Enum)
-
-		alpha2 = Alpha2(friction, gauss, ig)
-
-		for i in 1:numnodes
-			for j in 1:numnodes
-				Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*alpha2*basis[i]*basis[j]
-				Ke.values[2*i  ,2*j  ] += gauss.weights[ig]*Jdet*alpha2*basis[i]*basis[j]
+	phi=GetGroundedPortion(element, xyz_list)
+
+	if(phi>0)
+		basis = Vector{Float64}(undef,numnodes)
+		friction = CoreFriction(element)
+
+		#Start integrating
+		gauss = GaussTria(2)
+		for ig in 1:gauss.numgauss
+
+			Jdet = JacobianDeterminant(xyz_list, gauss)
+			NodalFunctions(element, basis, gauss, ig, P1Enum)
+
+			alpha2 = Alpha2(friction, gauss, ig)
+
+			for i in 1:numnodes
+				for j in 1:numnodes
+					Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*phi*alpha2*basis[i]*basis[j]
+					Ke.values[2*i  ,2*j  ] += gauss.weights[ig]*Jdet*phi*alpha2*basis[i]*basis[j]
+				end
 			end
 		end
 	end
-	
 
 	return Ke
Index: /issm/trunk-jpl/src/jl/solve/elements.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26704)
+++ /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26705)
@@ -89,4 +89,16 @@
 
 end # }}}
+function GetInputListOnVertices!(element::Tria, vector::Vector{Float64}, enum::IssmEnum) # {{{
+
+	#Get Input first 
+	input = GetInput(element, enum)
+
+	#Get value at each vertex (i.e. P1 Nodes)
+	gauss=GaussTria(P1Enum)
+	for i in 1:gauss.numgauss
+		vector[i] = GetInputValue(input, gauss, i)
+	end
+
+end # }}}
 function FindParam(element::Tria,enum::IssmEnum) # {{{
 
@@ -107,4 +119,45 @@
 
 end # }}}
+function GetGroundedPortion(element::Tria, xyz_list::Matrix{Float64}) #{{{
+
+	level = Vector{Float64}(undef,3)
+	GetInputListOnVertices!(element, level, MaskOceanLevelsetEnum)
+
+	#Be sure that values are not zero
+	epsilon = 1.e-15
+	for i in 1:3
+		if(level[i]==0.) level[i]=level[i]+epsilon end
+	end
+
+	if level[1]>0 && level[2]>0 && level[3]>0
+      #Completely grounded
+		phi = 1.0
+   elseif level[1]<0 && level[2]<0 && level[3]<0
+      #Completely floating
+      phi = 0.0
+   else
+		#Partially floating,
+		if(level[1]*level[2]>0) #Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+			s1=level[3]/(level[3]-level[2]);
+			s2=level[3]/(level[3]-level[1]);
+		elseif(level[2]*level[3]>0) #Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+			s1=level[1]/(level[1]-level[2]);
+			s2=level[1]/(level[1]-level[3]);
+		elseif(level[1]*level[3]>0) #Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+			s1=level[2]/(level[2]-level[1]);
+			s2=level[2]/(level[2]-level[3]);
+		else
+			error("not supposed to be here...")
+		end
+
+		if(level[1]*level[2]*level[3]>0)
+			phi = s1*s2
+		else
+			phi = (1-s1*s2)
+		end
+	end
+
+	return phi
+end#}}}
 
 #Finite Element stuff
Index: /issm/trunk-jpl/src/jl/solve/friction.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/friction.jl	(revision 26704)
+++ /issm/trunk-jpl/src/jl/solve/friction.jl	(revision 26705)
@@ -1,4 +1,6 @@
 #Friction class definition
-struct Friction#{{{
+
+abstract type CoreFriction end
+struct CoreBuddFriction#{{{
 	H_input::ElementInput
 	b_input::ElementInput
@@ -8,19 +10,25 @@
 	g::Float64
 end# }}}
-function Friction(element::Tria) #{{{
+function CoreFriction(element::Tria) #{{{
 
-	H_input  = GetInput(element, ThicknessEnum)
-	b_input  = GetInput(element, BaseEnum)
-	c_input  = GetInput(element, FrictionCoefficientEnum)
+	frictionlaw = FindParam(element, FrictionLawEnum)
 
-	rho_ice   = FindParam(element, MaterialsRhoIceEnum)
-	rho_water = FindParam(element, MaterialsRhoSeawaterEnum)
-	g         = FindParam(element, ConstantsGEnum)
+	if frictionlaw==1
+		H_input  = GetInput(element, ThicknessEnum)
+		b_input  = GetInput(element, BaseEnum)
+		c_input  = GetInput(element, FrictionCoefficientEnum)
 
-	return Friction(H_input, b_input, c_input, rho_ice, rho_water, g)
+		rho_ice   = FindParam(element, MaterialsRhoIceEnum)
+		rho_water = FindParam(element, MaterialsRhoSeawaterEnum)
+		g         = FindParam(element, ConstantsGEnum)
+
+		return CoreBuddFriction(H_input, b_input, c_input, rho_ice, rho_water, g)
+	else
+		error("Friction ",typeof(md.friction)," not supported yet")
+	end
 end#}}}
 
 #vertices functions
-function Alpha2(friction::Friction, gauss::GaussTria, i::Int64) #{{{
+function Alpha2(friction::CoreBuddFriction, gauss::GaussTria, i::Int64) #{{{
 
 	#Get effective pressure