Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17652)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17653)
@@ -4320,5 +4320,5 @@
 	xDelete<int>(cs_list);
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH(Elements* elements,Parameters* parameters){/*{{{*/
+void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters){/*{{{*/
 
 	/*Intermediaries*/
@@ -4330,5 +4330,4 @@
 	IssmDouble *xyz_list = NULL;
 	IssmDouble  Jdet,r;
-	Gauss*      gauss = NULL;
 
 	parameters->FindParam(&r,AugmentedLagrangianREnum);
@@ -4391,5 +4390,5 @@
 		}
 
-		gauss=element->NewGauss(5);
+		Gauss* gauss=element->NewGauss(5);
 		for(int ig=gauss->begin();ig<gauss->end();ig++){
 			gauss->GaussPoint(ig);
@@ -4522,4 +4521,72 @@
 		}
 
+		/*Clean up*/
+		delete gauss;
+		xDelete<IssmDouble>(xyz_list);
+		xDelete<IssmDouble>(tbasis);
+		xDelete<IssmDouble>(Ke);
+		xDelete<IssmDouble>(pe_xx); xDelete<IssmDouble>(d_xx);
+		xDelete<IssmDouble>(pe_yy); xDelete<IssmDouble>(d_yy);
+		xDelete<IssmDouble>(pe_zz); xDelete<IssmDouble>(d_zz);
+		xDelete<IssmDouble>(pe_xy); xDelete<IssmDouble>(d_xy);
+		xDelete<IssmDouble>(pe_xz); xDelete<IssmDouble>(d_xz);
+		xDelete<IssmDouble>(pe_yz); xDelete<IssmDouble>(d_yz);
+	}
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters){/*{{{*/
+
+	/*Intermediaries*/
+	int         dim,tausize,meshtype;
+	IssmDouble  epsxx,epsyy,epszz,epsxy,epsxz,epsyz,D_scalar;
+	IssmDouble  d_xx,d_yy,d_zz,d_xy,d_xz,d_yz;
+	IssmDouble  sigmapxx,sigmapyy,sigmapzz,sigmapxy,sigmapxz,sigmapyz;
+	IssmDouble  dvx[3],dvy[3],dvz[3];
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  Jdet,r;
+
+	parameters->FindParam(&r,AugmentedLagrangianREnum);
+	parameters->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; tausize = 3; break;
+		case Mesh3DEnum:         dim = 3; tausize = 6; break;
+		case Mesh3DtetrasEnum:   dim = 3; tausize = 6; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+
+		/*Get inputs and parameters*/
+		element->GetVerticesCoordinates(&xyz_list);
+		Input* vx_input=element->GetInput(VxEnum);                 _assert_(vx_input);
+		Input* vy_input=element->GetInput(VyEnum);                 _assert_(vy_input);
+		Input* vz_input;
+		if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+
+		/*Get previous tau*/
+		Input* sigmapxx_input=element->GetInput(DeviatoricStressxxEnum); _assert_(sigmapxx_input);
+		Input* sigmapyy_input=element->GetInput(DeviatoricStressyyEnum); _assert_(sigmapyy_input);
+		Input* sigmapxy_input=element->GetInput(DeviatoricStressxyEnum); _assert_(sigmapxy_input);
+		Input* sigmapzz_input=NULL; Input* sigmapxz_input=NULL; Input* sigmapyz_input=NULL;
+		if(dim==3){
+			sigmapzz_input=element->GetInput(DeviatoricStresszzEnum); _assert_(sigmapzz_input);
+			sigmapxz_input=element->GetInput(DeviatoricStressxzEnum); _assert_(sigmapxz_input);
+			sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
+		}
+
+		/*Get NEW 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);
+		}
+
+		/*Fetch number of nodes and dof for this finite element*/
+		int tnumnodes = element->GetNumberOfVertices();      //Tensors, P1 DG
+
 		/*Update tau accordingly*/
 		IssmDouble* tau_xx = xNew<IssmDouble>(tnumnodes);
@@ -4529,6 +4596,5 @@
 		IssmDouble* tau_xz = NULL;
 		IssmDouble* tau_yz = NULL;
-		delete gauss;
-		gauss = element->NewGauss();
+		Gauss* gauss = element->NewGauss();
 		for(int ig=0;ig<tnumnodes;ig++){
 			gauss->GaussNode(P1DGEnum,ig);
@@ -4559,12 +4625,22 @@
 			}
 
+			/*Get new d*/
+			epsxx_input->GetInputValue(&d_xx,gauss);
+			epsyy_input->GetInputValue(&d_yy,gauss);
+			epsxy_input->GetInputValue(&d_xy,gauss);
+			if(dim==3){
+				epszz_input->GetInputValue(&d_zz,gauss);
+				epsxz_input->GetInputValue(&d_xz,gauss);
+				epsyz_input->GetInputValue(&d_yz,gauss);
+			}
+
 			/*Get d and update tau accordingly*/
-			tau_xx[ig] = sigmapxx + r*(epsxx - d_xx[ig]);
-			tau_yy[ig] = sigmapyy + r*(epsyy - d_yy[ig]);
-			tau_xy[ig] = sigmapxy + r*(epsxy - d_xy[ig]);
+			tau_xx[ig] = sigmapxx + r*(epsxx - d_xx);
+			tau_yy[ig] = sigmapyy + r*(epsyy - d_yy);
+			tau_xy[ig] = sigmapxy + r*(epsxy - d_xy);
 			if(dim==3){
-				tau_zz[ig] = sigmapzz + r*(epszz - d_zz[ig]);
-				tau_xz[ig] = sigmapxz + r*(epsxz - d_xz[ig]);
-				tau_yz[ig] = sigmapyz + r*(epsyz - d_yz[ig]);
+				tau_zz[ig] = sigmapzz + r*(epszz - d_zz);
+				tau_xz[ig] = sigmapxz + r*(epsxz - d_xz);
+				tau_yz[ig] = sigmapyz + r*(epsyz - d_yz);
 			}
 		}
@@ -4583,12 +4659,10 @@
 		delete gauss;
 		xDelete<IssmDouble>(xyz_list);
-		xDelete<IssmDouble>(tbasis);
-		xDelete<IssmDouble>(Ke);
-		xDelete<IssmDouble>(pe_xx); xDelete<IssmDouble>(d_xx);   xDelete<IssmDouble>(tau_xx);
-		xDelete<IssmDouble>(pe_yy); xDelete<IssmDouble>(d_yy);   xDelete<IssmDouble>(tau_yy);
-		xDelete<IssmDouble>(pe_zz); xDelete<IssmDouble>(d_zz);   xDelete<IssmDouble>(tau_zz);
-		xDelete<IssmDouble>(pe_xy); xDelete<IssmDouble>(d_xy);   xDelete<IssmDouble>(tau_xy);
-		xDelete<IssmDouble>(pe_xz); xDelete<IssmDouble>(d_xz);   xDelete<IssmDouble>(tau_xz);
-		xDelete<IssmDouble>(pe_yz); xDelete<IssmDouble>(d_yz);   xDelete<IssmDouble>(tau_yz);
+		xDelete<IssmDouble>(tau_xx);
+		xDelete<IssmDouble>(tau_yy);
+		xDelete<IssmDouble>(tau_zz);
+		xDelete<IssmDouble>(tau_xy);
+		xDelete<IssmDouble>(tau_xz);
+		xDelete<IssmDouble>(tau_yz);
 	}
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 17652)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 17653)
@@ -81,5 +81,6 @@
 		void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
 		void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
-		void InputUpdateFromSolutionFSXTH(Elements* elements,Parameters* parameters);
+		void InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters);
+		void InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters);
 		void InitializeXTH(Elements* elements,Parameters* parameters);
 		/*Coupling*/
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 17652)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 17653)
@@ -22,4 +22,7 @@
 	int  configuration_type,max_nonlinear_iterations;
 
+	IssmDouble r     = .6;
+	IssmDouble theta = 0.; // 0<theta<.5   -> .15<theta<.45
+
 	/*Create analysis*/
 	StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
@@ -30,4 +33,5 @@
 	femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
 	femmodel->parameters->SetParam(.6,AugmentedLagrangianREnum);
+	femmodel->parameters->SetParam(theta,AugmentedLagrangianThetaEnum);
 
 	/*Update constraints and initialize d and tau if necessary*/
@@ -45,4 +49,9 @@
 		delete ug_old;ug_old=ug;
 
+		/*Calculate d*/
+		if(theta){
+			analysis->InputUpdateFromSolutionFSXTH_d(femmodel->elements,femmodel->parameters);
+		}
+
 		/*Solve KU=F*/
 		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
@@ -57,5 +66,6 @@
 
 		/*Update d and tau accordingly*/
-		analysis->InputUpdateFromSolutionFSXTH(femmodel->elements,femmodel->parameters);
+		analysis->InputUpdateFromSolutionFSXTH_d(  femmodel->elements,femmodel->parameters);
+		analysis->InputUpdateFromSolutionFSXTH_tau(femmodel->elements,femmodel->parameters);
 
 		/*Check for convergence*/