Index: /issm/trunk-jpl/src/m/solvers/solver_newton.m
===================================================================
--- /issm/trunk-jpl/src/m/solvers/solver_newton.m	(revision 11336)
+++ /issm/trunk-jpl/src/m/solvers/solver_newton.m	(revision 11336)
@@ -0,0 +1,62 @@
+function femmodel=solver_newton(femmodel)
+%SOLVER_NEWTON - core solver of diagnostic run
+%
+%   Usage:
+%      [femmodel]=solver_newton(femmodel)
+
+	%Branch on partitioning schema requested
+	maxiter=femmodel.parameters.DiagnosticMaxiter;
+	configuration_type=femmodel.parameters.ConfigurationType;
+	[femmodel.nodes]=UpdateConstraints(femmodel.nodes,femmodel.constraints,femmodel.parameters);
+
+	%initialize solution vector
+	converged=0; count=1;
+
+	%Start non-linear iteration using input velocity: 
+	ug=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
+	uf=Reducevectorgtof( ug, femmodel.nodes,femmodel.parameters);
+
+	%Update the solution to make sure that vx and vxold are similar
+	[femmodel.elements femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);
+	[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug);
+
+	while(~converged),
+
+		%save pointer to old velocity
+		old_ug=ug;
+		old_uf=uf;
+
+		%Solver forward model
+		[K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+		ys=CreateNodalConstraints(femmodel.nodes,configuration_type);
+		p_f=Reduceload( p_f, K_fs, ys);
+		issmprintf(VerboseSolver(),'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
+		uf=Solver(K_ff,p_f,old_uf,df,femmodel.parameters);
+		ug=Mergesolutionfromftog( uf, ys, femmodel.nodes,femmodel.parameters); 
+		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug);
+
+		%Check convergence
+		converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);
+		if(converged==1) break; end
+		if(count>maxiter),
+			issmprintf(true,'%s%i%s','      maximum number of iterations ',maxiter,' exceeded');
+			break;
+		end
+
+		%Prepare next iteration using Newton's method
+		[K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+		ys=CreateNodalConstraints(femmodel.nodes,configuration_type);
+		p_f=Reduceload( p_f, K_fs, ys);
+		pJf=p_f-K_ff*uf;
+		Jff=CreateJacobianMatrix(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,kmax);
+		duf=Solver(Jff,pJf,[],[],femmodel.parameters);
+		ys=0*ys;
+		uf=uf+duf;
+		ug=Mergesolutionfromftog(uf,ys,femmodel.nodes,femmodel.parameters); 
+		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug);
+
+		%increase count
+		count=count+1;
+
+	end
+end