Index: /issm/trunk-jpl/src/m/classes/initialization.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.m	(revision 23738)
+++ /issm/trunk-jpl/src/m/classes/initialization.m	(revision 23739)
@@ -12,4 +12,5 @@
 		pressure      = NaN;
 		temperature   = NaN;
+		enthalpy      = NaN;
 		waterfraction = NaN;
 		sediment_head = NaN;
@@ -25,4 +26,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','layer',1);
@@ -110,4 +112,5 @@
 			fielddisplay(self,'pressure','pressure field [Pa]');
 			fielddisplay(self,'temperature','temperature [K]');
+			fielddisplay(self,'enthalpy','enthalpy [J]');
 			fielddisplay(self,'waterfraction','fraction of water in the ice');
 			fielddisplay(self,'sediment_head','sediment water head of subglacial system [m]');
@@ -126,4 +129,5 @@
 			WriteData(fid,prefix,'object',self,'fieldname','pressure','format','DoubleMat','mattype',1);
 			WriteData(fid,prefix,'object',self,'fieldname','temperature','format','DoubleMat','mattype',1);
+			WriteData(fid,prefix,'object',self,'fieldname','enthalpy','format','DoubleMat','mattype',1);
 			WriteData(fid,prefix,'object',self,'fieldname','waterfraction','format','DoubleMat','mattype',1);
 			WriteData(fid,prefix,'object',self,'fieldname','sediment_head','format','DoubleMat','mattype',1);
@@ -133,8 +137,13 @@
 
 			if md.thermal.isenthalpy,
-				tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
-				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);
+				if numel(self.enthalpy) <= 1,
+					%reconstruct enthalpy
+					tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
+					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);
+				else
+					enthalpy = self.enthalpy;
+				end
 				WriteData(fid,prefix,'data',enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy');
 			end
@@ -148,4 +157,5 @@
 			writejs1Darray(fid,[modelname '.initialization.pressure'],self.pressure);
 			writejs1Darray(fid,[modelname '.initialization.temperature'],self.temperature);
+			writejs1Darray(fid,[modelname '.initialization.enthalpy'],self.enthalpy);
 			writejs1Darray(fid,[modelname '.initialization.waterfraction'],self.waterfraction);
 			writejs1Darray(fid,[modelname '.initialization.sediment_head'],self.sediment_head);