Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16822)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16823)
@@ -8447,5 +8447,5 @@
 	IssmDouble  z_list[NUMVERTICES];
 	IssmDouble  slope[2];
-	IssmDouble  slope2,constant_part;
+	IssmDouble  slope2,constant_part,drag;
 	IssmDouble  rho_ice,gravity,n,B;
 	IssmDouble  ub,vb,z_g,surface,thickness;
@@ -8468,4 +8468,5 @@
 	Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
 	Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
+	Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
 	for(i=0;i<NUMVERTICES;i++)z_list[i]=xyz_list[i][2];
 
@@ -8512,7 +8513,15 @@
 		/*Deal with lower surface*/
 		if(IsOnBed()){
-			constant_part=-1.58*1.e-10*rho_ice*gravity*thickness;
-			ub=constant_part*slope[0];
-			vb=constant_part*slope[1];
+			drag_input->GetInputValue(&drag,gauss);
+
+			/*No ref...*/
+			ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
+			vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
+			///*Ritz et al. 1996*/
+			//ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
+			//vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
+			///*Rutt et al. 2009*/
+			//ub=-drag*rho_ice*gravity*thickness*slope[0];
+			//vb=-drag*rho_ice*gravity*thickness*slope[1];
 
 			pe->values[2*node0+0]+=ub/connectivity[0];
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16822)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16823)
@@ -4220,6 +4220,6 @@
 
 	/*Intermediaries */
-	IssmDouble constant_part,ub,vb;
-	IssmDouble rho_ice,gravity,n,B;
+	IssmDouble ub,vb;
+	IssmDouble rho_ice,gravity,n,B,drag;
 	IssmDouble slope2,thickness,connectivity;
 	IssmDouble slope[2];
@@ -4239,4 +4239,5 @@
 	Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
 	Input* thickness_input=inputs->GetInput(ThicknessEnum);  _assert_(thickness_input);
+	Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
 
 	/*Spawn 3 sing elements: */
@@ -4249,12 +4250,18 @@
 
 		thickness_input->GetInputValue(&thickness,gauss);
+		drag_input->GetInputValue(&drag,gauss);
 		slopex_input->GetInputValue(&slope[0],gauss);
 		slopey_input->GetInputValue(&slope[1],gauss);
-		slope2=pow(slope[0],2)+pow(slope[1],2);
-
-		constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
-
+		slope2=slope[0]*slope[0]+slope[1]*slope[1];
+
+		/*No ref...*/
 		ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
 		vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
+		///*Ritz et al. 1996*/
+		//ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
+		//vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
+		///*Rutt et al. 2009*/
+		//ub=-drag*rho_ice*gravity*thickness*slope[0];
+		//vb=-drag*rho_ice*gravity*thickness*slope[1];
 
 		pe->values[2*i+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity;