Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18347)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18348)
@@ -722,4 +722,20 @@
 	}
 	delete gauss;
+}
+/*}}}*/
+void       Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,int input_enum){/*{{{*/
+
+	Input* input = GetInput(input_enum);
+	if(!input){
+		_error_("Input "<<EnumToStringx(input_enum)<<" not found");
+	}
+
+	int numnodes = this->GetNumberOfNodes();
+	IssmDouble input_min = input->Min();
+	IssmDouble input_max = input->Max();
+	for(int i=0;i<numnodes;i++){
+		if(min[nodes[i]->Sid()]>input_min) min[nodes[i]->Sid()] = input_min;
+		if(max[nodes[i]->Sid()]<input_max) max[nodes[i]->Sid()] = input_max;
+	}
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18347)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18348)
@@ -85,4 +85,5 @@
 		void       GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
 		void       GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
+		void       GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,int input_enum);
 		void       GetInputValue(bool* pvalue,int enum_type);
 		void       GetInputValue(int* pvalue,int enum_type);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 18347)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 18348)
@@ -1125,4 +1125,37 @@
 
 }/*}}}*/
+void FemModel::GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,int input_enum){/*{{{*/
+
+	/*Get vector sizes for current configuration*/
+	int configuration_type;
+	this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	int numnodes = this->nodes->NumberOfNodes(configuration_type);
+
+	/*Initialize output vectors*/
+	IssmDouble* uLmin_local = xNew<IssmDouble>(numnodes);
+	IssmDouble* uLmax_local = xNew<IssmDouble>(numnodes);
+	IssmDouble* uLmin = xNew<IssmDouble>(numnodes);
+	IssmDouble* uLmax = xNew<IssmDouble>(numnodes);
+	for(int i=0;i<numnodes;i++){
+		uLmin_local[i] = +1.e+50;
+		uLmax_local[i] = -1.e+50;
+	}
+
+	for(int i=0;i<this->elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->GetInputLocalMinMaxOnNodes(uLmin_local,uLmax_local,input_enum);
+	}
+
+	/*Synchronize all CPUs*/
+	ISSM_MPI_Allreduce((void*)uLmin_local,(void*)uLmin,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
+	ISSM_MPI_Allreduce((void*)uLmax_local,(void*)uLmax,numnodes,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+	xDelete<IssmDouble>(uLmin_local);
+	xDelete<IssmDouble>(uLmax_local);
+
+	/*Assign output pointers: */
+	*pmin=uLmin;
+	*pmax=uLmax;
+
+}/*}}}*/
 void FemModel::IceVolumex(IssmDouble* pV){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 18347)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 18348)
@@ -61,4 +61,5 @@
 
 		/*Modules*/ 
+		void GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,int input_enum);
 		void MassFluxx(IssmDouble* presponse);
 		void MaxAbsVxx(IssmDouble* presponse);
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18347)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18348)
@@ -8,4 +8,109 @@
 #include "../modules/modules.h"
 
+#ifdef _HAVE_PETSC_
+void CreateDMatrix(Mat* pD,Mat K){/*{{{*/
+	/*Create D matrix such that:
+	 *
+	 * d_ij = max( -k_ij,0,-k_ji) off diagonal
+	 *
+	 * d_ii = - sum_{i!=j} d_ij for the diagonal
+	 *
+	 */
+
+	/*Intermediaries*/
+	int        ncols,ncols2,rstart,rend;
+	double     d,diagD;
+	Mat        D        = NULL;
+	Mat        K_transp = NULL;
+	int*       cols  = NULL;
+	int*       cols2 = NULL;
+	double*    vals  = NULL;
+	double*    vals2 = NULL;
+
+	/*First, we need to transpose K so that we access both k_ij and k_ji*/
+	MatTranspose(K,MAT_INITIAL_MATRIX,&K_transp);
+
+	/*Initialize output (D has the same non zero pattern as K)*/
+	MatDuplicate(K,MAT_SHARE_NONZERO_PATTERN,&D);
+
+	/*Go through the rows of K an K' and build D*/
+	MatGetOwnershipRange(K,&rstart,&rend);
+	for(int row=rstart; row<rend; row++){
+		diagD = 0.;
+		MatGetRow(K       ,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatGetRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+		_assert_(ncols==ncols2);
+		for(int j=0; j<ncols; j++) {
+			_assert_(cols[j]==cols2[j]);
+			d = max(max(-vals[j],-vals2[j]),0.);
+			MatSetValue(D,row,cols[j],(const double)d,INSERT_VALUES);
+			if(cols[j]!=row) diagD -= d;
+		}
+		MatSetValue(D,row,row,(const double)diagD,INSERT_VALUES);
+		MatRestoreRow(K       ,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatRestoreRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+	}
+	MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY);
+	MatAssemblyEnd(  D,MAT_FINAL_ASSEMBLY);
+
+	/*Clean up and assign output pointer*/
+	MatFree(&K_transp);
+	*pD = D;
+}/*}}}*/
+void CreateLHS(Mat* pLHS,IssmDouble* pdmax,Mat K,Mat D,Vec Ml,IssmDouble theta,IssmDouble deltat,FemModel* femmodel,int configuration_type){/*{{{*/
+	/*Create Left Hand side of Lower order solution
+	 *
+	 * LHS = [ML − theta*detlat *(K+D)^n+1]
+	 *
+	 */
+
+	/*Intermediaries*/
+	int        dof,ncols,ncols2,rstart,rend;
+	double     d,mi,dmax = 0.;
+	Mat        LHS   = NULL;
+	int*       cols  = NULL;
+	int*       cols2 = NULL;
+	double*    vals  = NULL;
+	double*    vals2 = NULL;
+
+	MatDuplicate(K,MAT_SHARE_NONZERO_PATTERN,&LHS);
+	MatGetOwnershipRange(K,&rstart,&rend);
+	for(int row=rstart; row<rend; row++){
+		MatGetRow(K,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatGetRow(D,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+		_assert_(ncols==ncols2);
+		for(int j=0; j<ncols; j++) {
+			_assert_(cols[j]==cols2[j]);
+			d = -theta*deltat*(vals[j] + vals2[j]);
+			if(cols[j]==row){
+				VecGetValues(Ml,1,(const int*)&cols[j],&mi);
+				d += mi;
+			}
+			if(fabs(d)>dmax) dmax = fabs(d);
+			MatSetValue(LHS,row,cols[j],(const double)d,INSERT_VALUES);
+		}
+		MatRestoreRow(K,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatRestoreRow(D,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+	}
+
+	/*Penalize Dirichlet boundary*/
+	dmax = dmax * 1.e+3;
+	for(int i=0;i<femmodel->constraints->Size();i++){
+		Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+		if(constraint->InAnalysis(configuration_type)){
+			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+			if(dof!=-1){
+				MatSetValue(LHS,dof,dof,(const double)dmax,INSERT_VALUES);
+			}
+		}
+	}
+	MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
+	MatAssemblyEnd(  LHS,MAT_FINAL_ASSEMBLY);
+
+	/*Clean up and assign output pointer*/
+	*pdmax = dmax;
+	*pLHS  = LHS;
+}/*}}}*/
+#endif
 void solutionsequence_fct(FemModel* femmodel){
 
@@ -18,9 +123,8 @@
 	Vector<IssmDouble>*  ys = NULL;
 
-	IssmDouble theta,deltat;
+	IssmDouble theta,deltat,dmax=0.;
 	int        dof,ncols,ncols2,ncols3,rstart,rend;
 	int        configuration_type,analysis_type;
-	double     d,diagD,mi;
-	double     dmax = 0.;
+	double     d,mi;
 	int*       cols  = NULL;
 	int*       cols2 = NULL;
@@ -45,7 +149,5 @@
 	analysis->FctKMatrix(&K,NULL,femmodel);
 
-	/*Create D Matrix*/
-	#ifdef _HAVE_PETSC_
-	Mat K_transp = NULL;
+	/*Convert matrices to PETSc matrices*/
 	Mat D_petsc  = NULL;
 	Mat LHS      = NULL;
@@ -53,28 +155,12 @@
 	Vec Ml_petsc = Ml->pvector->vector;
 	Mat Mc_petsc = Mc->pmatrix->matrix;
-	MatTranspose(K_petsc,MAT_INITIAL_MATRIX,&K_transp);
-	MatDuplicate(K_petsc,MAT_SHARE_NONZERO_PATTERN,&D_petsc);
-	MatGetOwnershipRange(K_transp,&rstart,&rend);
-	for(int row=rstart; row<rend; row++){
-		diagD = 0.;
-		MatGetRow(K_petsc ,row,&ncols, (const int**)&cols, (const double**)&vals);
-		MatGetRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-		_assert_(ncols==ncols2);
-		for(int j=0; j<ncols; j++) {
-			_assert_(cols[j]==cols2[j]);
-			d = max(max(-vals[j],-vals2[j]),0.);
-			MatSetValues(D_petsc,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
-			if(cols[j]!=row) diagD -= d;
-		}
-		MatSetValues(D_petsc,1,&row,1,&row,(const double*)&diagD,INSERT_VALUES);
-		MatRestoreRow(K_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
-		MatRestoreRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-	}
-	MatAssemblyBegin(D_petsc,MAT_FINAL_ASSEMBLY);
-	MatAssemblyEnd(  D_petsc,MAT_FINAL_ASSEMBLY);
-	MatFree(&K_transp);
+
+	/*Create D Matrix*/
+	#ifdef _HAVE_PETSC_
+	CreateDMatrix(&D_petsc,K_petsc);
 
 	/*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
 	MatDuplicate(K_petsc,MAT_SHARE_NONZERO_PATTERN,&LHS);
+	MatGetOwnershipRange(K_petsc,&rstart,&rend);
 	for(int row=rstart; row<rend; row++){
 		MatGetRow(K_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);