Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17628)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17629)
@@ -3552,5 +3552,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(Dstar);
-	xDelete<IssmDouble>(d);
+xDelete<IssmDouble>(d);
 	xDelete<IssmDouble>(D);
 	xDelete<IssmDouble>(tau);
@@ -4244,4 +4244,5 @@
 	int         dim,tausize,meshtype;
 	IssmDouble  epsxx,epsyy,epszz,epsxy,epsxz,epsyz,D_scalar;
+	IssmDouble  epsxx_old,epsyy_old,epszz_old,epsxy_old,epsxz_old,epsyz_old;
 	IssmDouble  sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz;
 	IssmDouble  dvx[3],dvy[3],dvz[3],B,n;
@@ -4287,4 +4288,15 @@
 		}
 
+		/*Get previous d*/
+		Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
+		Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
+		Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
+		Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
+		if(dim==3){
+			epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
+			epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
+			epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
+		}
+
 		/*Get tau*/
 		Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
@@ -4312,4 +4324,14 @@
 				sigmapxz_input->GetInputValue(&sigmapxz,gauss);
 				sigmapyz_input->GetInputValue(&sigmapyz,gauss);
+			}
+
+			/*Get previous d*/
+			epsxx_input->GetInputValue(&epsxx_old,gauss);
+			epsyy_input->GetInputValue(&epsyy_old,gauss);
+			epsxy_input->GetInputValue(&epsxy_old,gauss);
+			if(dim==3){
+				epszz_input->GetInputValue(&epszz_old,gauss);
+				epsxz_input->GetInputValue(&epsxz_old,gauss);
+				epsyz_input->GetInputValue(&epsyz_old,gauss);
 			}
 
@@ -4333,5 +4355,5 @@
 			B_input->GetInputValue(&B,gauss);
 			n_input->GetInputValue(&n,gauss);
-			coef1 = (B/2.)*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0
+			coef1 = B*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 = 2 * B/(2* (1/sqrt(2)  )^(n-1)/n ) 
 			coef2 = r;
 			if(dim==2){
@@ -4353,5 +4375,12 @@
 			}
 			IssmDouble dnorm;
-			NewtonSolveDnorm(&dnorm,coef1,coef2,coef3,n);
+			if(dim==2){
+				dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + 2.*epsxy_old*epsxy_old );
+			}
+			else{
+				dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old 
+							+2.*(epsxy_old*epsxy_old + epsxz_old*epsxz_old + epsyz_old*epsyz_old));
+			}
+			NewtonSolveDnorm(&dnorm,coef1,coef2,coef3,n,dnorm);
 
 			/*Create Ke*/
@@ -4412,13 +4441,4 @@
 		IssmDouble* tau_xz = NULL;
 		IssmDouble* tau_yz = NULL;
-		Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
-		Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
-		Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
-		Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
-		if(dim==3){
-			epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
-			epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
-			epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
-		}
 		delete gauss;
 		gauss = element->NewGauss();
Index: /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp	(revision 17628)
+++ /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp	(revision 17629)
@@ -3,5 +3,5 @@
 #include "../Exceptions/exceptions.h"
 
-int NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n){
+int NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n,IssmDouble dnorm){
 	/* solve the following equation using Newton-Raphson
 	 *
@@ -26,5 +26,5 @@
 
 	/*Initial guess*/
-	IssmDouble y1 = 0.;
+	IssmDouble y1 = log10(dnorm);
 
 	while(true){
Index: /issm/trunk-jpl/src/c/shared/Numerics/numerics.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 17628)
+++ /issm/trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 17629)
@@ -38,5 +38,5 @@
 int         cubic(IssmDouble a, IssmDouble b, IssmDouble c, IssmDouble d, double X[3], int *num);
 
-int         NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n);
+int         NewtonSolveDnorm(IssmDouble* pdnorm,IssmDouble c1,IssmDouble c2,IssmDouble c3,IssmDouble n,IssmDouble dnorm);
 
 #endif //ifndef _NUMERICS_H_