Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 6023)
+++ /issm/trunk/src/c/Makefile.am	(revision 6024)
@@ -1149,4 +1149,5 @@
 					./solvers/solver_adjoint_linear.cpp\
 					./solvers/solver_diagnostic_nonlinear.cpp\
+					./solvers/solver_stokescoupling_nonlinear.cpp\
 					./solvers/solver_thermal_nonlinear.cpp\
 					./modules/Bamgx/Bamgx.cpp\
Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 6023)
+++ /issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 6024)
@@ -63,5 +63,5 @@
 	}
 
-	if (ismacayealpattyn){
+	if (ismacayealpattyn^isstokes){
 		
 		if(verbose)_printf_("%s\n"," computing horizontal velocities...");
@@ -70,6 +70,11 @@
 	}
 	
-	
-	if (dim==3 & ~isstokes){
+	if (ismacayealpattyn && isstokes){
+
+		if(verbose)_printf_("%s\n"," computing coupling macayealpattyn and stokes velocities and pressure ...");
+		solver_stokescoupling_nonlinear(femmodel,conserve_loads);
+	}
+
+	if (dim==3 & (ishutter || ismacayealpattyn)){
 
 		if(verbose)_printf_("%s\n"," computing vertical velocities...");
@@ -78,10 +83,4 @@
 	}
 
-	if (isstokes){
-
-		if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
-		femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
-		solver_diagnostic_nonlinear(femmodel,conserve_loads);
-	}
 
 	if(solution_type==DiagnosticSolutionEnum && !control_analysis){
Index: /issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp	(revision 6024)
+++ /issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp	(revision 6024)
@@ -0,0 +1,109 @@
+/*!\file: solver_stokescoupling_nonlinear.cpp
+ * \brief: core of the coupling between stokes and macayealpattyn
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../modules/modules.h"
+#include "../solutions/solutions.h"
+#include "./solvers.h"
+
+void solver_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads){
+
+	/*intermediary: */
+	Mat Kgg_horiz = NULL, Kff_horiz = NULL, Kfs_horiz   = NULL;
+	Vec ug_horiz  = NULL, uf_horiz  = NULL, old_ug_horiz= NULL, old_uf_horiz = NULL;
+	Vec pg_horiz  = NULL, pf_horiz  = NULL;
+	Mat Kgg_vert  = NULL, Kff_vert  = NULL, Kfs_vert    = NULL;
+	Vec ug_vert   = NULL, uf_vert   = NULL;
+	Vec pg_vert   = NULL, pf_vert   = NULL;
+	int converged;
+	int constraints_converged;
+	int num_unstable_constraints;
+	int count;
+
+	/*parameters:*/
+	int  verbose=0;
+	bool kffpartitioning=false;
+	int min_mechanical_constraints;
+	int max_nonlinear_iterations;
+
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
+	femmodel->parameters->FindParam(&kffpartitioning,KffEnum);
+	femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum);
+	femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum);
+	
+	count=1;
+	converged=0;
+
+	/*First get ug_horiz:*/
+	femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
+	GetSolutionFromInputsx(&ug_horiz, femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters);
+	Reducevectorgtofx(&uf_horiz, ug_horiz, femmodel->nodesets,femmodel->parameters);
+
+
+	for(;;){
+
+
+		/*First diagnostic horiz:*/
+		femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
+		
+		//Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
+		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz);
+
+		//save pointer to old velocity
+		VecFree(&old_ug_horiz);old_ug_horiz=ug_horiz;
+		VecFree(&old_uf_horiz);old_uf_horiz=uf_horiz;
+
+		if(kffpartitioning){
+			SystemMatricesx(NULL,&Kff_horiz, &Kfs_horiz, NULL,&pf_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			Reduceloadx(pf_horiz, Kfs_horiz, femmodel->ys,femmodel->parameters); MatFree(&Kfs_horiz);
+		}
+		else{
+			SystemMatricesx(&Kgg_horiz, NULL, NULL, &pg_horiz,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			Reducematrixfromgtofx(&Kff_horiz,&Kfs_horiz,Kgg_horiz,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg_horiz);
+			Reduceloadfromgtofx(&pf_horiz, pg_horiz, Kfs_horiz, femmodel->ys, femmodel->nodesets,femmodel->parameters); VecFree(&pg_horiz); MatFree(&Kfs_horiz);
+		}
+		
+		Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, femmodel->parameters);
+
+		Mergesolutionfromftogx(&ug_horiz, uf_horiz,femmodel->ys,femmodel->nodesets,femmodel->parameters);
+
+		InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz);
+
+		convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); MatFree(&Kff_horiz);VecFree(&pf_horiz);
+
+		/*Second compute vertical velocity: */
+		femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum);
+		if(kffpartitioning){
+			SystemMatricesx(NULL,&Kff_vert, &Kfs_vert, NULL,&pf_vert, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			Reduceloadx(pf_vert, Kfs_vert, femmodel->ys,femmodel->parameters); MatFree(&Kfs_vert);
+		}
+		else{
+			SystemMatricesx(&Kgg_vert, NULL, NULL, &pg_vert,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			Reducematrixfromgtofx(&Kff_vert,&Kfs_vert,Kgg_vert,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg_vert);
+			Reduceloadfromgtofx(&pf_vert, pg_vert, Kfs_vert, femmodel->ys, femmodel->nodesets,femmodel->parameters);VecFree(&pg_vert); MatFree(&Kfs_vert);
+		}
+
+		Solverx(&uf_vert, Kff_vert, pf_vert, NULL, femmodel->parameters); MatFree(&Kff_vert); VecFree(&pf_vert);
+		Mergesolutionfromftogx(&ug_vert, uf_vert,femmodel->ys,femmodel->nodesets,femmodel->parameters);VecFree(&uf_vert);
+
+		InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_vert); VecFree(&ug_vert); VecFree(&uf_vert);
+
+		/*Increase count: */
+		count++;
+		if(converged==1)break;
+		if(count>=max_nonlinear_iterations){
+			_printf_("   maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); 
+			break;
+		}
+	}
+
+	/*clean-up*/
+	VecFree(&uf_horiz);
+	VecFree(&ug_horiz);
+	VecFree(&old_uf_horiz);
+	VecFree(&old_ug_horiz);
+}
Index: /issm/trunk/src/c/solvers/solvers.h
===================================================================
--- /issm/trunk/src/c/solvers/solvers.h	(revision 6023)
+++ /issm/trunk/src/c/solvers/solvers.h	(revision 6024)
@@ -14,4 +14,5 @@
 void solver_thermal_nonlinear(FemModel* femmodel);
 void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads);
+void solver_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solver_linear(FemModel* femmodel);
 void solver_adjoint_linear(FemModel* femmodel);