Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 24119)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 24120)
@@ -376,16 +376,15 @@
 				break;
 			case 2:
+				/*Streamline upwinding*/
+				element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+				vxaverage_input->GetInputAverage(&vx);
 				if(dim==1){
 					vel=fabs(vx)+1.e-8;
-					D[0]=h/(2*vel)*vx*vx;
 				}
 				else{
-					/*Streamline upwinding*/
+					vyaverage_input->GetInputAverage(&vy);
 					vel=sqrt(vx*vx+vy*vy)+1.e-8;
-					D[0*dim+0]=h/(2*vel)*vx*vx;
-					D[1*dim+0]=h/(2*vel)*vy*vx;
-					D[0*dim+1]=h/(2*vel)*vx*vy;
-					D[1*dim+1]=h/(2*vel)*vy*vy;
 				}
+				tau=h/(2*vel);
 				break;
 			case 5:
@@ -403,5 +402,6 @@
 				_error_("Stabilization "<<stabilization<<" not supported yet");
 		}
-		if(stabilization==1 || stabilization==2){
+		if(stabilization==1){
+			/*SSA*/
 			if(dim==1) D[0]=D_scalar*D[0];
 			else{
@@ -416,4 +416,12 @@
 						Bprime,dim,numnodes,0,
 						&Ke->values[0],1);
+		}
+		if(stabilization==2){
+			/*Streamline upwind*/
+			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+i]+vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
+				}
+			}
 		}
 		if(stabilization==5){/*{{{*/