Index: /issm/trunk/src/m/solutions/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/diagnostic_core.m	(revision 5819)
+++ /issm/trunk/src/m/solutions/diagnostic_core.m	(revision 5820)
@@ -47,5 +47,5 @@
 	end
 			
-	if ismacayealpattyn,
+	if xor(ismacayealpattyn,isstokes), %if macayealpattyn and stokes, we must use stokescoupling
 
 		displaystring(verbose,'\n%s',['computing horizontal velocities...']);
@@ -54,17 +54,16 @@
 	end
 	
+	if ismacayealpattyn & isstokes,
 
-	if dim==3 & ~isstokes,
+		displaystring(verbose,'\n%s',['computing coupling velocities for macayealpattyn and stokes...']);
+		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
+		femmodel=solver_stokescoupling_nonlinear(femmodel,conserve_loads); 
+	end
+
+	if dim==3 & (ismacayealpattyn | ishutter),
 	
 		displaystring(verbose,'\n%s',['computing vertical velocities...']);
 		femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
 		femmodel=solver_linear(femmodel);
-	end
-
-	if isstokes,
-
-		displaystring(verbose,'\n%s',['computing stokes velocities and pressure...']);
-		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
-		femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads); 
 	end
 
Index: /issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m	(revision 5819)
+++ /issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m	(revision 5820)
@@ -10,27 +10,20 @@
 	%First get ug=ug_horiz+ug_vert
 	femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
-	ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
+	ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
 
 	femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
-	ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
-
-	ug=ug_horiz+ug_vert;
+	ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
 
 	while(~converged),
 
-		old_ug=ug;
-		
 		%First compute the horizontal velocity
 		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
 
-		%keep a copy of loads for now
-		loads=femmodel.loads;
-
 		%Start non-linear iteration using input velocity: 
-		ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
+		ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
 		uf_horiz=Reducevectorgtof( ug_horiz, femmodel.nodesets,femmodel.parameters);
 
 		%Update the solution to make sure that vx and vxold are similar
-		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug_horiz);
+		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);
 
 		%save pointer to old velocity
@@ -38,5 +31,5 @@
 		old_uf_horiz=uf_horiz;
 
-		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
+		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 
 		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); 
@@ -47,8 +40,8 @@
 		ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters); 
 
-		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug_horiz);
+		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);
 
 		%Figure out if convergence have been reached
-		converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);
+		converged=convergence(K_ff,p_f,uf_horiz,old_uf_horiz,femmodel.parameters);
 
 		%Then compute vertical velocity
@@ -66,8 +59,4 @@
 		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert);
 
-		%Finally sum the two velocities and check the convergence
-		ug=ug_horiz+ug_vert;
-		converged=convergence([],[],ug,old_ug,femmodel.parameters);
-
 	end