Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10448)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10449)
@@ -4603,4 +4603,6 @@
 	/*material parameters: */
 	double mu_water=0.001787;  //unit= [N s/m2]
+	double VelocityFactor;  // This factor represents the number 12 in laminar flow velocity which can vary by differnt hydrology.CR
+	double n_man,CR;
 	double w;
 	double rho_ice, rho_water, g;
@@ -4614,4 +4616,6 @@
 	rho_water=matpar->GetRhoWater();
 	g=matpar->GetG();
+	CR=matpar->GetHydrologyCR(); // To have Lebrocq equavalent equation: CR=0.01,n_man=0.02
+	n_man=matpar->GetHydrologyN(); 
 	Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
 	Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
@@ -4619,4 +4623,7 @@
 	Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum);         _assert_(bedslopey_input);
 	Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);     _assert_(watercolumn_input);
+
+	/* compute VelocityFactor */
+	VelocityFactor= n_man*pow(CR,2)*rho_water*g/mu_water;
 
 	gauss=new GaussTria();
@@ -4630,8 +4637,9 @@
 
 		/* Water velocity x and y components */
-		vx[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
-		vy[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
-		//vx[iv]=0.001;
-		//vy[iv]=0;
+	//	vx[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
+	//	vy[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
+	
+		vx[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
+		vy[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
 	}