Index: /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp	(revision 21144)
+++ /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp	(revision 21145)
@@ -8,4 +8,6 @@
 #define OMEGA 0.001    // parameter controlling transition to nonlinear resistance in basal system (dimensionless)
 #define NU    1.787e-6 //kinematic water viscosity m^2/s
+#define CT    7.5e-8  // Clapeyron slope (K/Pa) 
+#define CW    4.22e3   // specific heat capacity of water (J/kg/K)
 
 /*Model processing*/
@@ -176,4 +178,7 @@
 	/*Get conductivity from inputs*/
 	IssmDouble conductivity = GetConductivity(element);
+//if(element->Id()==1){
+//	printf("Conductivity at CreateKMatrix: %g \n",conductivity);
+//}
 
 	/* Start  looping on the number of gaussian points: */
@@ -208,4 +213,5 @@
 	IssmDouble  lr,br,vx,vy,beta;
 	IssmDouble  alpha2,frictionheat;
+   IssmDouble  PMPheat,dpressure_water[2],dbed[2];	
 	IssmDouble* xyz_list = NULL;
 
@@ -238,4 +244,7 @@
 	/*Get conductivity from inputs*/
 	IssmDouble conductivity = GetConductivity(element);
+//if(element->Id()==1){
+//	printf("Conductivity in CreatePVector: %g \n",conductivity);
+//}
 
 	/*Build friction element, needed later: */
@@ -251,4 +260,5 @@
 		geothermalflux_input->GetInputValue(&G,gauss);
 		base_input->GetInputValue(&bed,gauss);
+		base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss);
 		thickness_input->GetInputValue(&thickness,gauss);
 		gap_input->GetInputValue(&gap,gauss);
@@ -283,7 +293,12 @@
 		if(pressure_water>pressure_ice) pressure_water = pressure_ice;
 
-		meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
+		/*Compute change in sensible heat due to changes in pressure melting point*/
+   	dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
+		dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
+		PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
+
+		meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
 		_assert_(meltrate>0.);
-
+		
 		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
 		 (
@@ -315,4 +330,7 @@
 	IssmDouble* xyz_list = NULL;
 
+	/*Get gravity from parameters*/
+	   IssmDouble  g = element->GetMaterialParameter(ConstantsGEnum);
+
 	/*Fetch number of nodes for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
@@ -323,4 +341,5 @@
 
 	/*Get thickness and base on nodes to apply cap on water head*/
+   IssmDouble* eff_pressure = xNew<IssmDouble>(numnodes);
 	IssmDouble* thickness = xNew<IssmDouble>(numnodes);
 	IssmDouble* bed       = xNew<IssmDouble>(numnodes);
@@ -344,4 +363,7 @@
 		}
 
+		/*Calculate effective pressure*/
+		eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]);
+	
 		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
 		if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
@@ -350,4 +372,5 @@
 	/*Add input to the element: */
 	element->AddInput(HydrologyHeadEnum,values,element->GetElementType());
+   element->AddInput(EffectivePressureEnum,eff_pressure,P1Enum);
 
 	/*Update reynolds number according to new solution*/
@@ -356,4 +379,8 @@
 	head_input->GetInputDerivativeAverageValue(&dh[0],xyz_list);
 	IssmDouble conductivity = GetConductivity(element);
+//if(element->Id()==1){
+//	printf("Conductivity in UpdateInputFromSolution: %g \n",conductivity);
+//}
+
 	IssmDouble reynolds = conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1])/(2.*NU);
 	element->AddInput(HydrologyReynoldsEnum,&reynolds,P0Enum);
@@ -365,4 +392,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<int>(doflist);
+	xDelete<IssmDouble>(eff_pressure);
 }/*}}}*/
 void           HydrologySommersAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
@@ -385,5 +413,5 @@
 	reynolds_input->GetInputAverage(&reynolds);
 	gap_input->GetInputAverage(&gap);
-
+	
 	/*Compute conductivity*/
 	IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds));
@@ -414,5 +442,5 @@
 	IssmDouble  alpha2,frictionheat;
 	IssmDouble* xyz_list = NULL;
-
+   IssmDouble  dpressure_water[2],dbed[2],PMPheat;
 
 	/*Retrieve all inputs and parameters*/
@@ -438,5 +466,7 @@
 	/*Get conductivity from inputs*/
 	IssmDouble conductivity = GetConductivity(element);
-
+//if(element->Id()==1){
+//	printf("Conductivity at gap update: %g \n",conductivity);
+//}
 	/*Build friction element, needed later: */
 	Friction* friction=new Friction(element,2);
@@ -454,4 +484,5 @@
 		geothermalflux_input->GetInputValue(&G,gauss);
 		base_input->GetInputValue(&bed,gauss);
+		base_input->GetInputDerivativeValue(&dbed[0],xyz_list,gauss);
 		thickness_input->GetInputValue(&thickness,gauss);
 		gap_input->GetInputValue(&gap,gauss);
@@ -486,6 +517,10 @@
 		if(pressure_water>pressure_ice) pressure_water = pressure_ice;
       
-
-		meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1]));
+      /* Compute change in sensible heat due to changes in pressure melting point*/
+	   dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
+		dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
+		PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
+	
+		meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
 		_assert_(meltrate>0.);
 
Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 21144)
+++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 21145)
@@ -88,6 +88,6 @@
 		if(save_results){
 			if(VerboseSolution()) _printf0_("   saving results \n");
-			int outputs[2] = {HydrologyHeadEnum,HydrologyGapHeightEnum};
-			femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+			int outputs[3] = {HydrologyHeadEnum,HydrologyGapHeightEnum,EffectivePressureEnum};
+			femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
 			
 			/*unload results*/