Index: /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 23965)
+++ /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 23966)
@@ -61,4 +61,5 @@
 
 	/*Intermediaries */
+	const IssmPDouble yts = 365*24*3600.;
 	int        domaintype,dim;
 	IssmDouble Jdet,D_scalar,onboundary;
@@ -104,40 +105,59 @@
 		GetBprime(Bprime,element,dim,xyz_list,gauss);
 
-		D_scalar=gauss->weight*Jdet;
-
 		/*Get velocity*/
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
-		bc_input->GetInputValue(&onboundary,gauss);
-		if(onboundary>0.){
-			/*We do not want to advect garbage, make sure only diffusion is applied on boundary*/
-			vx = 0.; vy = 0.;
-		}
-
-		/*Diffusion */
-		if(sqrt(vx*vx+vy*vy)<1000./31536000.){
-			IssmPDouble kappa = -10.;
+
+		if(false){
+			/*Streamline diffusion*/
+			vel = sqrt(vx*vx+vy*vy);
+			if(vel<10./yts){
+				vx = 0.; vy = 0.;
+				vel = 30./yts*500000.;
+			}
+
 			for(int i=0;i<numnodes;i++){
 				for(int j=0;j<numnodes;j++){
-					Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
+					Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+								(vx*dbasis[0*numnodes+i] + vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j])
+								+ vel/500000.*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]));
 				}
 			}
 		}
-
-		/*Advection: */
-		for(int i=0;i<numnodes;i++){
-			for(int j=0;j<numnodes;j++){
-				Ke->values[i*numnodes+j] += (D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]))*1e-2;
-			}
-		}
-
-		/*Artificial diffusivity*/
-		vel=sqrt(vx*vx + vy*vy)+1.e-14;
-		D[0][0]=D_scalar*h/(2.*vel)*fabs(vx*vx);  D[0][1]=D_scalar*h/(2.*vel)*fabs(vx*vy);
-		D[1][0]=D_scalar*h/(2.*vel)*fabs(vy*vx);  D[1][1]=D_scalar*h/(2.*vel)*fabs(vy*vy);
-		TripleMultiply(Bprime,dim,numnodes,1,
-					&D[0][0],2,2,0,
-					Bprime,dim,numnodes,0,
-					&Ke->values[0],1);
+		else{
+			D_scalar=gauss->weight*Jdet;
+
+			bc_input->GetInputValue(&onboundary,gauss);
+			if(onboundary>0.){
+				/*We do not want to advect garbage, make sure only diffusion is applied on boundary*/
+				vx = 0.; vy = 0.;
+			}
+
+			/*Diffusion */
+			if(sqrt(vx*vx+vy*vy)<1000./31536000.){
+				IssmPDouble kappa = -10.;
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
+					}
+				}
+			}
+
+			/*Advection: */
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += (D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]))*1e-2;
+				}
+			}
+
+			/*Artificial diffusivity*/
+			vel=sqrt(vx*vx + vy*vy)+1.e-14;
+			D[0][0]=D_scalar*h/(2.*vel)*fabs(vx*vx);  D[0][1]=D_scalar*h/(2.*vel)*fabs(vx*vy);
+			D[1][0]=D_scalar*h/(2.*vel)*fabs(vy*vx);  D[1][1]=D_scalar*h/(2.*vel)*fabs(vy*vy);
+			TripleMultiply(Bprime,dim,numnodes,1,
+						&D[0][0],2,2,0,
+						Bprime,dim,numnodes,0,
+						&Ke->values[0],1);
+		}
 	}