Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17367)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17368)
@@ -98,5 +98,5 @@
 	IssmDouble Jdet, dt, D_scalar;
 	IssmDouble h,hx,hy,hz;
-	IssmDouble vel,v[dim];
+	IssmDouble vel,v[dim],w[dim];
 	IssmDouble calvingrate, c[dim];
 	IssmDouble dlsf[dim], norm_dlsf, normal[dim];
@@ -153,12 +153,19 @@
 		norm_dlsf=sqrt(norm_dlsf);
 
-		if(norm_dlsf>1.e-10)
-			for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
+		if(norm_dlsf>1.e-10){
+			for(i=0;i<dim;i++) normal[i]=dlsf[i]/norm_dlsf;
+			for(i=0;i<dim;i++) c[i]=calvingrate*normal[i];
+		}
 		else
 			for(i=0;i<dim;i++) c[i]=0.;
+		
+		for(i=0;i<dim;i++) w[i]=v[i]-c[i];
 
 		for(row=0;row<dim;row++)
 			for(col=0;col<dim;col++)
-				D[row][col]=((row==col)?D_scalar*(v[row]-c[row]):0.);
+				if(row==col)
+					D[row][col]=D_scalar*w[row];
+				else
+				   D[row][col]=0.;
 
 		TripleMultiply(B,dim,numnodes,1,
@@ -179,10 +186,12 @@
 				/* Artificial Diffusion */
 				element->ElementSizes(&hx,&hy,&hz);
-				h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); //FIXME: is this correct?
-
-				kappa=h*vel/2.; //FIXME: insert suitable value for kappa
+				h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); 
+				kappa=h*vel/2.;
 				for(row=0;row<dim;row++)
 					for(col=0;col<dim;col++)
-						D[row][col]=((row==col)?D_scalar*kappa:0.);
+					if(row==col)
+						D[row][col]=D_scalar*kappa;
+					else
+						D[row][col]=0.;
 
 				TripleMultiply(Bprime,dim,numnodes,1,