Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Material.h	(revision 21699)
+++ /issm/trunk-jpl/src/c/classes/Materials/Material.h	(revision 21700)
@@ -53,4 +53,7 @@
 		virtual void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
 		virtual void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
+		virtual void       ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
+		virtual void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
+		virtual void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
 
 };
Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp	(revision 21699)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp	(revision 21700)
@@ -160,4 +160,5 @@
 /*}}}*/
 IssmDouble Matestar::GetB(){/*{{{*/
+
 	/*Output*/
 	IssmDouble B;
@@ -168,4 +169,5 @@
 /*}}}*/
 IssmDouble Matestar::GetBbar(){/*{{{*/
+
 	/*Output*/
 	IssmDouble Bbar;
@@ -193,4 +195,13 @@
 }
 /*}}}*/
+IssmDouble Matestar::GetEcbar(){/*{{{*/
+
+	/*Output*/
+	IssmDouble Ecbar;
+
+	element->inputs->GetInputAverage(&Ecbar,MaterialsRheologyEcbarEnum);
+	return Ecbar;
+}
+/*}}}*/
 IssmDouble Matestar::GetEs(){/*{{{*/
 
@@ -202,4 +213,13 @@
 }
 /*}}}*/
+IssmDouble Matestar::GetEsbar(){/*{{{*/
+
+	/*Output*/
+	IssmDouble Esbar;
+
+	element->inputs->GetInputAverage(&Esbar,MaterialsRheologyEsbarEnum);
+	return Esbar;
+}
+/*}}}*/
 IssmDouble Matestar::GetN(){/*{{{*/
 
@@ -211,75 +231,32 @@
 void  Matestar::GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
 	_error_("not implemented yet");
-	/*From a string tensor and a material object, return viscosity, using Glen's flow law.
-								(1-D) B
-	  viscosity= -------------------------
-						  2 eps_eff ^[(n-1)/n]
-
-	  where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
-	  and n the flow law exponent.
-
-	  If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we 
-	  return 10^14, initial viscosity.
-	  */
+}
+/*}}}*/
+void  Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+void  Matestar::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}
+/*}}}*/
+IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
 
 	/*output: */
 	IssmDouble viscosity;
 
-	/*Intermediary: */
-	IssmDouble B,D=0.,n;
-
-	/*Get B and n*/
-	B=GetB(); _assert_(B>0.);
-	n=GetN(); _assert_(n>0.);
-
-	if (n==1.){
-		/*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
-		viscosity=(1.-D)*B/2.;
-	}
-	else{
-
-		/*if no strain rate, return maximum viscosity*/
-		if(eps_eff==0.){
-			viscosity = 1.e+14/2.;
-			//viscosity = B;
-			//viscosity=2.5*pow(10.,17);
-		}
-
-		else{
-			viscosity=(1.-D)*B/(2.*pow(eps_eff,(n-1.)/n));
-		}
-	}
-
-	/*Checks in debugging mode*/
-	if(viscosity<=0) _error_("Negative viscosity");
-
-	/*Return: */
-	*pviscosity=viscosity;
-}
-/*}}}*/
-void  Matestar::GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matestar::GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matestar::GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-void  Matestar::GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}
-/*}}}*/
-IssmDouble Matestar::GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble viscosity;
+	/*Intermediaries*/
 	IssmDouble epseff,epsprime_norm;
 	IssmDouble lambdas;
 	IssmDouble vmag,dvmag[3];
-	IssmDouble B,Ec,Es,E;
+	IssmDouble B,Ec,Es,E,n;
 
 	/*Calculate velocity magnitude and its derivative*/
@@ -300,7 +277,15 @@
 
 	/*Get B and enhancement*/
-	B=GetB(); _assert_(B>0.);
-	Ec=GetEc(); _assert_(Ec>=0.);
-	Es=GetEs(); _assert_(Es>=0.);
+	n=GetN(); _assert_(n>0.);
+	if (isdepthaveraged==0.){
+		B=GetB(); _assert_(B>0.);
+		Ec=GetEc(); _assert_(Ec>=0.);
+		Es=GetEs(); _assert_(Es>=0.);
+	}
+	else{
+		B=GetBbar(); _assert_(B>0.);
+		Ec=GetEcbar(); _assert_(Ec>=0.);
+		Es=GetEsbar(); _assert_(Es>=0.);
+	}
 
 	/*Get total enhancement factor E(lambdas)*/
@@ -315,5 +300,5 @@
 	}
 	else{
-		viscosity = B/(2.*pow(E*epseff*epseff,1./3.));
+		viscosity = B/(2.*pow(E,1./n)*pow(epseff,2./n));
 	}
 
@@ -322,23 +307,53 @@
 }
 /*}}}*/
+IssmDouble Matestar::GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble dmudB;
+	IssmDouble epseff,epsprime_norm;
+	IssmDouble lambdas;
+	IssmDouble vmag,dvmag[3];
+	IssmDouble Ec,Es,E;
+
+	/*Calculate velocity magnitude and its derivative*/
+	vmag = sqrt(vx*vx+vy*vy+vz*vz);
+	if(vmag<1e-12){
+		dvmag[0]=0;
+		dvmag[1]=0;
+		dvmag[2]=0;
+	}
+	else{
+		dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
+		dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
+		dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
+	}
+
+	EstarStrainrateQuantities(&epseff,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
+	lambdas=EstarLambdaS(epseff,epsprime_norm);
+
+	/*Get enhancement*/
+	if (isdepthaveraged==0.){
+		Ec=GetEc(); _assert_(Ec>=0.);
+		Es=GetEs(); _assert_(Es>=0.);
+	}
+	else{
+		Ec=GetEcbar(); _assert_(Ec>=0.);
+		Es=GetEsbar(); _assert_(Es>=0.);
+	}
+
+	/*Get total enhancement factor E(lambdas)*/
+	E = Ec + (Es-Ec)*lambdas*lambdas; _assert_(E>0.);
+
+	/*Compute dmudB*/
+	if(epseff==0.) dmudB = 0.;
+	else           dmudB = 1./(2.*pow(E,1./3.)*pow(epseff,2./3.));
+
+	/*Assign output*/
+	return dmudB;
+
+}
+/*}}}*/
 void  Matestar::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){/*{{{*/
-	/*output: */
-	IssmDouble dmudB;
-
-	/*Intermediary: */
-	IssmDouble E=1.,n;
-
-	n=GetN(); _assert_(n>0.);
-	if(n==1.){
-		/*Linear Viscous behavior (Newtonian fluid) dmudB=B/2E: */
-		dmudB=1./(2.*E);
-	}
-	else{
-		if(eps_eff==0.) dmudB = 0.;
-		else            dmudB = 1./(2.*pow(E*eps_eff*eps_eff,1./3.));
-	}
-
-	/*Return: */
-	*pdmudB=dmudB;
+	 _error_("not implemented yet");
 }
 /*}}}*/
@@ -393,10 +408,10 @@
 }
 /*}}}*/
-void  Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+void  Matestar::ViscosityBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
 
 	/*Intermediaries*/
 	IssmDouble vx,vy,vz;
 	IssmDouble dvx[3],dvy[3],dvz[3];
-	IssmDouble B,Ec,Es;
+	bool isdepthaveraged=0.;
 
 	/*Get velocity derivatives in all directions*/
@@ -418,16 +433,14 @@
 	}
 
-	/*Compute viscosity*/
-	*pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
-}
-/*}}}*/
-void  Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
-}/*}}}*/
-void  Matestar::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+	/*Compute dmudB*/
+	*pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
+}
+/*}}}*/
+void  Matestar::ViscosityBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
 
 	/*Intermediaries*/
 	IssmDouble vx,vy,vz;
 	IssmDouble dvx[3],dvy[3],dvz[3];
+	bool isdepthaveraged=0.;
 
 	/*Get velocity derivatives in all directions*/
@@ -450,17 +463,11 @@
 
 	/*Compute viscosity*/
-	*pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
-}/*}}}*/
-void  Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
-void  Matestar::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
-void  Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+	*pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
+}/*}}}*/
+void  Matestar::ViscosityBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
 	/*Intermediaries*/
 	IssmDouble vx,vy,vz;
 	IssmDouble dvx[3],dvy[3],dvz[3];
-	IssmDouble B,Ec,Es;
+	bool isdepthaveraged=1.;
 
 	/*Get velocity derivatives in all directions*/
@@ -486,5 +493,102 @@
 
 	/*Compute viscosity*/
-	*pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
+	*pdmudB=GetViscosity_BGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
+}/*}}}*/
+void  Matestar::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble vx,vy,vz;
+	IssmDouble dvx[3],dvy[3],dvz[3];
+	bool isdepthaveraged=0.;
+
+	/*Get velocity derivatives in all directions*/
+	_assert_(dim>1);
+	_assert_(vx_input);
+	vx_input->GetInputValue(&vx,gauss);
+	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+	_assert_(vy_input);
+	vy_input->GetInputValue(&vy,gauss);
+	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+	if(dim==3){
+		_assert_(vz_input);
+		vz_input->GetInputValue(&vz,gauss);
+		vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+	}
+	else{
+		vz = 0.;
+		dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
+	}
+
+	/*Compute viscosity*/
+	*pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
+}
+/*}}}*/
+void  Matestar::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
+}/*}}}*/
+void  Matestar::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble vx,vy,vz;
+	IssmDouble dvx[3],dvy[3],dvz[3];
+	bool isdepthaveraged=0.;
+
+	/*Get velocity derivatives in all directions*/
+	_assert_(dim==2 || dim==3);
+	_assert_(vx_input);
+	vx_input->GetInputValue(&vx,gauss);
+	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+	if(dim==3){
+		_assert_(vy_input);
+		vy_input->GetInputValue(&vy,gauss);
+		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+	}
+	else{
+		dvx[2] = 0.;
+		vy = 0.;
+		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
+	}
+	vz = 0.;
+	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
+
+	/*Compute viscosity*/
+	*pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
+}/*}}}*/
+void  Matestar::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void  Matestar::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void  Matestar::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble vx,vy,vz;
+	IssmDouble dvx[3],dvy[3],dvz[3];
+	bool isdepthaveraged=1.;
+
+	/*Get velocity derivatives in all directions*/
+	_assert_(dim==1 || dim==2);
+	_assert_(vx_input);
+	vx_input->GetInputValue(&vx,gauss);
+	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+	if(dim==2){
+		_assert_(vy_input);
+		vy_input->GetInputValue(&vy,gauss);
+		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+	}
+	else{
+		dvx[1] = 0.;
+		dvx[2] = 0.;
+		vy = 0.;
+		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
+	}
+	dvx[2] = 0.;
+	dvy[2] = 0.;
+	vz = 0.;
+	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
+
+	/*Compute viscosity*/
+	*pviscosity=GetViscosityGeneral(vx,vy,vz,&dvx[0],&dvy[0],&dvz[0],isdepthaveraged);
 }/*}}}*/
 void  Matestar::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matestar.h	(revision 21699)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.h	(revision 21700)
@@ -70,5 +70,7 @@
 		IssmDouble GetDbar();
 		IssmDouble GetEc();
+		IssmDouble GetEcbar();
 		IssmDouble GetEs();
+		IssmDouble GetEsbar();
 		IssmDouble GetN();
 		bool       IsDamage();
@@ -84,6 +86,10 @@
 		void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
 		void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
+		void       ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
+		void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
+		void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
 		/*}}}*/
-		IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
+		IssmDouble GetViscosityGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged);
+		IssmDouble GetViscosity_BGeneral(IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,bool isdepthaveraged);
 };
 
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 21699)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 21700)
@@ -88,4 +88,7 @@
 		void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
 		void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
+		void       ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
+		void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
+		void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
 		/*}}}*/
 };
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 21699)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 21700)
@@ -125,4 +125,7 @@
 		void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
 		void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");};
+		void       ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
+		void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
+		void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
 		/*}}}*/
 		/*Numerics: {{{*/