Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 21793)
+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 21794)
@@ -469,7 +469,8 @@
 
 	/*Intermediaries */
 	IssmDouble  Jdet,dt;
-	IssmDouble  ms,mb,gmb,fmb,thickness,phi;
+	IssmDouble  ms,mb,gmb,fmb,thickness;
+	IssmDouble  gllevelset,phi=1.;
 	IssmDouble* xyz_list = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -484,7 +485,7 @@
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
 	Input* gmb_input           = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
 	Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
-	Input* groundedice_input   = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(groundedice_input);
+	Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
 	Input* ms_input            = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
 	Input* thickness_input     = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
 
@@ -499,7 +500,7 @@
 		ms_input->GetInputValue(&ms,gauss);
 		gmb_input->GetInputValue(&gmb,gauss);
 		fmb_input->GetInputValue(&fmb,gauss);
-		groundedice_input->GetInputValue(&phi,gauss);
+		gllevelset_input->GetInputValue(&phi,gauss);
 		thickness_input->GetInputValue(&thickness,gauss);
 		if(phi>0.) mb=gmb;
 		else mb=fmb;
@@ -519,9 +520,13 @@
 	if(!element->IsIceInElement()) return NULL;
 
 	/*Intermediaries */
-	IssmDouble  Jdet,dt;
-	IssmDouble  ms,mb,gmb,fmb,thickness,phi;
+	int         migration_style,point1;
+	bool        mainlyfloating;
+	IssmDouble  Jdet,fraction1,fraction2;
+	IssmDouble  dt,gllevelset;
+	IssmDouble  ms,mb,gmb,fmb,thickness,phi=1.;
 	IssmDouble* xyz_list = NULL;
+	Gauss*      gauss     = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
@@ -536,11 +541,21 @@
 	Input* gmb_input           = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
 	Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
 	Input* ms_input            = element->GetInput(SmbMassBalanceEnum);          _assert_(ms_input);
-	Input* groundedice_input   = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(groundedice_input);
+	Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(gllevelset_input);
 	Input* thickness_input     = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
 
+	/*Recover portion of element that is grounded*/
+	phi=element->GetGroundedPortion(xyz_list);
+	if(migration_style==SubelementMigration2Enum){
+		gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+	   gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,3);
+	}
+	else{
+		gauss = element->NewGauss(3);
+	}
+
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
@@ -550,10 +565,24 @@
 		ms_input->GetInputValue(&ms,gauss);
 		gmb_input->GetInputValue(&gmb,gauss);
 		fmb_input->GetInputValue(&fmb,gauss);
-		groundedice_input->GetInputValue(&phi,gauss);
+		gllevelset_input->GetInputValue(&gllevelset,gauss);
 		thickness_input->GetInputValue(&thickness,gauss);
-		if(phi>0) mb=gmb;
-		else mb=fmb;
+		if(migration_style==SubelementMigrationEnum){
+			if (phi<0.00000001) mb=gmb;
+			else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
+		}
+		else if(migration_style==SubelementMigration2Enum){
+			if(gllevelset>0.) mb=gmb;
+			else mb=fmb;
+		}
+		else if(migration_style==SubelementMigration3Enum){
+			if (phi<0.00000001) mb=gmb;
+			else mb=gmb;
+		}
+		else{
+			if(gllevelset>0) mb=gmb;
+			else mb=fmb;
+		}
 
 		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
 	}