Index: ../trunk-jpl/src/c/classes/Loads/Numericalflux.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Numericalflux.cpp	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Numericalflux.cpp	(revision 18926)
@@ -130,43 +130,6 @@
 /*}}}*/
 
 /*Object virtual functions definitions:*/
-void Numericalflux::Echo(void){/*{{{*/
-	_printf_("Numericalflux:\n");
-	_printf_("   id: " << id << "\n");
-	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
-	_printf_("   flux_type: " << this->flux_type<< "\n");
-	hnodes->Echo();
-	hvertices->Echo();
-	helement->Echo();
-	_printf_("   parameters: " << parameters << "\n");
-}
-/*}}}*/
-void Numericalflux::DeepEcho(void){/*{{{*/
-
-	_printf_("Numericalflux:\n");
-	_printf_("   id: " << id << "\n");
-	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
-	_printf_("   flux_type: " << this->flux_type<< "\n");
-	hnodes->DeepEcho();
-	hvertices->DeepEcho();
-	helement->DeepEcho();
-	_printf_("   parameters\n");
-	if(parameters)
-	 parameters->DeepEcho();
-	else
-	 _printf_("      NULL\n");
-}		
-/*}}}*/
-int Numericalflux::Id(void){/*{{{*/
-	return id;
-}
-/*}}}*/
-int Numericalflux::ObjectEnum(void){/*{{{*/
-
-	return NumericalfluxEnum;
-
-}
-/*}}}*/
 Object* Numericalflux::copy() {/*{{{*/
 
 	Numericalflux* numericalflux=NULL;
@@ -194,7 +157,44 @@
 	return numericalflux;
 }
 /*}}}*/
+void    Numericalflux::DeepEcho(void){/*{{{*/
 
+	_printf_("Numericalflux:\n");
+	_printf_("   id: " << id << "\n");
+	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
+	_printf_("   flux_type: " << this->flux_type<< "\n");
+	hnodes->DeepEcho();
+	hvertices->DeepEcho();
+	helement->DeepEcho();
+	_printf_("   parameters\n");
+	if(parameters)
+	 parameters->DeepEcho();
+	else
+	 _printf_("      NULL\n");
+}		
+/*}}}*/
+void    Numericalflux::Echo(void){/*{{{*/
+	_printf_("Numericalflux:\n");
+	_printf_("   id: " << id << "\n");
+	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
+	_printf_("   flux_type: " << this->flux_type<< "\n");
+	hnodes->Echo();
+	hvertices->Echo();
+	helement->Echo();
+	_printf_("   parameters: " << parameters << "\n");
+}
+/*}}}*/
+int     Numericalflux::Id(void){/*{{{*/
+	return id;
+}
+/*}}}*/
+int     Numericalflux::ObjectEnum(void){/*{{{*/
+
+	return NumericalfluxEnum;
+
+}
+/*}}}*/
+
 /*Load virtual functions definitions:*/
 void  Numericalflux::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
 
@@ -213,24 +213,6 @@
 	this->parameters=parametersin;
 }
 /*}}}*/
-void  Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
-
-}
-/*}}}*/
-void  Numericalflux::ResetHooks(){/*{{{*/
-
-	this->nodes=NULL;
-	this->vertices=NULL;
-	this->element=NULL;
-	this->parameters=NULL;
-
-	/*Get Element type*/
-	this->hnodes->reset();
-	this->hvertices->reset();
-	this->helement->reset();
-
-}
-/*}}}*/
 void  Numericalflux::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
 
 	/*recover some parameters*/
@@ -290,34 +272,34 @@
 
 }
 /*}}}*/
-void  Numericalflux::GetNodesSidList(int* sidlist){/*{{{*/
+void  Numericalflux::GetNodesLidList(int* lidlist){/*{{{*/
 
-	_assert_(sidlist);
+	_assert_(lidlist);
 	_assert_(nodes);
 
 	switch(this->flux_type){
 		case InternalEnum:
-			for(int i=0;i<NUMNODES_INTERNAL;i++) sidlist[i]=nodes[i]->Sid();
+			for(int i=0;i<NUMNODES_INTERNAL;i++) lidlist[i]=nodes[i]->Lid();
 			return;
 		case BoundaryEnum:
-			for(int i=0;i<NUMNODES_BOUNDARY;i++) sidlist[i]=nodes[i]->Sid();
+			for(int i=0;i<NUMNODES_BOUNDARY;i++) lidlist[i]=nodes[i]->Lid();
 			return;
 		default:
 			_error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
 	}
 }
 /*}}}*/
-void  Numericalflux::GetNodesLidList(int* lidlist){/*{{{*/
+void  Numericalflux::GetNodesSidList(int* sidlist){/*{{{*/
 
-	_assert_(lidlist);
+	_assert_(sidlist);
 	_assert_(nodes);
 
 	switch(this->flux_type){
 		case InternalEnum:
-			for(int i=0;i<NUMNODES_INTERNAL;i++) lidlist[i]=nodes[i]->Lid();
+			for(int i=0;i<NUMNODES_INTERNAL;i++) sidlist[i]=nodes[i]->Sid();
 			return;
 		case BoundaryEnum:
-			for(int i=0;i<NUMNODES_BOUNDARY;i++) lidlist[i]=nodes[i]->Lid();
+			for(int i=0;i<NUMNODES_BOUNDARY;i++) sidlist[i]=nodes[i]->Sid();
 			return;
 		default:
 			_error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
@@ -337,6 +319,11 @@
 
 }
 /*}}}*/
+bool  Numericalflux::InAnalysis(int in_analysis_type){/*{{{*/
+	if (in_analysis_type==this->analysis_type) return true;
+	else return false;
+}
+/*}}}*/
 bool  Numericalflux::IsPenalty(void){/*{{{*/
 	return false;
 }
@@ -355,11 +342,24 @@
 
 }
 /*}}}*/
-bool  Numericalflux::InAnalysis(int in_analysis_type){/*{{{*/
-	if (in_analysis_type==this->analysis_type) return true;
-	else return false;
+void  Numericalflux::ResetHooks(){/*{{{*/
+
+	this->nodes=NULL;
+	this->vertices=NULL;
+	this->element=NULL;
+	this->parameters=NULL;
+
+	/*Get Element type*/
+	this->hnodes->reset();
+	this->hvertices->reset();
+	this->helement->reset();
+
 }
 /*}}}*/
+void  Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
+
+}
+/*}}}*/
 void  Numericalflux::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
 
 	/*Output */
@@ -416,90 +416,52 @@
 /*}}}*/
 
 /*Numericalflux management*/
-ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/
+ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){/*{{{*/
 
 	switch(this->flux_type){
 		case InternalEnum:
-			return CreateKMatrixMasstransportInternal();
+			return CreateKMatrixAdjointBalancethicknessInternal();
 		case BoundaryEnum:
-			return CreateKMatrixMasstransportBoundary();
+			return CreateKMatrixAdjointBalancethicknessBoundary();
 		default:
 			_error_("type not supported yet");
 	}
 }
 /*}}}*/
-ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/
+ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){/*{{{*/
 
-	/* constants*/
-	const int numdof=NDOF1*NUMNODES_INTERNAL;
+	ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary();
+	if(Ke) Ke->Transpose();
+	return Ke;
+}
+/*}}}*/
+ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){/*{{{*/
 
-	/* Intermediaries*/
-	int        i,j,ig,index1,index2;
-	IssmDouble     DL1,DL2,Jdet,dt,vx,vy,UdotN;
-	IssmDouble     xyz_list[NUMVERTICES][3];
-	IssmDouble     normal[2];
-	IssmDouble     B[numdof];
-	IssmDouble     Bprime[numdof];
-	IssmDouble     Ke_g1[numdof][numdof];
-	IssmDouble     Ke_g2[numdof][numdof];
-	GaussTria *gauss;
+	ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal();
+	if (Ke) Ke->Transpose();
+	return Ke;
+}
+/*}}}*/
+ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){/*{{{*/
 
-	/*Initialize Element matrix and return if necessary*/
-	Tria*  tria=(Tria*)element;
-	if(!tria->IsIceInElement()) return NULL;
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);
-
-	/*Retrieve all inputs and parameters*/
-	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-	Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
-	Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
-	GetNormal(&normal[0],xyz_list);
-
-	/* Start  looping on the number of gaussian points: */
-	index1=tria->GetNodeIndex(nodes[0]);
-	index2=tria->GetNodeIndex(nodes[1]);
-	gauss=new GaussTria(index1,index2,2);
-	for(ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
-		tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
-
-		vxaverage_input->GetInputValue(&vx,gauss);
-		vyaverage_input->GetInputValue(&vy,gauss);
-		UdotN=vx*normal[0]+vy*normal[1];
-		tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
-		DL1=gauss->weight*Jdet*dt*UdotN/2;
-		DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
-
-		TripleMultiply(&B[0],1,numdof,1,
-					&DL1,1,1,0,
-					&Bprime[0],1,numdof,0,
-					&Ke_g1[0][0],0);
-		TripleMultiply(&B[0],1,numdof,1,
-					&DL2,1,1,0,
-					&B[0],1,numdof,0,
-					&Ke_g2[0][0],0);
-
-		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
-		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
+	switch(this->flux_type){
+		case InternalEnum:
+			return CreateKMatrixBalancethicknessInternal();
+		case BoundaryEnum:
+			return CreateKMatrixBalancethicknessBoundary();
+		default:
+			_error_("type not supported yet");
 	}
-
-	/*Clean up and return*/
-	delete gauss;
-	return Ke;
 }
 /*}}}*/
-ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
+ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/
 
 	/* constants*/
 	const int numdof=NDOF1*NUMNODES_BOUNDARY;
 
 	/* Intermediaries*/
 	int        i,j,ig,index1,index2;
-	IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN;
+	IssmDouble     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN;
 	IssmDouble     xyz_list[NUMVERTICES][3];
 	IssmDouble     normal[2];
 	IssmDouble     L[numdof];
@@ -513,9 +475,8 @@
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-	Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
-	Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
+	Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
+	Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
 	GetNormal(&normal[0],xyz_list);
 
 	/*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
@@ -548,7 +509,7 @@
 		vyaverage_input->GetInputValue(&vy,gauss);
 		UdotN=vx*normal[0]+vy*normal[1];
 		tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
-		DL=gauss->weight*Jdet*dt*UdotN;
+		DL=gauss->weight*Jdet*UdotN;
 
 		TripleMultiply(&L[0],1,numdof,1,
 					&DL,1,1,0,
@@ -563,18 +524,6 @@
 	return Ke;
 }
 /*}}}*/
-ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){/*{{{*/
-
-	switch(this->flux_type){
-		case InternalEnum:
-			return CreateKMatrixBalancethicknessInternal();
-		case BoundaryEnum:
-			return CreateKMatrixBalancethicknessBoundary();
-		default:
-			_error_("type not supported yet");
-	}
-}
-/*}}}*/
 ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/
 
 	/* constants*/
@@ -638,14 +587,26 @@
 	return Ke;
 }
 /*}}}*/
-ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/
+ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/
 
+	switch(this->flux_type){
+		case InternalEnum:
+			return CreateKMatrixMasstransportInternal();
+		case BoundaryEnum:
+			return CreateKMatrixMasstransportBoundary();
+		default:
+			_error_("type not supported yet");
+	}
+}
+/*}}}*/
+ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
+
 	/* constants*/
 	const int numdof=NDOF1*NUMNODES_BOUNDARY;
 
 	/* Intermediaries*/
 	int        i,j,ig,index1,index2;
-	IssmDouble     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN;
+	IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN;
 	IssmDouble     xyz_list[NUMVERTICES][3];
 	IssmDouble     normal[2];
 	IssmDouble     L[numdof];
@@ -659,8 +620,9 @@
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
-	Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
+	parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
+	Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
 	GetNormal(&normal[0],xyz_list);
 
 	/*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
@@ -693,7 +655,7 @@
 		vyaverage_input->GetInputValue(&vy,gauss);
 		UdotN=vx*normal[0]+vy*normal[1];
 		tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
-		DL=gauss->weight*Jdet*UdotN;
+		DL=gauss->weight*Jdet*dt*UdotN;
 
 		TripleMultiply(&L[0],1,numdof,1,
 					&DL,1,1,0,
@@ -708,62 +670,99 @@
 	return Ke;
 }
 /*}}}*/
-ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){/*{{{*/
+ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/
 
-	switch(this->flux_type){
-		case InternalEnum:
-			return CreateKMatrixAdjointBalancethicknessInternal();
-		case BoundaryEnum:
-			return CreateKMatrixAdjointBalancethicknessBoundary();
-		default:
-			_error_("type not supported yet");
+	/* constants*/
+	const int numdof=NDOF1*NUMNODES_INTERNAL;
+
+	/* Intermediaries*/
+	int        i,j,ig,index1,index2;
+	IssmDouble     DL1,DL2,Jdet,dt,vx,vy,UdotN;
+	IssmDouble     xyz_list[NUMVERTICES][3];
+	IssmDouble     normal[2];
+	IssmDouble     B[numdof];
+	IssmDouble     Bprime[numdof];
+	IssmDouble     Ke_g1[numdof][numdof];
+	IssmDouble     Ke_g2[numdof][numdof];
+	GaussTria *gauss;
+
+	/*Initialize Element matrix and return if necessary*/
+	Tria*  tria=(Tria*)element;
+	if(!tria->IsIceInElement()) return NULL;
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
+	Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
+	GetNormal(&normal[0],xyz_list);
+
+	/* Start  looping on the number of gaussian points: */
+	index1=tria->GetNodeIndex(nodes[0]);
+	index2=tria->GetNodeIndex(nodes[1]);
+	gauss=new GaussTria(index1,index2,2);
+	for(ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement());
+		tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement());
+
+		vxaverage_input->GetInputValue(&vx,gauss);
+		vyaverage_input->GetInputValue(&vy,gauss);
+		UdotN=vx*normal[0]+vy*normal[1];
+		tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
+		DL1=gauss->weight*Jdet*dt*UdotN/2;
+		DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
+
+		TripleMultiply(&B[0],1,numdof,1,
+					&DL1,1,1,0,
+					&Bprime[0],1,numdof,0,
+					&Ke_g1[0][0],0);
+		TripleMultiply(&B[0],1,numdof,1,
+					&DL2,1,1,0,
+					&B[0],1,numdof,0,
+					&Ke_g2[0][0],0);
+
+		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
+		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
 	}
-}
-/*}}}*/
-ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){/*{{{*/
 
-	ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal();
-	if (Ke) Ke->Transpose();
+	/*Clean up and return*/
+	delete gauss;
 	return Ke;
 }
 /*}}}*/
-ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){/*{{{*/
+ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/
 
-	ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary();
-	if(Ke) Ke->Transpose();
-	return Ke;
+	/*No PVector for the Adjoint*/
+	return NULL;
 }
 /*}}}*/
-ElementVector* Numericalflux::CreatePVectorMasstransport(void){/*{{{*/
+ElementVector* Numericalflux::CreatePVectorBalancethickness(void){/*{{{*/
 
 	switch(this->flux_type){
 		case InternalEnum:
-			return CreatePVectorMasstransportInternal();
+			return CreatePVectorBalancethicknessInternal();
 		case BoundaryEnum:
-			return CreatePVectorMasstransportBoundary();
+			return CreatePVectorBalancethicknessBoundary();
 		default:
 			_error_("type not supported yet");
 	}
 }
 /*}}}*/
-ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){/*{{{*/
+ElementVector* Numericalflux::CreatePVectorBalancethicknessBoundary(void){/*{{{*/
 
-	/*Nothing added to PVector*/
-	return NULL;
-
-}
-/*}}}*/
-ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/
-
 	/* constants*/
 	const int numdof=NDOF1*NUMNODES_BOUNDARY;
 
 	/* Intermediaries*/
 	int        i,ig,index1,index2;
-	IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness;
-	IssmDouble     xyz_list[NUMVERTICES][3];
-	IssmDouble     normal[2];
-	IssmDouble     L[numdof];
+	IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble normal[2];
+	IssmDouble L[numdof];
 	GaussTria *gauss;
 
 	/*Initialize Load Vector and return if necessary*/
@@ -773,10 +772,9 @@
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-	Input* vxaverage_input   =tria->inputs->GetInput(VxEnum);                     _assert_(vxaverage_input); 
-	Input* vyaverage_input   =tria->inputs->GetInput(VyEnum);                     _assert_(vyaverage_input);
-	Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
+	Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); 
+	Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
+	Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
 	GetNormal(&normal[0],xyz_list);
 
 	/*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
@@ -788,7 +786,6 @@
 	vxaverage_input->GetInputValue(&mean_vx,gauss);
 	vyaverage_input->GetInputValue(&mean_vy,gauss);
 	delete gauss;
-
 	UdotN=mean_vx*normal[0]+mean_vy*normal[1];
 	if (UdotN>0){
 		return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
@@ -807,12 +804,11 @@
 
 		vxaverage_input->GetInputValue(&vx,gauss);
 		vyaverage_input->GetInputValue(&vy,gauss);
-		spcthickness_input->GetInputValue(&thickness,gauss);
-		if(xIsNan<IssmDouble>(thickness)) _error_("Cannot weakly apply constraint because NaN was provided");
+		thickness_input->GetInputValue(&thickness,gauss);
 
 		UdotN=vx*normal[0]+vy*normal[1];
 		tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
-		DL= - gauss->weight*Jdet*dt*UdotN*thickness;
+		DL= - gauss->weight*Jdet*UdotN*thickness;
 
 		for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
 	}
@@ -822,36 +818,36 @@
 	return pe;
 }
 /*}}}*/
-ElementVector* Numericalflux::CreatePVectorBalancethickness(void){/*{{{*/
+ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){/*{{{*/
 
+	/*Nothing added to PVector*/
+	return NULL;
+
+}
+/*}}}*/
+ElementVector* Numericalflux::CreatePVectorMasstransport(void){/*{{{*/
+
 	switch(this->flux_type){
 		case InternalEnum:
-			return CreatePVectorBalancethicknessInternal();
+			return CreatePVectorMasstransportInternal();
 		case BoundaryEnum:
-			return CreatePVectorBalancethicknessBoundary();
+			return CreatePVectorMasstransportBoundary();
 		default:
 			_error_("type not supported yet");
 	}
 }
 /*}}}*/
-ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){/*{{{*/
+ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/
 
-	/*Nothing added to PVector*/
-	return NULL;
-
-}
-/*}}}*/
-ElementVector* Numericalflux::CreatePVectorBalancethicknessBoundary(void){/*{{{*/
-
 	/* constants*/
 	const int numdof=NDOF1*NUMNODES_BOUNDARY;
 
 	/* Intermediaries*/
 	int        i,ig,index1,index2;
-	IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble normal[2];
-	IssmDouble L[numdof];
+	IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness;
+	IssmDouble     xyz_list[NUMVERTICES][3];
+	IssmDouble     normal[2];
+	IssmDouble     L[numdof];
 	GaussTria *gauss;
 
 	/*Initialize Load Vector and return if necessary*/
@@ -861,9 +857,10 @@
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); 
-	Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
-	Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* vxaverage_input   =tria->inputs->GetInput(VxEnum);                     _assert_(vxaverage_input); 
+	Input* vyaverage_input   =tria->inputs->GetInput(VyEnum);                     _assert_(vyaverage_input);
+	Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
 	GetNormal(&normal[0],xyz_list);
 
 	/*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
@@ -875,6 +872,7 @@
 	vxaverage_input->GetInputValue(&mean_vx,gauss);
 	vyaverage_input->GetInputValue(&mean_vy,gauss);
 	delete gauss;
+
 	UdotN=mean_vx*normal[0]+mean_vy*normal[1];
 	if (UdotN>0){
 		return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
@@ -893,11 +891,12 @@
 
 		vxaverage_input->GetInputValue(&vx,gauss);
 		vyaverage_input->GetInputValue(&vy,gauss);
-		thickness_input->GetInputValue(&thickness,gauss);
+		spcthickness_input->GetInputValue(&thickness,gauss);
+		if(xIsNan<IssmDouble>(thickness)) _error_("Cannot weakly apply constraint because NaN was provided");
 
 		UdotN=vx*normal[0]+vy*normal[1];
 		tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
-		DL= - gauss->weight*Jdet*UdotN*thickness;
+		DL= - gauss->weight*Jdet*dt*UdotN*thickness;
 
 		for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
 	}
@@ -907,13 +906,14 @@
 	return pe;
 }
 /*}}}*/
-ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/
+ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){/*{{{*/
 
-	/*No PVector for the Adjoint*/
+	/*Nothing added to PVector*/
 	return NULL;
+
 }
 /*}}}*/
-void Numericalflux:: GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){/*{{{*/
+void           Numericalflux::GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){/*{{{*/
 
 	/*Build unit outward pointing vector*/
 	IssmDouble vector[2];
Index: ../trunk-jpl/src/c/classes/Loads/Penpair.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Penpair.cpp	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Penpair.cpp	(revision 18926)
@@ -44,33 +44,6 @@
 /*}}}*/
 
 /*Object virtual functions definitions:*/
-void Penpair::Echo(void){/*{{{*/
-
-	_printf_("Penpair:\n");
-	_printf_("   id: " << id << "\n");
-	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
-	hnodes->Echo();
-
-	return;
-}
-/*}}}*/
-void Penpair::DeepEcho(void){/*{{{*/
-
-	_printf_("Penpair:\n");
-	_printf_("   id: " << id << "\n");
-	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
-	hnodes->DeepEcho();
-
-	return;
-}		
-/*}}}*/
-int  Penpair::Id(void){ return id; }/*{{{*/
-/*}}}*/
-int  Penpair::ObjectEnum(void){/*{{{*/
-
-	return PenpairEnum;
-}
-/*}}}*/
 Object* Penpair::copy() {/*{{{*/
 
 	Penpair* penpair=NULL;
@@ -92,7 +65,34 @@
 
 }
 /*}}}*/
+void    Penpair::DeepEcho(void){/*{{{*/
 
+	_printf_("Penpair:\n");
+	_printf_("   id: " << id << "\n");
+	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
+	hnodes->DeepEcho();
+
+	return;
+}		
+/*}}}*/
+void    Penpair::Echo(void){/*{{{*/
+
+	_printf_("Penpair:\n");
+	_printf_("   id: " << id << "\n");
+	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
+	hnodes->Echo();
+
+	return;
+}
+/*}}}*/
+int     Penpair::Id(void){ return id; }/*{{{*/
+/*}}}*/
+int     Penpair::ObjectEnum(void){/*{{{*/
+
+	return PenpairEnum;
+}
+/*}}}*/
+
 /*Load virtual functions definitions:*/
 void  Penpair::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
 
@@ -108,20 +108,10 @@
 
 }
 /*}}}*/
-void  Penpair::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
-
+void  Penpair::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){/*{{{*/
+	this->CreateKMatrix(Jff,NULL);
 }
 /*}}}*/
-void  Penpair::ResetHooks(){/*{{{*/
-
-	this->nodes=NULL;
-	this->parameters=NULL;
-
-	/*Get Element type*/
-	this->hnodes->reset();
-
-}
-/*}}}*/
 void  Penpair::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
 	/*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: */
@@ -136,10 +126,6 @@
 
 }
 /*}}}*/
-void  Penpair::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){/*{{{*/
-	this->CreateKMatrix(Jff,NULL);
-}
-/*}}}*/
 void  Penpair::GetNodesSidList(int* sidlist){/*{{{*/
 
 	_assert_(sidlist);
@@ -161,10 +147,19 @@
 	return NUMVERTICES;
 }
 /*}}}*/
+bool  Penpair::InAnalysis(int in_analysis_type){/*{{{*/
+	if (in_analysis_type==this->analysis_type)return true;
+	else return false;
+}
+/*}}}*/
 bool  Penpair::IsPenalty(void){/*{{{*/
 	return true;
 }
 /*}}}*/
+void  Penpair::PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){/*{{{*/
+	this->PenaltyCreateKMatrix(Jff,NULL,kmax);
+}
+/*}}}*/
 void  Penpair::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
 
 	/*Retrieve parameters: */
@@ -195,13 +190,18 @@
 	return;
 }
 /*}}}*/
-void  Penpair::PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){/*{{{*/
-	this->PenaltyCreateKMatrix(Jff,NULL,kmax);
+void  Penpair::ResetHooks(){/*{{{*/
+
+	this->nodes=NULL;
+	this->parameters=NULL;
+
+	/*Get Element type*/
+	this->hnodes->reset();
+
 }
 /*}}}*/
-bool  Penpair::InAnalysis(int in_analysis_type){/*{{{*/
-	if (in_analysis_type==this->analysis_type)return true;
-	else return false;
+void  Penpair::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
+
 }
 /*}}}*/
 void  Penpair::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
@@ -278,6 +278,58 @@
 /*}}}*/
 
 /*Penpair management:*/
+ElementMatrix* Penpair::PenaltyCreateKMatrixMasstransport(IssmDouble kmax){/*{{{*/
+
+	const int numdof=NUMVERTICES*NDOF1;
+	IssmDouble penalty_factor;
+
+	/*Initialize Element vector and return if necessary*/
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
+
+	/*recover parameters: */
+	parameters->FindParam(&penalty_factor,MasstransportPenaltyFactorEnum);
+
+	//Create elementary matrix: add penalty to 
+	Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_factor);
+	Ke->values[0*numdof+1]=-kmax*pow(10.,penalty_factor);
+	Ke->values[1*numdof+0]=-kmax*pow(10.,penalty_factor);
+	Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_factor);
+
+	/*Clean up and return*/
+	return Ke;
+}
+/*}}}*/
+ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax){/*{{{*/
+
+	const int  numdof=NUMVERTICES*NDOF3;
+	IssmDouble penalty_offset;
+
+	/*Initialize Element vector and return if necessary*/
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
+
+	/*recover parameters: */
+	parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
+
+	//Create elementary matrix: add penalty to 
+	Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset);
+	Ke->values[0*numdof+3]=-kmax*pow(10.,penalty_offset);
+	Ke->values[3*numdof+0]=-kmax*pow(10.,penalty_offset);
+	Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset);
+
+	Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset);
+	Ke->values[1*numdof+4]=-kmax*pow(10.,penalty_offset);
+	Ke->values[4*numdof+1]=-kmax*pow(10.,penalty_offset);
+	Ke->values[4*numdof+4]=+kmax*pow(10.,penalty_offset);
+
+	Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset);
+	Ke->values[2*numdof+5]=-kmax*pow(10.,penalty_offset);
+	Ke->values[5*numdof+2]=-kmax*pow(10.,penalty_offset);
+	Ke->values[5*numdof+5]=+kmax*pow(10.,penalty_offset);
+
+	/*Clean up and return*/
+	return Ke;
+}
+/*}}}*/
 ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax){/*{{{*/
 
 	int    approximation0=nodes[0]->GetApproximation();
@@ -337,55 +389,3 @@
 	return Ke;
 }
 /*}}}*/
-ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax){/*{{{*/
-
-	const int  numdof=NUMVERTICES*NDOF3;
-	IssmDouble penalty_offset;
-
-	/*Initialize Element vector and return if necessary*/
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
-
-	/*recover parameters: */
-	parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
-
-	//Create elementary matrix: add penalty to 
-	Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset);
-	Ke->values[0*numdof+3]=-kmax*pow(10.,penalty_offset);
-	Ke->values[3*numdof+0]=-kmax*pow(10.,penalty_offset);
-	Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset);
-
-	Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset);
-	Ke->values[1*numdof+4]=-kmax*pow(10.,penalty_offset);
-	Ke->values[4*numdof+1]=-kmax*pow(10.,penalty_offset);
-	Ke->values[4*numdof+4]=+kmax*pow(10.,penalty_offset);
-
-	Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset);
-	Ke->values[2*numdof+5]=-kmax*pow(10.,penalty_offset);
-	Ke->values[5*numdof+2]=-kmax*pow(10.,penalty_offset);
-	Ke->values[5*numdof+5]=+kmax*pow(10.,penalty_offset);
-
-	/*Clean up and return*/
-	return Ke;
-}
-/*}}}*/
-ElementMatrix* Penpair::PenaltyCreateKMatrixMasstransport(IssmDouble kmax){/*{{{*/
-
-	const int numdof=NUMVERTICES*NDOF1;
-	IssmDouble penalty_factor;
-
-	/*Initialize Element vector and return if necessary*/
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
-
-	/*recover parameters: */
-	parameters->FindParam(&penalty_factor,MasstransportPenaltyFactorEnum);
-
-	//Create elementary matrix: add penalty to 
-	Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_factor);
-	Ke->values[0*numdof+1]=-kmax*pow(10.,penalty_factor);
-	Ke->values[1*numdof+0]=-kmax*pow(10.,penalty_factor);
-	Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_factor);
-
-	/*Clean up and return*/
-	return Ke;
-}
-/*}}}*/
Index: ../trunk-jpl/src/c/classes/Loads/Friction.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Friction.h	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Friction.h	(revision 18926)
@@ -28,14 +28,14 @@
 		~Friction();
 
 		void  Echo(void);
+		void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
 		void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
-		void  GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss);
-		void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
 		void  GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
 		void  GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss);
+		void  GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss);
 		void  GetAlpha2WaterLayer(IssmDouble* palpha2,Gauss* gauss);
+		void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
 		void  GetAlpha2WeertmanTemp(IssmDouble* palpha2,Gauss* gauss);
-		void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
 };
 
 #endif  /* _FRICTION_H_ */
Index: ../trunk-jpl/src/c/classes/Loads/Loads.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Loads.cpp	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Loads.cpp	(revision 18926)
@@ -47,20 +47,6 @@
 
 }
 /*}}}*/
-void Loads::ResetHooks(){/*{{{*/
-
-	vector<Object*>::iterator object;
-	Load* load=NULL;
-
-	for ( object=objects.begin() ; object < objects.end(); object++ ){
-
-		load=xDynamicCast<Load*>((*object));
-		load->ResetHooks();
-
-	}
-
-}
-/*}}}*/
 bool Loads::IsPenalty(int analysis_type){/*{{{*/
 
 	int ispenalty=0;
@@ -85,7 +71,7 @@
 	 return false;
 }
 /*}}}*/
-int Loads::MaxNumNodes(int analysis_type){/*{{{*/
+int  Loads::MaxNumNodes(int analysis_type){/*{{{*/
 
 	int max=0;
 	int allmax;
@@ -108,7 +94,7 @@
 	return max;
 }
 /*}}}*/
-int Loads::NumberOfLoads(void){/*{{{*/
+int  Loads::NumberOfLoads(void){/*{{{*/
 
 	int localloads;
 	int numberofloads;
@@ -123,7 +109,7 @@
 	return numberofloads;
 }
 /*}}}*/
-int Loads::NumberOfLoads(int analysis_type){/*{{{*/
+int  Loads::NumberOfLoads(int analysis_type){/*{{{*/
 
 	int localloads = 0;
 	int numberofloads;
@@ -144,25 +130,18 @@
 	return numberofloads;
 }
 /*}}}*/
-int Loads::Size(void){/*{{{*/
+void Loads::ResetHooks(){/*{{{*/
 
-	return this->DataSet::Size();
-}
-/*}}}*/
-int Loads::Size(int analysis_type){/*{{{*/
+	vector<Object*>::iterator object;
+	Load* load=NULL;
 
-	int localloads = 0;
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
 
-	/*Get number of local loads*/
-	for(int i=0;i<this->Size();i++){
+		load=xDynamicCast<Load*>((*object));
+		load->ResetHooks();
 
-		Load* load=xDynamicCast<Load*>(this->GetObjectByOffset(i));
-
-		/*Check that this load corresponds to our analysis currently being carried out: */
-		if (load->InAnalysis(analysis_type)) localloads++;
 	}
 
-	return localloads;
 }
 /*}}}*/
 void Loads::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){/*{{{*/
@@ -179,3 +158,24 @@
 
 }
 /*}}}*/
+int  Loads::Size(void){/*{{{*/
+
+	return this->DataSet::Size();
+}
+/*}}}*/
+int  Loads::Size(int analysis_type){/*{{{*/
+
+	int localloads = 0;
+
+	/*Get number of local loads*/
+	for(int i=0;i<this->Size();i++){
+
+		Load* load=xDynamicCast<Load*>(this->GetObjectByOffset(i));
+
+		/*Check that this load corresponds to our analysis currently being carried out: */
+		if (load->InAnalysis(analysis_type)) localloads++;
+	}
+
+	return localloads;
+}
+/*}}}*/
Index: ../trunk-jpl/src/c/classes/Loads/Riftfront.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Riftfront.h	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Riftfront.h	(revision 18926)
@@ -56,44 +56,44 @@
 		~Riftfront();
 		/*}}}*/
 		/*Object virtual functions definitions:{{{ */
-		void  Echo();
-		void  DeepEcho();
-		int   Id(); 
-		int   ObjectEnum();
-		Object* copy();
+		Object*  copy();
+		void     DeepEcho();
+		void     Echo();
+		int      Id(); 
+		int      ObjectEnum();
 		/*}}}*/
 		/*Update virtual functions resolution: {{{*/
-		void    InputUpdateFromVector(IssmDouble* vector, int name, int type);
-		void    InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows,int ncols, int name, int type){_error_("Not implemented yet!");}
-		void    InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");}
 		void    InputUpdateFromConstant(IssmDouble constant, int name);
 		void    InputUpdateFromConstant(int constant, int name){_error_("Not implemented yet!");}
 		void    InputUpdateFromConstant(bool constant, int name);
 		void    InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
+		void    InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows,int ncols, int name, int type){_error_("Not implemented yet!");}
+		void    InputUpdateFromVector(IssmDouble* vector, int name, int type);
+		void    InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");}
 		/*}}}*/
 		/*Load virtual functions definitions: {{{*/
 		void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
-		void  ResetHooks();
-		void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
+		void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
 		void  CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
 		void  CreatePVector(Vector<IssmDouble>* pf);
-		void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
+		void  GetNodesLidList(int* lidlist);
 		void  GetNodesSidList(int* sidlist);
-		void  GetNodesLidList(int* lidlist);
 		int   GetNumberOfNodes(void);
+		bool  InAnalysis(int analysis_type);
 		bool  IsPenalty(void);
 		void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};
 		void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);
 		void  PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);
+		void  ResetHooks();
+		void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
 		void  SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
-		bool  InAnalysis(int analysis_type);
 		/*}}}*/
 		/*Riftfront specific routines: {{{*/
+		int            Constrain(int* punstable);
+		void           FreezeConstraints(void);
+		bool           IsFrozen(void);
 		ElementMatrix* PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax);
 		ElementVector* PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax);
-		int   Constrain(int* punstable);
-		void  FreezeConstraints(void);
-		bool  IsFrozen(void);
 		/*}}}*/
 };
 #endif  /* _RIFTFRONT_H_ */
Index: ../trunk-jpl/src/c/classes/Loads/Numericalflux.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Numericalflux.h	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Numericalflux.h	(revision 18926)
@@ -39,56 +39,56 @@
 		~Numericalflux();
 		/*}}}*/
 		/*Object virtual functions definitions:{{{ */
+		Object *copy();
+		void    DeepEcho();
 		void    Echo();
-		void    DeepEcho();
 		int     Id();
 		int     ObjectEnum();
-		Object *copy();
 		/*}}}*/
 		/*Update virtual functions resolution: {{{*/
-		void InputUpdateFromVector(IssmDouble* vector, int name, int type){/*Do nothing*/}
-		void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*Do nothing*/}
-		void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*Do nothing*/}
 		void InputUpdateFromConstant(IssmDouble constant, int name){/*Do nothing*/};
 		void InputUpdateFromConstant(int constant, int name){/*Do nothing*/};
 		void InputUpdateFromConstant(bool constant, int name){_error_("Not implemented yet!");}
 		void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
+		void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*Do nothing*/}
+		void InputUpdateFromVector(IssmDouble* vector, int name, int type){/*Do nothing*/}
+		void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*Do nothing*/}
 		/*}}}*/
 		/*Load virtual functions definitions: {{{*/
 		void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
-		void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
-		void ResetHooks();
+		void CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
 		void CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
 		void CreatePVector(Vector<IssmDouble>* pf);
+		void GetNodesLidList(int* lidlist);
 		void GetNodesSidList(int* sidlist);
-		void GetNodesLidList(int* lidlist);
 		int  GetNumberOfNodes(void);
-		void CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
+		bool InAnalysis(int analysis_type);
 		bool IsPenalty(void);
 		void PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};
 		void PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);
 		void PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);
+		void ResetHooks();
 		void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
-		bool InAnalysis(int analysis_type);
+		void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
 		/*}}}*/
 		/*Numericalflux management:{{{*/
-		void           GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
+		ElementMatrix* CreateKMatrixAdjointBalancethickness(void);
+		ElementMatrix* CreateKMatrixAdjointBalancethicknessBoundary(void);
+		ElementMatrix* CreateKMatrixAdjointBalancethicknessInternal(void);
+		ElementMatrix* CreateKMatrixBalancethickness(void);
+		ElementMatrix* CreateKMatrixBalancethicknessBoundary(void);
+		ElementMatrix* CreateKMatrixBalancethicknessInternal(void);
 		ElementMatrix* CreateKMatrixMasstransport(void);
+		ElementMatrix* CreateKMatrixMasstransportBoundary(void);
 		ElementMatrix* CreateKMatrixMasstransportInternal(void);
-		ElementMatrix* CreateKMatrixMasstransportBoundary(void);
-		ElementMatrix* CreateKMatrixBalancethickness(void);
-		ElementMatrix* CreateKMatrixBalancethicknessInternal(void);
-		ElementMatrix* CreateKMatrixBalancethicknessBoundary(void);
-		ElementMatrix* CreateKMatrixAdjointBalancethickness(void);
-		ElementMatrix* CreateKMatrixAdjointBalancethicknessInternal(void);
-		ElementMatrix* CreateKMatrixAdjointBalancethicknessBoundary(void);
+		ElementVector* CreatePVectorAdjointBalancethickness(void);
+		ElementVector* CreatePVectorBalancethickness(void);
+		ElementVector* CreatePVectorBalancethicknessBoundary(void);
+		ElementVector* CreatePVectorBalancethicknessInternal(void);
 		ElementVector* CreatePVectorMasstransport(void);
+		ElementVector* CreatePVectorMasstransportBoundary(void);
 		ElementVector* CreatePVectorMasstransportInternal(void);
-		ElementVector* CreatePVectorMasstransportBoundary(void);
-		ElementVector* CreatePVectorBalancethickness(void);
-		ElementVector* CreatePVectorBalancethicknessInternal(void);
-		ElementVector* CreatePVectorBalancethicknessBoundary(void);
-		ElementVector* CreatePVectorAdjointBalancethickness(void);
+		void           GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
 		/*}}}*/
 
 };
Index: ../trunk-jpl/src/c/classes/Loads/Penpair.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Penpair.h	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Penpair.h	(revision 18926)
@@ -30,43 +30,43 @@
 		~Penpair();
 		/*}}}*/
 		/*Object virtual functions definitions:{{{ */
-		void  Echo();
-		void  DeepEcho();
-		int   Id(); 
-		int   ObjectEnum();
 		Object* copy();
+		void     DeepEcho();
+		void     Echo();
+		int      Id(); 
+		int      ObjectEnum();
 		/*}}}*/
 		/*Update virtual functions resolution: {{{*/
-		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
-		void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols,int name, int type){_error_("Not implemented yet!");}
-		void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");}
 		void  InputUpdateFromConstant(IssmDouble constant, int name);
 		void  InputUpdateFromConstant(int constant, int name);
 		void  InputUpdateFromConstant(bool constant, int name);
 		void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
+		void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols,int name, int type){_error_("Not implemented yet!");}
+		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
+		void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");}
 		/*}}}*/
 			/*Load virtual functions definitions: {{{*/
 		void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
-		void  ResetHooks();
-		void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
+		void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
 		void  CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
 		void  CreatePVector(Vector<IssmDouble>* pf);
-		void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
+		void  GetNodesLidList(int* lidlist);
 		void  GetNodesSidList(int* sidlist);
-		void  GetNodesLidList(int* lidlist);
 		int   GetNumberOfNodes(void);
+		bool  InAnalysis(int analysis_type);
 		bool  IsPenalty(void);
+		void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax);
 		void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff,Matrix<IssmDouble>* Kfs,IssmDouble kmax);
 		void  PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);
-		void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax);
+		void  ResetHooks();
+		void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
 		void  SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
-		bool  InAnalysis(int analysis_type);
 		/*}}}*/
 			/*Penpair management: {{{*/
+		ElementMatrix* PenaltyCreateKMatrixMasstransport(IssmDouble kmax);
+		ElementMatrix* PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax);
 		ElementMatrix* PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax);
 		ElementMatrix* PenaltyCreateKMatrixStressbalanceSSAHO(IssmDouble kmax);
-		ElementMatrix* PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax);
-		ElementMatrix* PenaltyCreateKMatrixMasstransport(IssmDouble kmax);
 		/*}}}*/
 };
 
Index: ../trunk-jpl/src/c/classes/Loads/Load.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Load.h	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Load.h	(revision 18926)
@@ -25,19 +25,19 @@
 	public: 
 		virtual       ~Load(){};
 		virtual void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
-		virtual void  ResetHooks()=0;
-		virtual bool  IsPenalty(void)=0;
-		virtual int   GetNumberOfNodes(void)=0;
-		virtual void  GetNodesSidList(int* sidlist)=0;
-		virtual void  GetNodesLidList(int* lidlist)=0;
-		virtual void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
+		virtual void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
 		virtual void  CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs)=0;
 		virtual void  CreatePVector(Vector<IssmDouble>* pf)=0;
-		virtual void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
+		virtual void  GetNodesLidList(int* lidlist)=0;
+		virtual void  GetNodesSidList(int* sidlist)=0;
+		virtual int   GetNumberOfNodes(void)=0;
+		virtual bool  InAnalysis(int analysis_type)=0;
+		virtual bool  IsPenalty(void)=0;
 		virtual void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax)=0;
 		virtual void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs, IssmDouble kmax)=0;
 		virtual void  PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax)=0;
-		virtual bool  InAnalysis(int analysis_type)=0;
+		virtual void  ResetHooks()=0;
+		virtual void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
 		virtual void  SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
 };
 #endif
Index: ../trunk-jpl/src/c/classes/Loads/Loads.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Loads.h	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Loads.h	(revision 18926)
@@ -23,11 +23,11 @@
 
 		/*numerics*/
 		void  Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
-		void  ResetHooks();
 		bool  IsPenalty(int analysis);
 		int   MaxNumNodes(int analysis);
 		int   NumberOfLoads(void);
 		int   NumberOfLoads(int analysis);
+		void  ResetHooks();
 		void  SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
 		int   Size(int analysis);
 		int   Size(void);
Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 18926)
@@ -80,31 +80,6 @@
 /*}}}*/
 
 /*Object virtual functions definitions:*/
-void Pengrid::Echo(void){/*{{{*/
-	this->DeepEcho();
-}
-/*}}}*/
-void Pengrid::DeepEcho(void){/*{{{*/
-
-	_printf_("Pengrid:\n");
-	_printf_("   id: " << id << "\n");
-	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
-	hnode->DeepEcho();
-	helement->DeepEcho();
-	hmatpar->DeepEcho();
-	_printf_("   active " << this->active << "\n");
-	_printf_("   zigzag_counter " << this->zigzag_counter << "\n");
-	_printf_("   parameters\n");
-	parameters->DeepEcho();
-}
-/*}}}*/
-int    Pengrid::Id(void){ return id; }/*{{{*/
-/*}}}*/
-int Pengrid::ObjectEnum(void){/*{{{*/
-
-	return PengridEnum;
-}
-/*}}}*/
 Object* Pengrid::copy() {/*{{{*/
 
 	Pengrid* pengrid=NULL;
@@ -136,7 +111,32 @@
 
 }
 /*}}}*/
+void    Pengrid::DeepEcho(void){/*{{{*/
 
+	_printf_("Pengrid:\n");
+	_printf_("   id: " << id << "\n");
+	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
+	hnode->DeepEcho();
+	helement->DeepEcho();
+	hmatpar->DeepEcho();
+	_printf_("   active " << this->active << "\n");
+	_printf_("   zigzag_counter " << this->zigzag_counter << "\n");
+	_printf_("   parameters\n");
+	parameters->DeepEcho();
+}
+/*}}}*/
+void    Pengrid::Echo(void){/*{{{*/
+	this->DeepEcho();
+}
+/*}}}*/
+int     Pengrid::Id(void){ return id; }/*{{{*/
+/*}}}*/
+int     Pengrid::ObjectEnum(void){/*{{{*/
+
+	return PengridEnum;
+}
+/*}}}*/
+
 /*Load virtual functions definitions:*/
 void  Pengrid::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
 
@@ -155,24 +155,6 @@
 	this->parameters=parametersin;
 }
 /*}}}*/
-void  Pengrid::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
-
-}
-/*}}}*/
-void  Pengrid::ResetHooks(){/*{{{*/
-
-	this->node=NULL;
-	this->element=NULL;
-	this->matpar=NULL;
-	this->parameters=NULL;
-
-	/*Get Element type*/
-	this->hnode->reset();
-	this->helement->reset();
-	this->hmatpar->reset();
-
-}
-/*}}}*/
 void  Pengrid::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
 
 	/*No loads applied, do nothing: */
@@ -202,20 +184,20 @@
 
 }
 /*}}}*/
-void  Pengrid::GetNodesSidList(int* sidlist){/*{{{*/
+void  Pengrid::GetNodesLidList(int* lidlist){/*{{{*/
 
-	_assert_(sidlist);
+	_assert_(lidlist);
 	_assert_(node);
 
-	sidlist[0]=node->Sid();
+	lidlist[0]=node->Lid();
 }
 /*}}}*/
-void  Pengrid::GetNodesLidList(int* lidlist){/*{{{*/
+void  Pengrid::GetNodesSidList(int* sidlist){/*{{{*/
 
-	_assert_(lidlist);
+	_assert_(sidlist);
 	_assert_(node);
 
-	lidlist[0]=node->Lid();
+	sidlist[0]=node->Sid();
 }
 /*}}}*/
 int   Pengrid::GetNumberOfNodes(void){/*{{{*/
@@ -223,6 +205,15 @@
 	return NUMVERTICES;
 }
 /*}}}*/
+bool  Pengrid::InAnalysis(int in_analysis_type){/*{{{*/
+	if (in_analysis_type==this->analysis_type)return true;
+	else return false;
+}
+/*}}}*/
+bool  Pengrid::IsPenalty(void){/*{{{*/
+	return true;
+}
+/*}}}*/
 void  Pengrid::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
 
 	/*Retrieve parameters: */
@@ -281,13 +272,22 @@
 	}
 }
 /*}}}*/
-bool  Pengrid::InAnalysis(int in_analysis_type){/*{{{*/
-	if (in_analysis_type==this->analysis_type)return true;
-	else return false;
+void  Pengrid::ResetHooks(){/*{{{*/
+
+	this->node=NULL;
+	this->element=NULL;
+	this->matpar=NULL;
+	this->parameters=NULL;
+
+	/*Get Element type*/
+	this->hnode->reset();
+	this->helement->reset();
+	this->hmatpar->reset();
+
 }
 /*}}}*/
-bool  Pengrid::IsPenalty(void){/*{{{*/
-	return true;
+void  Pengrid::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
+
 }
 /*}}}*/
 void  Pengrid::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
@@ -342,18 +342,6 @@
 /*}}}*/
 
 /*Update virtual functions definitions:*/
-void  Pengrid::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
-	/*Nothing updated yet*/
-}
-/*}}}*/
-void  Pengrid::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
-	/*Nothing updated yet*/
-}
-/*}}}*/
-void  Pengrid::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
-	/*Nothing updated yet*/
-}
-/*}}}*/
 void  Pengrid::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
 	/*Nothing*/
 }
@@ -373,9 +361,21 @@
 	}
 }
 /*}}}*/
+void  Pengrid::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+void  Pengrid::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+void  Pengrid::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
 
 /*Pengrid management:*/
-void  Pengrid::ConstraintActivate(int* punstable){/*{{{*/
+void           Pengrid::ConstraintActivate(int* punstable){/*{{{*/
 
 	int analysis_type;
 
@@ -403,8 +403,68 @@
 
 }
 /*}}}*/
-void  Pengrid::ConstraintActivateThermal(int* punstable){/*{{{*/
+void           Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){/*{{{*/
 
+	//   The penalty is stable if it doesn't change during two consecutive iterations.   
+	int        unstable=0;
+	int        new_active;
+	int        penalty_lock;
+	IssmDouble pressure;
+	IssmDouble h;
+	IssmDouble h_max;	
+	HydrologyDCInefficientAnalysis* inefanalysis = NULL;
+
+	/*check that pengrid is not a clone (penalty to be added only once)*/
+	if(node->IsClone()){
+		unstable=0;
+		*punstable=unstable;
+		return;
+	}
+	if(!element->IsOnBase()){
+		unstable=0;
+		active=0;
+		*punstable=unstable;
+		return;
+	}
+
+	/*Get sediment water head h*/
+	inefanalysis = new HydrologyDCInefficientAnalysis();
+	element->GetInputValue(&h,node,SedimentHeadEnum);
+	inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node);
+	parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
+
+	if (h>h_max){
+	 new_active=1;
+	}
+	else{
+	 new_active=0;
+	}
+
+	if(this->active==new_active){
+		unstable=0;
+	}
+	else{
+		unstable=1;
+		if(penalty_lock)zigzag_counter++;
+	}
+
+	/*If penalty keeps zigzagging more than penalty_lock times: */
+	if(penalty_lock){
+		if(zigzag_counter>penalty_lock){
+			unstable=0;
+			active=1;
+		}
+	}
+	/*Set penalty flag*/
+	this->active=new_active;
+
+	/*Assign output pointers:*/
+	delete inefanalysis;
+	*punstable=unstable;
+}
+/*}}}*/
+void           Pengrid::ConstraintActivateThermal(int* punstable){/*{{{*/
+
 	//   The penalty is stable if it doesn't change during to successive iterations.   
 	IssmDouble pressure;
 	IssmDouble temperature;
@@ -466,6 +526,38 @@
 	*punstable=unstable;
 }
 /*}}}*/
+ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){/*{{{*/
+
+	IssmDouble moulin_load,dt;
+
+	/*Initialize Element matrix*/
+	ElementVector* pe=new ElementVector(&node,1,this->parameters);
+
+	this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
+	parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+	if(dt!=0.0) pe->values[0]=moulin_load*dt;
+
+	/*Clean up and return*/
+	return pe;
+ }
+/*}}}*/
+ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){/*{{{*/
+	IssmDouble    penalty_factor;
+
+	/*Retrieve parameters*/
+	parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
+
+	/*Initialize Element matrix and return if necessary*/
+	if(!this->active) return NULL;
+	ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
+
+	Ke->values[0]=kmax*pow(10.,penalty_factor);
+
+	/*Clean up and return*/
+	return Ke;
+}
+/*}}}*/
 ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(IssmDouble kmax){/*{{{*/
 
 	IssmDouble pressure,temperature,t_pmp;
@@ -511,6 +603,30 @@
 	return Ke;
 }
 /*}}}*/
+ElementVector* Pengrid::PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax){/*{{{*/
+
+	IssmDouble h_max;
+	IssmDouble penalty_factor;
+	HydrologyDCInefficientAnalysis* inefanalysis = NULL;
+
+	/*Initialize Element matrix and return if necessary*/
+	if(!this->active) return NULL;
+	ElementVector* pe=new ElementVector(&node,1,this->parameters);
+	inefanalysis = new HydrologyDCInefficientAnalysis();
+
+	/*Retrieve parameters*/
+	parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
+
+	/*Get h_max and compute penalty*/
+	inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node);
+
+	pe->values[0]=kmax*pow(10.,penalty_factor)*h_max;
+
+	/*Clean up and return*/
+	delete inefanalysis;
+	return pe;
+}
+/*}}}*/
 ElementVector* Pengrid::PenaltyCreatePVectorMelting(IssmDouble kmax){/*{{{*/
 
 	IssmDouble pressure;
@@ -577,123 +693,7 @@
 	return pe;
 }
 /*}}}*/
-void  Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){/*{{{*/
-
-	//   The penalty is stable if it doesn't change during two consecutive iterations.   
-	int        unstable=0;
-	int        new_active;
-	int        penalty_lock;
-	IssmDouble pressure;
-	IssmDouble h;
-	IssmDouble h_max;	
-	HydrologyDCInefficientAnalysis* inefanalysis = NULL;
-
-	/*check that pengrid is not a clone (penalty to be added only once)*/
-	if(node->IsClone()){
-		unstable=0;
-		*punstable=unstable;
-		return;
-	}
-	if(!element->IsOnBase()){
-		unstable=0;
-		active=0;
-		*punstable=unstable;
-		return;
-	}
-
-	/*Get sediment water head h*/
-	inefanalysis = new HydrologyDCInefficientAnalysis();
-	element->GetInputValue(&h,node,SedimentHeadEnum);
-	inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node);
-	parameters->FindParam(&penalty_lock,HydrologydcPenaltyLockEnum);
-
-	if (h>h_max){
-	 new_active=1;
-	}
-	else{
-	 new_active=0;
-	}
-
-	if(this->active==new_active){
-		unstable=0;
-	}
-	else{
-		unstable=1;
-		if(penalty_lock)zigzag_counter++;
-	}
-
-	/*If penalty keeps zigzagging more than penalty_lock times: */
-	if(penalty_lock){
-		if(zigzag_counter>penalty_lock){
-			unstable=0;
-			active=1;
-		}
-	}
-	/*Set penalty flag*/
-	this->active=new_active;
-
-	/*Assign output pointers:*/
-	delete inefanalysis;
-	*punstable=unstable;
-}
-/*}}}*/
-ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){/*{{{*/
-	IssmDouble    penalty_factor;
-
-	/*Retrieve parameters*/
-	parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
-
-	/*Initialize Element matrix and return if necessary*/
-	if(!this->active) return NULL;
-	ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
-
-	Ke->values[0]=kmax*pow(10.,penalty_factor);
-
-	/*Clean up and return*/
-	return Ke;
-}
-/*}}}*/
-ElementVector* Pengrid::PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax){/*{{{*/
-
-	IssmDouble h_max;
-	IssmDouble penalty_factor;
-	HydrologyDCInefficientAnalysis* inefanalysis = NULL;
-
-	/*Initialize Element matrix and return if necessary*/
-	if(!this->active) return NULL;
-	ElementVector* pe=new ElementVector(&node,1,this->parameters);
-	inefanalysis = new HydrologyDCInefficientAnalysis();
-
-	/*Retrieve parameters*/
-	parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
-
-	/*Get h_max and compute penalty*/
-	inefanalysis->GetHydrologyDCInefficientHmax(&h_max,element,node);
-
-	pe->values[0]=kmax*pow(10.,penalty_factor)*h_max;
-
-	/*Clean up and return*/
-	delete inefanalysis;
-	return pe;
-}
-/*}}}*/
-ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){/*{{{*/
-
-	IssmDouble moulin_load,dt;
-
-	/*Initialize Element matrix*/
-	ElementVector* pe=new ElementVector(&node,1,this->parameters);
-
-	this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
-	parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-
-	if(dt!=0.0) pe->values[0]=moulin_load*dt;
-
-	/*Clean up and return*/
-	return pe;
- }
-/*}}}*/
-void  Pengrid::ResetConstraint(void){/*{{{*/
+void           Pengrid::ResetConstraint(void){/*{{{*/
 	active         = 0;
 	zigzag_counter = 0;
 }
Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 18926)
@@ -39,45 +39,22 @@
 	_printf_("   dim: " << this->dim<< "\n");
 }
 /*}}}*/
-void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
 
-	switch(this->law){
-		case 1:
-			GetAlpha2Viscous(palpha2,gauss);
-			break;
-		case 2:
-			GetAlpha2Weertman(palpha2,gauss);
-			break;
-		case 3:
-			GetAlpha2Hydro(palpha2,gauss);
-			break;
-		case 4:
-			GetAlpha2Temp(palpha2,gauss);
-			break;
-		case 5:
-			GetAlpha2WaterLayer(palpha2,gauss);
-			break;
-		case 6:
-			GetAlpha2WeertmanTemp(palpha2,gauss);
-			break;
-	  default:
-			_error_("not supported");
-	}
+	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
+	 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
+	 * alpha_complement= Neff ^r * vel ^s*/
 
-}/*}}}*/
-void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+	if(this->law!=1)_error_("not supported");
 
-	/*This routine calculates the basal friction coefficient 
-	  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
-
 	/*diverse: */
 	IssmDouble  r,s;
-	IssmDouble  drag_p, drag_q;
+	IssmDouble  vx,vy,vz,vmag;
+	IssmDouble  drag_p,drag_q;
 	IssmDouble  Neff;
-	IssmDouble  thickness,bed;
-	IssmDouble  vx,vy,vz,vmag;
 	IssmDouble  drag_coefficient;
-	IssmDouble  alpha2;
+	IssmDouble  bed,thickness;
+	IssmDouble  alpha_complement;
 
 	/*Recover parameters: */
 	element->GetInputValue(&drag_p,FrictionPEnum);
@@ -97,6 +74,7 @@
 	Neff=gravity*(rho_ice*thickness+rho_water*bed);
 	if(Neff<0)Neff=0;
 
+	//We need the velocity magnitude to evaluate the basal stress:
 	switch(dim){
 		case 1:
 			element->GetInputValue(&vx,gauss,VxEnum);
@@ -118,53 +96,38 @@
 	}
 
 	/*Check to prevent dividing by zero if vmag==0*/
-	if(vmag==0. && (s-1.)<0.) alpha2=0.;
-	else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
-	_assert_(!xIsNan<IssmDouble>(alpha2));
+	if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
+	else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan<IssmDouble>(alpha_complement));
 
 	/*Assign output pointers:*/
-	*palpha2=alpha2;
-}/*}}}*/
-void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+	*palpha_complement=alpha_complement;
+}
+/*}}}*/
+void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 
-	/*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */
-
-	/*diverse: */
-	IssmDouble  C,m;
-	IssmDouble  vx,vy,vz,vmag;
-	IssmDouble  alpha2;
-
-	/*Recover parameters: */
-	element->GetInputValue(&C,gauss,FrictionCEnum);
-	element->GetInputValue(&m,FrictionMEnum);
-
-	switch(dim){
+	switch(this->law){
 		case 1:
-			element->GetInputValue(&vx,gauss,VxEnum);
-			vmag=sqrt(vx*vx);
+			GetAlpha2Viscous(palpha2,gauss);
 			break;
 		case 2:
-			element->GetInputValue(&vx,gauss,VxEnum);
-			element->GetInputValue(&vy,gauss,VyEnum);
-			vmag=sqrt(vx*vx+vy*vy);
+			GetAlpha2Weertman(palpha2,gauss);
 			break;
 		case 3:
-			element->GetInputValue(&vx,gauss,VxEnum);
-			element->GetInputValue(&vy,gauss,VyEnum);
-			element->GetInputValue(&vz,gauss,VzEnum);
-			vmag=sqrt(vx*vx+vy*vy+vz*vz);
+			GetAlpha2Hydro(palpha2,gauss);
 			break;
-		default:
+		case 4:
+			GetAlpha2Temp(palpha2,gauss);
+			break;
+		case 5:
+			GetAlpha2WaterLayer(palpha2,gauss);
+			break;
+		case 6:
+			GetAlpha2WeertmanTemp(palpha2,gauss);
+			break;
+	  default:
 			_error_("not supported");
 	}
 
-	/*Check to prevent dividing by zero if vmag==0*/
-	if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
-	else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
-	_assert_(!xIsNan<IssmDouble>(alpha2));
-
-	/*Assign output pointers:*/
-	*palpha2=alpha2;
 }/*}}}*/
 void Friction::GetAlpha2Hydro(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 
@@ -259,7 +222,7 @@
 	/*Assign output pointers:*/
 	*palpha2=alpha2;
 }/*}}}*/
-void Friction::GetAlpha2WaterLayer(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 
 	/*This routine calculates the basal friction coefficient 
 	  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
@@ -270,7 +233,7 @@
 	IssmDouble  Neff;
 	IssmDouble  thickness,bed;
 	IssmDouble  vx,vy,vz,vmag;
-	IssmDouble  drag_coefficient,water_layer;
+	IssmDouble  drag_coefficient;
 	IssmDouble  alpha2;
 
 	/*Recover parameters: */
@@ -279,7 +242,6 @@
 	element->GetInputValue(&thickness, gauss,ThicknessEnum);
 	element->GetInputValue(&bed, gauss,BaseEnum);
 	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
-	element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum);
 	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
 	IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
 	IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
@@ -289,9 +251,8 @@
 	s=1./drag_p;
 
 	//From bed and thickness, compute effective pressure when drag is viscous:
-	if(bed>0) bed=0;
-	Neff=gravity*(rho_ice*thickness+rho_water*(bed-water_layer));
-	if(Neff<0) Neff=0;
+	Neff=gravity*(rho_ice*thickness+rho_water*bed);
+	if(Neff<0)Neff=0;
 
 	switch(dim){
 		case 1:
@@ -321,49 +282,19 @@
 	/*Assign output pointers:*/
 	*palpha2=alpha2;
 }/*}}}*/
-void Friction::GetAlpha2WeertmanTemp(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
-	/*Here, we want to parameterize the friction as a function of temperature
-	 *
-	 * alpha2 = alpha2_weertman * 1/f(T)
-	 *
-	 * where f(T) = exp((T-Tpmp)/gamma)
-	 */
+void Friction::GetAlpha2WaterLayer(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 
-	/*Intermediaries: */
-	IssmDouble  f,T,pressure,Tpmp,gamma;
-	IssmDouble  alpha2;
+	/*This routine calculates the basal friction coefficient 
+	  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
 
-	/*Get viscous part*/
-	this->GetAlpha2Weertman(&alpha2,gauss);
-
-	/*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
-	element->GetInputValue(&T,gauss,TemperatureEnum);
-	element->GetInputValue(&pressure,gauss,PressureEnum);
-	Tpmp = element->TMeltingPoint(pressure);
-
-	/*Compute scaling parameter*/
-	element->parameters->FindParam(&gamma,FrictionGammaEnum);
-	alpha2 = alpha2 / exp((T-Tpmp)/gamma);
-
-	/*Assign output pointers:*/
-	*palpha2=alpha2;
-}/*}}}*/
-void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
-
-	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
-	 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
-	 * alpha_complement= Neff ^r * vel ^s*/
-
-	if(this->law!=1)_error_("not supported");
-
 	/*diverse: */
 	IssmDouble  r,s;
+	IssmDouble  drag_p, drag_q;
+	IssmDouble  Neff;
+	IssmDouble  thickness,bed;
 	IssmDouble  vx,vy,vz,vmag;
-	IssmDouble  drag_p,drag_q;
-	IssmDouble  Neff;
-	IssmDouble  drag_coefficient;
-	IssmDouble  bed,thickness;
-	IssmDouble  alpha_complement;
+	IssmDouble  drag_coefficient,water_layer;
+	IssmDouble  alpha2;
 
 	/*Recover parameters: */
 	element->GetInputValue(&drag_p,FrictionPEnum);
@@ -371,6 +302,7 @@
 	element->GetInputValue(&thickness, gauss,ThicknessEnum);
 	element->GetInputValue(&bed, gauss,BaseEnum);
 	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
+	element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum);
 	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
 	IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
 	IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
@@ -380,10 +312,10 @@
 	s=1./drag_p;
 
 	//From bed and thickness, compute effective pressure when drag is viscous:
-	Neff=gravity*(rho_ice*thickness+rho_water*bed);
-	if(Neff<0)Neff=0;
+	if(bed>0) bed=0;
+	Neff=gravity*(rho_ice*thickness+rho_water*(bed-water_layer));
+	if(Neff<0) Neff=0;
 
-	//We need the velocity magnitude to evaluate the basal stress:
 	switch(dim){
 		case 1:
 			element->GetInputValue(&vx,gauss,VxEnum);
@@ -405,10 +337,78 @@
 	}
 
 	/*Check to prevent dividing by zero if vmag==0*/
-	if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
-	else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan<IssmDouble>(alpha_complement));
+	if(vmag==0. && (s-1.)<0.) alpha2=0.;
+	else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
+	_assert_(!xIsNan<IssmDouble>(alpha2));
 
 	/*Assign output pointers:*/
-	*palpha_complement=alpha_complement;
-}
-/*}}}*/
+	*palpha2=alpha2;
+}/*}}}*/
+void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+
+	/*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */
+
+	/*diverse: */
+	IssmDouble  C,m;
+	IssmDouble  vx,vy,vz,vmag;
+	IssmDouble  alpha2;
+
+	/*Recover parameters: */
+	element->GetInputValue(&C,gauss,FrictionCEnum);
+	element->GetInputValue(&m,FrictionMEnum);
+
+	switch(dim){
+		case 1:
+			element->GetInputValue(&vx,gauss,VxEnum);
+			vmag=sqrt(vx*vx);
+			break;
+		case 2:
+			element->GetInputValue(&vx,gauss,VxEnum);
+			element->GetInputValue(&vy,gauss,VyEnum);
+			vmag=sqrt(vx*vx+vy*vy);
+			break;
+		case 3:
+			element->GetInputValue(&vx,gauss,VxEnum);
+			element->GetInputValue(&vy,gauss,VyEnum);
+			element->GetInputValue(&vz,gauss,VzEnum);
+			vmag=sqrt(vx*vx+vy*vy+vz*vz);
+			break;
+		default:
+			_error_("not supported");
+	}
+
+	/*Check to prevent dividing by zero if vmag==0*/
+	if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
+	else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));
+	_assert_(!xIsNan<IssmDouble>(alpha2));
+
+	/*Assign output pointers:*/
+	*palpha2=alpha2;
+}/*}}}*/
+void Friction::GetAlpha2WeertmanTemp(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+	/*Here, we want to parameterize the friction as a function of temperature
+	 *
+	 * alpha2 = alpha2_weertman * 1/f(T)
+	 *
+	 * where f(T) = exp((T-Tpmp)/gamma)
+	 */
+
+	/*Intermediaries: */
+	IssmDouble  f,T,pressure,Tpmp,gamma;
+	IssmDouble  alpha2;
+
+	/*Get viscous part*/
+	this->GetAlpha2Weertman(&alpha2,gauss);
+
+	/*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
+	element->GetInputValue(&T,gauss,TemperatureEnum);
+	element->GetInputValue(&pressure,gauss,PressureEnum);
+	Tpmp = element->TMeltingPoint(pressure);
+
+	/*Compute scaling parameter*/
+	element->parameters->FindParam(&gamma,FrictionGammaEnum);
+	alpha2 = alpha2 / exp((T-Tpmp)/gamma);
+
+	/*Assign output pointers:*/
+	*palpha2=alpha2;
+}/*}}}*/
Index: ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp	(revision 18925)
+++ ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp	(revision 18926)
@@ -107,52 +107,6 @@
 /*}}}*/
 
 /*Object virtual functions definitions:*/
-void Riftfront::Echo(void){/*{{{*/
-
-	_printf_("Riftfront:\n");
-	_printf_("   id: " << id << "\n");
-	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
-	_printf_("   hnodes: " << hnodes << "\n");
-	_printf_("   helements: " << helements << "\n");
-	_printf_("   hmatpar: " << hmatpar << "\n");
-	_printf_("   parameters: " << parameters << "\n");
-	_printf_("   internal parameters: \n");
-	_printf_("   normal: " << normal[0] << "|" << normal[1] << "\n");
-	_printf_("   length: " << length << "\n");
-	_printf_("   penalty_lock: " << penalty_lock << "\n");
-	_printf_("   active: " <<(active ? "true":"false") << "\n");
-	_printf_("   counter: " << counter << "\n");
-	_printf_("   prestable: " << (prestable ? "true":"false") << "\n");
-	_printf_("   material_converged: " << (material_converged ? "true":"false") << "\n");
-	_printf_("   fill: " << fill << "\n");
-	_printf_("   friction: " << friction << "\n");
-	_printf_("   fraction: " << fraction << "\n");
-	_printf_("   fractionincrement: " << fractionincrement << "\n");
-	_printf_("   state: " << state << "\n");
-	_printf_("   frozen: " << (frozen ? "true":"false") << "\n");
-
-}
-/*}}}*/
-void Riftfront::DeepEcho(void){/*{{{*/
-
-	_printf_("Riftfront:\n");
-	_printf_("   id: " << id << "\n");
-	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
-	hnodes->DeepEcho();
-	helements->DeepEcho();
-	hmatpar->DeepEcho();
-	_printf_("   parameters\n");
-	if(parameters)parameters->DeepEcho();
-}
-/*}}}*/
-int    Riftfront::Id(void){ return id; }/*{{{*/
-/*}}}*/
-int Riftfront::ObjectEnum(void){/*{{{*/
-
-	return RiftfrontEnum;
-
-}
-/*}}}*/
 Object* Riftfront::copy() {/*{{{*/
 
 	Riftfront* riftfront=NULL;
@@ -198,7 +152,53 @@
 
 }
 /*}}}*/
+void    Riftfront::DeepEcho(void){/*{{{*/
 
+	_printf_("Riftfront:\n");
+	_printf_("   id: " << id << "\n");
+	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
+	hnodes->DeepEcho();
+	helements->DeepEcho();
+	hmatpar->DeepEcho();
+	_printf_("   parameters\n");
+	if(parameters)parameters->DeepEcho();
+}
+/*}}}*/
+void    Riftfront::Echo(void){/*{{{*/
+
+	_printf_("Riftfront:\n");
+	_printf_("   id: " << id << "\n");
+	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
+	_printf_("   hnodes: " << hnodes << "\n");
+	_printf_("   helements: " << helements << "\n");
+	_printf_("   hmatpar: " << hmatpar << "\n");
+	_printf_("   parameters: " << parameters << "\n");
+	_printf_("   internal parameters: \n");
+	_printf_("   normal: " << normal[0] << "|" << normal[1] << "\n");
+	_printf_("   length: " << length << "\n");
+	_printf_("   penalty_lock: " << penalty_lock << "\n");
+	_printf_("   active: " <<(active ? "true":"false") << "\n");
+	_printf_("   counter: " << counter << "\n");
+	_printf_("   prestable: " << (prestable ? "true":"false") << "\n");
+	_printf_("   material_converged: " << (material_converged ? "true":"false") << "\n");
+	_printf_("   fill: " << fill << "\n");
+	_printf_("   friction: " << friction << "\n");
+	_printf_("   fraction: " << fraction << "\n");
+	_printf_("   fractionincrement: " << fractionincrement << "\n");
+	_printf_("   state: " << state << "\n");
+	_printf_("   frozen: " << (frozen ? "true":"false") << "\n");
+
+}
+/*}}}*/
+int     Riftfront::Id(void){ return id; }/*{{{*/
+/*}}}*/
+int     Riftfront::ObjectEnum(void){/*{{{*/
+
+	return RiftfrontEnum;
+
+}
+/*}}}*/
+
 /*Update virtual functions definitions:*/
 void  Riftfront::InputUpdateFromConstant(bool constant,int name){/*{{{*/
 }
@@ -233,28 +233,46 @@
 
 }
 /*}}}*/
-void  Riftfront::ResetHooks(){/*{{{*/
+void  Riftfront::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
+	/*do nothing: */
+	return;
+}
+/*}}}*/
+void  Riftfront::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
+	/*do nothing: */
+	return;
+}
+/*}}}*/
+void  Riftfront::GetNodesLidList(int* lidlist){/*{{{*/
 
-	this->nodes=NULL;
-	this->elements=NULL;
-	this->matpar=NULL;
-	this->parameters=NULL;
+	_assert_(lidlist);
+	_assert_(nodes);
 
-	/*Get Element type*/
-	this->hnodes->reset();
-	this->helements->reset();
-	this->hmatpar->reset();
+	for(int i=0;i<NUMVERTICES;i++) lidlist[i]=nodes[i]->Lid();
+}
+/*}}}*/
+void  Riftfront::GetNodesSidList(int* sidlist){/*{{{*/
 
+	_assert_(sidlist);
+	_assert_(nodes);
+
+	for(int i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->Sid();
 }
 /*}}}*/
+int   Riftfront::GetNumberOfNodes(void){/*{{{*/
+
+	return NUMVERTICES;
+}
+/*}}}*/
+bool  Riftfront::InAnalysis(int in_analysis_type){/*{{{*/
+	if (in_analysis_type==this->analysis_type) return true;
+	else return false;
+}
+/*}}}*/
 bool  Riftfront::IsPenalty(void){/*{{{*/
 	return true;
 }
 /*}}}*/
-void  Riftfront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
-
-}
-/*}}}*/
 void  Riftfront::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
 
 	/*Retrieve parameters: */
@@ -305,42 +323,24 @@
 	}
 }
 /*}}}*/
-void  Riftfront::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
-	/*do nothing: */
-	return;
-}
-/*}}}*/
-void  Riftfront::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
-	/*do nothing: */
-	return;
-}
-/*}}}*/
-void  Riftfront::GetNodesSidList(int* sidlist){/*{{{*/
+void  Riftfront::ResetHooks(){/*{{{*/
 
-	_assert_(sidlist);
-	_assert_(nodes);
+	this->nodes=NULL;
+	this->elements=NULL;
+	this->matpar=NULL;
+	this->parameters=NULL;
 
-	for(int i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->Sid();
-}
-/*}}}*/
-void  Riftfront::GetNodesLidList(int* lidlist){/*{{{*/
+	/*Get Element type*/
+	this->hnodes->reset();
+	this->helements->reset();
+	this->hmatpar->reset();
 
-	_assert_(lidlist);
-	_assert_(nodes);
-
-	for(int i=0;i<NUMVERTICES;i++) lidlist[i]=nodes[i]->Lid();
 }
 /*}}}*/
-int   Riftfront::GetNumberOfNodes(void){/*{{{*/
+void  Riftfront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
 
-	return NUMVERTICES;
 }
 /*}}}*/
-bool  Riftfront::InAnalysis(int in_analysis_type){/*{{{*/
-	if (in_analysis_type==this->analysis_type) return true;
-	else return false;
-}
-/*}}}*/
 void  Riftfront::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
 
 	/*Output */
