Index: /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp	(revision 16017)
+++ /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp	(revision 16018)
@@ -19,10 +19,10 @@
 	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
-	if(VerboseSolution()) _printf0_("call computational core:\n");
-	//femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum);
-	//solutionsequence_linear(femmodel);
-	//femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum);
-	//solutionsequence_linear(femmodel);
-	surfaceslope_core(femmodel);
+	if(VerboseSolution()) _printf0_("computing smoothed slopes:\n");
+	femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum);
+	solutionsequence_linear(femmodel);
+	femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum);
+	solutionsequence_linear(femmodel);
+	//surfaceslope_core(femmodel);
 
 	if(VerboseSolution()) _printf0_("call computational core:\n");
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16017)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16018)
@@ -7231,6 +7231,7 @@
 	IssmDouble h,gamma,thickness;
 	IssmDouble hnx,hny,dhnx[2],dhny[2];
-	int        i,j;
-	IssmDouble D_scalar;
+	IssmDouble vel,vx,vy,dvxdx,dvydy;
+	IssmDouble D[2][2];
+	int        stabilization=1;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -7239,4 +7240,5 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
+	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
@@ -7295,4 +7297,33 @@
 							);
 			}
+		}
+
+		vx=hnx;
+		vy=hny;
+		if(stabilization==2){
+			/*Streamline upwinding*/
+			vel=sqrt(vx*vx+vy*vy)+1.e-8;
+			D[0][0]=h/(2*vel)*vx*vx;
+			D[1][0]=h/(2*vel)*vy*vx;
+			D[0][1]=h/(2*vel)*vx*vy;
+			D[1][1]=h/(2*vel)*vy*vy;
+		}
+		else if(stabilization==1){
+			/*MacAyeal*/
+			D[0][0]=h/2.0*fabs(vx);
+			D[0][1]=0.;
+			D[1][0]=0.;
+			D[1][1]=h/2.0*fabs(vy);
+		}
+		if(stabilization==1 || stabilization==2){
+			GetBprimeMasstransport(B,&xyz_list[0][0],gauss);
+			D[0][0]=gauss->weight*Jdettria*D[0][0];
+			D[1][0]=gauss->weight*Jdettria*D[1][0];
+			D[0][1]=gauss->weight*Jdettria*D[0][1];
+			D[1][1]=gauss->weight*Jdettria*D[1][1];
+			TripleMultiply(B,2,numnodes,1,
+						&D[0][0],2,2,0,
+						B,2,numnodes,0,
+						&Ke->values[0],1);
 		}
 	}
@@ -7306,4 +7337,5 @@
 	xDelete<IssmDouble>(HNx);
 	xDelete<IssmDouble>(HNy);
+	xDelete<IssmDouble>(B);
 	delete gauss;
 	return Ke;
@@ -7317,5 +7349,5 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble D[2][2];
-	IssmDouble l=12.;
+	IssmDouble l=15.;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -7339,4 +7371,5 @@
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		thickness_input->GetInputValue(&thickness,gauss);
+		if(thickness<50.) thickness=50.;
 
 		D_scalar=gauss->weight*Jdet*(l*thickness)*(l*thickness);