Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16544)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16545)
@@ -3824,5 +3824,4 @@
 
 		/*Conduction: */  
-		/*Need to change that depending on enthalpy value -> cold or temperate ice: */  
 		GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
 
@@ -3873,5 +3872,5 @@
 		}
 
-		/*Artifficial diffusivity*/
+		/*Artificial diffusivity*/
 		if(stabilization==1){
 			/*Build K: */
@@ -3879,7 +3878,9 @@
 			vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
 			h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
+
 			K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
 			K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
 			K[2][0]=h/(2*vel)*vz*vx;  K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz;
+
 			D_scalar_stab=gauss->weight*Jdet;
 			if(reCast<bool,IssmDouble>(dt)) D_scalar_stab=D_scalar_stab*dt;
@@ -4070,5 +4071,4 @@
 
 		/*Advection: */
-
 		GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); 
 		GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 
@@ -4960,6 +4960,7 @@
 	IssmDouble  vx[NUMVERTICES],vy[NUMVERTICES],vz[NUMVERTICES];
 	IssmDouble  geothermalflux[NUMVERTICES];
-	IssmDouble  dt;
-	IssmDouble  meltingrate_enthalpy;
+	IssmDouble  dt, yts;
+	IssmDouble  melting_overshoot,meltingrate_enthalpy;
+	IssmDouble  lambda,heating;
 	Friction   *friction  = NULL;
 
@@ -4998,6 +4999,9 @@
 		checkpositivethickness=true;
 
+		_assert_(watercolumn[iv]>=0.);
+
 		/*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
 		meltingrate_enthalpy=0.;
+		heating=0.;
 		if((watercolumn[iv]>0.) && (enthalpy[iv]<matpar->PureIceEnthalpy(pressure[iv]))){
 			/*ensure that no ice is at T<Tm(p), if water layer present*/
@@ -5032,14 +5036,28 @@
 			matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]);
 			// -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66
-			meltingrate_enthalpy=(heatflux+basalfriction+geothermalflux[iv])/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent 
-		}
+			heating=(heatflux+basalfriction+geothermalflux[iv]);
+			meltingrate_enthalpy=heating/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent 
+		}
+
 		/*Update water column, basal meltingrate*/
-		//basalmeltingrate[iv]+=meltingrate_enthalpy;
-		basalmeltingrate[iv]=meltingrate_enthalpy;
 		this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-		if(reCast<bool,IssmDouble>(dt))
-			watercolumn[iv]+=dt*meltingrate_enthalpy; 
-		else
+		if(reCast<bool,IssmDouble>(dt)){
+			if(watercolumn[iv]+meltingrate_enthalpy*dt<0.){				
+				melting_overshoot=watercolumn[iv]+meltingrate_enthalpy*dt;
+				lambda=melting_overshoot/(meltingrate_enthalpy*dt); _assert_(lambda>0); _assert_(lambda<1);
+				basalmeltingrate[iv]=(1.-lambda)*meltingrate_enthalpy;
+				watercolumn[iv]=0.;
+				yts=365*24*60*60;
+				enthalpy[iv]+=dt/yts*lambda*heating;
+			}
+			else{
+				basalmeltingrate[iv]=meltingrate_enthalpy;
+				watercolumn[iv]+=dt*meltingrate_enthalpy; 
+			}
+		}
+		else{
+			basalmeltingrate[iv]=meltingrate_enthalpy;
 			watercolumn[iv]+=meltingrate_enthalpy;
+		}	  
 	}  
 	/*feed updated variables back into model*/
@@ -5072,6 +5090,6 @@
 	pentabase=this->GetBasalElement();
 	
-	GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
-	GetInputListOnVertices(&pressure[0],PressureEnum);
+	this->GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
+	this->GetInputListOnVertices(&pressure[0],PressureEnum);
 	pentabase->GetInputListOnVertices(&watercolumnbase[0], WatercolumnEnum);
 
@@ -5089,4 +5107,5 @@
 		watercolumnbase[iv%NUMVERTICES2D]+=dwc;
 	}
+
 	/*feed updated results back into model*/
 	this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
@@ -5094,5 +5113,4 @@
 	pentabase->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumnbase,P1Enum));
 
-	delete pentabase;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16544)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16545)
@@ -377,10 +377,15 @@
 /*FUNCTION Matpar::GetEnthalpyDiffusionParameter{{{*/
 IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
-	if(enthalpy<PureIceEnthalpy(pressure)){
+	/*if (enthalpy<=PureIceEnthalpy(pressure))
 		return thermalconductivity/heatcapacity;
-	}
-	else{
-		return 1.045*1.e-4; // K0=1.045*1e-4 from Aschwanden 2012. TODO: fetch K0 from model
-	}
+	else
+		return thermalconductivity/heatcapacity*pow(10.,-5);
+	*/
+
+	IssmDouble eps=0.05*heatcapacity;
+	IssmDouble hpmp=PureIceEnthalpy(pressure);
+	IssmDouble kappa_c=thermalconductivity/heatcapacity;
+	IssmDouble kappa_t=thermalconductivity/heatcapacity*pow(10.,-1);
+	return 1./(1.+exp(-(enthalpy-(hpmp))/eps))*(kappa_t-kappa_c) + kappa_c;
 }
 /*}}}*/