Index: ../trunk-jpl/src/m/classes/model.py
===================================================================
--- ../trunk-jpl/src/m/classes/model.py	(revision 18577)
+++ ../trunk-jpl/src/m/classes/model.py	(revision 18578)
@@ -29,6 +29,7 @@
 from stressbalance import stressbalance
 from groundingline import groundingline
 from hydrologyshreve import hydrologyshreve
+from hydrologydc import hydrologydc
 from masstransport import masstransport
 from thermal import thermal
 from steadystate import steadystate
@@ -79,10 +80,10 @@
 		self.cluster          = generic()
 
 		self.balancethickness = balancethickness()
-		self.stressbalance       = stressbalance()
+		self.stressbalance    = stressbalance()
 		self.groundingline    = groundingline()
 		self.hydrology        = hydrologyshreve()
-		self.masstransport       = masstransport()
+		self.masstransport    = masstransport()
 		self.thermal          = thermal()
 		self.steadystate      = steadystate()
 		self.transient        = transient()
@@ -127,7 +128,7 @@
 		        'thermal',\
 		        'steadystate',\
 		        'transient',\
-				  'gia',\
+                        'gia',\
 		        'autodiff',\
 		        'flaim',\
 		        'inversion',\
@@ -614,8 +615,13 @@
 			md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node')
 		if not numpy.any(numpy.isnan(md.initialization.watercolumn)):
 			md.initialization.watercolumn=project3d(md,'vector',md.initialization.watercolumn,'type','node')
+                if not numpy.any(numpy.isnan(md.initialization.sediment_head)):
+                        md.initialization.sediment_head=project3d(md,'vector',md.initialization.sediment_head,'type','node','layer',1)
+                if not numpy.any(numpy.isnan(md.initialization.epl_head)):
+                        md.initialization.epl_head=project3d(md,'vector',md.initialization.epl_head,'type','node','layer',1)
+                if not numpy.any(numpy.isnan(md.initialization.epl_thickness)):
+                        md.initialization.epl_thickness=project3d(md,'vector',md.initialization.epl_thickness,'type','node','layer',1)
 
-
 		#bedinfo and surface info
 		md.mesh.vertexonbase=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
 		md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
@@ -653,6 +659,14 @@
 		md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node')
 		md.masstransport.calvingrate=project3d(md,'vector',md.masstransport.calvingrate,'type','node')
 
+                # Hydrologydc variables
+                if hasattr(md.hydrology,'hydrologydc'):
+                        md.hydrology.spcsediment_head=project3d(md,'vector',md.hydrology.spcsediment_head,'type','node','layer',1)
+                        md.hydrology.spcepl_head=project3d(md,'vector',md.hydrology.spcepl_head,'type','node','layer',1)
+                        md.hydrology.mask_eplactive_node=project3d(md,'vector',md.hydrology.mask_eplactive_node,'type','node','layer',1)
+                        md.hydrology.sediment_transmitivity=project3d(md,'vector',md.hydrology.sediment_transmitivity,'type','node','layer',1)
+                        md.hydrology.basal_moulin_input=project3d(md,'vector',md.hydrology.basal_moulin_input,'type','node','layer',1)
+
 		#connectivity
 		md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1))
 		md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1
Index: ../trunk-jpl/src/m/classes/model.m
===================================================================
--- ../trunk-jpl/src/m/classes/model.m	(revision 18577)
+++ ../trunk-jpl/src/m/classes/model.m	(revision 18578)
@@ -210,9 +210,11 @@
 			if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end;
 			if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end;
 			if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end;
-			if ~isnan(md.initialization.temperature),	md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end;
+			if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end;
 			if ~isnan(md.initialization.pressure),md.initialization.pressure=project2d(md,md.initialization.pressure,1);end;
-
+			if ~isnan(md.initialization.sediment_head),md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);end;
+			if ~isnan(md.initialization.epl_head),md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);end;
+			if ~isnan(md.initialization.epl_thickness),md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);end;
 			%gia
 			if ~isnan(md.gia.mantle_viscosity), md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1); end
 			if ~isnan(md.gia.lithosphere_thickness), md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1); end
