Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 19112)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 19113)
@@ -10,40 +10,62 @@
 
 	/*Intermediary*/
-	int count;
-	IssmDouble yts;
-
-	/*Fetch parameters: */
-	iomodel->Constant(&yts,ConstantsYtsEnum);
+	bool        isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
+	int         Mz,Nz;
+	IssmDouble *spcvz = NULL;
 
 	/*return if not 3d mesh*/
 	if(iomodel->domaintype!=Domain3DEnum) return;
 
-	/*Fetch data: */
-	iomodel->FetchData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum);
-
-	/*Initialize counter*/
-	count=0;
-
-	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
-	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(FlowequationBorderFSEnum)[i])){
-				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceVerticalAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for FS
-				count++;
-			}
-			else if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvzEnum)[i])){
-				constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,
-								iomodel->Data(StressbalanceSpcvzEnum)[i],StressbalanceVerticalAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
-				count++;
-
-			}
-		} 
-	}
-
-	/*Free data: */
-	iomodel->DeleteData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum);
+	/*Fetch parameters: */
+	iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
+	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
+	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
+	iomodel->Constant(&isHO,FlowequationIsHOEnum);
+	iomodel->Constant(&isFS,FlowequationIsFSEnum);
+
+	/*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, use P1 elements only*/
+	if(!iscoupling){
+		IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvzEnum,StressbalanceVerticalAnalysisEnum,P1Enum,0);
+	}
+	else{
+		/*Fetch data: */
+		iomodel->FetchData(1,FlowequationBorderFSEnum);
+		/*Fetch Spc*/
+		iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
+		if(Nz>1) _error_("not supported yet (needs to be coded)");
+
+		/*Initialize counter*/
+		int count=0;
+
+		/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+		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(FlowequationBorderFSEnum)[i])){
+					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceVerticalAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for FS
+					count++;
+				}
+				else if (!xIsNan<IssmDouble>(spcvz[i])){
+					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,
+									spcvz[i],StressbalanceVerticalAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+					count++;
+
+				}
+			} 
+		}
+
+		/*Free data: */
+		iomodel->DeleteData(1,FlowequationBorderFSEnum);
+		iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
+	}
 
 }/*}}}*/