Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 21478)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 21479)
@@ -347,5 +347,4 @@
 		basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		basalelement ->NodalFunctions(basis,gauss);
-
 		epl_storing	= EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input);
 		/*Loading term*/
@@ -570,11 +569,9 @@
 				/*Compute first the effective pressure in the EPL*/
 				EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
-				//				if(EPL_N<0.0)EPL_N=0.0;
+				if(EPL_N<0.0)EPL_N=0.0;
 				/*Get then the square of the gradient of EPL heads*/
 				EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i])+(epl_slopeY[i]*epl_slopeY[i]);
 				/*And proceed to the real thing*/
-				thickness[i] = old_thickness[i]*
-					(2.0+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-(2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))))/
-					(2.0-((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2+(2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n))));
+				thickness[i] = old_thickness[i]/(1.0-((rho_water*gravity*epl_conductivity*EPLgrad2*dt)/(rho_ice*latentheat))+((2.0*A*dt*pow(EPL_N,n))/(pow(n,n))));
 				/*Take care of otherthikening*/
 				if(thickness[i]>max_thick){
@@ -735,18 +732,18 @@
 	storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity));
 
-	//porosity for unconfined region
-	if (water_sheet<=0.9*epl_thickness){
-		epl_storing=epl_porosity;
-	}
-	//continuity ramp
-	else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){
-		epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing;
-	}
-	//storing coefficient for confined
-	else{
-		epl_storing=storing;
-	}
-	//return storing;
- 	return epl_storing;
+	/* //porosity for unconfined region */
+	/* if (water_sheet<=0.9*epl_thickness){ */
+	/* 	epl_storing=epl_porosity; */
+	/* } */
+	/* //continuity ramp */
+	/* else if((water_sheet<epl_thickness) && (water_sheet>0.9*epl_thickness)){ */
+	/* 	epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */
+	/* } */
+	/* //storing coefficient for confined */
+	/* else{ */
+	/* 	epl_storing=storing; */
+	/* } */
+	return storing;
+ 	//return epl_storing;
 }/*}}}*/
 
@@ -760,9 +757,8 @@
 	base_input->GetInputValue(&base_elev,gauss);
 
-	//water_sheet=max(0.0,(prestep_head-base_elev));
-	water_sheet=prestep_head-base_elev;
+	water_sheet=max(0.0,(prestep_head-base_elev));
+	//water_sheet=prestep_head-base_elev;
 	
 	//epl_transmitivity=epl_conductivity*epl_thickness;
-	//if (element->Sid()==8)printf("water sheet is %f for Elt %i\n",water_sheet,element->Sid());
 	epl_transmitivity=epl_conductivity*min(water_sheet,epl_thickness);
 	return epl_transmitivity;