Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24489)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 24490)
@@ -2761,4 +2761,67 @@
 }
 /*}}}*/
+void FemModel::ThicknessAverage(){/*{{{*/
+   
+	int elementswidth                   = this->GetElementsWidth();//just 2D mesh, tria elements
+   int numberofvertices                = this->vertices->NumberOfVertices();//total number of vertices
+
+   IssmDouble weight                   = 0.;
+   IssmDouble* totalweight             = NULL;
+	IssmDouble* Hserial						= NULL;
+   IssmDouble* H                       = xNew<IssmDouble>(elementswidth);
+   Vector<IssmDouble>* vecH				= new Vector<IssmDouble>(numberofvertices);
+   Vector<IssmDouble>* vectotalweight  = new Vector<IssmDouble>(numberofvertices);
+
+   for(int i=0;i<this->elements->Size();i++){
+      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+     
+		/*check if there is ice in this element*/
+		if(!element->IsIceInElement()) continue;  
+	
+		/*get H on the vertices*/
+		element->GetInputListOnVertices(H,ThicknessEnum);
+
+      /*weight to calculate the smoothed H*/
+      weight=1.;//simple average
+
+		/*add in the serial vector*/
+      vecH->SetValue(element->vertices[0]->Sid(),weight*H[0],ADD_VAL);
+      vecH->SetValue(element->vertices[1]->Sid(),weight*H[1],ADD_VAL);
+      vecH->SetValue(element->vertices[2]->Sid(),weight*H[2],ADD_VAL);
+      /*total weight*/
+      vectotalweight->SetValue(element->vertices[0]->Sid(),weight,ADD_VAL);
+      vectotalweight->SetValue(element->vertices[1]->Sid(),weight,ADD_VAL);
+      vectotalweight->SetValue(element->vertices[2]->Sid(),weight,ADD_VAL);
+   }
+
+   /*Assemble and serialize*/
+   vecH->Assemble();
+   vectotalweight->Assemble();
+   Hserial=vecH->ToMPISerial();
+   totalweight=vectotalweight->ToMPISerial();
+
+   /*Divide for the total weight*/
+   for(int i=0;i<numberofvertices;i++){
+      _assert_(totalweight[i]>0);
+      Hserial[i]=Hserial[i]/totalweight[i];
+   }
+
+   /*Set element inputs*/
+   for(int i=0;i<this->elements->Size();i++){
+      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+		H[0]=Hserial[element->vertices[0]->Sid()];
+		H[1]=Hserial[element->vertices[1]->Sid()];
+		H[2]=Hserial[element->vertices[2]->Sid()];
+		element->AddInput2(ThicknessEnum,H,P1Enum);
+	}
+
+ 	/*Cleanup*/
+   delete vecH;
+   delete vectotalweight;
+   xDelete<IssmDouble>(H);
+   xDelete<IssmDouble>(Hserial);
+   xDelete<IssmDouble>(totalweight);
+}
+/*}}}*/
 void FemModel::ThicknessPositivex(IssmDouble* pJ){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 24489)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 24490)
@@ -151,4 +151,5 @@
 		void OmegaAbsGradientx( IssmDouble* pJ);
 		void EtaDiffx( IssmDouble* pJ);
+		void ThicknessAverage();
 		void ThicknessAbsGradientx( IssmDouble* pJ);
 		void ThicknessPositivex(IssmDouble* pJ);
Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp	(revision 24489)
+++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp	(revision 24490)
@@ -70,4 +70,9 @@
 		else{
 			solutionsequence_linear(femmodel);
+			// ThicknessAverage: method not totally tested
+			//if(stabilization==3){
+			//	if(VerboseSolution()) _printf0_("   call thickness average core\n");
+			//	femmodel->ThicknessAverage();
+			//}
 		}
 		femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);