Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp	(revision 14604)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp	(revision 14605)
@@ -826,4 +826,106 @@
 	/*Assign output pointers:*/
 	*pdoflist=doflist;
+}
+/*}}}*/
+/*FUNCTION Penta::GetGroundedPortion{{{*/
+IssmDouble Penta::GetGroundedPortion(IssmDouble* xyz_list){
+	/*Computeportion of the element that is grounded*/ 
+
+	bool               mainlyfloating = true;
+	const IssmPDouble  epsilon= 1.e-15;
+	IssmDouble         phi,s1,s2,area_init,area_grounded;
+	IssmDouble         gl[3];
+	IssmDouble         xyz_bis[3][3];
+
+	/*Recover parameters and values*/
+	GetInputListOnVertices(&gl[0],GLlevelsetEnum);
+	
+	/*Be sure that values are not zero*/
+	if(gl[0]==0) gl[0]=gl[0]+epsilon;
+	if(gl[1]==0) gl[1]=gl[1]+epsilon;
+	if(gl[2]==0) gl[2]=gl[2]+epsilon;
+
+	/*Check that not all nodes are grounded or floating*/
+	if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
+		phi=1;
+	}
+	else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
+		phi=0;
+	}
+	else{
+		/*Figure out if two nodes are floating or grounded*/
+		if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=false;
+
+		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+			/*Coordinates of point 2: same as initial point 2*/
+			xyz_bis[2][0]=*(xyz_list+3*2+0);
+			xyz_bis[2][1]=*(xyz_list+3*2+1);
+			xyz_bis[2][2]=*(xyz_list+3*2+2);
+
+			/*Portion of the segments*/
+			s1=gl[2]/(gl[2]-gl[1]);
+			s2=gl[2]/(gl[2]-gl[0]);
+
+			/*New point 1*/
+			xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0));
+			xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1));
+			xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2));
+
+			/*New point 0*/
+			xyz_bis[0][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0));
+			xyz_bis[0][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1));
+			xyz_bis[0][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2));
+		}
+		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+			/*Coordinates of point 0: same as initial point 2*/
+			xyz_bis[0][0]=*(xyz_list+3*0+0);
+			xyz_bis[0][1]=*(xyz_list+3*0+1);
+			xyz_bis[0][2]=*(xyz_list+3*0+2);
+
+			/*Portion of the segments*/
+			s1=gl[0]/(gl[0]-gl[1]);
+			s2=gl[0]/(gl[0]-gl[2]);
+
+			/*New point 1*/
+			xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0));
+			xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1));
+			xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2));
+
+			/*New point 2*/
+			xyz_bis[2][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0));
+			xyz_bis[2][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1));
+			xyz_bis[2][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*0+2));
+		}
+		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+			/*Coordinates of point 1: same as initial point 2*/
+			xyz_bis[1][0]=*(xyz_list+3*1+0);
+			xyz_bis[1][1]=*(xyz_list+3*1+1);
+			xyz_bis[1][2]=*(xyz_list+3*1+2);
+
+			/*Portion of the segments*/
+			s1=gl[1]/(gl[1]-gl[0]);
+			s2=gl[1]/(gl[1]-gl[2]);
+
+			/*New point 0*/
+			xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0));
+			xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1));
+			xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2));
+
+			/*New point 2*/
+			xyz_bis[2][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0));
+			xyz_bis[2][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1));
+			xyz_bis[2][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2));
+		}
+
+		/*Compute fraction of grounded element*/
+		GetJacobianDeterminant(&area_init, xyz_list,NULL);
+		GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
+		if(mainlyfloating==true) area_grounded=area_init-area_grounded;
+		phi=area_grounded/area_init;
+	}
+
+	if(phi>1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
+
+	return phi;
 }
 /*}}}*/
@@ -2053,4 +2155,5 @@
 				name==WaterfractionEnum||
 				name==FrictionCoefficientEnum ||
+				name==GLlevelsetEnum ||
 				name==GradientEnum ||
 				name==OldGradientEnum  ||
@@ -2181,7 +2284,8 @@
 	int     i,migration_style;
 	bool    floatingelement = false;
+	bool    groundedelement = false;
 	IssmDouble  bed_hydro,yts,gl_melting_rate;
 	IssmDouble  rho_water,rho_ice,density;
-	IssmDouble  melting[NUMVERTICES];
+	IssmDouble  melting[NUMVERTICES],phi[NUMVERTICES];
 	IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
 
@@ -2196,4 +2300,5 @@
 	GetInputListOnVertices(&b[0],BedEnum);
 	GetInputListOnVertices(&ba[0],BathymetryEnum);
+	if(migration_style==SubelementMigrationEnum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
 	rho_water=matpar->GetRhoWater();
 	rho_ice=matpar->GetRhoIce();
@@ -2217,5 +2322,5 @@
 			if (bed_hydro>ba[i]){
 				/*Unground only if the element is connected to the ice shelf*/
-				if(migration_style==AgressiveMigrationEnum){
+				if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum){
 					s[i]=(1-density)*h[i];
 					b[i]=-density*h[i];
@@ -2229,13 +2334,27 @@
 					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
 				}
+				else{
+					if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
+				}
 			}
 		}
 	}
 
-	/*If at least one vertex is now floating, the element is now floating*/
-	for(i=0;i<NUMVERTICES;i++){
-		if(nodes[i]->IsFloating()){
-			floatingelement=true;
-			break;
+	/*SubelementMigrationEnum: if one grounded, all grounded*/
+	if(migration_style==SubelementMigrationEnum){
+		for(i=0;i<NUMVERTICES;i++){
+			if(nodes[i]->IsGrounded()){
+				groundedelement=true;
+				break;
+			}
+		}
+		floatingelement=!groundedelement;
+	}
+	else{
+		for(i=0;i<NUMVERTICES;i++){
+			if(nodes[i]->IsFloating()){
+				floatingelement=true;
+				break;
+			}
 		}
 	}
@@ -2251,4 +2370,11 @@
 	this->inputs->AddInput(new PentaP1Input(BedEnum,&b[0]));
    this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
+
+	/*Recalculate phi*/
+	if(migration_style==SubelementMigrationEnum){
+		for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
+		this->inputs->AddInput(new PentaP1Input(GLlevelsetEnum,&phi[0]));
+		this->InputExtrude(GLlevelsetEnum,ElementEnum);
+	}
 
 	/*Extrude inputs*/
@@ -6767,8 +6893,9 @@
 	/*Intermediaries */
 	int       i,j;
-	int       analysis_type;
+	int       analysis_type,migration_style;
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
 	IssmDouble slope_magnitude,alpha2,Jdet;
+	IssmDouble phi=1.0;
 	IssmDouble slope[3]={0.0,0.0,0.0};
 	IssmDouble MAXSLOPE=.06; // 6 %
@@ -6797,4 +6924,7 @@
 	friction=new Friction("2d",inputs,matpar,analysis_type);
 
+	/*Recover portion of element that is grounded*/
+	if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list[0][0]);
+
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussPenta(0,1,2,2);
@@ -6809,4 +6939,5 @@
 		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); 
 		slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
+		if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
 
 		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h	(revision 14604)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h	(revision 14605)
@@ -181,16 +181,17 @@
 		ElementVector* CreatePVectorPrognostic(void);
 		ElementVector* CreatePVectorSlope(void);
-		void	  GetDofList(int** pdoflist,int approximation_enum,int setenum);
-		void	  GetVertexPidList(int* doflist);
-		void    GetVertexSidList(int* sidlist);
-		void    GetConnectivityList(int* connectivity);
-		int     GetElementType(void);
-		void    GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
-		void    GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
-		void    GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
-		void    GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
-		void	  GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
-		void	  GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solutiong);
-		IssmDouble  GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
+		void	         GetDofList(int** pdoflist,int approximation_enum,int setenum);
+		void	         GetVertexPidList(int* doflist);
+		void           GetVertexSidList(int* sidlist);
+		void           GetConnectivityList(int* connectivity);
+		IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
+		int            GetElementType(void);
+		void           GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
+		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
+		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
+		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
+		void	         GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
+		void	         GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solutiong);
+		IssmDouble     GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
 		void    GetStrainRate3dPattyn(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
 		void    GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
Index: /issm/trunk-jpl/src/m/classes/groundingline.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/groundingline.m	(revision 14604)
+++ /issm/trunk-jpl/src/m/classes/groundingline.m	(revision 14605)
@@ -44,7 +44,4 @@
 					md = checkmessage(md,['bathymetry superior to bed on floating ice!']);
 				end
-				if strcmp(obj.migration,'SubelementMigration') & md.mesh.dimension==3,
-					md = checkmessage(md,['SubelementMigration only implemented in 2d!']);
-				end
 			end
 
Index: /issm/trunk-jpl/src/m/classes/groundingline.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/groundingline.py	(revision 14604)
+++ /issm/trunk-jpl/src/m/classes/groundingline.py	(revision 14605)
@@ -59,7 +59,4 @@
 			if any(md.geometry.bathymetry[pos]-md.geometry.bed[pos]>10**-9):
 				md.checkmessage("bathymetry superior to bed on floating ice!")
-			if strcmp(self.migration,'SubelementMigration'):
-				if md.mesh.dimension==3:
-					md.checkmessage("SubelementMigration only implemented in 2d!")
 
 		return md