Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 17014)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 17015)
@@ -873,8 +873,5 @@
 	IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum);
 	IssmDouble rho_ice    = element->GetMaterialParameter(MaterialsRhoIceEnum);
-	//Input* watercolumn_input      = inputs->GetInput(WatercolumnEnum);                 _assert_(watercolumn_input);
 	Input* enthalpy_input         = element->GetInput(EnthalpyEnum);                    _assert_(enthalpy_input);
-	//  Input* pressure_input         = inputs->GetInput(PressureEnum);                    _assert_(pressure_input);
-	//	Input* basalmeltingrate_input = inputs->GetInput(BasalforcingsMeltingRateEnum);    _assert_(basalmeltingrate_input);
 	Input* geothermalflux_input   = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
 	Input* vx_input               = element->GetInput(VxEnum);                          _assert_(vx_input);
@@ -909,5 +906,5 @@
 		
 		bool checkpositivethickness=true;
-		_assert_(watercolumn>=0.);
+		_assert_(watercolumn[vertexdown]>=0.);
 
 		/*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
@@ -947,14 +944,14 @@
 			element->EnthalpyToThermal(&temperature,&waterfraction,enthalpy[vertexdown],pressure[vertexdown]);
 			geothermalflux_input->GetInputValue(&geothermalflux,gauss);
-			// -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66
+			/* -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66*/
 			heating[is]=(heatflux+basalfriction+geothermalflux);
 			meltingrate_enthalpy[is]=heating[is]/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent //????
 		}
 	}
-	// enthalpy might have been changed, update 
+	/* enthalpy might have been changed, update */
 	//element->AddInput(EnthalpyEnum,enthalpy,P1Enum);
 
 	/******** DRAINAGE *****************************************/
-	IssmDouble* drainrate_column = xNew<IssmDouble>(numsegments); //TODO: xDelete?
+	IssmDouble* drainrate_column  = xNew<IssmDouble>(numsegments); //TODO: xDelete?
 	IssmDouble* drainrate_element = xNew<IssmDouble>(numsegments);
 	for(is=0;is<numsegments;is++)	drainrate_column[is]=0.;
@@ -968,5 +965,5 @@
 		elementi=elementi->GetUpperElement();			
 	}
-	// add drained water to melting rate
+	/* add drained water to melting rate*/
 	for(is=0;is<numsegments;is++) meltingrate_enthalpy[is]+=drainrate_column[is];
 
@@ -976,5 +973,5 @@
 		vertexdown = pairindices[is*2+0];
 		vertexup   = pairindices[is*2+1];
-		if(reCast<bool,IssmDouble>(dt)){
+		if(dt!=0.){
 			if(watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt<0.){	// prevent too much freeze on			
 				melting_overshoot=watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt;
@@ -982,5 +979,5 @@
 				basalmeltingrate[vertexdown]=(1.-lambda)*meltingrate_enthalpy[is];
 				watercolumn[vertexdown]=0.;
-				yts=365*24*60*60;
+				yts=365.*24.*60.*60.;
 				enthalpy[vertexdown]+=dt/yts*lambda*heating[is];
 			}