Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 23142)
+++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 23143)
@@ -676,10 +676,24 @@
 }/*}}}*/
 void Friction::GetAlpha2PISM(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
-	/*Here, we want to parameterize the friction as
-	 *
-	 * ... more comments to come
+	/*Here, we want to parameterize the friction using a pseudoplastic friction law,
+	 * computing the basal shear stress as
+	 *
+	 * alpha2 = tau_c (u_b/(abs(u_b)^(1-q)*u_0^q))
+	 *
+	 * The yield stress tau_c is a function of the effective pressure N
+	 * using a Mohr-Coloumb criterion, so that
+	 * tau_c = tan(phi)*N,
+	 * where phi is the till friction angle, representing sediment strength
+	 *
+	 * The effective pressure is given by Eq. (5) in Aschwanden et al. 2016:
+	 *
+	 * N = delta * P0 * 10^((e_0/Cc)(1-(W/Wmax)))
+	 *
+	 * W is calculated by a non-conserving hydrology model in HydrologyPismAnalysis.cpp
+	 *
+	 * see Aschwanden et al. 2016 and Bueler and Brown, 2009 for more details
 	 */
 
-	/*compute ice pressure P0*/
+	/*compute ice overburden pressure P0*/
 	IssmDouble thickness,base,P0;
 	element->GetInputValue(&thickness, gauss,ThicknessEnum);
@@ -697,9 +711,29 @@
 	element->GetInputValue(&W,gauss,WatercolumnEnum);
 	element->GetInputValue(&Wmax,gauss,FrictionWatercolumnMaxEnum);
+
+// 	/*Check that water column height is within 0 and upper bound, correct if needed*/
+// 	// if watercolumn height is higher than the maximum allowed height, set height to upper bound
+// 	if(W>Wmax){
+// 		W=Wmax;
+// 	}
+// 	// if watercolumn height is negative (shouldn't happen), set it to 0
+// 	if(W<0){
+// 		W=0;
+// 	}
+// 	// if watercolumn height is within 0 and upper bound, nothing to be done
+// 	else{
+// 		//do nothing
+// 	}
+
 	N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax));
 
+	/*Get till friction angles, defined by user [deg]*/
+	IssmDouble phi,pi;
+	element->GetInputValue(&phi,gauss,FrictionTillFrictionAngleEnum);
+
+	/*Convert till friction angle from user-defined deg to rad, which Matlab uses*/
+	phi = phi*pi/180.;
+
 	/*Compute yield stress following a Mohr-Colomb criterion*/
-	IssmDouble phi;
-	element->GetInputValue(&phi,gauss,FrictionTillFrictionAngleEnum);
 	IssmDouble tau_c = N*tan(phi);