Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 17443)
+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 17444)
@@ -123,10 +123,10 @@
 
 	iomodel->FetchData(1,MeshVertexonbedEnum);
 
-	//create penalties for nodes: no node can have a temperature over the melting point
+	//create penalties for nodes: no node can have water above the max
 	CreateSingleNodeToElementConnectivity(iomodel);
 	for(int i=0;i<iomodel->numberofvertices;i++){
-		if (iomodel->meshtype==Mesh3DEnum){
+		if (iomodel->meshtype!=Mesh3DEnum){
 			/*keep only this partition's nodes:*/
 			if(iomodel->my_vertices[i]){
 				loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
@@ -279,7 +279,7 @@
 	/*Intermediaries */
 	bool       active_element,isefficientlayer;
 	IssmDouble dt,scalar;
-	IssmDouble moulin_load,water_head;
+	IssmDouble water_head;
 	IssmDouble water_load,transfer;
 	IssmDouble Jdet;
 
@@ -306,7 +306,6 @@
 	Input* thickness_input   = basalelement->GetInput(ThicknessEnum);
 	Input* bed_input         = basalelement->GetInput(BedEnum);
 	Input* water_input       = basalelement->GetInput(BasalforcingsMeltingRateEnum);    _assert_(water_input);
-	Input* moulin_input      = basalelement->GetInput(HydrologydcBasalMoulinInputEnum); _assert_(moulin_input);
 	if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);             _assert_(old_wh_input);}
 
 	IssmDouble sediment_storing    = SedimentStoring(basalelement);
@@ -326,10 +325,8 @@
 
 		/*Loading term*/
 		water_input->GetInputValue(&water_load,gauss);
-		moulin_input->GetInputValue(&moulin_load,gauss);
-	
+
 		scalar = Jdet*gauss->weight*(water_load);
-		scalar = scalar + Jdet* moulin_load;
 		if(dt!=0.) scalar = scalar*dt;
 		for(int i=0;i<numnodes;i++){
 			pe->values[i]+=scalar*basis[i];
Index: ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 17443)
+++ ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 17444)
@@ -75,7 +75,11 @@
 		ElementVector* pe = analysis->CreatePVector(element);
 		element->ReduceMatrices(Ke,pe);
 		if(Ke) Ke->AddToGlobal(Kff,Kfs);
-		if(pe) pe->AddToGlobal(pf);
+		if(pe){
+			pe->AddToGlobal(pf);
+			/* printf("-------------------USUAL \n"); */
+			/* pf->Echo(); */
+		}
 		delete Ke;
 		delete pe;
 	}
@@ -86,6 +90,8 @@
 		if(load->InAnalysis(configuration_type)){
 			load->CreateKMatrix(Kff,Kfs);
 			load->CreatePVector(pf);
+			/* printf("-------------------LOADING \n"); */
+			/* pf->Echo(); */
 		}
 	}
 
@@ -96,6 +102,8 @@
 			if(load->InAnalysis(configuration_type)){
 				load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
 				load->PenaltyCreatePVector(pf,kmax);
+				/* printf("-------------------PENALTY \n"); */
+				/* pf->Echo(); */
 			}
 		}
 	}
Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 17443)
+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 17444)
@@ -180,9 +180,24 @@
 /*FUNCTION Pengrid::CreatePVector {{{*/
 void  Pengrid::CreatePVector(Vector<IssmDouble>* pf){
 
-	/*No loads applied, do nothing: */
-	return;
+	ElementVector* pe=NULL;
+	int analysis_type;
+	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
+	switch(analysis_type){
+		
+	case HydrologyDCInefficientAnalysisEnum:
+		pe = CreatePVectorHydrologyDCInefficient();
+		break;
+	default:
+		/*No loads applied, do nothing: */
+		return;
+	}
+	if(pe){
+		pe->AddToGlobal(pf);
+		delete pe;
+	}
+
 }
 /*}}}*/
 /*FUNCTION Pengrid::GetNodesSidList{{{*/
@@ -792,6 +807,23 @@
 	return pe;
 }
 /*}}}*/
+/*FUNCTION Pengrid::CreatePVectorHydrologyDCInefficient {{{*/
+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;
+ }
+/*}}}*/
 /*FUNCTION Pengrid::ResetConstraint {{{*/
 void  Pengrid::ResetConstraint(void){
 	active         = 0;
Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 17443)
+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 17444)
@@ -93,6 +93,7 @@
 		ElementVector* PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax);
 		void           ConstraintActivateHydrologyDCInefficient(int* punstable);
 		void  ConstraintActivate(int* punstable);
+		ElementVector* CreatePVectorHydrologyDCInefficient(void);
 		void  ResetConstraint(void);
 		/*}}}*/