Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp	(revision 17881)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp	(revision 17882)
@@ -331,6 +331,4 @@
 	IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
 	IssmDouble  g         = element->GetMaterialParameter(ConstantsGEnum);
-	IssmDouble  CR        = element->GetMaterialParameter(HydrologyshreveCREnum);
-	IssmDouble  n_man     = element->GetMaterialParameter(HydrologyshreveNEnum);
 	IssmDouble  mu_water  = element->GetMaterialParameter(MaterialsMuWaterEnum);
 	Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
@@ -339,7 +337,4 @@
 	Input* bedslopey_input     = element->GetInput(BedSlopeYEnum);     _assert_(bedslopey_input);
 	Input* watercolumn_input   = element->GetInput(WatercolumnEnum);   _assert_(watercolumn_input);
-
-	/* compute VelocityFactor */
-	IssmDouble VelocityFactor = n_man*CR*CR*rho_water*g/mu_water;
 
 	/*Fetch number of vertices and allocate output*/
@@ -358,8 +353,6 @@
 
 		/* Water velocity x and y components */
-	//	vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
-	//	vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
-		vx[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
-		vy[iv]= - w*w/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
+		vx[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
+		vy[iv]= - w*w/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
 	}
 
@@ -373,2 +366,25 @@
 	xDelete<IssmDouble>(vy);
 }/*}}}*/
+
+
+
+/*Needed changes to switch to the Johnson formulation*//*{{{*/
+/*All the changes are to be done in the velocity computation.
+	The new velocity needs some new parameter that should be introduce in the hydrologyshreve class:
+	'p' and 'q' which are the exponent of the Manning formula for laminar (p=2,q=1) or turbulent (p=2/3,q=1/2) flow
+	'R' the hydraulic radius
+	'n' the manning roughness coeficient
+
+	With these, the velocity reads ;
+
+	v= - (1/n)* pow(R,p)*pow((grad phi(rho_water*g)),q)
+
+	you should also redefine the water pressure potential 'phi' with respect to the effective pressure deffinition given in Johson:
+	phi=(rho_ice*g*( surface + ((rho_water/rho_ice)-1)*base - k_n*((thickness* grad(base))/omega) )
+
+	where 
+	'omega' is the fractional area of the base occupied by the water film
+	'k_n' is a parameter
+	This last equation derives from the effective pressure definition developped in Alley 1989
+*/
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 17881)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 17882)
@@ -60,13 +60,5 @@
 			_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
 	}
-
-	if(hydrology_model==HydrologyshreveEnum){
-		iomodel->Constant(&this->hydro_CR,HydrologyshreveCREnum);
-		iomodel->Constant(&this->hydro_kn,HydrologyshreveKnEnum);
-		iomodel->Constant(&this->hydro_n,HydrologyshreveNEnum);
-		iomodel->Constant(&this->hydro_p,HydrologyshrevePEnum);
-		iomodel->Constant(&this->hydro_q,HydrologyshreveQEnum);
-	}
-	else if(hydrology_model==HydrologydcEnum){
+	if(hydrology_model==HydrologydcEnum){
 		iomodel->Constant(&this->sediment_compressibility,HydrologydcSedimentCompressibilityEnum);
 		iomodel->Constant(&this->sediment_porosity,HydrologydcSedimentPorosityEnum);