Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 22235)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 22236)
@@ -240,4 +240,5 @@
 		/*Diffusivity*/
 		D_scalar=epl_transmitivity*gauss->weight*Jdet;
+		//D_scalar=gauss->weight*Jdet;
 		if(dt!=0.) D_scalar=D_scalar*dt;
 		D[0][0]=D_scalar;
@@ -253,4 +254,5 @@
 			basalelement->NodalFunctions(&basis[0],gauss);
 			D_scalar=epl_storing*gauss->weight*Jdet;
+			//D_scalar=(epl_storing/epl_transmitivity)*gauss->weight*Jdet;
 			TripleMultiply(basis,numnodes,1,0,
 						&D_scalar,1,1,0,
@@ -262,4 +264,5 @@
 			transfer=GetHydrologyKMatrixTransfer(basalelement);
 			D_scalar=dt*transfer*gauss->weight*Jdet;
+			//D_scalar=dt*(transfer/epl_transmitivity)*gauss->weight*Jdet;
 			TripleMultiply(basis,numnodes,1,0,
 										 &D_scalar,1,1,0,
@@ -312,5 +315,5 @@
 	IssmDouble dt,scalar,water_head;
 	IssmDouble water_load,transfer;
-	IssmDouble epl_storing;
+	IssmDouble epl_storing,epl_transmitivity;
 	IssmDouble Jdet;
 	IssmDouble residual,connectivity;
@@ -348,7 +351,10 @@
 		basalelement ->NodalFunctions(basis,gauss);
 		epl_storing	= EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
+		epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
+
 		/*Loading term*/
 		water_input->GetInputValue(&water_load,gauss);
 		scalar = Jdet*gauss->weight*(water_load);
+		//scalar = Jdet*gauss->weight*(water_load)/epl_transmitivity;
 		if(dt!=0.) scalar = scalar*dt;
 		for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
@@ -360,4 +366,5 @@
 			transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input);
 			scalar = Jdet*gauss->weight*((water_head*epl_storing)+(dt*transfer));
+			//scalar = Jdet*gauss->weight*((water_head*epl_storing)+(dt*transfer))/epl_transmitivity;
 			for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
 		}
@@ -369,7 +376,9 @@
 	for(int iv=0;iv<numvertices;iv++){
 		gauss->GaussVertex(iv);
+		epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
 		connectivity = IssmDouble(basalelement->VertexConnectivity(iv));
 		residual_input->GetInputValue(&residual,gauss);
 		pe->values[iv]+=residual/connectivity;
+		//pe->values[iv]+=residual/(epl_transmitivity*connectivity);
 	}
 	/*Clean up and return*/
@@ -415,4 +424,6 @@
 	/*Fetch dof list and allocate solution vector*/
 	IssmDouble* eplHeads    = xNew<IssmDouble>(numnodes);
+	IssmDouble* basevalue    = xNew<IssmDouble>(numnodes);
+	IssmDouble* initvalue    = xNew<IssmDouble>(numnodes);
 	basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 
@@ -422,4 +433,5 @@
 	/*Use the dof list to index into the solution vector: */
 	/*If the EPL is not active we revert to the bedrock elevation*/
+
 	if(active_element){
 		for(int i=0;i<numnodes;i++){
@@ -430,6 +442,9 @@
 	}
 	else{
-		basalelement->GetInputListOnVertices(&eplHeads[0],BaseEnum);
+		//tradeof between keeping initial condition and not getting to far from head at deactivation
+		basalelement->GetInputListOnVertices(&basevalue[0],BaseEnum);
+		basalelement->GetInputListOnVertices(&initvalue[0],EplHeadEnum);
 		for(int i=0;i<numnodes;i++){
+			eplHeads[i]=max(basevalue[i],initvalue[i]);
 			if(xIsNan<IssmDouble>(eplHeads[i])) _error_("NaN found in solution vector");
 			if(xIsInf<IssmDouble>(eplHeads[i])) _error_("Inf found in solution vector");
@@ -440,4 +455,6 @@
 	/*Free ressources:*/
 	xDelete<IssmDouble>(eplHeads);
+	xDelete<IssmDouble>(basevalue);
+	xDelete<IssmDouble>(initvalue);
 	xDelete<int>(doflist);
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
@@ -464,5 +481,5 @@
 	case 1:
 		element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
-		transfer=leakage;
+		transfer=+leakage;
 		break;
 	default:
@@ -489,5 +506,5 @@
 		sed_head_input->GetInputValue(&sediment_head,gauss);
 		element->FindParam(&leakage,HydrologydcLeakageFactorEnum);
-		transfer=sediment_head*leakage;
+		transfer=+sediment_head*leakage;
 		break;
 	default:
@@ -558,5 +575,5 @@
 		element->GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
 		element->GetInputListOnVertices(&bed[0],BaseEnum);
-		
+
 		if(!active_element){
 			/*Keeping thickness to initial value if EPL is not active*/
@@ -731,5 +748,4 @@
 	base_input->GetInputValue(&base_elev,gauss);
 	water_sheet=max(0.0,(prestep_head-base_elev));
-
 	storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));
 
@@ -760,5 +776,4 @@
 
 	water_sheet=max(0.0,(prestep_head-base_elev));
-	
 	epl_transmitivity=epl_conductivity*epl_thickness;
 	//epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness)));