Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 17165)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 17166)
@@ -821,25 +821,39 @@
 
 /*Modules*/
-/*{{{*/
-void EnthalpyAnalysis::PostProcessing(FemModel* femmodel){
+void EnthalpyAnalysis::PostProcessing(FemModel* femmodel){/*{{{*/
+
 	/*Intermediaries*/
-	int solution_type;
-
-	/*Compute basal melting rates: */
-	for(int i=0;i<femmodel->elements->Size();i++){
-		Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		ComputeBasalMeltingrate(element);
-	}
-
-	/*Update basal dirichlet BCs for enthalpy in transient runs: */
-	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-	if(solution_type==TransientSolutionEnum){
-		for(int i=0;i<femmodel->elements->Size();i++){
+	int solution_type, i;
+	bool computebasalmeltingrates=true;
+	bool isdrainage=true;
+	bool updatebasalconstraints=true;
+
+	if(isdrainage){
+		/*Drain excess water fraction in ice column: */
+		for(i=0;i<femmodel->elements->Size();i++){
 			Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-			UpdateBasalConstraints(element);
-		}
-	}
-}/*}}}*/
-
+			DrainWaterfractionIcecolumn(element);
+		}
+	}
+
+	if(computebasalmeltingrates){
+		/*Compute basal melting rates: */
+		for(i=0;i<femmodel->elements->Size();i++){
+			Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+			ComputeBasalMeltingrate(element);
+		}
+	}
+
+	if(updatebasalconstraints){
+		/*Update basal dirichlet BCs for enthalpy in transient runs: */
+		femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+		if(solution_type==TransientSolutionEnum){
+			for(i=0;i<femmodel->elements->Size();i++){
+				Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+				UpdateBasalConstraints(element);
+			}
+		}
+	}
+}/*}}}*/
 void EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
 	/*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
@@ -946,21 +960,4 @@
 		}
 	}
-
-	/******** DRAINAGE *****************************************/
-	IssmDouble* drainrate_column  = xNew<IssmDouble>(numsegments);
-	IssmDouble* drainrate_element = xNew<IssmDouble>(numsegments);
-	for(is=0;is<numsegments;is++)	drainrate_column[is]=0.;
-	Element* elementi = element;
-	for(;;){
-		for(is=0;is<numsegments;is++)	drainrate_element[is]=0.;
-		DrainWaterfraction(elementi,drainrate_element); // TODO: make sure every vertex is only drained once
-		for(is=0;is<numsegments;is++)	drainrate_column[is]+=drainrate_element[is];
-
-		if(elementi->IsOnSurface()) break;
-		elementi=elementi->GetUpperElement();			
-	}
-	/* add drained water to melting rate*/
-	for(is=0;is<numsegments;is++) meltingrate_enthalpy[is]+=drainrate_column[is];
-
 	/******** UPDATE MELTINGRATES AND WATERCOLUMN **************/
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
@@ -970,10 +967,8 @@
 		if(dt!=0.){
 			if(watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt<0.){	// prevent too much freeze on			
-				melting_overshoot=watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt;
-				lambda=melting_overshoot/(meltingrate_enthalpy[is]*dt); _assert_(lambda>0); _assert_(lambda<1);
-				basalmeltingrate[vertexdown]=(1.-lambda)*meltingrate_enthalpy[is];
+				lambda = -watercolumn[vertexdown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
 				watercolumn[vertexdown]=0.;
-				yts=365.*24.*60.*60.;
-				enthalpy[vertexdown]+=dt/yts*lambda*heating[is];
+				basalmeltingrate[vertexdown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
+				enthalpy[vertexdown]+=(1.-lambda)*meltingrate_enthalpy[is]*dt*latentheat; // use rest of energy to cool down base
 			}
 			else{
@@ -1004,9 +999,44 @@
 	xDelete<IssmDouble>(meltingrate_enthalpy);
 	xDelete<IssmDouble>(heating);
+	xDelete<IssmDouble>(xyz_list);
+}/*}}}*/
+void EnthalpyAnalysis::DrainWaterfractionIcecolumn(Element* element){/*{{{*/
+
+	/* Only drain waterfraction of ice column from element at base*/
+	if(!element->IsOnBed()) return; //FIXME: allow freeze on for floating elements
+
+	/* Intermediaries*/
+	int is, numvertices, numsegments;
+	int *pairindices   = NULL;
+
+	numvertices=element->GetNumberOfVertices();
+	element->VerticalSegmentIndices(&pairindices,&numsegments);
+
+	IssmDouble* watercolumn = xNew<IssmDouble>(numvertices);
+	IssmDouble* drainrate_column  = xNew<IssmDouble>(numsegments);
+	IssmDouble* drainrate_element = xNew<IssmDouble>(numsegments);
+
+	element->GetInputListOnVertices(watercolumn,WatercolumnEnum);
+
+	for(is=0;is<numsegments;is++)	drainrate_column[is]=0.;
+	Element* elementi = element;
+	for(;;){
+		for(is=0;is<numsegments;is++)	drainrate_element[is]=0.;
+		DrainWaterfraction(elementi,drainrate_element); // TODO: make sure every vertex is only drained once
+		for(is=0;is<numsegments;is++)	drainrate_column[is]+=drainrate_element[is];
+
+		if(elementi->IsOnSurface()) break;
+		elementi=elementi->GetUpperElement();			
+	}
+	/* add drained water to water column*/
+	for(is=0;is<numsegments;is++) watercolumn[is]+=drainrate_column[is];
+	/* Feed updated water column back into model */
+	element->AddInput(WatercolumnEnum,watercolumn,P1Enum);
+
+	delete pairindices;
 	xDelete<IssmDouble>(drainrate_column);
 	xDelete<IssmDouble>(drainrate_element);
-	xDelete<IssmDouble>(xyz_list);
-}/*}}}*/
-
+	xDelete<IssmDouble>(watercolumn);
+}/*}}}*/
 void EnthalpyAnalysis::DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element){/*{{{*/
 
@@ -1056,4 +1086,5 @@
 		height_element=fabs(xyz_list[vertexup*3+2]-xyz_list[vertexdown*3+2]);
 		pdrainrate_element[is]=(deltawaterfractions[vertexdown]+deltawaterfractions[vertexup])/2.*height_element; // return water equivalent of drainage
+		_assert_(pdrainrate_element[is]>=0.);
 	}
 
@@ -1066,5 +1097,4 @@
 	xDelete<IssmDouble>(deltawaterfractions);
 }/*}}}*/
-
 void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
 
@@ -1120,11 +1150,9 @@
 			setspc = false;
 
-
-
 		node=element->GetNode(indices[i]);
 		if (setspc) 
-			node->ApplyConstraint(1,h_pmp); /*apply spc*/ //nodes[indices[i]]->ApplyConstraint(1,h_pmp);
+			node->ApplyConstraint(1,h_pmp); /*apply spc*/ 
 		else			
-			node->DofInFSet(0); /*remove spc*/ //nodes[indices[i]]->DofInFSet(0);
+			node->DofInFSet(0); /*remove spc*/ 
 	}
 
Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 17165)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 17166)
@@ -40,4 +40,5 @@
 		static void PostProcessing(FemModel* femmodel);
 		static void ComputeBasalMeltingrate(Element* element);
+		static void DrainWaterfractionIcecolumn(Element* element);
 		static void DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element);
 		static void UpdateBasalConstraints(Element* element);