Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 170)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 171)
@@ -341,7 +341,10 @@
 	/* material data: */
 	double viscosity; //viscosity
+	double oldviscosity; //viscosity
+	double newviscosity; //viscosity
 
 	/* strain rate: */
 	double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+	double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
 
 	/* matrices: */
@@ -367,4 +370,6 @@
 	double* velocity_param=NULL;
 	double  vxvy_list[numgrids][2];
+	double* oldvelocity_param=NULL;
+	double  oldvxvy_list[numgrids][2];
 	double  thickness;
 
@@ -403,4 +408,5 @@
 		/*recover extra inputs from users, at current convergence iteration: */
 		velocity_param=ParameterInputsRecover(inputs,"velocity"); 
+		oldvelocity_param=ParameterInputsRecover(inputs,"old_velocity"); 
 
 		/* Get node coordinates and dof list: */
@@ -427,4 +433,16 @@
 		}
 			
+		/*Initialize oldvxvy_list: */
+		for(i=0;i<numgrids;i++){
+			if(oldvelocity_param){
+				oldvxvy_list[i][0]=oldvelocity_param[doflist[i*numberofdofspernode+0]];
+				oldvxvy_list[i][1]=oldvelocity_param[doflist[i*numberofdofspernode+1]];
+			}
+			else{
+				oldvxvy_list[i][0]=0;
+				oldvxvy_list[i][1]=0;
+			}
+		}
+
 		#ifdef _DEBUGELEMENTS_
 		if(my_rank==RANK && id==ELID){ 
@@ -486,7 +504,9 @@
 				/*Get strain rate from velocity: */
 				GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
+				GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
 			
 				/*Get viscosity: */
 				matice->GetViscosity3d(&viscosity, &epsilon[0]);
+				matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
 
 				/*Get B and Bprime matrices: */
@@ -499,5 +519,7 @@
 				/*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
 				  onto this scalar matrix, so that we win some computational time: */
-				D_scalar=viscosity*gauss_weight*Jdet;
+				
+				newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+				D_scalar=newviscosity*gauss_weight*Jdet;
 				for (i=0;i<5;i++){
 					D[i][i]=D_scalar;
@@ -510,5 +532,5 @@
 						&Ke_gg_gaussian[0][0],0);
 
-				/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
+				/* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
 				for( i=0; i<numdofs; i++){
 					for(j=0;j<numdofs;j++){
@@ -594,9 +616,9 @@
 	/*If this element is on the bed rock, spawn a Tria element out of the bed triangle, and use it 
 	 * to create the stifness matrix for the base velocity: */
-	if(onbed){
-		tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 for the base
-		tria->CreateKMatrixDiagnosticBaseVert(Kgg,inputs, analysis_type);
-		delete tria;
-	}
+	//if(onbed){
+//		tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 for the base
+//		tria->CreateKMatrixDiagnosticBaseVert(Kgg,inputs, analysis_type);
+//		delete tria;
+//	}
 
 	/*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness 
@@ -1678,7 +1700,9 @@
 	/*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the 
 	 *diagnostic base vertical stifness: */
-	tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
-	tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type);
-	delete tria;
+	if(onbed){
+		tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
+		tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type);
+		delete tria;
+	}
 
 	/*Now, handle the standard penta element: */
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 170)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 171)
@@ -663,11 +663,5 @@
 		GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
 
-		if (velocity_param){
-			DL_scalar=alpha2*gauss_weight*Jdet;
-		}
-		else{
-			DL_scalar=0;
-		}
-			
+		DL_scalar=alpha2*gauss_weight*Jdet;
 		for (i=0;i<2;i++){
 			DL[i][i]=DL_scalar;
@@ -2126,4 +2120,5 @@
 					&Ke_gg_gaussian[0][0],0);
 
+		
 		/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
 		for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
@@ -2355,6 +2350,6 @@
 		/*Get bed slope: */
 		GetParameterDerivativeValue(&slope[0], &b[0],&xyz_list[0][0], gauss_l1l2l3);
-		dbdx=slope[0];
-		dbdy=slope[1];
+		dbdx=-slope[0];
+		dbdy=-slope[1];
 
 		/* Get Jacobian determinant: */