Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 19290)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 19291)
@@ -9,38 +9,57 @@
 void StressbalanceSIAAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
+	/*Intermediaries*/
+	bool       isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
+
 	/*Fetch parameters: */
-	bool isSIA;
 	iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
+	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
+	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
+	iomodel->Constant(&isHO,FlowequationIsHOEnum);
+	iomodel->Constant(&isFS,FlowequationIsFSEnum);
 
 	/*Now, is the flag isSIA on? otherwise, do nothing: */
 	if (!isSIA) return;
 
-	/*Fetch data: */
-	iomodel->FetchData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);
-
-	/*Initialize conunter*/
-	int count=0;
-
-	/*vx and vy are spc'd if we are not on nodeonSIA: */
-	for(int i=0;i<iomodel->numberofvertices;i++){
-		/*keep only this partition's nodes:*/
-		if((iomodel->my_vertices[i])){
-			if (reCast<int,IssmDouble>(iomodel->Data(FlowequationVertexEquationEnum)[i])!=SIAApproximationEnum){
-
-				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceSIAAnalysisEnum));
-				count++;
-
-				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceSIAAnalysisEnum));
-				count++;
-			}
-			else{
-				if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvxEnum)[i])){
-					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,iomodel->Data(StressbalanceSpcvxEnum)[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+	/*Do we have coupling*/
+	if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
+	 iscoupling = true;
+	else
+	 iscoupling = false;
+
+	/*If no coupling, call Regular IoModelToConstraintsx, else, OLD stuff, keep for now*/
+	if(!iscoupling){
+		IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvxEnum,StressbalanceSIAAnalysisEnum,P1Enum,0);
+		IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvyEnum,StressbalanceSIAAnalysisEnum,P1Enum,1);
+	}
+	else{
+		/*Fetch data: */
+		iomodel->FetchData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);
+
+		/*Initialize conunter*/
+		int count=0;
+
+		/*vx and vy are spc'd if we are not on nodeonSIA: */
+		for(int i=0;i<iomodel->numberofvertices;i++){
+			/*keep only this partition's nodes:*/
+			if((iomodel->my_vertices[i])){
+				if (reCast<int,IssmDouble>(iomodel->Data(FlowequationVertexEquationEnum)[i])!=SIAApproximationEnum){
+
+					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceSIAAnalysisEnum));
+					count++;
+
+					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceSIAAnalysisEnum));
 					count++;
 				}
-
-				if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvyEnum)[i])){
-					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->Data(StressbalanceSpcvyEnum)[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
-					count++;
+				else{
+					if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvxEnum)[i])){
+						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,iomodel->Data(StressbalanceSpcvxEnum)[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+						count++;
+					}
+
+					if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvyEnum)[i])){
+						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->Data(StressbalanceSpcvyEnum)[i],StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
+						count++;
+					}
 				}
 			}
