Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 23984)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 23985)
@@ -285,4 +285,5 @@
 	IssmDouble Jdet,D_scalar,dt,h;
 	IssmDouble vel,vx,vy,dvxdx,dvydy;
+	IssmDouble lambda_x,lambda_y,tau;
 	IssmDouble dvx[2],dvy[2];
 	IssmDouble* xyz_list = NULL;
@@ -389,4 +390,12 @@
 				}
 				break;
+			case 5:
+				/*SUPG*/
+				if(dim!=2) _error_("Stabilization "<<stabilization<<" not supported yet for dim != 2");
+				vel=sqrt(vx*vx+vy*vy)+1.e-8;
+				tau=0.25;
+				lambda_x=tau*h*vx/vel;
+				lambda_y=tau*h*vy/vel;
+				break;
 			default:
 				_error_("Stabilization "<<stabilization<<" not supported yet");
@@ -405,4 +414,32 @@
 						Bprime,dim,numnodes,0,
 						&Ke->values[0],1);
+		}
+		if(stabilization==5){
+			 /*Mass matrix - part 2*/
+			D_scalar=gauss->weight*Jdet;
+			D[0*dim+0]=D_scalar*lambda_x;
+			D[1*dim+1]=D_scalar*lambda_y;
+			TripleMultiply(Bprime,dim,numnodes,1,
+								D,dim,dim,0,
+								B,dim,numnodes,0,
+								&Ke->values[0],1);
+			
+			/*Advection matrix - part 2, A*/
+			D_scalar=dt*gauss->weight*Jdet;
+			D[0*dim+0]=D_scalar*dvxdx*lambda_x;
+			D[1*dim+1]=D_scalar*dvydy*lambda_y;
+			TripleMultiply(Bprime,dim,numnodes,1,
+						D,dim,dim,0,
+						B,dim,numnodes,0,
+						&Ke->values[0],1);
+			
+			/*Advection matrix - part 2, B*/
+			D[0*dim+0]=D_scalar*vx*lambda_x;
+			D[1*dim+1]=D_scalar*vy*lambda_y;
+			TripleMultiply(Bprime,dim,numnodes,1,
+						D,dim,dim,0,
+						Bprime,dim,numnodes,0,
+						&Ke->values[0],1);
+		
 		}
 	}
@@ -516,4 +553,5 @@
 
 	/*Intermediaries */
+	int			stabilization,dim,domaintype;
 	int         melt_style,point1;
 	bool        mainlyfloating;
@@ -521,8 +559,18 @@
 	IssmDouble  Jdet,dt;
 	IssmDouble  ms,mb,gmb,fmb,thickness;
+	IssmDouble  vx,vy,vel,lambda_x,lambda_y,h,tau;
 	IssmDouble  gllevelset,phi=1.;
 	IssmDouble* xyz_list = NULL;
 	Gauss*      gauss     = NULL;
 
+	/*Get problem dimension*/
+	element->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DverticalEnum:   dim = 1; break;
+		case Domain2DhorizontalEnum: dim = 2; break;
+		case Domain3DEnum:           dim = 2; break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
+	  
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
@@ -531,4 +579,5 @@
 	ElementVector* pe    = element->NewElementVector();
 	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis= xNew<IssmDouble>(dim*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -536,4 +585,5 @@
 	element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	element->FindParam(&stabilization,MasstransportStabilizationEnum);
 	Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
 	Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
@@ -541,5 +591,9 @@
 	Input* ms_input         = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
 	Input* thickness_input  = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
-
+	Input* vxaverage_input  = element->GetInput(VxAverageEnum);										_assert_(vxaverage_input);
+	Input* vyaverage_input  = element->GetInput(VyAverageEnum);										_assert_(vyaverage_input);
+
+	h=element->CharacteristicLength();
+	
 	/*Recover portion of element that is grounded*/
 	phi=element->GetGroundedPortion(xyz_list);
@@ -584,4 +638,23 @@
 
 		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
+	
+		if(stabilization==5){ //SUPG
+			/*Prepare coefficients*/
+			vxaverage_input->GetInputValue(&vx,gauss);
+			vxaverage_input->GetInputValue(&vy,gauss);
+			vel = sqrt(vx*vx+vy*vy)+1.e-8;
+			tau=0.25;
+			lambda_x = tau*h*vx/vel;
+			lambda_y = tau*h*vy/vel;
+	
+			/*Get nodal derivatives*/
+			element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			
+			/*Force vector - part 2*/
+			for(int i=0;i<numnodes;i++){
+				pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(lambda_x*dbasis[numnodes*0+i]+lambda_y*dbasis[numnodes*1+i]);
+			}
+		}
+	
 	}
 
@@ -589,4 +662,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
 	delete gauss;
 	return pe;
Index: /issm/trunk-jpl/src/m/classes/masstransport.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/masstransport.m	(revision 23984)
+++ /issm/trunk-jpl/src/m/classes/masstransport.m	(revision 23985)
@@ -92,5 +92,5 @@
 			md = checkfield(md,'fieldname','masstransport.isfreesurface','values',[0 1]);
 			md = checkfield(md,'fieldname','masstransport.hydrostatic_adjustment','values',{'Absolute' 'Incremental'});
-			md = checkfield(md,'fieldname','masstransport.stabilization','values',[0 1 2 3 4]);
+			md = checkfield(md,'fieldname','masstransport.stabilization','values',[0 1 2 3 4 5]);
 			md = checkfield(md,'fieldname','masstransport.min_thickness','>',0);
 			md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1);
@@ -103,5 +103,5 @@
 			fielddisplay(self,'min_thickness','minimum ice thickness allowed [m]');
 			fielddisplay(self,'hydrostatic_adjustment','adjustment of ice shelves surface and bed elevations: ''Incremental'' or ''Absolute'' ');
-			fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: streamline upwinding, 3: discontinuous Galerkin, 4: Flux Correction Transport');
+			fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: streamline upwinding, 3: discontinuous Galerkin, 4: Flux Correction Transport, 5: streamline upwind Petrov-Galerkin (SUPG)');
 
 			disp(sprintf('\n      %s','Penalty options:'));