Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 15385)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 15386)
@@ -75,13 +75,9 @@
 		virtual void   GetVectorFromResults(Vector<IssmDouble>* vector,int id,int enum_in,int interp)=0;
 		virtual void   InputArtificialNoise(int enum_type,IssmDouble min,IssmDouble max)=0;
-		virtual void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
 		virtual int*   GetHorizontalNeighboorSids(void)=0;
 		virtual IssmDouble TimeAdapt()=0;
-		virtual void   MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding)=0;
-		virtual void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
 		virtual void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm)=0;
 		virtual void   Delta18oParameterization(void)=0;
 		virtual void   SmbGradients()=0;
-		virtual int    UpdatePotentialUngrounding(IssmDouble* potential_sheet_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf)=0;
 		virtual void UpdateConstraints(void)=0; 
 		virtual void   ResetCoordinateSystem()=0;
@@ -137,4 +133,15 @@
 		virtual void HydrologyEPLGetActive(Vector<IssmDouble>* active)=0;
 		#endif
+		
+		#ifdef _HAVE_GROUNDINGLINE_
+		virtual void   MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding)=0;
+		virtual void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
+		virtual int    UpdatePotentialUngrounding(IssmDouble* potential_sheet_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf)=0;
+		#endif
+
+		#ifdef _HAVE_WRAPPERS_
+		virtual void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
+		#endif
+
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15385)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15386)
@@ -142,9 +142,4 @@
 
 /*Other*/
-/*FUNCTION Penta::AverageOntoPartition {{{*/
-void  Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
-	_error_("Not supported yet!");
-}
-/*}}}*/
 /*FUNCTION Penta::BedNormal {{{*/
 void Penta::BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]){
@@ -2318,111 +2313,4 @@
 
 }/*}}}*/
-/*FUNCTION Penta::MigrateGroundingLine{{{*/
-void  Penta::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
-
-	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],phi[NUMVERTICES];
-	IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
-
-	if(!IsOnBed()) return;
-
-	/*Recover info at the vertices: */
-	parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
-	parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
-	parameters->FindParam(&yts,ConstantsYtsEnum);
-	GetInputListOnVertices(&h[0],ThicknessEnum);
-	GetInputListOnVertices(&s[0],SurfaceEnum);
-	GetInputListOnVertices(&b[0],BedEnum);
-	GetInputListOnVertices(&ba[0],BathymetryEnum);
-	if(migration_style==SubelementMigrationEnum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
-	rho_water=matpar->GetRhoWater();
-	rho_ice=matpar->GetRhoIce();
-	density=rho_ice/rho_water;
-
-	/*go through vertices, and update inputs, considering them to be PentaVertex type: */
-	for(i=0;i<NUMVERTICES;i++){
-		/*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
-		if(reCast<bool,IssmDouble>(old_floating_ice[nodes[i]->Sid()])){
-			if(b[i]<=ba[i]){ 
-				b[i]=ba[i];
-				s[i]=b[i]+h[i];
-				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
-				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
-			}
-		}
-		/*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
-		/*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
-		else{
-			bed_hydro=-density*h[i];
-			if (bed_hydro>ba[i]){
-				/*Unground only if the element is connected to the ice shelf*/
-				if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum){
-					s[i]=(1-density)*h[i];
-					b[i]=-density*h[i];
-					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
-					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
-				}
-				else if(migration_style==SoftMigrationEnum && reCast<int,IssmDouble>(sheet_ungrounding[nodes[i]->Sid()])){
-					s[i]=(1-density)*h[i];
-					b[i]=-density*h[i];
-					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
-					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
-				}
-				else{
-					if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
-				}
-			}
-		}
-	}
-
-	/*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;
-			}
-		}
-	}
-
-   /*Add basal melting rate if element just ungrounded*/
-	if(!this->IsFloating() && floatingelement==true){
-		for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
-		this->inputs->AddInput(new PentaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
-	} 
-
-	/*Update inputs*/
-	this->inputs->AddInput(new PentaP1Input(SurfaceEnum,&s[0]));
-	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*/
-	this->InputExtrude(SurfaceEnum,ElementEnum);
-	this->InputExtrude(BedEnum,ElementEnum);
-	this->InputExtrude(MaskElementonfloatingiceEnum,ElementEnum);
-	this->InputExtrude(MaskVertexonfloatingiceEnum,NodeEnum);
-	this->InputExtrude(MaskVertexongroundediceEnum,NodeEnum);
-}
-/*}}}*/
 /*FUNCTION Penta::MinEdgeLength{{{*/
 IssmDouble Penta::MinEdgeLength(IssmDouble xyz_list[6][3]){
@@ -2581,32 +2469,4 @@
 	/*clean-up*/
 	delete gauss;
-}
-/*}}}*/
-/*FUNCTION Penta::PotentialUngrounding{{{*/
-void  Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
-
-	int     i;
-	IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
-	IssmDouble  bed_hydro;
-	IssmDouble  rho_water,rho_ice,density;
-
-	/*material parameters: */
-	rho_water=matpar->GetRhoWater();
-	rho_ice=matpar->GetRhoIce();
-	density=rho_ice/rho_water;
-	GetInputListOnVertices(&h[0],ThicknessEnum);
-	GetInputListOnVertices(&ba[0],BathymetryEnum);
-
-	/*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
-	for(i=0;i<NUMVERTICES;i++){
-		/*Find if grounded vertices want to start floating*/
-		if (!nodes[i]->IsFloating()){
-			bed_hydro=-density*h[i];
-			if (bed_hydro>ba[i]){
-				/*Vertex that could potentially unground, flag it*/
-				potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
-			}
-		}
-	}
 }
 /*}}}*/
@@ -3206,24 +3066,4 @@
 }
 /*}}}*/
-/*FUNCTION Penta::UpdatePotentialUngrounding{{{*/
-int Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
-
-	int i;
-	int nflipped=0;
-
-	/*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
-	for(i=0;i<NUMVERTICES;i++){
-		if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
-			vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
-
-			/*If node was not on ice shelf, we flipped*/
-			if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
-				nflipped++;
-			}
-		}
-	}
-	return nflipped;
-}
-/*}}}*/
 /*FUNCTION Penta::UpdateConstraints{{{*/
 void Penta::UpdateConstraints(void){
@@ -9425,5 +9265,5 @@
 	xDelete<int>(doflist);
 }
-/*{{{*/
+/*}}}*/
 /*FUNCTION Penta::HydrologyEPLGetActive {{{*/
 void Penta::HydrologyEPLGetActive(Vector<IssmDouble>* active_vec){
@@ -9516,2 +9356,168 @@
 /*}}}*/
 #endif
+
+#ifdef _HAVE_GROUNDINGLINE_
+/*FUNCTION Penta::MigrateGroundingLine{{{*/
+void  Penta::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
+
+	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],phi[NUMVERTICES];
+	IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
+
+	if(!IsOnBed()) return;
+
+	/*Recover info at the vertices: */
+	parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
+	parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
+	parameters->FindParam(&yts,ConstantsYtsEnum);
+	GetInputListOnVertices(&h[0],ThicknessEnum);
+	GetInputListOnVertices(&s[0],SurfaceEnum);
+	GetInputListOnVertices(&b[0],BedEnum);
+	GetInputListOnVertices(&ba[0],BathymetryEnum);
+	if(migration_style==SubelementMigrationEnum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
+	rho_water=matpar->GetRhoWater();
+	rho_ice=matpar->GetRhoIce();
+	density=rho_ice/rho_water;
+
+	/*go through vertices, and update inputs, considering them to be PentaVertex type: */
+	for(i=0;i<NUMVERTICES;i++){
+		/*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
+		if(reCast<bool,IssmDouble>(old_floating_ice[nodes[i]->Sid()])){
+			if(b[i]<=ba[i]){ 
+				b[i]=ba[i];
+				s[i]=b[i]+h[i];
+				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
+				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
+			}
+		}
+		/*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
+		/*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
+		else{
+			bed_hydro=-density*h[i];
+			if (bed_hydro>ba[i]){
+				/*Unground only if the element is connected to the ice shelf*/
+				if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum){
+					s[i]=(1-density)*h[i];
+					b[i]=-density*h[i];
+					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
+					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
+				}
+				else if(migration_style==SoftMigrationEnum && reCast<int,IssmDouble>(sheet_ungrounding[nodes[i]->Sid()])){
+					s[i]=(1-density)*h[i];
+					b[i]=-density*h[i];
+					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
+					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
+				}
+				else{
+					if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
+				}
+			}
+		}
+	}
+
+	/*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;
+			}
+		}
+	}
+
+   /*Add basal melting rate if element just ungrounded*/
+	if(!this->IsFloating() && floatingelement==true){
+		for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
+		this->inputs->AddInput(new PentaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
+	} 
+
+	/*Update inputs*/
+	this->inputs->AddInput(new PentaP1Input(SurfaceEnum,&s[0]));
+	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*/
+	this->InputExtrude(SurfaceEnum,ElementEnum);
+	this->InputExtrude(BedEnum,ElementEnum);
+	this->InputExtrude(MaskElementonfloatingiceEnum,ElementEnum);
+	this->InputExtrude(MaskVertexonfloatingiceEnum,NodeEnum);
+	this->InputExtrude(MaskVertexongroundediceEnum,NodeEnum);
+}
+/*}}}*/
+#endif
+
+#ifdef _HAVE_WRAPPERS_
+/*FUNCTION Penta::PotentialUngrounding{{{*/
+void  Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
+
+	int     i;
+	IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
+	IssmDouble  bed_hydro;
+	IssmDouble  rho_water,rho_ice,density;
+
+	/*material parameters: */
+	rho_water=matpar->GetRhoWater();
+	rho_ice=matpar->GetRhoIce();
+	density=rho_ice/rho_water;
+	GetInputListOnVertices(&h[0],ThicknessEnum);
+	GetInputListOnVertices(&ba[0],BathymetryEnum);
+
+	/*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
+	for(i=0;i<NUMVERTICES;i++){
+		/*Find if grounded vertices want to start floating*/
+		if (!nodes[i]->IsFloating()){
+			bed_hydro=-density*h[i];
+			if (bed_hydro>ba[i]){
+				/*Vertex that could potentially unground, flag it*/
+				potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
+			}
+		}
+	}
+}
+/*}}}*/
+/*FUNCTION Penta::AverageOntoPartition {{{*/
+void  Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
+	_error_("Not supported yet!");
+}
+/*}}}*/
+/*FUNCTION Penta::UpdatePotentialUngrounding{{{*/
+int Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
+
+	int i;
+	int nflipped=0;
+
+	/*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
+	for(i=0;i<NUMVERTICES;i++){
+		if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
+			vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
+
+			/*If node was not on ice shelf, we flipped*/
+			if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
+				nflipped++;
+			}
+		}
+	}
+	return nflipped;
+}
+/*}}}*/
+#endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 15385)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 15386)
@@ -75,5 +75,4 @@
 		/*}}}*/
 		/*Element virtual functions definitions: {{{*/
-		void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
 		void   BasalFrictionCreateInput(void);
 		void   ComputeBasalStress(Vector<IssmDouble>* sigma_b);
@@ -106,6 +105,4 @@
 
 		void   InputToResult(int enum_type,int step,IssmDouble time);
-		void   MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding);
-		void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
 		void   RequestedOutput(int output_enum,int step,IssmDouble time);
 		void   ListResultsInfo(int** results_enums,int** results_size,IssmDouble** results_times,int** results_steps,int* num_results);
@@ -117,5 +114,4 @@
 		IssmDouble SurfaceArea(void);
 		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
-		int    UpdatePotentialUngrounding(IssmDouble* potential_sheet_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
 		int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
 		IssmDouble TimeAdapt();
@@ -173,4 +169,15 @@
 		void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
 		#endif
+		
+		#ifdef _HAVE_GROUNDINGLINE_
+		void   MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding);
+		void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
+		int    UpdatePotentialUngrounding(IssmDouble* potential_sheet_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
+		#endif
+
+		#ifdef _HAVE_WRAPPERS_
+		void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
+		#endif
+
 		/*}}}*/
 		/*Penta specific routines:{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15385)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15386)
@@ -126,43 +126,4 @@
 
 /*Other*/
-/*FUNCTION Tria::AverageOntoPartition {{{*/
-void  Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
-
-	bool       already = false;
-	int        i,j;
-	int        partition[NUMVERTICES];
-	int        offsetsid[NUMVERTICES];
-	int        offsetdof[NUMVERTICES];
-	IssmDouble area;
-	IssmDouble mean;
-
-	/*First, get the area: */
-	area=this->GetArea();
-
-	/*Figure out the average for this element: */
-	this->GetVertexSidList(&offsetsid[0]);
-	this->GetVertexPidList(&offsetdof[0]);
-	mean=0;
-	for(i=0;i<NUMVERTICES;i++){
-		partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
-		mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
-	}
-
-	/*Add contribution: */
-	for(i=0;i<NUMVERTICES;i++){
-		already=false;
-		for(j=0;j<i;j++){
-			if (partition[i]==partition[j]){
-				already=true;
-				break;
-			}
-		}
-		if(!already){
-			partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
-			partition_areas->SetValue(partition[i],area,ADD_VAL);
-		};
-	}
-}
-/*}}}*/
 /*FUNCTION Tria::SetwiseNodeConnectivity{{{*/
 void Tria::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int set1_enum,int set2_enum){
@@ -282,254 +243,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixMelting {{{*/
-ElementMatrix* Tria::CreateKMatrixMelting(void){
-
-	/*Constants*/
-	const int  numdof=NUMVERTICES*NDOF1;
-
-	/*Intermediaries */
-	IssmDouble heatcapacity,latentheat;
-	IssmDouble Jdet,D_scalar;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble basis[3];
-	GaussTria *gauss=NULL;
-
-	/*Initialize Element matrix*/
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
-
-	/*Retrieve all inputs and parameters*/
-	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	latentheat=matpar->GetLatentHeat();
-	heatcapacity=matpar->GetHeatCapacity();
-
-	/* Start looping on the number of gauss  (nodes on the bedrock) */
-	gauss=new GaussTria(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		GetNodalFunctions(&basis[0],gauss);
-		GetJacobianDeterminant(&Jdet, &xyz_list[0][0], gauss);
-
-		D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
-
-		TripleMultiply(&basis[0],numdof,1,0,
-					&D_scalar,1,1,0,
-					&basis[0],1,numdof,0,
-					&Ke->values[0],1);
-	}
-
-	/*Clean up and return*/
-	delete gauss;
-	return Ke;
-}
-/*}}}*/
-/*FUNCTION Tria::CreateKMatrixPrognostic {{{*/
-ElementMatrix* Tria::CreateKMatrixPrognostic(void){
-
-	switch(GetElementType()){
-		case P1Enum:
-			return CreateKMatrixPrognostic_CG();
-		case P1DGEnum:
-			return CreateKMatrixPrognostic_DG();
-		default:
-			_error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION Tria::CreateKMatrixPrognostic_CG {{{*/
-ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){
-
-	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	/*Intermediaries */
-	int        stabilization;
-	int        dim;
-	IssmDouble Jdettria,DL_scalar,dt,h;
-	IssmDouble vel,vx,vy,dvxdx,dvydy;
-	IssmDouble dvx[2],dvy[2];
-	IssmDouble v_gauss[2]={0.0};
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble basis[NUMVERTICES];
-	IssmDouble B[2][NUMVERTICES];
-	IssmDouble Bprime[2][NUMVERTICES];
-	IssmDouble K[2][2]                        ={0.0};
-	IssmDouble KDL[2][2]                      ={0.0};
-	IssmDouble DL[2][2]                        ={0.0};
-	IssmDouble DLprime[2][2]                   ={0.0};
-
-	/*Initialize Element matrix*/
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
-
-	/*Retrieve all inputs and parameters*/
-	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-	this->parameters->FindParam(&dim,MeshDimensionEnum);
-	this->parameters->FindParam(&stabilization,PrognosticStabilizationEnum);
-	Input* vxaverage_input=NULL;
-	Input* vyaverage_input=NULL;
-	if(dim==2){
-		vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
-		vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
-	}
-	else{
-		vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
-		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
-	}
-	h=sqrt(2*this->GetArea());
-
-	/* Start  looping on the number of gaussian points: */
-	GaussTria *gauss=new GaussTria(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetNodalFunctions(&basis[0],gauss);
-
-		vxaverage_input->GetInputValue(&vx,gauss);
-		vyaverage_input->GetInputValue(&vy,gauss);
-		vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
-		vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
-
-		DL_scalar=gauss->weight*Jdettria;
-
-		TripleMultiply(&basis[0],1,numdof,1,
-					&DL_scalar,1,1,0,
-					&basis[0],1,numdof,0,
-					&Ke->values[0],1);
-
-		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
-		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
-
-		dvxdx=dvx[0];
-		dvydy=dvy[1];
-		DL_scalar=dt*gauss->weight*Jdettria;
-
-		DL[0][0]=DL_scalar*dvxdx;
-		DL[1][1]=DL_scalar*dvydy;
-		DLprime[0][0]=DL_scalar*vx;
-		DLprime[1][1]=DL_scalar*vy;
-
-		TripleMultiply( &B[0][0],2,numdof,1,
-					&DL[0][0],2,2,0,
-					&B[0][0],2,numdof,0,
-					&Ke->values[0],1);
-
-		TripleMultiply( &B[0][0],2,numdof,1,
-					&DLprime[0][0],2,2,0,
-					&Bprime[0][0],2,numdof,0,
-					&Ke->values[0],1);
-
-		if(stabilization==2){
-			/*Streamline upwinding*/
-			vel=sqrt(vx*vx+vy*vy)+1.e-8;
-			K[0][0]=h/(2*vel)*vx*vx;
-			K[1][0]=h/(2*vel)*vy*vx;
-			K[0][1]=h/(2*vel)*vx*vy;
-			K[1][1]=h/(2*vel)*vy*vy;
-		}
-		else if(stabilization==1){
-			/*MacAyeal*/
-			vxaverage_input->GetInputAverage(&vx);
-			vyaverage_input->GetInputAverage(&vy);
-			K[0][0]=h/2.0*fabs(vx);
-			K[0][1]=0.;
-			K[1][0]=0.;
-			K[1][1]=h/2.0*fabs(vy);
-		}
-		if(stabilization==1 || stabilization==2){
-			KDL[0][0]=DL_scalar*K[0][0];
-			KDL[1][0]=DL_scalar*K[1][0];
-			KDL[0][1]=DL_scalar*K[0][1];
-			KDL[1][1]=DL_scalar*K[1][1];
-			TripleMultiply( &Bprime[0][0],2,numdof,1,
-						&KDL[0][0],2,2,0,
-						&Bprime[0][0],2,numdof,0,
-						&Ke->values[0],1);
-		}
-	}
-
-	/*Clean up and return*/
-	delete gauss;
-	return Ke;
-}
-/*}}}*/
-/*FUNCTION Tria::CreateKMatrixPrognostic_DG {{{*/
-ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){
-
-	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	/*Intermediaries */
-	int        dim;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble Jdettria,dt,vx,vy;
-	IssmDouble basis[NUMVERTICES];
-	IssmDouble B[2][NUMVERTICES];
-	IssmDouble Bprime[2][NUMVERTICES];
-	IssmDouble DL[2][2]={0.0};
-	IssmDouble DLprime[2][2]={0.0};
-	IssmDouble DL_scalar;
-	GaussTria  *gauss=NULL;
-
-	/*Initialize Element matrix*/
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
-
-	/*Retrieve all inputs and parameters*/
-	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-	this->parameters->FindParam(&dim,MeshDimensionEnum);
-	Input* vxaverage_input=NULL;
-	Input* vyaverage_input=NULL;
-	if(dim==2){
-		vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
-		vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
-	}
-	else{
-		vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
-		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
-	}
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		vxaverage_input->GetInputValue(&vx,gauss);
-		vyaverage_input->GetInputValue(&vy,gauss);
-
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetNodalFunctions(&basis[0],gauss);
-
-		DL_scalar=gauss->weight*Jdettria;
-
-		TripleMultiply(&basis[0],1,numdof,1,
-					&DL_scalar,1,1,0,
-					&basis[0],1,numdof,0,
-					&Ke->values[0],1);
-
-		/*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
-		GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
-		GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
-
-		DL_scalar=-dt*gauss->weight*Jdettria;
-
-		DLprime[0][0]=DL_scalar*vx;
-		DLprime[1][1]=DL_scalar*vy;
-
-		TripleMultiply( &B[0][0],2,numdof,1,
-					&DLprime[0][0],2,2,0,
-					&Bprime[0][0],2,numdof,0,
-					&Ke->values[0],1);
-	}
-
-	/*Clean up and return*/
-	delete gauss;
-	return Ke;
-}
-/*}}}*/
 /*FUNCTION Tria::CreateMassMatrix {{{*/
 ElementMatrix* Tria::CreateMassMatrix(void){
@@ -641,111 +352,4 @@
 		delete pe;
 	}
-}
-/*}}}*/
-/*FUNCTION Tria::CreatePVectorPrognostic{{{*/
-ElementVector* Tria::CreatePVectorPrognostic(void){
-
-	switch(GetElementType()){
-		case P1Enum:
-			return CreatePVectorPrognostic_CG();
-		case P1DGEnum:
-			return CreatePVectorPrognostic_DG();
-		default:
-			_error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
-	}
-}
-/*}}}*/
-/*FUNCTION Tria::CreatePVectorPrognostic_CG {{{*/
-ElementVector* Tria::CreatePVectorPrognostic_CG(void){
-
-	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	/*Intermediaries */
-	IssmDouble Jdettria,dt;
-	IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble basis[NUMVERTICES];
-	GaussTria* gauss=NULL;
-
-	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
-
-	/*Retrieve all inputs and parameters*/
-	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-	Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
-	Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
-	Input* basal_melting_correction_input=inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
-	Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
-
-	/*Initialize basal_melting_correction_g to 0, do not forget!:*/
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetNodalFunctions(&basis[0],gauss);
-
-		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
-		basal_melting_input->GetInputValue(&basal_melting_g,gauss);
-		thickness_input->GetInputValue(&thickness_g,gauss);
-		if(basal_melting_correction_input)
-		 basal_melting_correction_input->GetInputValue(&basal_melting_correction_g,gauss);
-		else
-		 basal_melting_correction_g=0.;
-
-		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
-	}
-
-	/*Clean up and return*/
-	delete gauss;
-	return pe;
-}
-/*}}}*/
-/*FUNCTION Tria::CreatePVectorPrognostic_DG {{{*/
-ElementVector* Tria::CreatePVectorPrognostic_DG(void){
-
-	/*Constants*/
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	/*Intermediaries */
-	IssmDouble Jdettria,dt;
-	IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble basis[NUMVERTICES];
-	GaussTria* gauss=NULL;
-
-	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
-
-	/*Retrieve all inputs and parameters*/
-	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
-	Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
-	Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetNodalFunctions(&basis[0],gauss);
-
-		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
-		basal_melting_input->GetInputValue(&basal_melting_g,gauss);
-		thickness_input->GetInputValue(&thickness_g,gauss);
-
-		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
-	}
-
-	/*Clean up and return*/
-	delete gauss;
-	return pe;
 }
 /*}}}*/
@@ -2272,102 +1876,4 @@
 
 }/*}}}*/
-/*FUNCTION Tria::MigrateGroundingLine{{{*/
-void  Tria::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
-
-	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],phi[NUMVERTICES];;
-	IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
-
-	/*Recover info at the vertices: */
-	parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
-	parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
-	parameters->FindParam(&yts,ConstantsYtsEnum);
-	GetInputListOnVertices(&h[0],ThicknessEnum);
-	GetInputListOnVertices(&s[0],SurfaceEnum);
-	GetInputListOnVertices(&b[0],BedEnum);
-	GetInputListOnVertices(&ba[0],BathymetryEnum);
-	if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
-	rho_water=matpar->GetRhoWater();
-	rho_ice=matpar->GetRhoIce();
-	density=rho_ice/rho_water;
-
-	/*go through vertices, and update inputs, considering them to be TriaVertex type: */
-	for(i=0;i<NUMVERTICES;i++){
-		/*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
-		if(reCast<bool>(old_floating_ice[nodes[i]->Sid()])){
-			if(b[i]<=ba[i]){ 
-				b[i]=ba[i];
-				s[i]=b[i]+h[i];
-				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
-				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
-			}
-		}
-		/*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
-		/*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
-		else{
-			bed_hydro=-density*h[i];
-			if (bed_hydro>ba[i]){
-				/*Unground only if the element is connected to the ice shelf*/
-				if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
-					s[i]=(1-density)*h[i];
-					b[i]=-density*h[i];
-					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
-					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
-				}
-				else if(migration_style==SoftMigrationEnum && reCast<bool>(sheet_ungrounding[nodes[i]->Sid()])){
-					s[i]=(1-density)*h[i];
-					b[i]=-density*h[i];
-					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
-					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
-				}
-				else{
-					if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
-				}
-			}
-		}
-	}
-
-	/*SubelementMigrationEnum: if one grounded, all grounded*/
-	if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
-		for(i=0;i<NUMVERTICES;i++){
-			if(nodes[i]->IsGrounded()){
-				groundedelement=true;
-				break;
-			}
-		}
-		floatingelement=!groundedelement;
-	}
-	else{
-		/*Otherwise: if one floating, all floating*/
-		for(i=0;i<NUMVERTICES;i++){
-			if(nodes[i]->IsFloating()){
-				floatingelement=true;
-				break;
-			}
-		}
-	}
-
-   /*Add basal melting rate if element just ungrounded*/
-	if(!this->IsFloating() && floatingelement==true){
-		for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
-		this->inputs->AddInput(new TriaInput(BasalforcingsMeltingRateEnum,&melting[0],P1Enum));
-	} 
-
-	/*Update inputs*/
-   this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
-	this->inputs->AddInput(new TriaInput(SurfaceEnum,&s[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(BedEnum,&b[0],P1Enum));
-
-	/*Recalculate phi*/
-	if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
-		for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
-		this->inputs->AddInput(new TriaInput(GLlevelsetEnum,&phi[0],P1Enum));
-	}
-}
-/*}}}*/
 /*FUNCTION Tria::NodalValue {{{*/
 int    Tria::NodalValue(IssmDouble* pvalue, int index, int natureofdataenum){
@@ -2452,32 +1958,4 @@
 	*pnumvertices=NUMVERTICES;
 	*pnumnodes=numnodes;
-}
-/*}}}*/
-/*FUNCTION Tria::PotentialUngrounding{{{*/
-void  Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
-
-	int     i;
-	IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
-	IssmDouble  bed_hydro;
-	IssmDouble  rho_water,rho_ice,density;
-
-	/*material parameters: */
-	rho_water=matpar->GetRhoWater();
-	rho_ice=matpar->GetRhoIce();
-	density=rho_ice/rho_water;
-	GetInputListOnVertices(&h[0],ThicknessEnum);
-	GetInputListOnVertices(&ba[0],BathymetryEnum);
-
-	/*go through vertices, and figure out which ones are grounded and want to unground: */
-	for(i=0;i<NUMVERTICES;i++){
-		/*Find if grounded vertices want to start floating*/
-		if (!nodes[i]->IsFloating()){
-			bed_hydro=-density*h[i];
-			if (bed_hydro>ba[i]){
-				/*Vertex that could potentially unground, flag it*/
-				potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
-			}
-		}
-	}
 }
 /*}}}*/
@@ -2834,24 +2312,4 @@
 	if(IsOnWater()) return;
 
-}
-/*}}}*/
-/*FUNCTION Tria::UpdatePotentialUngrounding{{{*/
-int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
-
-	int i;
-	int nflipped=0;
-
-	/*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
-	for(i=0;i<3;i++){
-		if (reCast<bool>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
-			vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
-
-			/*If node was not on ice shelf, we flipped*/
-			if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
-				nflipped++;
-			}
-		}
-	}
-	return nflipped;
 }
 /*}}}*/
@@ -5852,4 +5310,50 @@
 #endif
 
+#ifdef _HAVE_THERMAL_
+/*FUNCTION Tria::CreateKMatrixMelting {{{*/
+ElementMatrix* Tria::CreateKMatrixMelting(void){
+
+	/*Constants*/
+	const int  numdof=NUMVERTICES*NDOF1;
+
+	/*Intermediaries */
+	IssmDouble heatcapacity,latentheat;
+	IssmDouble Jdet,D_scalar;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble basis[3];
+	GaussTria *gauss=NULL;
+
+	/*Initialize Element matrix*/
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	latentheat=matpar->GetLatentHeat();
+	heatcapacity=matpar->GetHeatCapacity();
+
+	/* Start looping on the number of gauss  (nodes on the bedrock) */
+	gauss=new GaussTria(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetNodalFunctions(&basis[0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0], gauss);
+
+		D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
+
+		TripleMultiply(&basis[0],numdof,1,0,
+					&D_scalar,1,1,0,
+					&basis[0],1,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
+#endif
+
 #ifdef _HAVE_HYDROLOGY_
 /*FUNCTION Tria::AllActive{{{*/
@@ -6632,4 +6136,321 @@
 #endif
 
+#ifdef _HAVE_PROGNOSTIC_
+/*FUNCTION Tria::CreateKMatrixPrognostic {{{*/
+ElementMatrix* Tria::CreateKMatrixPrognostic(void){
+
+	switch(GetElementType()){
+		case P1Enum:
+			return CreateKMatrixPrognostic_CG();
+		case P1DGEnum:
+			return CreateKMatrixPrognostic_DG();
+		default:
+			_error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
+	}
+}
+/*}}}*/
+/*FUNCTION Tria::CreateKMatrixPrognostic_CG {{{*/
+ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){
+
+	/*Constants*/
+	const int    numdof=NDOF1*NUMVERTICES;
+
+	/*Intermediaries */
+	int        stabilization;
+	int        dim;
+	IssmDouble Jdettria,DL_scalar,dt,h;
+	IssmDouble vel,vx,vy,dvxdx,dvydy;
+	IssmDouble dvx[2],dvy[2];
+	IssmDouble v_gauss[2]={0.0};
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble basis[NUMVERTICES];
+	IssmDouble B[2][NUMVERTICES];
+	IssmDouble Bprime[2][NUMVERTICES];
+	IssmDouble K[2][2]                        ={0.0};
+	IssmDouble KDL[2][2]                      ={0.0};
+	IssmDouble DL[2][2]                        ={0.0};
+	IssmDouble DLprime[2][2]                   ={0.0};
+
+	/*Initialize Element matrix*/
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	this->parameters->FindParam(&dim,MeshDimensionEnum);
+	this->parameters->FindParam(&stabilization,PrognosticStabilizationEnum);
+	Input* vxaverage_input=NULL;
+	Input* vyaverage_input=NULL;
+	if(dim==2){
+		vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
+		vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
+	}
+	else{
+		vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
+		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+	}
+	h=sqrt(2*this->GetArea());
+
+	/* Start  looping on the number of gaussian points: */
+	GaussTria *gauss=new GaussTria(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetNodalFunctions(&basis[0],gauss);
+
+		vxaverage_input->GetInputValue(&vx,gauss);
+		vyaverage_input->GetInputValue(&vy,gauss);
+		vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
+		vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
+
+		DL_scalar=gauss->weight*Jdettria;
+
+		TripleMultiply(&basis[0],1,numdof,1,
+					&DL_scalar,1,1,0,
+					&basis[0],1,numdof,0,
+					&Ke->values[0],1);
+
+		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
+		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
+
+		dvxdx=dvx[0];
+		dvydy=dvy[1];
+		DL_scalar=dt*gauss->weight*Jdettria;
+
+		DL[0][0]=DL_scalar*dvxdx;
+		DL[1][1]=DL_scalar*dvydy;
+		DLprime[0][0]=DL_scalar*vx;
+		DLprime[1][1]=DL_scalar*vy;
+
+		TripleMultiply( &B[0][0],2,numdof,1,
+					&DL[0][0],2,2,0,
+					&B[0][0],2,numdof,0,
+					&Ke->values[0],1);
+
+		TripleMultiply( &B[0][0],2,numdof,1,
+					&DLprime[0][0],2,2,0,
+					&Bprime[0][0],2,numdof,0,
+					&Ke->values[0],1);
+
+		if(stabilization==2){
+			/*Streamline upwinding*/
+			vel=sqrt(vx*vx+vy*vy)+1.e-8;
+			K[0][0]=h/(2*vel)*vx*vx;
+			K[1][0]=h/(2*vel)*vy*vx;
+			K[0][1]=h/(2*vel)*vx*vy;
+			K[1][1]=h/(2*vel)*vy*vy;
+		}
+		else if(stabilization==1){
+			/*MacAyeal*/
+			vxaverage_input->GetInputAverage(&vx);
+			vyaverage_input->GetInputAverage(&vy);
+			K[0][0]=h/2.0*fabs(vx);
+			K[0][1]=0.;
+			K[1][0]=0.;
+			K[1][1]=h/2.0*fabs(vy);
+		}
+		if(stabilization==1 || stabilization==2){
+			KDL[0][0]=DL_scalar*K[0][0];
+			KDL[1][0]=DL_scalar*K[1][0];
+			KDL[0][1]=DL_scalar*K[0][1];
+			KDL[1][1]=DL_scalar*K[1][1];
+			TripleMultiply( &Bprime[0][0],2,numdof,1,
+						&KDL[0][0],2,2,0,
+						&Bprime[0][0],2,numdof,0,
+						&Ke->values[0],1);
+		}
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Tria::CreateKMatrixPrognostic_DG {{{*/
+ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){
+
+	/*Constants*/
+	const int    numdof=NDOF1*NUMVERTICES;
+
+	/*Intermediaries */
+	int        dim;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble Jdettria,dt,vx,vy;
+	IssmDouble basis[NUMVERTICES];
+	IssmDouble B[2][NUMVERTICES];
+	IssmDouble Bprime[2][NUMVERTICES];
+	IssmDouble DL[2][2]={0.0};
+	IssmDouble DLprime[2][2]={0.0};
+	IssmDouble DL_scalar;
+	GaussTria  *gauss=NULL;
+
+	/*Initialize Element matrix*/
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	this->parameters->FindParam(&dim,MeshDimensionEnum);
+	Input* vxaverage_input=NULL;
+	Input* vyaverage_input=NULL;
+	if(dim==2){
+		vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
+		vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
+	}
+	else{
+		vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
+		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		vxaverage_input->GetInputValue(&vx,gauss);
+		vyaverage_input->GetInputValue(&vy,gauss);
+
+		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetNodalFunctions(&basis[0],gauss);
+
+		DL_scalar=gauss->weight*Jdettria;
+
+		TripleMultiply(&basis[0],1,numdof,1,
+					&DL_scalar,1,1,0,
+					&basis[0],1,numdof,0,
+					&Ke->values[0],1);
+
+		/*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
+		GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
+		GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
+
+		DL_scalar=-dt*gauss->weight*Jdettria;
+
+		DLprime[0][0]=DL_scalar*vx;
+		DLprime[1][1]=DL_scalar*vy;
+
+		TripleMultiply( &B[0][0],2,numdof,1,
+					&DLprime[0][0],2,2,0,
+					&Bprime[0][0],2,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorPrognostic{{{*/
+ElementVector* Tria::CreatePVectorPrognostic(void){
+
+	switch(GetElementType()){
+		case P1Enum:
+			return CreatePVectorPrognostic_CG();
+		case P1DGEnum:
+			return CreatePVectorPrognostic_DG();
+		default:
+			_error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
+	}
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorPrognostic_CG {{{*/
+ElementVector* Tria::CreatePVectorPrognostic_CG(void){
+
+	/*Constants*/
+	const int    numdof=NDOF1*NUMVERTICES;
+
+	/*Intermediaries */
+	IssmDouble Jdettria,dt;
+	IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble basis[NUMVERTICES];
+	GaussTria* gauss=NULL;
+
+	/*Initialize Element vector*/
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
+	Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
+	Input* basal_melting_correction_input=inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
+
+	/*Initialize basal_melting_correction_g to 0, do not forget!:*/
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetNodalFunctions(&basis[0],gauss);
+
+		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
+		basal_melting_input->GetInputValue(&basal_melting_g,gauss);
+		thickness_input->GetInputValue(&thickness_g,gauss);
+		if(basal_melting_correction_input)
+		 basal_melting_correction_input->GetInputValue(&basal_melting_correction_g,gauss);
+		else
+		 basal_melting_correction_g=0.;
+
+		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return pe;
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorPrognostic_DG {{{*/
+ElementVector* Tria::CreatePVectorPrognostic_DG(void){
+
+	/*Constants*/
+	const int    numdof=NDOF1*NUMVERTICES;
+
+	/*Intermediaries */
+	IssmDouble Jdettria,dt;
+	IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble basis[NUMVERTICES];
+	GaussTria* gauss=NULL;
+
+	/*Initialize Element vector*/
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
+
+	/*Retrieve all inputs and parameters*/
+	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
+	Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
+	Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetNodalFunctions(&basis[0],gauss);
+
+		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
+		basal_melting_input->GetInputValue(&basal_melting_g,gauss);
+		thickness_input->GetInputValue(&thickness_g,gauss);
+
+		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return pe;
+}
+/*}}}*/
+#endif
+
 #ifdef _HAVE_DAKOTA_
 /*FUNCTION Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);{{{*/
@@ -7058,2 +6879,193 @@
 /*}}}*/
 #endif
+
+#ifdef _HAVE_GROUNDINGLINE_
+/*FUNCTION Tria::MigrateGroundingLine{{{*/
+void  Tria::MigrateGroundingLine(IssmDouble* old_floating_ice,IssmDouble* sheet_ungrounding){
+
+	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],phi[NUMVERTICES];;
+	IssmDouble  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
+
+	/*Recover info at the vertices: */
+	parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
+	parameters->FindParam(&gl_melting_rate,GroundinglineMeltingRateEnum);
+	parameters->FindParam(&yts,ConstantsYtsEnum);
+	GetInputListOnVertices(&h[0],ThicknessEnum);
+	GetInputListOnVertices(&s[0],SurfaceEnum);
+	GetInputListOnVertices(&b[0],BedEnum);
+	GetInputListOnVertices(&ba[0],BathymetryEnum);
+	if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum) GetInputListOnVertices(&phi[0],GLlevelsetEnum);
+	rho_water=matpar->GetRhoWater();
+	rho_ice=matpar->GetRhoIce();
+	density=rho_ice/rho_water;
+
+	/*go through vertices, and update inputs, considering them to be TriaVertex type: */
+	for(i=0;i<NUMVERTICES;i++){
+		/*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
+		if(reCast<bool>(old_floating_ice[nodes[i]->Sid()])){
+			if(b[i]<=ba[i]){ 
+				b[i]=ba[i];
+				s[i]=b[i]+h[i];
+				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
+				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
+			}
+		}
+		/*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
+		/*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
+		else{
+			bed_hydro=-density*h[i];
+			if (bed_hydro>ba[i]){
+				/*Unground only if the element is connected to the ice shelf*/
+				if(migration_style==AgressiveMigrationEnum || migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
+					s[i]=(1-density)*h[i];
+					b[i]=-density*h[i];
+					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
+					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
+				}
+				else if(migration_style==SoftMigrationEnum && reCast<bool>(sheet_ungrounding[nodes[i]->Sid()])){
+					s[i]=(1-density)*h[i];
+					b[i]=-density*h[i];
+					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
+					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
+				}
+				else{
+					if(migration_style!=SoftMigrationEnum) _error_("Error: migration should be Aggressive, Soft or Subelement");
+				}
+			}
+		}
+	}
+
+	/*SubelementMigrationEnum: if one grounded, all grounded*/
+	if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
+		for(i=0;i<NUMVERTICES;i++){
+			if(nodes[i]->IsGrounded()){
+				groundedelement=true;
+				break;
+			}
+		}
+		floatingelement=!groundedelement;
+	}
+	else{
+		/*Otherwise: if one floating, all floating*/
+		for(i=0;i<NUMVERTICES;i++){
+			if(nodes[i]->IsFloating()){
+				floatingelement=true;
+				break;
+			}
+		}
+	}
+
+   /*Add basal melting rate if element just ungrounded*/
+	if(!this->IsFloating() && floatingelement==true){
+		for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
+		this->inputs->AddInput(new TriaInput(BasalforcingsMeltingRateEnum,&melting[0],P1Enum));
+	} 
+
+	/*Update inputs*/
+   this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
+	this->inputs->AddInput(new TriaInput(SurfaceEnum,&s[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(BedEnum,&b[0],P1Enum));
+
+	/*Recalculate phi*/
+	if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
+		for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
+		this->inputs->AddInput(new TriaInput(GLlevelsetEnum,&phi[0],P1Enum));
+	}
+}
+/*}}}*/
+/*FUNCTION Tria::PotentialUngrounding{{{*/
+void  Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){
+
+	int     i;
+	IssmDouble  h[NUMVERTICES],ba[NUMVERTICES];
+	IssmDouble  bed_hydro;
+	IssmDouble  rho_water,rho_ice,density;
+
+	/*material parameters: */
+	rho_water=matpar->GetRhoWater();
+	rho_ice=matpar->GetRhoIce();
+	density=rho_ice/rho_water;
+	GetInputListOnVertices(&h[0],ThicknessEnum);
+	GetInputListOnVertices(&ba[0],BathymetryEnum);
+
+	/*go through vertices, and figure out which ones are grounded and want to unground: */
+	for(i=0;i<NUMVERTICES;i++){
+		/*Find if grounded vertices want to start floating*/
+		if (!nodes[i]->IsFloating()){
+			bed_hydro=-density*h[i];
+			if (bed_hydro>ba[i]){
+				/*Vertex that could potentially unground, flag it*/
+				potential_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
+			}
+		}
+	}
+}
+/*}}}*/
+/*FUNCTION Tria::UpdatePotentialUngrounding{{{*/
+int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){
+
+	int i;
+	int nflipped=0;
+
+	/*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
+	for(i=0;i<3;i++){
+		if (reCast<bool>(vertices_potentially_ungrounding[nodes[i]->Sid()])){
+			vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
+
+			/*If node was not on ice shelf, we flipped*/
+			if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
+				nflipped++;
+			}
+		}
+	}
+	return nflipped;
+}
+/*}}}*/
+#endif
+
+#ifdef _HAVE_WRAPPERS_
+/*FUNCTION Tria::AverageOntoPartition {{{*/
+void  Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
+
+	bool       already = false;
+	int        i,j;
+	int        partition[NUMVERTICES];
+	int        offsetsid[NUMVERTICES];
+	int        offsetdof[NUMVERTICES];
+	IssmDouble area;
+	IssmDouble mean;
+
+	/*First, get the area: */
+	area=this->GetArea();
+
+	/*Figure out the average for this element: */
+	this->GetVertexSidList(&offsetsid[0]);
+	this->GetVertexPidList(&offsetdof[0]);
+	mean=0;
+	for(i=0;i<NUMVERTICES;i++){
+		partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
+		mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
+	}
+
+	/*Add contribution: */
+	for(i=0;i<NUMVERTICES;i++){
+		already=false;
+		for(j=0;j<i;j++){
+			if (partition[i]==partition[j]){
+				already=true;
+				break;
+			}
+		}
+		if(!already){
+			partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
+			partition_areas->SetValue(partition[i],area,ADD_VAL);
+		};
+	}
+}
+/*}}}*/
+#endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15385)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15386)
@@ -72,5 +72,4 @@
 		/*}}}*/
 		/*Element virtual functions definitions: {{{*/
-		void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
 		void   ComputeBasalStress(Vector<IssmDouble>* sigma_b);
 		void   ComputeStrainRate(Vector<IssmDouble>* eps);
@@ -105,7 +104,5 @@
 		void   DeleteResults(void);
 		void   MaterialUpdateFromTemperature(void){_error_("not implemented yet");};
-		void   MigrateGroundingLine(IssmDouble* oldfloating,IssmDouble* sheet_ungrounding);
 		int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
-		void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
 		void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
 		void   RequestedOutput(int output_enum,int step,IssmDouble time);
@@ -117,5 +114,4 @@
 		IssmDouble SurfaceArea(void);
 		void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);
-		int    UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
 		IssmDouble TimeAdapt();
 		int*   GetHorizontalNeighboorSids(void);
@@ -175,4 +171,13 @@
 		IssmDouble SurfaceAverageVelMisfit(int weight_index);
 		void   InputControlUpdate(IssmDouble scalar,bool save_parameter);
+		#endif
+		#ifdef _HAVE_GROUNDINGLINE_
+		void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
+		void   MigrateGroundingLine(IssmDouble* oldfloating,IssmDouble* sheet_ungrounding);
+		int    UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
+		#endif
+
+		#ifdef _HAVE_WRAPPERS_
+		void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
 		#endif