Index: /issm/trunk-jpl/src/m/classes/initialization.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.py	(revision 23811)
+++ /issm/trunk-jpl/src/m/classes/initialization.py	(revision 23812)
@@ -20,4 +20,5 @@
 		self.vz            = float('NaN')
 		self.vel           = float('NaN')
+                self.enthalpy      = float('NaN')
 		self.pressure      = float('NaN')
 		self.temperature   = float('NaN')
@@ -40,4 +41,5 @@
 		string="%s\n%s"%(string,fielddisplay(self,'pressure','pressure [Pa]'))
 		string="%s\n%s"%(string,fielddisplay(self,'temperature','temperature [K]'))
+                string="%s\n%s"%(string,fielddisplay(self,'enthalpy','enthalpy [J]'))
 		string="%s\n%s"%(string,fielddisplay(self,'waterfraction','fraction of water in the ice'))
 		string="%s\n%s"%(string,fielddisplay(self,'watercolumn','thickness of subglacial water [m]'))
@@ -54,4 +56,5 @@
 		self.vel=project3d(md,'vector',self.vel,'type','node')
 		self.temperature=project3d(md,'vector',self.temperature,'type','node')
+                self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node')
 		self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node')
 		self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node')
@@ -133,9 +136,11 @@
 		
 		if md.thermal.isenthalpy:
+                    if (np.size(self.enthalpy)<=1):
 			tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
 			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,)
-			WriteData(fid,prefix,'data',enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy');
+			self.enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
+			self.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',self.enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy');
 
 	# }}}