Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 17443)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 17444)
@@ -124,8 +124,8 @@
 	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]){
@@ -280,5 +280,5 @@
 	bool       active_element,isefficientlayer;
 	IssmDouble dt,scalar;
-	IssmDouble moulin_load,water_head;
+	IssmDouble water_head;
 	IssmDouble water_load,transfer;
 	IssmDouble Jdet;
@@ -307,5 +307,4 @@
 	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);}
 
@@ -327,8 +326,6 @@
 		/*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++){
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 17443)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 17444)
@@ -181,6 +181,21 @@
 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;
+	}
 
 }
@@ -793,4 +808,21 @@
 }
 /*}}}*/
+/*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){
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 17443)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 17444)
@@ -94,4 +94,5 @@
 		void           ConstraintActivateHydrologyDCInefficient(int* punstable);
 		void  ConstraintActivate(int* punstable);
+		ElementVector* CreatePVectorHydrologyDCInefficient(void);
 		void  ResetConstraint(void);
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 17443)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 17444)
@@ -76,5 +76,9 @@
 		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;
@@ -87,4 +91,6 @@
 			load->CreateKMatrix(Kff,Kfs);
 			load->CreatePVector(pf);
+			/* printf("-------------------LOADING \n"); */
+			/* pf->Echo(); */
 		}
 	}
@@ -97,4 +103,6 @@
 				load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
 				load->PenaltyCreatePVector(pf,kmax);
+				/* printf("-------------------PENALTY \n"); */
+				/* pf->Echo(); */
 			}
 		}