Index: /issm/trunk/src/m/solvers/solver_couplinstokes_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_couplinstokes_nonlinear.m	(revision 5811)
+++ /issm/trunk/src/m/solvers/solver_couplinstokes_nonlinear.m	(revision 5811)
@@ -0,0 +1,74 @@
+function femmodel=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
+%SOLVER_COUPLINGSTOKES_NONLINEAR - core solver of coupling run
+%
+%   Usage:
+%      [femmodel]=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
+
+	%initialize solution vector
+	converged=0; count=1;
+
+	%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);
+
+	femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
+	ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
+
+	ug=ug_horiz+ug_vert;
+
+	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);
+		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);
+
+		%save pointer to old velocity
+		old_ug_horiz=ug_horiz;
+		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_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); 
+		p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,femmodel.parameters);
+
+		uf_horiz=Solver(K_ff,p_f,old_uf_horiz,femmodel.parameters);
+
+		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);
+
+		%Figure out if convergence have been reached
+		converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);
+
+		%Then compute vertical velocity
+		femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
+		[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); 
+		p_f = Reduceloadfromgtof( p_g,  K_fs, femmodel.ys, femmodel.nodesets,femmodel.parameters);
+
+		displaystring(femmodel.parameters.Verbose>1,'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
+		uf_vert=Solver(K_ff,p_f,[],femmodel.parameters);
+
+		ug_vert= Mergesolutionfromftog( uf_vert, femmodel.ys, femmodel.nodesets,femmodel.parameters); 
+
+		[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
+
+end