Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 23999)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 24000)
@@ -285,5 +285,5 @@
 	IssmDouble Jdet,D_scalar,dt,h;
 	IssmDouble vel,vx,vy,dvxdx,dvydy;
-	IssmDouble lambda_x,lambda_y,tau;
+	IssmDouble xi,tau;
 	IssmDouble dvx[2],dvy[2];
 	IssmDouble* xyz_list = NULL;
@@ -304,4 +304,5 @@
 	ElementMatrix* Ke     = element->NewElementMatrix();
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*		dbasis = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
@@ -367,5 +368,5 @@
 		switch(stabilization){
 			case 0:
-				/*Nothing to be onde*/
+				/*Nothing to be done*/
 				break;
 			case 1:
@@ -393,8 +394,11 @@
 				/*SUPG*/
 				if(dim!=2) _error_("Stabilization "<<stabilization<<" not supported yet for dim != 2");
+				element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+				vxaverage_input->GetInputAverage(&vx);
+				vyaverage_input->GetInputAverage(&vy);
 				vel=sqrt(vx*vx+vy*vy)+1.e-8;
-				tau=0.25;
-				lambda_x=tau*h*vx/vel;
-				lambda_y=tau*h*vy/vel;
+				//xi=0.3130;
+				xi=1;
+				tau=xi*h/(2*vel);
 				break;
 			default:
@@ -415,32 +419,44 @@
 						&Ke->values[0],1);
 		}
-		if(stabilization==5){
+		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);
+			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]);	
+				}
+			}
+			/*Mass matrix - part 3*/
+			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]*(basis[i]*dvxdx+basis[i]*dvydy);	
+				}
+			}
 			
 			/*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);
-			
+			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]);
+            }
+         }
+			/*Advection matrix - part 3, 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])*(basis[i]*dvxdx+basis[i]*dvydy);
+				}
+         }
+
 			/*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);
-		
-		}
+			for(int i=0;i<numnodes;i++){
+            for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
+				}
+         }
+			/*Advection matrix - part 3, B*/
+			for(int i=0;i<numnodes;i++){
+            for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(basis[i]*dvxdx+basis[i]*dvydy);	
+				}
+			}
+		}/*}}}*/
 	}
 
@@ -448,4 +464,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
 	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(Bprime);
@@ -559,5 +576,6 @@
 	IssmDouble  Jdet,dt;
 	IssmDouble  ms,mb,gmb,fmb,thickness;
-	IssmDouble  vx,vy,vel,lambda_x,lambda_y,h,tau;
+	IssmDouble  vx,vy,vel,dvxdx,dvydy,xi,h,tau;
+	IssmDouble  dvx[2],dvy[2];
 	IssmDouble  gllevelset,phi=1.;
 	IssmDouble* xyz_list = NULL;
@@ -640,18 +658,23 @@
 	
 		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);
+			vxaverage_input->GetInputAverage(&vx);
+			vyaverage_input->GetInputAverage(&vy);
+			vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+			vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+			vel=sqrt(vx*vx+vy*vy)+1.e-8;
+			dvxdx=dvx[0];
+			dvydy=dvy[1];
+			//xi=0.3130;
+			xi=1;
+			tau=xi*h/(2*vel);
 			
 			/*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]);
+				pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(tau*vx*dbasis[0*numnodes+i]+tau*vy*dbasis[1*numnodes+i]);
+			}
+			/*Force vector - part 3*/
+			for(int i=0;i<numnodes;i++){
+				pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(tau*basis[i]*dvxdx+tau*basis[i]*dvydy);
 			}
 		}