Index: ../trunk-jpl/src/c/classes/Loads/Friction.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Friction.h	(revision 18769)
+++ ../trunk-jpl/src/c/classes/Loads/Friction.h	(revision 18770)
@@ -33,6 +33,7 @@
 		void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
 		void  GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
 		void  GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss);
+		void  GetAlpha2WaterLayer(IssmDouble* palpha2,Gauss* gauss);
 		void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
 };
 
Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 18769)
+++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 18770)
@@ -54,6 +54,9 @@
 		case 4:
 			GetAlpha2Temp(palpha2,gauss);
 			break;
+		case 5:
+			GetAlpha2WaterLayer(palpha2,gauss);
+			break;
 	  default:
 			_error_("not supported");
 	}
@@ -249,6 +252,66 @@
 	/*Assign output pointers:*/
 	*palpha2=alpha2;
 }/*}}}*/
+void Friction::GetAlpha2WaterLayer(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+
+	/*This routine calculates the basal friction coefficient 
+	  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
+
+	/*diverse: */
+	IssmDouble  r,s;
+	IssmDouble  drag_p, drag_q;
+	IssmDouble  Neff;
+	IssmDouble  thickness,bed;
+	IssmDouble  vx,vy,vz,vmag;
+	IssmDouble  drag_coefficient;
+	IssmDouble  alpha2;
+
+	/*Recover parameters: */
+	element->GetInputValue(&drag_p,FrictionPEnum);
+	element->GetInputValue(&drag_q,FrictionQEnum);
+	element->GetInputValue(&thickness, gauss,ThicknessEnum);
+	element->GetInputValue(&bed, gauss,BaseEnum);
+	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
+	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
+	IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
+
+	//compute r and q coefficients: */
+	r=drag_q/drag_p;
+	s=1./drag_p;
+
+	//From bed and thickness, compute effective pressure when drag is viscous:
+	Neff=gravity*(rho_ice*thickness+rho_water*bed);
+	if(Neff<0)Neff=0;
+
+	switch(dim){
+		case 1:
+			element->GetInputValue(&vx,gauss,VxEnum);
+			vmag=sqrt(vx*vx);
+			break;
+		case 2:
+			element->GetInputValue(&vx,gauss,VxEnum);
+			element->GetInputValue(&vy,gauss,VyEnum);
+			vmag=sqrt(vx*vx+vy*vy);
+			break;
+		case 3:
+			element->GetInputValue(&vx,gauss,VxEnum);
+			element->GetInputValue(&vy,gauss,VyEnum);
+			element->GetInputValue(&vz,gauss,VzEnum);
+			vmag=sqrt(vx*vx+vy*vy+vz*vz);
+			break;
+		default:
+			_error_("not supported");
+	}
+
+	/*Check to prevent dividing by zero if vmag==0*/
+	if(vmag==0. && (s-1.)<0.) alpha2=0.;
+	else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
+	_assert_(!xIsNan<IssmDouble>(alpha2));
+
+	/*Assign output pointers:*/
+	*palpha2=alpha2;
+}/*}}}*/
 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
 
 	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
Index: ../trunk-jpl/src/m/classes/frictionwaterlayer.m
===================================================================
--- ../trunk-jpl/src/m/classes/frictionwaterlayer.m	(revision 0)
+++ ../trunk-jpl/src/m/classes/frictionwaterlayer.m	(revision 18770)
@@ -0,0 +1,54 @@
+%FRICTIONWATERLAYER class definition
+%
+%   Usage:
+%      frictionwaterlayer=frictionwaterlayer();
+
+classdef frictiontemp
+	properties (SetAccess=public) 
+		coefficient = NaN;
+		p           = NaN;
+		q           = NaN;
+		water_layer = NaN;
+	end
+	methods
+		function obj = frictionwaterlayer(varargin) % {{{
+			switch nargin
+				case 0
+					obj=setdefaultparameters(obj);
+				case 1
+					obj=structtoobj(frictiontemp(),varargin{1});
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function obj = setdefaultparameters(obj) % {{{
+
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+
+			%Early return
+			if ~ismember(StressbalanceAnalysisEnum(),analyses) & ~ismember(ThermalAnalysisEnum(),analyses), return; end
+
+			md = checkfield(md,'fieldname','friction.coefficient','forcing',1,'NaN',1);
+			md = checkfield(md,'fieldname','friction.q','NaN',1,'size',[md.mesh.numberofelements 1]);
+			md = checkfield(md,'fieldname','friction.p','NaN',1,'size',[md.mesh.numberofelements 1]);
+			md = checkfield(md,'fieldname','thermal.spctemperature','forcing',1,'>=',0.);
+
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('Basal shear stress parameters: tau_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b * 1/f(T)\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*(bed+water_layer), r=q/p and s=1/p)'));
+			fielddisplay(obj,'coefficient','frictiontemp coefficient [SI]');
+			fielddisplay(obj,'p','p exponent');
+			fielddisplay(obj,'q','q exponent');
+			fielddisplay(obj,'water_layer','water thickness at the base of the ice (m)');
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+
+			WriteData(fid,'enum',FrictionLawEnum,'data',5,'format','Integer');
+			WriteData(fid,'class','friction','object',obj,'fieldname','coefficient','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
+			WriteData(fid,'class','friction','object',obj,'fieldname','p','format','DoubleMat','mattype',2);
+			WriteData(fid,'class','friction','object',obj,'fieldname','q','format','DoubleMat','mattype',2);
+			WriteData(fid,'class','friction','object',obj,'fieldname','water_layer','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
+		end % }}}
+	end
+end
