function femmodel=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
%SOLVER_COUPLINGSTOKES_NONLINEAR - core solver of coupling run
%
%   Usage:
%      [femmodel]=solver_couplingstokes_nonlinear(femmodel,conserve_loads)

	%Get parameters
	kffpartitioning=femmodel.parameters.Kff;

	%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, femmodel.loads, femmodel.materials, femmodel.parameters);

	femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
	ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);

	while(~converged),

		%First compute the horizontal velocity
		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);

		%Start non-linear iteration using input velocity: 
		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,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);

		%save pointer to old velocity
		old_ug_horiz=ug_horiz;
		old_uf_horiz=uf_horiz;

		if kffpartitioning,
			[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
			p_f = Reduceload( p_f, K_fs, femmodel.ys,femmodel.parameters);
		else
			[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);
		end

		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,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);

		%Figure out if convergence have been reached
		converged=convergence(K_ff,p_f,uf_horiz,old_uf_horiz,femmodel.parameters);
		converged=1;

%		%Then compute vertical velocity
%		femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
%	if kffpartitioning,
%		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
%		p_f = Reduceload( p_f, K_fs, femmodel.ys,femmodel.parameters);
%	else
%		[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);
%	end
%
%		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);

	end
end
