Index: /issm/trunk-jpl/src/m/classes/initialization.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.m	(revision 27504)
+++ /issm/trunk-jpl/src/m/classes/initialization.m	(revision 27505)
@@ -26,4 +26,5 @@
 		sample              = NaN;
 		debris              = NaN;
+		age                 = NaN;
 	end
 	methods
@@ -127,4 +128,9 @@
 				if ~isnan(md.initialization.debris)
 					md = checkfield(md,'fieldname','initialization.debris','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
+				end
+			end
+			if ismember('AgeAnalysis',analyses),
+				if ~isnan(md.initialization.age)
+					md = checkfield(md,'fieldname','initialization.age','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
 				end
 			end
@@ -152,4 +158,5 @@
 			fielddisplay(self,'str','Steric sea level.');
 			fielddisplay(self,'debris','Surface debris layer [m]');
+			fielddisplay(self,'age','Initial age [yr]');
 		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
@@ -174,5 +181,5 @@
 			WriteData(fid,prefix,'object',self,'fieldname','hydraulic_potential','format','DoubleMat','mattype',1);
 			WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1);
-			WriteData(fid,prefix,'object',self,'fieldname','debris','format','DoubleMat','mattype',1);
+			WriteData(fid,prefix,'object',self,'fieldname','debris','format','DoubleMat','mattype',1,'scale',yts);
 
 			if md.thermal.isenthalpy,
@@ -206,4 +213,5 @@
 			self.str=project3d(md,'vector',self.str,'type','node','layer',1);
 			self.str=project3d(md,'vector',self.debris,'type','node','layer',1);
+			self.str=project3d(md,'vector',self.age,'type','node','layer',1);
 
 			%Lithostatic pressure by default
@@ -228,4 +236,5 @@
 			writejs1Darray(fid,[modelname '.initialization.sample'],self.sample);
 			writejs1Darray(fid,[modelname '.initialization.debris'],self.debris);
+			writejs1Darray(fid,[modelname '.initialization.age'],self.age);
 
 		end % }}}
Index: /issm/trunk-jpl/src/m/classes/initialization.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.py	(revision 27504)
+++ /issm/trunk-jpl/src/m/classes/initialization.py	(revision 27505)
@@ -15,24 +15,25 @@
 
     def __init__(self):  #{{{
-        self.vx = np.nan
-        self.vy = np.nan
-        self.vz = np.nan
-        self.vel = np.nan
-        self.pressure = np.nan
-        self.temperature = np.nan
-        self.enthalpy = np.nan
-        self.waterfraction = np.nan
-        self.sediment_head = np.nan
-        self.epl_head = np.nan
-        self.epl_thickness = np.nan
-        self.watercolumn = np.nan
+        self.vx                  = np.nan
+        self.vy                  = np.nan
+        self.vz                  = np.nan
+        self.vel                 = np.nan
+        self.pressure            = np.nan
+        self.temperature         = np.nan
+        self.enthalpy            = np.nan
+        self.waterfraction       = np.nan
+        self.sediment_head       = np.nan
+        self.epl_head            = np.nan
+        self.epl_thickness       = np.nan
+        self.watercolumn         = np.nan
         self.hydraulic_potential = np.nan
-        self.channelarea = np.nan
-        self.sealevel = np.nan
-        self.bottompressure = np.nan
-        self.dsl = np.nan
-        self.str = np.nan
-        self.sample = np.nan
-        self.debris = np.nan
+        self.channelarea         = np.nan
+        self.sealevel            = np.nan
+        self.bottompressure      = np.nan
+        self.dsl                 = np.nan
+        self.str                 = np.nan
+        self.sample              = np.nan
+        self.debris              = np.nan
+        self.age                 = np.nan
 
         self.setdefaultparameters()
@@ -57,4 +58,5 @@
         s += '{}\n'.format(fielddisplay(self, 'sample', 'Realization of a Gaussian random field'))
         s += '{}\n'.format(fielddisplay(self, 'debris', 'Surface debris layer [m]'))
+        s += '{}\n'.format(fielddisplay(self, 'age', 'initial ice age [yr]'))
         return s
     # }}}
@@ -127,4 +129,8 @@
                 if (solution == 'TransientSolution' and md.transient.ishydrology) or solution == 'HydrologySolution':
                     md = checkfield(md, 'fieldname', 'initialization.debris', 'NaN', 1,'Inf', 1, 'size', [md.mesh.numberofvertices])
+        if 'AgeAnalysis' in analyses:
+            if not np.isnan(md.initialization.age):
+                if (solution == 'TransientSolution' and md.transient.ishydrology) or solution == 'HydrologySolution':
+                    md = checkfield(md, 'fieldname', 'initialization.age', 'NaN', 1,'Inf', 1, 'size', [md.mesh.numberofvertices])
         return md
     # }}}
@@ -151,4 +157,5 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'sample', 'format', 'DoubleMat', 'mattype', 1)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'debris', 'format', 'DoubleMat', 'mattype', 1)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'age', 'format', 'DoubleMat', 'mattype', 1, 'scale', yts)
 
         if md.thermal.isenthalpy:
@@ -180,4 +187,5 @@
         self.str = project3d(md, 'vector', self.str, 'type', 'node', 'layer', 1)
         self.debris = project3d(md, 'vector', self.debris, 'type', 'node', 'layer', 1)
+        self.age = project3d(md, 'vector', self.age, 'type', 'node', 'layer', 1)
 
         # Lithostatic pressure by default