Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 21796)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 21797)
@@ -469,8 +469,12 @@
 
 	/*Intermediaries */
+	int         migration_style,point1;
+	bool        mainlyfloating;
+	IssmDouble  fraction1,fraction2;
 	IssmDouble  Jdet,dt;
 	IssmDouble  ms,mb,gmb,fmb,thickness;
 	IssmDouble  gllevelset,phi=1.;
 	IssmDouble* xyz_list = NULL;
+	Gauss*      gauss     = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -490,6 +494,16 @@
 	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);
@@ -501,8 +515,23 @@
 		gmb_input->GetInputValue(&gmb,gauss);
 		fmb_input->GetInputValue(&fmb,gauss);
-		gllevelset_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];
@@ -521,11 +550,7 @@
 
 	/*Intermediaries */
-	int         migration_style,point1;
-	bool        mainlyfloating;
-	IssmDouble  Jdet,fraction1,fraction2;
-	IssmDouble  dt,gllevelset;
+	IssmDouble  Jdet,dt;
 	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*/
@@ -545,16 +570,6 @@
 	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);
@@ -566,22 +581,9 @@
 		gmb_input->GetInputValue(&gmb,gauss);
 		fmb_input->GetInputValue(&fmb,gauss);
-		gllevelset_input->GetInputValue(&gllevelset,gauss);
+		gllevelset_input->GetInputValue(&phi,gauss);
 		thickness_input->GetInputValue(&thickness,gauss);
-		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;
-		}
+
+		if(phi>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];