Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp	(revision 11326)
+++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp	(revision 11327)
@@ -10,10 +10,11 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 
-void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){
+void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,double kmax){
 	
-	int      connectivity;
+	int      i,connectivity;
 	int      numberofdofspernode;
 	int      fsize,configuration_type;
 	Element *element = NULL;
+	Load    *load    = NULL;
 	Mat      Jff     = NULL;
 
@@ -31,7 +32,12 @@
 	
 	/*Create and assemble matrix*/
-	for(int i=0;i<elements->Size();i++){
+	for(i=0;i<elements->Size();i++){
 		element=(Element*)elements->GetObjectByOffset(i);
 		element->CreateJacobianMatrix(Jff);
+	}
+	for (i=0;i<loads->Size();i++){
+		load=(Load*)loads->GetObjectByOffset(i);
+		if(load->InAnalysis(configuration_type)) load->CreateJacobianMatrix(Jff);
+		if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax);
 	}
 	MatAssemblyBegin(Jff,MAT_FINAL_ASSEMBLY);
Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h	(revision 11326)
+++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h	(revision 11327)
@@ -10,5 +10,5 @@
 
 /* local prototypes: */
-void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters);
+void CreateJacobianMatrixx(Mat* pJff,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,double kmax);
 
 #endif  /* _CREATEJACOBIANMATRIXX_H */
Index: /issm/trunk-jpl/src/c/objects/Loads/Icefront.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Icefront.h	(revision 11326)
+++ /issm/trunk-jpl/src/c/objects/Loads/Icefront.h	(revision 11327)
@@ -76,4 +76,6 @@
 		void  CreateKMatrix(Mat Kff, Mat Kfs);
 		void  CreatePVector(Vec pf);
+		void  CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
+		void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");};
 		void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
 		void  PenaltyCreatePVector(Vec pf, double kmax);
Index: /issm/trunk-jpl/src/c/objects/Loads/Load.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Load.h	(revision 11326)
+++ /issm/trunk-jpl/src/c/objects/Loads/Load.h	(revision 11327)
@@ -29,4 +29,6 @@
 		virtual void  CreateKMatrix(Mat Kff, Mat Kfs)=0;
 		virtual void  CreatePVector(Vec pf)=0;
+		virtual void  CreateJacobianMatrix(Mat Jff)=0;
+		virtual void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax)=0;
 		virtual void  PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax)=0;
 		virtual void  PenaltyCreatePVector(Vec pf, double kmax)=0;
Index: /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.h	(revision 11326)
+++ /issm/trunk-jpl/src/c/objects/Loads/Numericalflux.h	(revision 11327)
@@ -72,4 +72,6 @@
 		void  CreateKMatrix(Mat Kff, Mat Kfs);
 		void  CreatePVector(Vec pf);
+		void  CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
+		void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");};
 		void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
 		void  PenaltyCreatePVector(Vec pf, double kmax);
Index: /issm/trunk-jpl/src/c/objects/Loads/Pengrid.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Pengrid.h	(revision 11326)
+++ /issm/trunk-jpl/src/c/objects/Loads/Pengrid.h	(revision 11327)
@@ -77,4 +77,6 @@
 		void  CreateKMatrix(Mat Kff, Mat Kfs);
 		void  CreatePVector(Vec pf);
+		void  CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
+		void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");};
 		void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
 		void  PenaltyCreatePVector(Vec pf, double kmax);
Index: /issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp	(revision 11326)
+++ /issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp	(revision 11327)
@@ -209,4 +209,12 @@
 void  Penpair::CreatePVector(Vec pf){
 
+	/*No loads applied, do nothing: */
+	return;
+
+}
+/*}}}1*/
+/*FUNCTION Penpair::CreateJacobianMatrix{{{1*/
+void  Penpair::CreateJacobianMatrix(Mat Jff){
+	/*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/
 	/*No loads applied, do nothing: */
 	return;
@@ -244,4 +252,9 @@
 	/*No loads applied, do nothing: */
 	return;
+}
+/*}}}1*/
+/*FUNCTION Penpair::PenaltyCreateJacobianMatrix{{{1*/
+void  Penpair::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){
+	this->PenaltyCreateKMatrix(Jff,NULL,kmax);
 }
 /*}}}1*/
Index: /issm/trunk-jpl/src/c/objects/Loads/Penpair.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Penpair.h	(revision 11326)
+++ /issm/trunk-jpl/src/c/objects/Loads/Penpair.h	(revision 11327)
@@ -64,6 +64,8 @@
 		void  CreateKMatrix(Mat Kff, Mat Kfs);
 		void  CreatePVector(Vec pf);
-		void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
+		void  CreateJacobianMatrix(Mat Jff);
+		void  PenaltyCreateKMatrix(Mat Kff,Mat Kfs,double kmax);
 		void  PenaltyCreatePVector(Vec pf, double kmax);
+		void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax);
 		bool  InAnalysis(int analysis_type);
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/objects/Loads/Riftfront.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Loads/Riftfront.h	(revision 11326)
+++ /issm/trunk-jpl/src/c/objects/Loads/Riftfront.h	(revision 11327)
@@ -84,4 +84,6 @@
 		void  CreateKMatrix(Mat Kff, Mat Kfs);
 		void  CreatePVector(Vec pf);
+		void  CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
+		void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");};
 		void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
 		void  PenaltyCreatePVector(Vec pf, double kmax);
Index: /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp	(revision 11326)
+++ /issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp	(revision 11327)
@@ -251,4 +251,10 @@
 	double* localvalues=NULL;
 
+	/*If Kfs is not provided, call the other function*/
+	if(!Kfs){
+		this->AddToGlobal(Kff);
+		return;
+	}
+
 	/*In debugging mode, check consistency (no NaN, and values not too big)*/
 	this->CheckConsistency();
@@ -300,4 +306,7 @@
 	double* localvalues=NULL;
 
+	/*Check that Jff is not NULL*/
+	_assert_(Jff); 
+
 	/*In debugging mode, check consistency (no NaN, and values not too big)*/
 	this->CheckConsistency();
Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11326)
+++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 11327)
@@ -14,4 +14,8 @@
 
 	/*intermediary: */
+	bool   converged;
+	int    num_unstable_constraints;
+	int    count;
+	double kmax;
 	Mat Kff = NULL, Kfs    = NULL, Jff = NULL;
 	Vec ug  = NULL, old_ug = NULL;
@@ -20,8 +24,4 @@
 	Vec df  = NULL;
 	Vec ys  = NULL;
-	
-	bool converged;
-	int num_unstable_constraints;
-	int count;
 
 	/*parameters:*/
@@ -69,6 +69,6 @@
 
 		/*Prepare next iteration using Newton's method*/
-		CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-		SystemMatricesx(&Kff,NULL,&pf,NULL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		SystemMatricesx(&Kff,NULL,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		VecDuplicate(pf,&pJf);
@@ -78,5 +78,5 @@
 		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters);
 		MatFree(&Jff);            VecFree(&pJf);
-		VecAXPY(uf,+1.,duf);      VecFree(&duf);
+		VecAXPY(uf,1.,duf);      VecFree(&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);