Changeset 23739

Timestamp:

02/19/19 11:06:40 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: adding enthalpy to temperature for Josh that has non P1 elements

File:

• issm/trunk-jpl/src/m/classes/initialization.m (modified) (6 diffs)

issm/trunk-jpl/src/m/classes/initialization.m

@@ -12 +12 @@
                 pressure      = NaN;
                 temperature   = NaN;
+                enthalpy      = NaN;
                 waterfraction = NaN;
                 sediment_head = NaN;
@@ -25 +26 @@
                         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 +112 @@
                         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 +129 @@
                         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 +137 @@
 
                         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 +157 @@
                         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);