Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 26419)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 26420)
@@ -215,8 +215,4 @@
 	IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
 	IssmDouble*    dbasis   = xNew<IssmDouble>(2*numnodes);
-	IssmDouble*    Bprime = NULL;
-	if(stabilization==2){
-		Bprime   = xNew<IssmDouble>(dim*numnodes);
-	}
 
 	/*Retrieve all inputs and parameters*/
@@ -224,4 +220,6 @@
 	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
 	basalelement->FindParam(&migrationmax,MigrationMaxEnum);
+
+	h = element->CharacteristicLength();
 
 	Input* mf_vx_input        = NULL;
@@ -291,5 +289,5 @@
 		switch(stabilization){
 			case 0:
-				// no stabilization, do nothing
+				/*Nothing to be done*/
 				break;
 			case 1:
@@ -309,16 +307,42 @@
 				  {
 					/* Streamline Upwinding */
-					basalelement->ElementSizes(&hx,&hy,&hz);
-					h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
+					mf_vx_input->GetInputAverage(&w[0]);
+					mf_vy_input->GetInputAverage(&w[1]);
+					vel=sqrt(w[0]*w[0]+w[1]*w[1])+1.e-8;
+					IssmDouble tau=h/(2*vel);
 					for(int i=0;i<numnodes;i++){
 						for(int j=0;j<numnodes;j++){
-							Ke->values[i*numnodes+j] += D_scalar*h/(2.*vel)*(
-										dbasis[0*numnodes+i] *(w[0]*w[0]*dbasis[0*numnodes+j] + w[0]*w[1]*dbasis[1*numnodes+j]) +
-										dbasis[1*numnodes+i] *(w[1]*w[0]*dbasis[0*numnodes+j] + w[1]*w[1]*dbasis[1*numnodes+j]) 
-										);
+							Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(
+										w[0]*dbasis[0*numnodes+i]+w[1]*dbasis[1*numnodes+i])*(w[0]*dbasis[0*numnodes+j]+w[1]*dbasis[1*numnodes+j]);
 						}
 					}
 				  }
 				break;
+			case 5:{
+				/*SUPG*/
+				mf_vx_input->GetInputAverage(&w[0]);
+				mf_vy_input->GetInputAverage(&w[1]);
+				vel=sqrt(w[0]*w[0]+w[1]*w[1])+1.e-8;
+				IssmPDouble xi=1.;
+				IssmDouble  tau=xi*h/(2*vel);
+
+				IssmDouble vx,vy;
+				vx = w[0]; vy = w[1];
+				/*Mass matrix - part 2*/
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+					}
+				}
+
+				/*Advection matrix - part 2, A*/
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+					}
+				}
+
+				break;
+			}
 			default:
 				_error_("unknown type of stabilization in LevelsetAnalysis.cpp");
@@ -330,5 +354,4 @@
 	xDelete<IssmDouble>(basis);
 	xDelete<IssmDouble>(dbasis);
-	xDelete<IssmDouble>(Bprime);
 	delete gauss;
 	if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
@@ -341,5 +364,5 @@
 
 	/*Intermediaries */
-	int domaintype;
+	int         domaintype,stabilization;
 	IssmDouble  Jdet,dt;
 	IssmDouble  lsf;
@@ -348,35 +371,56 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = basalelement->GetNumberOfNodes();
+	basalelement->FindParam(&stabilization,MasstransportStabilizationEnum);
 
 	/*Initialize Element vector*/
 	ElementVector* pe = basalelement->NewElementVector();
-	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-
-	if(dt!=0.){
-		/*Initialize basis vector*/
-		IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-		/*Retrieve all inputs and parameters*/
-		basalelement->GetVerticesCoordinates(&xyz_list);
-		Input* levelset_input  = basalelement->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input);
-
-		/* Start  looping on the number of gaussian points: */
-		Gauss* gauss=basalelement->NewGauss(2);
-		while(gauss->next()){
-			basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
-			basalelement->NodalFunctions(basis,gauss);
-
-			/* old function value */
-			levelset_input->GetInputValue(&lsf,gauss);
-			for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
-		}
-
-		/*Clean up and return*/
-		xDelete<IssmDouble>(xyz_list);
-		xDelete<IssmDouble>(basis);
-		basalelement->FindParam(&domaintype,DomainTypeEnum);
-		if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
-		delete gauss;
-	}
+	basalelement->FindParam(&dt,TimesteppingTimeStepEnum); _assert_(dt>0.);
+
+	/*Initialize basis vector*/
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis = NULL;
+	if(stabilization==5) dbasis= xNew<IssmDouble>(2*numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	Input* levelset_input = basalelement->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input);
+	Input* mf_vx_input    = basalelement->GetInput(MovingFrontalVxEnum); _assert_(mf_vx_input);
+	Input* mf_vy_input    = basalelement->GetInput(MovingFrontalVyEnum); _assert_(mf_vy_input);
+
+	IssmDouble h=element->CharacteristicLength();
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=basalelement->NewGauss(2);
+	while(gauss->next()){
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
+
+		/* old function value */
+		levelset_input->GetInputValue(&lsf,gauss);
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
+
+		if(stabilization==5){ /*SUPG*/
+         IssmDouble vx,vy,vel;
+			basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			mf_vx_input->GetInputAverage(&vx);
+			mf_vy_input->GetInputAverage(&vy);
+			vel=sqrt(vx*vx+vy*vy)+1.e-8;
+			IssmPDouble xi=1;
+			IssmDouble  tau=xi*h/(2*vel);
+
+			/*Force vector - part 2*/
+			for(int i=0;i<numnodes;i++){
+				pe->values[i]+=Jdet*gauss->weight*lsf*(tau*vx*dbasis[0*numnodes+i]+tau*vy*dbasis[1*numnodes+i]);
+			}
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+   xDelete<IssmDouble>(dbasis);
+	basalelement->FindParam(&domaintype,DomainTypeEnum);
+	if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
+	delete gauss;
 
 	return pe;