Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 27003)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 27004)
@@ -10,4 +10,5 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+#include <math.h>
 
 void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
@@ -423,4 +424,71 @@
 				break;
 			}
+			case 6:{
+				/*SUPG*/
+				IssmDouble vx,vy;
+				mf_vx_input->GetInputAverage(&vx);
+				mf_vy_input->GetInputAverage(&vy);
+				vel=sqrt(vx*vx+vy*vy)+1.e-8;
+				IssmDouble ECN, K;
+				ECN = vel *dt /h;
+				K = 1./tanh(ECN) - 1./ECN;
+//				if (ECN<1e-6) K = ECN /3.0;
+
+				/*According to Hilmar, xi=K is too large*/
+				IssmPDouble xi=0.1*K;
+
+				IssmDouble  tau=xi*h/(2*vel);
+				Input* levelset_input = NULL;
+
+
+				IssmDouble kappa;
+				IssmDouble p=4, q=4;
+				IssmDouble phi[3];
+
+			   levelset_input=basalelement->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input);
+				levelset_input->GetInputValue(&phi[0], gauss);
+
+				IssmDouble dphidx=0., dphidy=0.;
+				IssmDouble nphi;
+
+				for(int i=0;i<numnodes;i++){
+					dphidx += phi[i]*dbasis[0*numnodes+i];
+					dphidy += phi[i]*dbasis[1*numnodes+i];
+				}
+				nphi = sqrt(dphidx*dphidx+dphidy*dphidy);
+			
+				if (nphi >= 1) {
+					kappa = 1 - 1.0/nphi;
+				}
+				else {
+					kappa = 0.5/M_PI *sin(2*M_PI*nphi)/nphi;
+				}
+
+				kappa = kappa * vel / h;
+
+				/*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]);
+					}
+				}
+				/*Add the pertubation term \nabla\cdot(\kappa*\nabla\phi)*/
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						for(int k=0;k<dim;k++){
+								Ke->values[i*numnodes+j]+= dt*gauss->weight*Jdet*kappa*dbasis[k*numnodes+j]*dbasis[k*numnodes+i];
+						}
+					}
+				}
+
+				break;
+			}
 			default:
 				_error_("unknown type of stabilization in LevelsetAnalysis.cpp");
@@ -458,5 +526,5 @@
 	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 	IssmDouble*    dbasis = NULL;
-	if(stabilization==5) dbasis= xNew<IssmDouble>(2*numnodes);
+	if((stabilization==5) |(stabilization == 6)) dbasis= xNew<IssmDouble>(2*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -479,5 +547,5 @@
 
 		if(stabilization==5){ /*SUPG*/
-         IssmDouble vx,vy,vel;
+			IssmDouble vx,vy,vel;
 			basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 			mf_vx_input->GetInputAverage(&vx);
@@ -492,4 +560,26 @@
 			}
 		}
+		else if (stabilization ==6) {
+			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;
+
+			IssmDouble ECN, K;
+			ECN = vel *dt /h;
+			K = 1./tanh(ECN) - 1./ECN;
+	//		if (ECN<1e-6) K = ECN /3.0;
+
+			/*According to Hilmar, xi=K is too large*/
+			IssmPDouble xi=0.1*K;
+
+			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]+vy*dbasis[1*numnodes+i]);
+			}
+		}
 	}
 
@@ -497,5 +587,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
-   xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(dbasis);
 	basalelement->FindParam(&domaintype,DomainTypeEnum);
 	if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};
Index: /issm/trunk-jpl/src/m/classes/levelset.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/levelset.m	(revision 27003)
+++ /issm/trunk-jpl/src/m/classes/levelset.m	(revision 27004)
@@ -62,5 +62,5 @@
 
 			md = checkfield(md,'fieldname','levelset.spclevelset','Inf',1,'timeseries',1);
-			md = checkfield(md,'fieldname','levelset.stabilization','values',[0 1 2 5]);
+			md = checkfield(md,'fieldname','levelset.stabilization','values',[0 1 2 5 6]);
 			md = checkfield(md,'fieldname','levelset.kill_icebergs','numel',1,'values',[0 1]);
 			md = checkfield(md,'fieldname','levelset.migration_max','numel',1,'NaN',1,'Inf',1,'>',0);