Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17658)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17659)
@@ -90,6 +90,7 @@
 
 	/*Intermediaries*/
-	int         numoutputs;
-	char**      requestedoutputs = NULL;
+	int     fe_FS;
+	int     numoutputs;
+	char**  requestedoutputs = NULL;
 
 	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
@@ -109,4 +110,11 @@
 	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum));
+
+	/*XTH parameters*/
+	iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
+	if(fe_FS==XTaylorHoodEnum){
+		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum));
+		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum));
+	}
 
 	/*Requested outputs*/
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 17658)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 17659)
@@ -19,9 +19,6 @@
 	Vector<IssmDouble>*  df     = NULL;
 	Vector<IssmDouble>*  ys     = NULL;
-	IssmDouble           eps_rel;
-	int  configuration_type,max_nonlinear_iterations;
-
-	IssmDouble r     = .6;
-	IssmDouble theta = 0.; // 0<theta<.5   -> .15<theta<.45
+	IssmDouble           eps_rel,r,theta; // 0<theta<.5   -> .15<theta<.45
+	int                  configuration_type,max_nonlinear_iterations;
 
 	/*Create analysis*/
@@ -32,6 +29,6 @@
 	femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
 	femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
-	femmodel->parameters->SetParam(.6,AugmentedLagrangianREnum);
-	femmodel->parameters->SetParam(theta,AugmentedLagrangianThetaEnum);
+	femmodel->parameters->FindParam(&r,AugmentedLagrangianREnum);
+	femmodel->parameters->FindParam(&theta,AugmentedLagrangianThetaEnum);
 
 	/*Update constraints and initialize d and tau if necessary*/
Index: /issm/trunk-jpl/src/m/classes/flowequation.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 17658)
+++ /issm/trunk-jpl/src/m/classes/flowequation.m	(revision 17659)
@@ -14,4 +14,6 @@
 		fe_HO            = '';
 		fe_FS            = '';
+		XTH_r            = 1.;
+		XTH_theta        = 0.;
 		vertex_equation  = NaN;
 		element_equation = NaN;
@@ -88,4 +90,6 @@
 				md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P1xP3','P2xP4'});
 				md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z'});
+				md = checkfield(md,'fieldname','flowequation.XTH_r','numel',[1],'>',0.);
+				md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',0.5);
 				md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices 1],'values',[0 1]);
 				md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices 1],'values',[0 1]);
@@ -143,4 +147,6 @@
 			WriteData(fid,'enum',FlowequationFeHOEnum() ,'data',StringToEnum(obj.fe_HO) ,'format','Integer');
 			WriteData(fid,'enum',FlowequationFeFSEnum() ,'data',StringToEnum(obj.fe_FS) ,'format','Integer');
+			WriteData(fid,'enum',AugmentedLagrangianREnum() ,'data',obj.XTH_r ,'format','Double');
+			WriteData(fid,'enum',AugmentedLagrangianThetaEnum() ,'data',obj.XTH_theta ,'format','Double');
 			WriteData(fid,'object',obj,'fieldname','borderSSA','format','DoubleMat','mattype',1);
 			WriteData(fid,'object',obj,'fieldname','borderHO','format','DoubleMat','mattype',1);
Index: /issm/trunk-jpl/src/m/classes/flowequation.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 17658)
+++ /issm/trunk-jpl/src/m/classes/flowequation.py	(revision 17659)
@@ -25,4 +25,6 @@
 		self.fe_HO            = ''
 		self.fe_FS            = ''
+		self.XTH_r            = 1.
+		self.XTH_theta        = 0.
 		self.vertex_equation  = float('NaN')
 		self.element_equation = float('NaN')
@@ -80,4 +82,6 @@
 			md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices],'values',[0,1])
 			md = checkfield(md,'fieldname','flowequation.borderFS','size',[md.mesh.numberofvertices],'values',[0,1])
+			md = checkfield(md,'fieldname','flowequation.XTH_r','numel',[1],'>',0.)
+			md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',.5)
 			if m.strcmp(md.mesh.meshtype(),'2Dhorizontal'):
 				md = checkfield(md,'fieldname','flowequation.vertex_equation','size',[md.mesh.numberofvertices],'values',[1,2])
@@ -107,4 +111,6 @@
 		WriteData(fid,'enum',FlowequationFeHOEnum() ,'data',StringToEnum(self.fe_HO)[0] ,'format','Integer')
 		WriteData(fid,'enum',FlowequationFeFSEnum() ,'data',StringToEnum(self.fe_FS)[0] ,'format','Integer')
+		WriteData(fid,'enum',AugmentedLagrangianREnum() ,'data',self.XTH_r ,'format','Double')
+		WriteData(fid,'enum',AugmentedLagrangianThetaEnum() ,'data',self.XTH_theta ,'format','Double')
 		WriteData(fid,'object',self,'fieldname','borderSSA','format','DoubleMat','mattype',1)
 		WriteData(fid,'object',self,'fieldname','borderHO','format','DoubleMat','mattype',1)