Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15898)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15899)
@@ -6493,8 +6493,7 @@
 	int        stabilization;
 	int        dim;
 	IssmDouble Jdettria,D_scalar,dt,h;
-	IssmDouble vel,vx,vy,dvxdx,dvydy;
-	IssmDouble dvx[2],dvy[2];
+	IssmDouble vel,vx,vy;
 	IssmDouble xyz_list[NUMVERTICES][3];
 
 	/*Fetch number of nodes for this finite element*/
@@ -6526,8 +6525,6 @@
 
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
-		vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
-		vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
 
 		D_scalar=gauss->weight*Jdettria;
 
@@ -6539,8 +6536,6 @@
 		GetBMasstransport(B,&xyz_list[0][0],gauss);
 		GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
 
-		dvxdx=dvx[0];
-		dvydy=dvy[1];
 		D_scalar=dt*gauss->weight*Jdettria;
 
 		D[0][0]=D_scalar*vx;
@@ -6596,8 +6591,7 @@
 	int        stabilization;
 	int        dim;
 	IssmDouble Jdettria,D_scalar,dt,h;
-	IssmDouble vel,vx,vy,dvxdx,dvydy;
-	IssmDouble dvx[2],dvy[2];
+	IssmDouble vel,vx,vy;
 	IssmDouble xyz_list[NUMVERTICES][3];
 
 	/*Fetch number of nodes for this finite element*/
@@ -6629,8 +6623,6 @@
 
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
-		vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
-		vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
 
 		D_scalar=gauss->weight*Jdettria;
 
@@ -6642,8 +6634,6 @@
 		GetBMasstransport(B,&xyz_list[0][0],gauss);
 		GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
 
-		dvxdx=dvx[0];
-		dvydy=dvy[1];
 		D_scalar=dt*gauss->weight*Jdettria;
 
 		D[0][0]=D_scalar*vx;
@@ -6819,7 +6809,7 @@
 	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 	Input* vz_input     = inputs->GetInput(VzEnum);                         _assert_(vz_input);
 	Input* ms_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
-	Input* surface_input= inputs->GetInput(ThicknessEnum);              _assert_(surface_input);
+	Input* surface_input= inputs->GetInput(SurfaceEnum);                    _assert_(surface_input);
 
 	/*Initialize mb_correction to 0, do not forget!:*/
 	/* Start  looping on the number of gaussian points: */
@@ -6831,6 +6821,7 @@
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
+		vz_input->GetInputValue(&vz,gauss);
 		ms_input->GetInputValue(&ms,gauss);
 		surface_input->GetInputValue(&surface,gauss);
 
@@ -6876,6 +6867,7 @@
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
+		vz_input->GetInputValue(&vz,gauss);
 		mb_input->GetInputValue(&mb,gauss);
 		bed_input->GetInputValue(&bed,gauss);
 		if(mb_correction_input)