Index: ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18544)
+++ ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 18545)
@@ -29,6 +29,7 @@
 	iomodel->FetchDataToInput(elements,SurfaceEnum);
 	iomodel->FetchDataToInput(elements,BaseEnum);
 	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
 	iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
 	iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
 	iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
@@ -104,7 +105,7 @@
 ElementVector* Balancethickness2Analysis::CreatePVector(Element* element){/*{{{*/
 
 	/*Intermediaries */
-	IssmDouble  dhdt,mb,ms,dD[2],db[2],Jdet;
+	IssmDouble  dhdt,mb,ms,D,db[2],Jdet;
 	IssmDouble* xyz_list = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -113,6 +114,7 @@
 	/*Initialize Element vector and other vectors*/
 	ElementVector* pe    = element->NewElementVector();
 	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	IssmDouble*   dbasis = xNew<IssmDouble>(2*numnodes);
 
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
@@ -129,17 +131,20 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 		ms_input->GetInputValue(&ms,gauss);
 		mb_input->GetInputValue(&mb,gauss);
 		dhdt_input->GetInputValue(&dhdt,gauss);
+		D_input->GetInputValue(&D,gauss);
 		bed_input->GetInputDerivativeValue(&db[0],xyz_list,gauss);
-		D_input->GetInputDerivativeValue(&dD[0],xyz_list,gauss);
-		db[0]=0.;
-		db[1]=0.;
+		//db[0]=0.;
+		//db[1]=0.;
 
-		/*Since grad(b) is constant div(D grad(b) ) = grad(D).grad(b)*/
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt +dD[0]*db[0]+dD[1]*db[1])*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(
+					(ms-mb-dhdt)*basis[i]
+					-D*db[0]*dbasis[0*numnodes+i] - D*db[1]*dbasis[1*numnodes+i]
+					);
 	}
 
 	/*Clean up and return*/
@@ -166,8 +171,8 @@
 void Balancethickness2Analysis::CreateDiffusionCoefficient(Element* element){/*{{{*/
 
 	/*Intermediaries */
-	IssmDouble       omega,h,mu0,ds[2],Cmu,B;
-	const int        n = 3.;
+	IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k;
+	const int        n = 3;
 	const IssmDouble Hstar = 500.;
 	const IssmDouble Lstar = 500.e+3;
 	IssmDouble *xyz_list  = NULL;
@@ -182,6 +187,7 @@
 	Input* thickness_input  = element->GetInput(ThicknessEnum);          _assert_(thickness_input);
 	Input* surface_input    = element->GetInput(SurfaceEnum);            _assert_(surface_input);
 	Input* B_input          = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
+	Input* k_input          = element->GetInput(FrictionCoefficientEnum);_assert_(k_input);
 	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
 	/*Calculate damage evolution source term: */
@@ -190,14 +196,15 @@
 		gauss->GaussNode(element->GetElementType(),i);
 		
 		B_input->GetInputValue(&B,gauss);
+		k_input->GetInputValue(&k,gauss);
 		thickness_input->GetInputValue(&h,gauss);
 		surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
 
-		mu0   = pow(2.,(1.-3*n)/(2.*n)) * B;
-		omega = pow(rhog,3) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
-		Cmu   = 0.;
+		mu0   = pow(2.,(1-3*n)/(2.*n)) * B;
+		Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
+		Cmu   = mu0/Hstar/(k*k);
 
-		D[i] = omega*(ds[0]*ds[0]+ds[1]*ds[1])*pow(h,4)*(1./5.*h - Cmu);
+		D[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)*pow(h,n+1)*(1./(n+2.)*h - Cmu);
 	}
 
 	/*Add input*/
