Index: /issm/trunk-jpl/src/m/classes/initialization.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.py	(revision 27174)
+++ /issm/trunk-jpl/src/m/classes/initialization.py	(revision 27175)
@@ -14,5 +14,5 @@
     """
 
-    def __init__(self): #{{{
+    def __init__(self):  #{{{
         self.vx = np.nan
         self.vy = np.nan
@@ -37,5 +37,6 @@
         self.setdefaultparameters()
     #}}}
-    def __repr__(self): #{{{
+
+    def __repr__(self):  #{{{
         s = '   initial field values:\n'
         s += '{}\n'.format(fielddisplay(self, 'vx', 'x component of velocity [m/yr]'))
@@ -53,11 +54,13 @@
         s += '{}\n'.format(fielddisplay(self, 'hydraulic_potential', 'Hydraulic potential (for GlaDS) [Pa]'))
         s += '{}\n'.format(fielddisplay(self, 'channelarea', 'subglaciale water channel area (for GlaDS) [m2]'))
-        s += '{}\n'.format(fielddisplay(self,'sample', 'Realization of a Gaussian random field'))
+        s += '{}\n'.format(fielddisplay(self, 'sample', 'Realization of a Gaussian random field'))
         return s
     #}}}
-    def setdefaultparameters(self): #{{{
+
+    def setdefaultparameters(self):  #{{{
         return
     #}}}
-    def checkconsistency(self, md, solution, analyses): #{{{
+
+    def checkconsistency(self, md, solution, analyses):  #{{{
         if 'StressbalanceAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isstressbalance:
             if not np.any(np.logical_or(np.isnan(md.initialization.vx), np.isnan(md.initialization.vy))):
@@ -93,9 +96,13 @@
                 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'):
+            if type(md.hydrology).__name__ == 'hydrologyshreve':
                 if (solution == 'TransientSolution' and md.transient.ishydrology) or solution == 'HydrologySolution':
                     md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
+        if 'HydrologyDCAnalysis' in analyses:
+            if type(md.hydrology).__name__ == 'hydrologydc':
+                if (solution == 'TransientSolution' and md.transient.ishydrology) or solution == 'HydrologySolution':
+                    md = checkfield(md, 'fieldname', 'initialization.sediment_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
         if 'HydrologyTwsAnalysis' in analyses:
-            if hasattr(md.hydrology, 'hydrologytws'):
+            if type(md.hydrology).__name__ == 'hydrologytws':
                 md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
         if 'SealevelchangeAnalysis' in analyses:
@@ -103,5 +110,5 @@
                 md = checkfield(md, 'fieldname', 'initialization.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
         if 'HydrologyGlaDSAnalysis' in analyses:
-            if hasattr(md.hydrology, 'hydrologyglads'):
+            if type(md.hydrology).__name__ == 'hydrologyglads':
                 md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
                 md = checkfield(md, 'fieldname', 'initialization.hydraulic_potential', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
@@ -112,4 +119,5 @@
         return md
     #}}}
+
     def marshall(self, prefix, md, fid): #{{{
         yts = md.constants.yts
@@ -119,5 +127,5 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'pressure', 'format', 'DoubleMat', 'mattype', 1)
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', 'format', 'DoubleMat', 'mattype', 1)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'str', 'format', 'DoubleMat', 'mattype', 1)
@@ -143,4 +151,5 @@
             WriteData(fid, prefix, 'data', self.enthalpy, 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.initialization.enthalpy')
     #}}}
+
     def extrude(self, md): #{{{
         self.vx = project3d(md, 'vector', self.vx, 'type', 'node')