Index: /issm/trunk-jpl/src/m/classes/initialization.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.m	(revision 21405)
+++ /issm/trunk-jpl/src/m/classes/initialization.m	(revision 21406)
@@ -76,4 +76,9 @@
 				md = checkfield(md,'fieldname','initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]);
 				md = checkfield(md,'fieldname','initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices 1]);
+				pos=find(md.initialization.waterfraction>0.);
+				if(~isempty(pos)),
+					md = checkfield(md,'fieldname', 'delta Tpmp', 'field',abs(md.initialization.temperature(pos)-(md.materials.meltingpoint-md.materials.beta*md.initialization.pressure(pos))),'<',1e-11,...
+					'message','set temperature to pressure melting point at locations with waterfraction>0');
+				end
 			end
 			if ismember('HydrologyShreveAnalysis',analyses),
@@ -129,5 +134,5 @@
 			if md.thermal.isenthalpy,
 				tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
-				pos  = find(md.initialization.temperature>=tpmp-1e-12);
+				pos  = find(md.initialization.waterfraction>0.);
 				enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
 				enthalpy(pos) = md.materials.heatcapacity*(tpmp(pos) - md.constants.referencetemperature) + md.materials.latentheat*md.initialization.waterfraction(pos);
Index: /issm/trunk-jpl/src/m/classes/initialization.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.py	(revision 21405)
+++ /issm/trunk-jpl/src/m/classes/initialization.py	(revision 21406)
@@ -100,4 +100,7 @@
 				md = checkfield(md,'fieldname','initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices])
 				md = checkfield(md,'fieldname','initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices])
+				pos = np.nonzero(md.initialization.waterfraction > 0.)[0]
+				if(pos.size):
+					md = checkfield(md,'fieldname', 'delta Tpmp', 'field', np.absolute(md.initialization.temperature[pos]-(md.materials.meltingpoint-md.materials.beta*md.initialization.pressure[pos])),'<',1e-11,	'message','set temperature to pressure melting point at locations with waterfraction>0');
 		if 'HydrologyShreveAnalysis' in analyses:
 			if hasattr(md.hydrology,'hydrologyshreve'):
@@ -131,7 +134,7 @@
 		if md.thermal.isenthalpy:
 			tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
-			pos  = np.nonzero(md.initialization.temperature > tpmp)[0]
+			pos  = np.nonzero(md.initialization.waterfraction > 0.)[0]
 			enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
-			enthalpy[pos] = md.materials.heatcapacity*tpmp[pos].reshape(-1,) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,)
+			enthalpy[pos] = md.materials.heatcapacity*(tpmp[pos].reshape(-1,) - md.constants.referencetemperature) + md.materials.latentheat*md.initialization.waterfraction[pos].reshape(-1,)
 			WriteData(fid,prefix,'data',enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy');