Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h	(revision 11321)
+++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h	(revision 11322)
@@ -28,4 +28,5 @@
 	DiagnosticAbstolEnum,
 	DiagnosticIcefrontEnum,
+	DiagnosticIsnewtonEnum,
 	DiagnosticMaxiterEnum,
 	DiagnosticNumRequestedOutputsEnum,
Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 11321)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 11322)
@@ -323,4 +323,6 @@
 					./modules/SystemMatricesx/SystemMatricesx.cpp\
 					./modules/SystemMatricesx/SystemMatricesx.h\
+					./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp\
+					./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h\
 					./modules/ConstraintsStatex/ConstraintsStatex.cpp\
 					./modules/ConstraintsStatex/ConstraintsStatex.h\
@@ -352,5 +354,6 @@
 					./solutions/CorePointerFromSolutionEnum.cpp\
 					./solvers/solver_linear.cpp\
-					./solvers/solver_nonlinear.cpp
+					./solvers/solver_nonlinear.cpp\
+					./solvers/solver_newton.cpp
 #}}}
 #DAKOTA sources  {{{1
Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp	(revision 11322)
+++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp	(revision 11322)
@@ -0,0 +1,42 @@
+/*!\file CreateJacobianMatrixx
+ * \brief: create system matrices (stiffness matrix, loads vector)
+ */
+
+#include "./CreateJacobianMatrixx.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+
+void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){
+	
+	int      connectivity;
+	int      numberofdofspernode;
+	int      fsize,configuration_type;
+	Element *element = NULL;
+	Mat      Jff     = NULL;
+
+	/*Checks*/
+	_assert_(nodes && elements);
+
+	/*Recover some parameters*/
+	parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
+	fsize=nodes->NumberOfDofs(configuration_type,FsetEnum);
+	numberofdofspernode=nodes->MaxNumDofs(configuration_type,GsetEnum);
+
+	/*Initialize Jacobian Matrix*/
+	Jff=NewMat(fsize,fsize,connectivity,numberofdofspernode);
+	
+	/*Create and assemble matrix*/
+	for(int i=0;i<elements->Size();i++){
+		element=(Element*)elements->GetObjectByOffset(i);
+		element->CreateJacobianMatrix(Jff);
+	}
+	MatAssemblyBegin(Jff,MAT_FINAL_ASSEMBLY);
+	MatAssemblyEnd(Jff,MAT_FINAL_ASSEMBLY);
+
+	/*Assign output pointer*/
+	*pJff=Jff;
+}
Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h	(revision 11322)
+++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h	(revision 11322)
@@ -0,0 +1,14 @@
+/*!\file:  CreateJacobianMatrixx.h
+ * \brief header file for degree of freedoms distribution routines.
+ */ 
+
+#ifndef _CREATEJACOBIANMATRIXX_H
+#define _CREATEJACOBIANMATRIXX_H
+
+#include "../../Container/Container.h"
+#include "../../objects/objects.h"
+
+/* local prototypes: */
+void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters);
+
+#endif  /* _CREATEJACOBIANMATRIXX_H */
Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp	(revision 11322)
@@ -32,4 +32,5 @@
 		case DiagnosticAbstolEnum : return "DiagnosticAbstol";
 		case DiagnosticIcefrontEnum : return "DiagnosticIcefront";
+		case DiagnosticIsnewtonEnum : return "DiagnosticIsnewton";
 		case DiagnosticMaxiterEnum : return "DiagnosticMaxiter";
 		case DiagnosticNumRequestedOutputsEnum : return "DiagnosticNumRequestedOutputs";
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 11322)
@@ -40,4 +40,5 @@
 	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticReltolEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticAbstolEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticIsnewtonEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticMaxiterEnum));
 	parameters->AddObject(iomodel->CopyConstantObject(SteadystateReltolEnum));
Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp	(revision 11322)
@@ -30,4 +30,5 @@
 	else if (strcmp(name,"DiagnosticAbstol")==0) return DiagnosticAbstolEnum;
 	else if (strcmp(name,"DiagnosticIcefront")==0) return DiagnosticIcefrontEnum;
+	else if (strcmp(name,"DiagnosticIsnewton")==0) return DiagnosticIsnewtonEnum;
 	else if (strcmp(name,"DiagnosticMaxiter")==0) return DiagnosticMaxiterEnum;
 	else if (strcmp(name,"DiagnosticNumRequestedOutputs")==0) return DiagnosticNumRequestedOutputsEnum;
Index: /issm/trunk-jpl/src/c/modules/modules.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/modules.h	(revision 11321)
+++ /issm/trunk-jpl/src/c/modules/modules.h	(revision 11322)
@@ -108,4 +108,5 @@
 #include "./SurfaceAreax/SurfaceAreax.h"
 #include "./SystemMatricesx/SystemMatricesx.h"
+#include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h"
 #include "./TimeAdaptx/TimeAdaptx.h"
 #include "./TriaSearchx/TriaSearchx.h"
Index: /issm/trunk-jpl/src/c/objects/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Element.h	(revision 11321)
+++ /issm/trunk-jpl/src/c/objects/Elements/Element.h	(revision 11322)
@@ -30,4 +30,5 @@
 		virtual void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df)=0;
 		virtual void   CreatePVector(Vec pf)=0;
+		virtual void   CreateJacobianMatrix(Mat Jff)=0;
 		virtual void   GetSolutionFromInputs(Vec solution)=0;
 		virtual int    GetNodeIndex(Node* node)=0;
Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11322)
@@ -767,4 +767,37 @@
 	/*clean up and return*/
 	return pe;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateJacobianMatrix{{{1*/
+void  Penta::CreateJacobianMatrix(Mat Jff){
+
+	/*retrieve parameters: */
+	ElementMatrix* Ke=NULL;
+	int analysis_type;
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+	/*Checks in debugging {{{2*/
+	_assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
+	/*}}}*/
+
+	/*Skip if water element*/
+	if(IsOnWater()) return;
+
+	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
+	switch(analysis_type){
+#ifdef _HAVE_DIAGNOSTIC_
+		case DiagnosticHorizAnalysisEnum:
+			Ke=CreateJacobianDiagnosticPattyn();
+			break;
+#endif
+		default:
+			_error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
+	}
+
+	/*Add to global matrix*/
+	if(Ke){
+		Ke->AddToGlobal(Jff);
+		delete Ke;
+	}
 }
 /*}}}*/
@@ -6224,4 +6257,71 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/
+ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
+
+	/*Constants*/
+	const int    numdof=NDOF2*NUMVERTICES;
+
+	/*Intermediaries */
+	int        i,j,ig;
+	int        approximation;
+	double     xyz_list[NUMVERTICES][3];
+	double     Jdet;
+	double     eps1dotdphii,eps1dotdphij;
+	double     eps2dotdphii,eps2dotdphij;
+	double     mu_prime;
+	double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+	double     eps1[3],eps2[3];
+	double     phi[NUMVERTICES];
+	double     dphi[3][NUMVERTICES];
+	GaussPenta *gauss=NULL;
+
+	/*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
+	ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
+
+	/*Retrieve all inputs and parameters*/
+	inputs->GetInputValue(&approximation,ApproximationEnum);
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussPenta(5,5);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
+
+		this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
+		eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
+		eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
+		eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
+
+		for(i=0;i<6;i++){
+			for(j=0;j<6;j++){
+				eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
+				eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
+				eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
+				eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
+
+				Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
+				Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
+				Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
+				Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
+			}
+		}
+	}
+
+	/*Transform Coordinate System*/
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+
+	/*Clean up and return*/
+	delete gauss;
+	return Ke;
+}
+/*}}}*/
 /*FUNCTION Penta::CreateKMatrixDiagnosticPattynFriction{{{1*/
 ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){
Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.h	(revision 11321)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h	(revision 11322)
@@ -88,4 +88,5 @@
 		void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df);
 		void   CreatePVector(Vec pf);
+		void   CreateJacobianMatrix(Mat Jff);
 		void   DeleteResults(void);
 		int    GetNodeIndex(Node* node);
@@ -234,4 +235,5 @@
 		ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
 		ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
+		ElementMatrix* CreateJacobianDiagnosticPattyn(void);
 		void           InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
 		void           InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp	(revision 11322)
@@ -1883,4 +1883,5 @@
 				name==InversionVyObsEnum ||
 				name==FrictionCoefficientEnum ||
+				name==MaterialsRheologyBbarEnum ||
 				name==GradientEnum ||
 				name==OldGradientEnum ||
@@ -4734,4 +4735,5 @@
 
 	int doflist1[NUMVERTICES];
+	Input *input=NULL;
 
 	/*Get out if this is not an element input*/
@@ -4742,6 +4744,10 @@
 
 	/*Get input (either in element or material)*/
-	Input* input=inputs->GetInput(control_enum);
-	if(!input) _error_("Input %s not found in element",EnumToStringx(control_enum));
+	if(control_enum==MaterialsRheologyBbarEnum){
+		input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
+	}
+	else{
+		input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
+	}
 
 	/*Check that it is a ControlInput*/
@@ -4772,5 +4778,4 @@
 	}
 	new_input = new TriaP1Input(control_enum,values);
-
 
 	if(control_enum==MaterialsRheologyBbarEnum){
Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Tria.h	(revision 11321)
+++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h	(revision 11322)
@@ -84,4 +84,5 @@
 		void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df);
 		void   CreatePVector(Vec pf);
+		void   CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
 		int    GetNodeIndex(Node* node);
 		int    Sid();
Index: /issm/trunk-jpl/src/c/objects/Materials/Matice.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Materials/Matice.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/objects/Materials/Matice.cpp	(revision 11322)
@@ -521,4 +521,38 @@
 	/*Return: */
 	*pviscosity_complement=viscosity_complement;
+}
+/*}}}*/
+/*FUNCTION Matice::GetViscosityDerivativeEpsSquare{{{1*/
+void  Matice::GetViscosityDerivativeEpsSquare(double* pmu_prime, double* epsilon){
+
+	/*output: */
+	double mu_prime;
+	double mu,n,eff2;
+
+	/*input strain rate: */
+	double exx,eyy,exy,exz,eyz;
+
+	/*Get visocisty and n*/
+	GetViscosity3d(&mu,epsilon);
+	n=GetN();
+
+	if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
+				(epsilon[3]==0) && (epsilon[4]==0)){
+		mu_prime=0.5*pow((double)10,(double)14);
+	}
+	else{
+		/*Retrive strain rate components: */
+		exx=epsilon[0];
+		eyy=epsilon[1];
+		exy=epsilon[2];
+		exz=epsilon[3];
+		eyz=epsilon[4];
+		eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
+
+		mu_prime=(1-n)/(2*n) * mu/eff2;
+	}
+
+	/*Assign output pointers:*/
+	*pmu_prime=mu_prime;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/objects/Materials/Matice.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Materials/Matice.h	(revision 11321)
+++ /issm/trunk-jpl/src/c/objects/Materials/Matice.h	(revision 11322)
@@ -68,4 +68,5 @@
 		void   GetViscosity3dStokes(double* pviscosity3d, double* epsilon);
 		void   GetViscosityComplement(double* pviscosity_complement, double* pepsilon);
+		void   GetViscosityDerivativeEpsSquare(double* pmu_prime, double* pepsilon);
 		double GetB();
 		double GetBbar();
Index: /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp	(revision 11322)
@@ -245,5 +245,5 @@
 
 /*ElementMatrix specific routines: */
-/*FUNCTION ElementMatrix::AddToGlobal{{{1*/
+/*FUNCTION ElementMatrix::AddToGlobal(Mat Kff, Mat Kfs){{{1*/
 void ElementMatrix::AddToGlobal(Mat Kff, Mat Kfs){
 
@@ -294,4 +294,38 @@
 }
 /*}}}*/
+/*FUNCTION ElementMatrix::AddToGlobal(Mat Jff){{{1*/
+void ElementMatrix::AddToGlobal(Mat Jff){
+
+	int i,j;
+	double* localvalues=NULL;
+
+	/*In debugging mode, check consistency (no NaN, and values not too big)*/
+	this->CheckConsistency();
+
+	if(this->dofsymmetrical){
+		/*only use row dofs to add values into global matrices: */
+
+		if(this->row_fsize){
+			/*first, retrieve values that are in the f-set from the g-set values matrix: */
+			localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double));
+			for(i=0;i<this->row_fsize;i++){
+				for(j=0;j<this->row_fsize;j++){
+					*(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]);
+				}
+			}
+			/*add local values into global  matrix, using the fglobaldoflist: */
+			MatSetValues(Jff,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES);
+
+			/*Free ressources:*/
+			xfree((void**)&localvalues);
+		}
+
+	}
+	else{
+		_error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!");
+	}
+
+}
+/*}}}*/
 /*FUNCTION ElementMatrix::CheckConsistency{{{1*/
 void ElementMatrix::CheckConsistency(void){
Index: /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.h	(revision 11321)
+++ /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.h	(revision 11322)
@@ -59,4 +59,5 @@
 		/*ElementMatrix specific routines {{{1*/
 		void AddToGlobal(Mat Kff, Mat Kfs);
+		void AddToGlobal(Mat Jff);
 		void Echo(void);
 		void CheckConsistency(void);
Index: /issm/trunk-jpl/src/c/solutions/controltao_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/solutions/controltao_core.cpp	(revision 11322)
@@ -25,5 +25,5 @@
 
 	/*TAO*/
-	int        ierr,numberofvertices;
+	int        ierr;
 	int        num_controls;
 	AppCtx     user;
@@ -45,17 +45,16 @@
 
 	/*Initialize TAO*/
+	_printf_(VerboseControl(),"%s\n","   Initializing the Toolkit for Advanced Optimization (TAO)");
 	TaoCreate(PETSC_COMM_WORLD,&tao);
 	PetscOptionsSetValue("-tao_monitor","");
 	TaoSetFromOptions(tao);
 	TaoSetType(tao,"tao_blmvm");
+	//TaoSetType(tao,"tao_cg");
+	//TaoSetType(tao,"tao_lmvm");
 
 	/*Prepare all TAO parameters*/
 	TaoSetMaximumFunctionEvaluations(tao,50);
 	TaoSetMaximumIterations(tao,10);
-	TaoSetTolerances(tao,10e-28,0.0000,0.0000,0.0000,10e-28);
-
-	numberofvertices = femmodel->vertices->NumberOfVertices();
-	XL=NewVec(numberofvertices); VecSet(XL,1.);
-	XU=NewVec(numberofvertices); VecSet(XU,200.);
+	TaoSetTolerances(tao,0.,0.,0.,0.,0.);
 
 	GetVectorFromControlInputsx(&X, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
@@ -71,4 +70,5 @@
 
 	/*Solver optimization problem*/
+	_printf_(VerboseControl(),"%s\n","   Starting optimization");
 	TaoSolve(tao);
 	TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
@@ -96,4 +96,5 @@
 
 	/*Set new variable*/
+	//VecView(X,PETSC_VIEWER_STDOUT_WORLD);
 	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
 
@@ -118,4 +119,5 @@
 	/*Compute objective function*/
 	CostFunctionx(fcn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	_printf_(true,"f(x)=%g\n",*fcn);
 
 	/*Compute gradient*/
@@ -123,4 +125,5 @@
 	VecCopy(gradient,G); VecFree(&gradient);
 	VecScale(G,-1.);
+	//VecView(G,PETSC_VIEWER_STDOUT_WORLD);
 
 	/*Clean-up and return*/
Index: /issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp	(revision 11322)
@@ -16,15 +16,16 @@
 
 	/*parameters: */
-	bool   dakota_analysis         = false;
-	int    dim                  = -1;
-	bool   ishutter             = false;
-	bool   ismacayealpattyn     = false;
-	bool   isstokes             = false;
-	bool   conserve_loads       = true;
-	bool   modify_loads         = true;
-	bool   control_analysis;
-	int    solution_type;
-	int    numoutputs         = 0;
-	int    *requested_outputs = NULL;
+	bool  dakota_analysis   = false;
+	int   dim               = -1;
+	bool  ishutter          = false;
+	bool  ismacayealpattyn  = false;
+	bool  isstokes          = false;
+	bool  isnewton          = false;
+	bool  conserve_loads    = true;
+	bool  modify_loads      = true;
+	bool  control_analysis;
+	int   solution_type;
+	int   numoutputs        = 0;
+	int  *requested_outputs = NULL;
 
 	/* recover parameters:*/
@@ -33,4 +34,5 @@
 	femmodel->parameters->FindParam(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
 	femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum);
+	femmodel->parameters->FindParam(&isnewton,DiagnosticIsnewtonEnum);
 	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
 	femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum);
@@ -72,5 +74,8 @@
 		_printf_(VerboseSolution(),"%s\n","   computing velocities");
 		femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
-		solver_nonlinear(femmodel,modify_loads); 
+		if(isnewton)
+		 solver_newton(femmodel); 
+		else
+		 solver_nonlinear(femmodel,conserve_loads); 
 	}
 
Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11322)
+++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11322)
@@ -0,0 +1,91 @@
+/*!\file: solver_nonlinear.cpp
+ * \brief: core of a non-linear solution, using fixed-point method 
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../io/io.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../modules/modules.h"
+#include "../solutions/solutions.h"
+#include "./solvers.h"
+
+void solver_newton(FemModel* femmodel){
+
+	/*intermediary: */
+	Mat Kff = NULL, Kfs    = NULL, Jff = NULL;
+	Vec ug  = NULL, old_ug = NULL;
+	Vec uf  = NULL, old_uf = NULL, duf = NULL;
+	Vec pf  = NULL, pJf    = NULL;
+	Vec df  = NULL;
+	Vec ys  = NULL;
+	
+	bool converged;
+	int num_unstable_constraints;
+	int count;
+
+	/*parameters:*/
+	int max_nonlinear_iterations;
+	int  configuration_type;
+
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&max_nonlinear_iterations,DiagnosticMaxiterEnum);
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
+
+	count=1;
+	converged=false;
+
+	/*Start non-linear iteration using input velocity: */
+	GetSolutionFromInputsx(&ug,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	Reducevectorgtofx(&uf,ug,femmodel->nodes,femmodel->parameters);
+
+	//Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
+	InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
+	InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
+
+	for(;;){
+
+		//save pointer to old velocity
+		VecFree(&old_ug);old_ug=ug;
+		VecFree(&old_uf);old_uf=uf;
+
+		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+		Reduceloadx(pf,Kfs,ys);MatFree(&Kfs);
+		Solverx(&uf,Kff,pf,old_uf,df,femmodel->parameters);
+
+		convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
+
+		CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);
+		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
+		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		VecDuplicate(pf,&pJf);
+		MatMultPatch(Kff,uf,pJf);
+		VecAXPY(pJf,-1.,pf);
+		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters);
+		MatFree(&Kff);VecFree(&pf); VecFree(&df);
+		VecAXPY(uf,-1.,duf);
+		Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
+
+		InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
+
+		/*Increase count: */
+		count++;
+		if(converged==true)break;
+		if(count>=max_nonlinear_iterations){
+			_printf_(true,"   maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); 
+			break;
+		}
+	}
+
+	_printf_(VerboseConvergence(),"\n   total number of iterations: %i\n",count-1);
+
+	/*clean-up*/
+	VecFree(&uf);
+	VecFree(&ug);
+	VecFree(&old_ug);
+	VecFree(&old_uf);
+}
Index: /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp	(revision 11321)
+++ /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp	(revision 11322)
@@ -89,4 +89,6 @@
 	}
 
+	_printf_(VerboseConvergence(),"\n   total number of iterations: %i\n",count-1);
+
 	/*clean-up*/
 	if(conserve_loads) delete loads;
Index: /issm/trunk-jpl/src/c/solvers/solvers.h
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solvers.h	(revision 11321)
+++ /issm/trunk-jpl/src/c/solvers/solvers.h	(revision 11322)
@@ -14,4 +14,5 @@
 void solver_thermal_nonlinear(FemModel* femmodel);
 void solver_nonlinear(FemModel* femmodel,bool conserve_loads);
+void solver_newton(FemModel* femmodel);
 void solver_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solver_linear(FemModel* femmodel);
Index: /issm/trunk-jpl/src/m/classes/diagnostic.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/diagnostic.m	(revision 11321)
+++ /issm/trunk-jpl/src/m/classes/diagnostic.m	(revision 11322)
@@ -12,4 +12,5 @@
 		reltol                   = 0;
 		abstol                   = 0;
+		isnewton                 = 0;
 		stokesreconditioning     = 0;
 		viscosity_overshoot      = 0;
@@ -82,4 +83,5 @@
 			checkfield(md,'diagnostic.reltol','size',[1 1]);
 			checkfield(md,'diagnostic.abstol','size',[1 1]);
+			checkfield(md,'diagnostic.isnewton','numel',1,'values',[0 1]);
 			checkfield(md,'diagnostic.stokesreconditioning','size',[1 1],'NaN',1);
 			checkfield(md,'diagnostic.viscosity_overshoot','size',[1 1],'NaN',1);
@@ -127,4 +129,5 @@
 			fielddisplay(obj,'reltol','velocity relative convergence criterion, NaN -> not applied');
 			fielddisplay(obj,'abstol','velocity absolute convergence criterion, NaN -> not applied');
+			fielddisplay(obj,'isnewton','Apply Newton''s method instead of a Picard fixed point method');
 			fielddisplay(obj,'maxiter','maximum number of nonlinear iterations');
 			fielddisplay(obj,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)');
@@ -158,4 +161,5 @@
 			WriteData(fid,'object',obj,'fieldname','reltol','format','Double');
 			WriteData(fid,'object',obj,'fieldname','abstol','format','Double');
+			WriteData(fid,'object',obj,'fieldname','isnewton','format','Boolean');
 			WriteData(fid,'object',obj,'fieldname','stokesreconditioning','format','Double');
 			WriteData(fid,'object',obj,'fieldname','viscosity_overshoot','format','Double');
Index: /issm/trunk-jpl/src/m/enum/DiagnosticIsnewtonEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/DiagnosticIsnewtonEnum.m	(revision 11322)
+++ /issm/trunk-jpl/src/m/enum/DiagnosticIsnewtonEnum.m	(revision 11322)
@@ -0,0 +1,11 @@
+function macro=DiagnosticIsnewtonEnum()
+%DIAGNOSTICISNEWTONENUM - Enum of DiagnosticIsnewton
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
+%            Please read src/c/EnumDefinitions/README for more information
+%
+%   Usage:
+%      macro=DiagnosticIsnewtonEnum()
+
+macro=StringToEnum('DiagnosticIsnewton');
Index: /issm/trunk-jpl/src/m/enum/PentaP1ElementResultEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/PentaP1ElementResultEnum.m	(revision 11322)
+++ /issm/trunk-jpl/src/m/enum/PentaP1ElementResultEnum.m	(revision 11322)
@@ -0,0 +1,11 @@
+function macro=PentaP1ElementResultEnum()
+%PENTAP1ELEMENTRESULTENUM - Enum of PentaP1ElementResult
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
+%            Please read src/c/EnumDefinitions/README for more information
+%
+%   Usage:
+%      macro=PentaP1ElementResultEnum()
+
+macro=StringToEnum('PentaP1ElementResult');
Index: /issm/trunk-jpl/src/m/enum/PentaP1InputEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/PentaP1InputEnum.m	(revision 11322)
+++ /issm/trunk-jpl/src/m/enum/PentaP1InputEnum.m	(revision 11322)
@@ -0,0 +1,11 @@
+function macro=PentaP1InputEnum()
+%PENTAP1INPUTENUM - Enum of PentaP1Input
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
+%            Please read src/c/EnumDefinitions/README for more information
+%
+%   Usage:
+%      macro=PentaP1InputEnum()
+
+macro=StringToEnum('PentaP1Input');
Index: sm/trunk-jpl/src/m/enum/PentaVertexElementResultEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/PentaVertexElementResultEnum.m	(revision 11321)
+++ 	(revision )
@@ -1,11 +1,0 @@
-function macro=PentaVertexElementResultEnum()
-%PENTAVERTEXELEMENTRESULTENUM - Enum of PentaVertexElementResult
-%
-%   WARNING: DO NOT MODIFY THIS FILE
-%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
-%            Please read src/c/EnumDefinitions/README for more information
-%
-%   Usage:
-%      macro=PentaVertexElementResultEnum()
-
-macro=StringToEnum('PentaVertexElementResult');
Index: sm/trunk-jpl/src/m/enum/PentaVertexInputEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/PentaVertexInputEnum.m	(revision 11321)
+++ 	(revision )
@@ -1,11 +1,0 @@
-function macro=PentaVertexInputEnum()
-%PENTAVERTEXINPUTENUM - Enum of PentaVertexInput
-%
-%   WARNING: DO NOT MODIFY THIS FILE
-%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
-%            Please read src/c/EnumDefinitions/README for more information
-%
-%   Usage:
-%      macro=PentaVertexInputEnum()
-
-macro=StringToEnum('PentaVertexInput');
Index: /issm/trunk-jpl/src/m/enum/TriaP1ElementResultEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/TriaP1ElementResultEnum.m	(revision 11322)
+++ /issm/trunk-jpl/src/m/enum/TriaP1ElementResultEnum.m	(revision 11322)
@@ -0,0 +1,11 @@
+function macro=TriaP1ElementResultEnum()
+%TRIAP1ELEMENTRESULTENUM - Enum of TriaP1ElementResult
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
+%            Please read src/c/EnumDefinitions/README for more information
+%
+%   Usage:
+%      macro=TriaP1ElementResultEnum()
+
+macro=StringToEnum('TriaP1ElementResult');
Index: /issm/trunk-jpl/src/m/enum/TriaP1InputEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/TriaP1InputEnum.m	(revision 11322)
+++ /issm/trunk-jpl/src/m/enum/TriaP1InputEnum.m	(revision 11322)
@@ -0,0 +1,11 @@
+function macro=TriaP1InputEnum()
+%TRIAP1INPUTENUM - Enum of TriaP1Input
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
+%            Please read src/c/EnumDefinitions/README for more information
+%
+%   Usage:
+%      macro=TriaP1InputEnum()
+
+macro=StringToEnum('TriaP1Input');
Index: sm/trunk-jpl/src/m/enum/TriaVertexElementResultEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/TriaVertexElementResultEnum.m	(revision 11321)
+++ 	(revision )
@@ -1,11 +1,0 @@
-function macro=TriaVertexElementResultEnum()
-%TRIAVERTEXELEMENTRESULTENUM - Enum of TriaVertexElementResult
-%
-%   WARNING: DO NOT MODIFY THIS FILE
-%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
-%            Please read src/c/EnumDefinitions/README for more information
-%
-%   Usage:
-%      macro=TriaVertexElementResultEnum()
-
-macro=StringToEnum('TriaVertexElementResult');
Index: sm/trunk-jpl/src/m/enum/TriaVertexInputEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/TriaVertexInputEnum.m	(revision 11321)
+++ 	(revision )
@@ -1,11 +1,0 @@
-function macro=TriaVertexInputEnum()
-%TRIAVERTEXINPUTENUM - Enum of TriaVertexInput
-%
-%   WARNING: DO NOT MODIFY THIS FILE
-%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
-%            Please read src/c/EnumDefinitions/README for more information
-%
-%   Usage:
-%      macro=TriaVertexInputEnum()
-
-macro=StringToEnum('TriaVertexInput');
Index: /issm/trunk-jpl/src/m/model/loadresultsfromdisk.m
===================================================================
--- /issm/trunk-jpl/src/m/model/loadresultsfromdisk.m	(revision 11321)
+++ /issm/trunk-jpl/src/m/model/loadresultsfromdisk.m	(revision 11322)
@@ -10,4 +10,7 @@
 	error('loadresultsfromdisk: error message.');
 end
+if ~exist(filename,'file'),
+	error(['binary file ' filename ' not found.']);
+end
 
 if ~md.qmu.isdakota,
@@ -20,4 +23,7 @@
 	%load results onto model
 	structure=parseresultsfromdisk(filename,~md.settings.io_gather);
+	if isempty(fieldnames(structure)),
+		error(['No result found in binary file ' filename '. Check for solution crash.']);
+	end
 	md.results.(structure(1).SolutionType)=structure;