Index: issm/trunk-jpl/test/NightlyRun/GetIds.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/GetIds.py	(revision 24568)
+++ issm/trunk-jpl/test/NightlyRun/GetIds.py	(revision 24569)
@@ -44,4 +44,6 @@
         if np.array([type(i) == int for i in ids_names]).all():
             ids = ids_names
+        if np.array([type(i) == np.int64 for i in ids_names]).all():
+            ids = ids_names
         elif np.array([type(i) == str for i in ids_names]).all():
             ids = np.concatenate([IdFromString(i) for i in ids_names])
Index: issm/trunk-jpl/test/NightlyRun/test900.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test900.m	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test900.m	(revision 24569)
@@ -0,0 +1,47 @@
+%Test Name:SquareNoDynUnConfinedHydroDC
+md=triangle(model(),'../Exp/Square.exp',100000.);
+md=setmask(md,'','');
+%reduced slab (20m long)
+md.mesh.x=md.mesh.x/5.0e4;
+md.mesh.y=md.mesh.y/5.0e4;
+md=parameterize(md,'../Par/SquareNoDyn.par');
+md.cluster=generic('name',oshostname(),'np',1);
+
+md.transient.ishydrology=1;
+md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
+
+%Hydro Model Parameters
+md.hydrology.isefficientlayer=0;
+md.hydrology.sedimentlimit_flag=0;
+md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1);
+md.hydrology.rel_tol=1.0e-6;
+md.hydrology.penalty_lock=0;
+md.hydrology.max_iter=200;
+md.hydrology.transfer_flag=0;
+md.hydrology.unconfined_flag=1;
+%Sediment
+md.hydrology.sediment_porosity=0.1;
+md.hydrology.sediment_thickness=10.0;
+md.hydrology.sediment_transmitivity=(1.0e-3*md.hydrology.sediment_thickness)*ones(md.mesh.numberofvertices,1);
+%init
+md.initialization.sediment_head=-5.0*ones(md.mesh.numberofvertices,1);
+%BC
+md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1);
+pos=find(md.mesh.x==0);
+md.hydrology.spcsediment_head(pos)=0.5;
+
+md.timestepping.time_step=5/md.constants.yts; %5s steppin
+md.settings.output_frequency=2;
+md.timestepping.final_time=300/md.constants.yts; %500s run
+
+md=solve(md,'Transient');
+
+%fields to track, results can also be found in
+%Wang 2009 Fig 6b (jouranl of Hydrology)
+field_names={'SedimentWaterHead1',...
+	     'SedimentWaterHead2'};
+field_tolerances={1e-13,...
+		  1e-13};
+field_values={md.results.TransientSolution(11).SedimentHead,...
+	      md.results.TransientSolution(31).SedimentHead};
Index: issm/trunk-jpl/test/NightlyRun/test900.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test900.py	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test900.py	(revision 24569)
@@ -0,0 +1,53 @@
+#Test Name:SquareNoDynUnConfinedHydroDC
+import numpy as np
+from model import *
+from setmask import *
+from triangle import triangle
+from parameterize import parameterize
+from solve import solve
+from socket import gethostname
+from generic import generic
+
+md = triangle(model(), '../Exp/Square.exp', 100000.)
+md = setmask(md, '', '')
+#reduced square (20m long)
+md.mesh.x = md.mesh.x / 5.0e4
+md.mesh.y = md.mesh.y / 5.0e4
+md = parameterize(md, '../Par/SquareNoDyn.py')
+md.cluster = generic('name', gethostname(), 'np', 1)
+
+md.transient.ishydrology = True
+md.hydrology = hydrologydc()
+md.hydrology = md.hydrology.initialize(md)
+
+
+#Hydro Model Parameters
+md.hydrology.isefficientlayer = 0
+md.hydrology.sedimentlimit_flag = 0
+md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices))
+md.hydrology.rel_tol = 1.0e-6
+md.hydrology.penalty_lock = 0
+md.hydrology.max_iter = 200
+md.hydrology.transfer_flag = 0
+md.hydrology.unconfined_flag = 1
+#Sediment
+md.hydrology.sediment_porosity = 0.1
+md.hydrology.sediment_thickness = 10.0
+md.hydrology.sediment_transmitivity = (1.0e-3 * md.hydrology.sediment_thickness) * np.ones((md.mesh.numberofvertices))
+#init
+md.initialization.sediment_head = -5.0 * np.ones((md.mesh.numberofvertices))
+#BC
+md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices))
+md.hydrology.spcsediment_head[np.where(md.mesh.x == 0)] = 0.5
+
+md.timestepping.time_step = 5 / md.constants.yts  #5s steppin
+md.settings.output_frequency = 2
+md.timestepping.final_time = 300 / md.constants.yts  #500s run
+
+md = solve(md, 'Transient')
+
+#fields to track, results can also be found in
+#Wang 2009 Fig 6b (journal of Hydrology)
+field_names = ['SedimentWaterHead1', 'SedimentWaterHead2']
+field_tolerances = [1e-13, 1e-13]
+field_values = [md.results.TransientSolution[10].SedimentHead, md.results.TransientSolution[30].SedimentHead]
Index: issm/trunk-jpl/test/NightlyRun/test901.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test901.m	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test901.m	(revision 24569)
@@ -0,0 +1,33 @@
+%Test Name: SquareNoDynHydrologyDCOneLayer
+md=triangle(model(),'../Exp/Square.exp',100000.);
+md=setmask(md,'','');
+md=parameterize(md,'../Par/IceCube.par');
+md.cluster=generic('name',oshostname(),'np',1);
+
+md.transient.ishydrology=1;
+md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
+
+md.hydrology.isefficientlayer=0;
+md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1);
+md.hydrology.sedimentlimit_flag=1;
+md.hydrology.sedimentlimit=8000.0;
+
+md.initialization.sediment_head=0.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1);
+pos=find(md.mesh.y==0);
+md.hydrology.spcsediment_head(pos)=0.0;
+md.basalforcings.groundedice_melting_rate = 2.0*ones(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate = 0.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.sediment_transmitivity=3*ones(md.mesh.numberofvertices,1);
+md.timestepping.time_step=0;
+md.timestepping.final_time=1.0;
+md=solve(md,'Hydrology');
+
+%Fields and tolerances to track changes
+%you can also compare with an analitic solution, but it is exact
+%only if no limits are applied
+%analitic=(md.mesh.y.^2-2*md.mesh.y*1.0e6)*(-2.0/(2*md.constants.yts*md.hydrology.sediment_transmitivity))
+field_names     ={'SedimentWaterHead','SedimentHeadResidual'};
+field_tolerances={1e-13, 3e-10};
+field_values={md.results.HydrologySolution.SedimentHead,md.results.HydrologySolution.SedimentHeadResidual};
Index: issm/trunk-jpl/test/NightlyRun/test901.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test901.py	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test901.py	(revision 24569)
@@ -0,0 +1,42 @@
+#Test Name: SquareNoDynHydrologyDCOneLayer
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import triangle
+from setmask import setmask
+from parameterize import parameterize
+from solve import solve
+from generic import generic
+
+md = triangle(model(), '../Exp/Square.exp', 100000.)
+md = setmask(md, '', '')
+md = parameterize(md, '../Par/SquareNoDyn.py')
+md.cluster = generic('name', gethostname(), 'np', 1)
+
+md.transient.ishydrology = True
+md.hydrology = hydrologydc()
+md.hydrology = md.hydrology.initialize(md)
+
+md.hydrology.isefficientlayer = 0
+md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices))
+md.hydrology.sedimentlimit_flag = 1
+md.hydrology.sedimentlimit = 8000.0
+
+md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices))
+md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices))
+pos = np.nonzero(md.mesh.y == 0.)[0]
+md.hydrology.spcsediment_head[pos] = 0.0
+md.basalforcings.groundedice_melting_rate = 2.0 * np.ones((md.mesh.numberofvertices))
+md.basalforcings.floatingice_melting_rate = 0.0 * np.ones((md.mesh.numberofvertices))
+md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices))
+md.timestepping.time_step = 0
+md.timestepping.final_time = 1.0
+md = solve(md, 'Hydrology')
+
+#Fields and tolerances to track changes
+#you can also compare with an analitic solution, but it is exact
+#only if no limits are applied
+#analitic=(md.mesh.y**2 - 2 * md.mesh.y * 1.0e6) * (-2.0 / (2 * md.constants.yts * md.hydrology.sediment_transmitivity))
+field_names = ['SedimentWaterHead', 'SedimentHeadResidual']
+field_tolerances = [1e-13, 3e-10]
+field_values = [md.results.HydrologySolution.SedimentHead, md.results.HydrologySolution.SedimentHeadResidual]
Index: issm/trunk-jpl/test/NightlyRun/test902.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test902.m	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test902.m	(revision 24569)
@@ -0,0 +1,75 @@
+%Test Name: SquareSheetHydrologyDCTwoLayers
+md=triangle(model(),'../Exp/Square.exp',100000.);
+md=setmask(md,'','');
+md=parameterize(md,'../Par/IceCube.par');
+md.cluster=generic('name',oshostname(),'np',1);
+
+md.transient.ishydrology=1;
+md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
+
+md.hydrology.isefficientlayer=1;
+md.hydrology.sedimentlimit_flag=1;
+md.hydrology.sedimentlimit=800.0;
+md.hydrology.transfer_flag = 0;
+md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1);
+md.initialization.sediment_head=0.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1);
+
+md.basalforcings.groundedice_melting_rate = 2.0*ones(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate = 0.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.sediment_transmitivity=3*ones(md.mesh.numberofvertices,1);
+
+md.initialization.epl_head=0.0*ones(md.mesh.numberofvertices,1);
+md.initialization.epl_thickness=1.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.spcepl_head=NaN*ones(md.mesh.numberofvertices,1);
+md.hydrology.mask_eplactive_node=0*ones(md.mesh.numberofvertices,1);
+md.hydrology.epl_conductivity=30;
+md.hydrology.epl_initial_thickness=1;
+md.hydrology.epl_colapse_thickness=1.0e-3;
+md.hydrology.epl_thick_comp=1;
+md.hydrology.epl_max_thickness=1;
+md.hydrology.steps_per_step=10;
+md.timestepping.time_step=2.0;
+md.timestepping.final_time=2.0;
+
+md=solve(md,'Transient');
+
+%re-run with no substeps
+mdfine=md;
+mdfine.results=struct();
+mdfine.hydrology.steps_per_step=1;
+mdfine.timestepping.time_step=0.2;
+mdfine=solve(mdfine,'Transient');
+
+%Fields and tolerances to track changes
+field_names     ={'SedimentWaterHead1','EplWaterHead1','SedimentHeadResidual1',...
+                  'SedimentWaterHead4','EplWaterHead4','SedimentHeadResidual4',...
+                  'SedimentWaterHead5','EplWaterHead5','SedimentHeadResidual5',...
+                  'SedimentWaterHead9','EplWaterHead9','SedimentHeadResidual9',...
+                  'EplWaterHead10', 'EplWaterHeadSubstep10', 'SedimentWaterHead10',...
+		  'SedimentWaterHeadSubstep10'};
+field_tolerances={...
+    1e-13, 1e-13, 1e-13,...
+    1e-13, 1e-13, 1e-13,...
+    1e-13, 5e-12, 1e-11,...
+    1e-13, 5e-12, 1e-11,...
+    1e-13, 1e-13, 1e-13,...
+    1e-13};
+field_values={mdfine.results.TransientSolution(1).SedimentHead, ...
+	      mdfine.results.TransientSolution(1).EplHead,...
+	      mdfine.results.TransientSolution(1).SedimentHeadResidual,...
+	      mdfine.results.TransientSolution(4).SedimentHead,...
+	      mdfine.results.TransientSolution(4).EplHead,...
+	      mdfine.results.TransientSolution(4).SedimentHeadResidual, ...
+	      mdfine.results.TransientSolution(5).SedimentHead,...
+	      mdfine.results.TransientSolution(5).EplHead,...
+	      mdfine.results.TransientSolution(5).SedimentHeadResidual, ...
+	      mdfine.results.TransientSolution(9).SedimentHead,...
+	      mdfine.results.TransientSolution(9).EplHead,...
+	      mdfine.results.TransientSolution(9).SedimentHeadResidual,...
+              md.results.TransientSolution(1).EplHead,...
+              md.results.TransientSolution(1).EplHeadSubstep,...
+              md.results.TransientSolution(1).SedimentHead,...
+              md.results.TransientSolution(1).SedimentHeadSubstep
+              };
Index: issm/trunk-jpl/test/NightlyRun/test902.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test902.py	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test902.py	(revision 24569)
@@ -0,0 +1,81 @@
+#Test Name: SquareNoDynHydrolgyDCTwoLayers
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import triangle
+from setmask import setmask
+from parameterize import parameterize
+from solve import solve
+from generic import generic
+
+md = triangle(model(), '../Exp/Square.exp', 100000.)
+md = setmask(md, '', '')
+md = parameterize(md, '../Par/SquareNoDyn.py')
+md.cluster = generic('name', gethostname(), 'np', 1)
+
+md.transient.ishydrology = True
+md.hydrology = hydrologydc()
+md.hydrology = md.hydrology.initialize(md)
+
+md.hydrology.isefficientlayer = 1
+md.hydrology.sedimentlimit_flag = 1
+md.hydrology.sedimentlimit = 800.0
+md.hydrology.transfer_flag = 0
+md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices))
+md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices))
+md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices))
+
+md.basalforcings.groundedice_melting_rate = 2.0 * np.ones((md.mesh.numberofvertices))
+md.basalforcings.floatingice_melting_rate = 0.0 * np.ones((md.mesh.numberofvertices))
+md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices))
+
+md.initialization.epl_head = np.zeros((md.mesh.numberofvertices))
+md.initialization.epl_thickness = np.ones((md.mesh.numberofvertices))
+md.hydrology.spcepl_head = np.nan * np.ones((md.mesh.numberofvertices))
+md.hydrology.mask_eplactive_node = np.zeros((md.mesh.numberofvertices))
+md.hydrology.epl_conductivity = 30
+md.hydrology.epl_initial_thickness = 1
+md.hydrology.epl_colapse_thickness = 1.0e-3
+md.hydrology.epl_thick_comp = 1
+md.hydrology.epl_max_thickness = 1
+md.hydrology.steps_per_step = 10
+md.timestepping.time_step = 2.0
+md.timestepping.final_time = 2.0
+
+md = solve(md, 'Transient')
+
+#re-run with no substeps
+mdfine = copy.deepcopy(md)
+mdfine.results = []
+mdfine.hydrology.steps_per_step = 1
+mdfine.timestepping.time_step = 0.2
+mdfine = solve(mdfine, 'Transient')
+
+field_names = ['SedimentWaterHead1', 'EplWaterHead1', 'SedimentHeadResidual1',
+               'SedimentWaterHead4', 'EplWaterHead4', 'SedimentHeadResidual4',
+               'SedimentWaterHead5', 'EplWaterHead5', 'SedimentHeadResidual5',
+               'SedimentWaterHead9', 'EplWaterHead9', 'SedimentHeadResidual9',
+               'EplWaterHead10', 'EplWaterHeadSubstep10', 'SedimentWaterHead10',
+               'SedimentWaterHeadSubstep10']
+field_tolerances = [1e-13, 1e-13, 1e-13,
+                    1e-13, 1e-13, 1e-13,
+                    1e-13, 5e-12, 1e-11,
+                    1e-13, 5e-12, 1e-11,
+                    1e-13, 1e-13, 1e-13,
+                    1e-13]
+field_values = [mdfine.results.TransientSolution[0].SedimentHead,
+                mdfine.results.TransientSolution[0].EplHead,
+                mdfine.results.TransientSolution[0].SedimentHeadResidual,
+                mdfine.results.TransientSolution[3].SedimentHead,
+                mdfine.results.TransientSolution[3].EplHead,
+                mdfine.results.TransientSolution[3].SedimentHeadResidual,
+                mdfine.results.TransientSolution[4].SedimentHead,
+                mdfine.results.TransientSolution[4].EplHead,
+                mdfine.results.TransientSolution[4].SedimentHeadResidual,
+                mdfine.results.TransientSolution[8].SedimentHead,
+                mdfine.results.TransientSolution[8].EplHead,
+                mdfine.results.TransientSolution[8].SedimentHeadResidual,
+                md.results.TransientSolution[-1].EplHead,
+                md.results.TransientSolution[-1].EplHeadSubstep,
+                md.results.TransientSolution[-1].SedimentHead,
+                md.results.TransientSolution[-1].SedimentHeadSubstep]
Index: issm/trunk-jpl/test/NightlyRun/test903.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test903.m	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test903.m	(revision 24569)
@@ -0,0 +1,35 @@
+%Test Name: SquareNoDynExtrudedHydrologyDCOneLayer
+md=triangle(model(),'../Exp/Square.exp',100000.);
+md=setmask(md,'','');
+md=parameterize(md,'../Par/IceCube.par');
+md.cluster=generic('name',oshostname(),'np',1);
+
+md.transient.ishydrology=1;
+md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
+
+md.hydrology.isefficientlayer=0;
+md.hydrology.sedimentlimit_flag=1;
+md.hydrology.sedimentlimit=8000.0;
+md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1);
+md.initialization.sediment_head=0.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1);
+pos=find(md.mesh.y==0);
+md.hydrology.spcsediment_head(pos)=0.0;
+
+md.basalforcings.groundedice_melting_rate = 2.0*ones(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate = 0.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.sediment_transmitivity= 3.0*ones(md.mesh.numberofvertices,1);
+
+md.timestepping.time_step=0;
+md.timestepping.final_time=1.0;
+md=extrude(md,3,1.1);
+md=solve(md,'Hydrology');
+
+%Fields and tolerances to track changes
+%you can also compare with an analitic solution, but it is exact
+%only if no limits are applied
+%analitic=(md.mesh.y.^2-2*md.mesh.y*1.0e6)*(-2.0/(2*md.constants.yts*md.hydrology.sediment_transmitivity))
+field_names     ={'SedimentWaterHead','SedimentHeadResidual'};
+field_tolerances={1e-13, 3e-10};
+field_values={md.results.HydrologySolution.SedimentHead,md.results.HydrologySolution.SedimentHeadResidual};
Index: issm/trunk-jpl/test/NightlyRun/test903.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test903.py	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test903.py	(revision 24569)
@@ -0,0 +1,43 @@
+#Test Name: SquareNoDynExtrudedHydrologyDCOneLayer
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import triangle
+from setmask import setmask
+from parameterize import parameterize
+from solve import solve
+from generic import generic
+
+md = triangle(model(), '../Exp/Square.exp', 100000.)
+md = setmask(md, '', '')
+md = parameterize(md, '../Par/SquareNoDyn.py')
+md.cluster = generic('name', gethostname(), 'np', 1)
+
+md.transient.ishydrology = True
+md.hydrology = hydrologydc()
+md.hydrology = md.hydrology.initialize(md)
+
+md.hydrology.isefficientlayer = 0
+md.hydrology.sedimentlimit_flag = 1
+md.hydrology.sedimentlimit = 8000.0
+md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices))
+md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices))
+md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices))
+md.hydrology.spcsediment_head[np.where(md.mesh.y == 0)] = 0.0
+
+md.basalforcings.groundedice_melting_rate = 2.0 * np.ones((md.mesh.numberofvertices))
+md.basalforcings.floatingice_melting_rate = 0.0 * np.ones((md.mesh.numberofvertices))
+md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices))
+
+md.timestepping.time_step = 0
+md.timestepping.final_time = 1.0
+md.extrude(3, 1.)
+md = solve(md, 'Hydrology')
+
+#Fields and tolerances to track changes
+#you can also compare with an analitic solution, but it is exact
+#only if no limits are applied
+#analitic=(md.mesh.y.^2 - 2 * md.mesh.y * 1.0e6) * (-2.0 / (2 * md.constants.yts * md.hydrology.sediment_transmitivity))
+field_names = ['SedimentWaterHead', 'SedimentHeadResidual']
+field_tolerances = [1e-13, 3e-10]
+field_values = [md.results.HydrologySolution.SedimentHead, md.results.HydrologySolution.SedimentHeadResidual]
Index: issm/trunk-jpl/test/NightlyRun/test904.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test904.m	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test904.m	(revision 24569)
@@ -0,0 +1,57 @@
+%Test Name: SquareNoDynExtrudedHydrologyDCTwoLayers
+md=triangle(model(),'../Exp/Square.exp',100000.);
+md=setmask(md,'','');
+md=parameterize(md,'../Par/IceCube.par');
+md.cluster=generic('name',oshostname(),'np',1);
+
+md.transient.ishydrology=1;
+md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
+
+md.hydrology.isefficientlayer=1;
+md.hydrology.sedimentlimit_flag=1;
+md.hydrology.transfer_flag = 0;
+md.hydrology.sedimentlimit=800.0;
+md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1);
+md.initialization.sediment_head=0.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1);
+md.basalforcings.groundedice_melting_rate = 2.0*ones(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate = 0.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.sediment_transmitivity=3*ones(md.mesh.numberofvertices,1);
+
+md.initialization.epl_head=0.0*ones(md.mesh.numberofvertices,1);
+md.initialization.epl_thickness=1.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.spcepl_head=NaN*ones(md.mesh.numberofvertices,1);
+md.hydrology.mask_eplactive_node=0*ones(md.mesh.numberofvertices,1);
+md.hydrology.epl_conductivity=30;
+md.hydrology.epl_initial_thickness=1;
+md.hydrology.epl_colapse_thickness=1.0e-3;
+md.hydrology.epl_thick_comp=1;
+md.hydrology.epl_max_thickness=1;
+md.timestepping.time_step=0.2;
+md.timestepping.final_time=2.0;
+
+md=extrude(md,3,1.);
+md=solve(md,'Transient');
+
+%Fields and tolerances to track changes
+field_names     ={'SedimentWaterHead1','EplWaterHead1','SedimentHeadResidual1',...
+		  'SedimentWaterHead4','EplWaterHead4','SedimentHeadResidual4',...
+		  'SedimentWaterHead5','EplWaterHead5','SedimentHeadResidual5',...
+		  'SedimentWaterHead9','EplWaterHead9','SedimentHeadResidual9'};
+field_tolerances={1e-13, 1e-13, 1e-13,...
+		  1e-13, 1e-13, 1e-13,...
+		  1e-13, 5e-12, 2e-11,...
+		  1e-13, 5e-12, 2e-11};
+field_values={md.results.TransientSolution(1).SedimentHead, ...
+	      md.results.TransientSolution(1).EplHead,...
+	      md.results.TransientSolution(1).SedimentHeadResidual,...
+	      md.results.TransientSolution(4).SedimentHead,...
+	      md.results.TransientSolution(4).EplHead,...
+	      md.results.TransientSolution(4).SedimentHeadResidual, ...
+	      md.results.TransientSolution(5).SedimentHead,...
+	      md.results.TransientSolution(5).EplHead,...
+	      md.results.TransientSolution(5).SedimentHeadResidual, ...
+	      md.results.TransientSolution(9).SedimentHead,...
+	      md.results.TransientSolution(9).EplHead,...
+	      md.results.TransientSolution(9).SedimentHeadResidual};
Index: issm/trunk-jpl/test/NightlyRun/test904.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test904.py	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test904.py	(revision 24569)
@@ -0,0 +1,66 @@
+#Test Name: SquareNoDynExtrudedHydrologyDCTwoLayers
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import triangle
+from setmask import setmask
+from parameterize import parameterize
+from solve import solve
+from generic import generic
+
+md = triangle(model(), '../Exp/Square.exp', 100000.)
+md = setmask(md, '', '')
+md = parameterize(md, '../Par/SquareNoDyn.py')
+md.cluster = generic('name', gethostname(), 'np', 1)
+
+md.transient.ishydrology = True
+md.hydrology = hydrologydc()
+md.hydrology = md.hydrology.initialize(md)
+
+md.hydrology.isefficientlayer = 1
+md.hydrology.sedimentlimit_flag = 1
+md.hydrology.transfer_flag = 0
+md.hydrology.sedimentlimit = 800.0
+md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices))
+md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices))
+md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices))
+md.basalforcings.groundedice_melting_rate = 2.0 * np.ones((md.mesh.numberofvertices))
+md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices))
+md.hydrology.sediment_transmitivity = 3 * np.ones((md.mesh.numberofvertices))
+
+md.initialization.epl_head = np.zeros((md.mesh.numberofvertices))
+md.initialization.epl_thickness = np.ones((md.mesh.numberofvertices))
+md.hydrology.spcepl_head = np.nan * np.ones((md.mesh.numberofvertices))
+md.hydrology.mask_eplactive_node = np.zeros((md.mesh.numberofvertices))
+md.hydrology.epl_conductivity = 30
+md.hydrology.epl_initial_thickness = 1
+md.hydrology.epl_colapse_thickness = 1.0e-3
+md.hydrology.epl_thick_comp = 1
+md.hydrology.epl_max_thickness = 1
+md.timestepping.time_step = 0.2
+md.timestepping.final_time = 2.0
+
+md.extrude(3, 1.)
+md = solve(md, 'Transient')
+
+#Fields and tolerances to track changes
+field_names = ['SedimentWaterHead1', 'EplWaterHead1', 'SedimentHeadResidual1',
+               'SedimentWaterHead4', 'EplWaterHead4', 'SedimentHeadResidual4',
+               'SedimentWaterHead5', 'EplWaterHead5', 'SedimentHeadResidual5',
+               'SedimentWaterHead9', 'EplWaterHead9', 'SedimentHeadResidual9']
+field_tolerances = [1e-13, 1e-13, 1e-13,
+                    1e-13, 1e-13, 1e-13,
+                    1e-13, 5e-12, 2e-11,
+                    1e-13, 5e-12, 2e-11]
+field_values = [md.results.TransientSolution[0].SedimentHead,
+                md.results.TransientSolution[0].EplHead,
+                md.results.TransientSolution[0].SedimentHeadResidual,
+                md.results.TransientSolution[3].SedimentHead,
+                md.results.TransientSolution[3].EplHead,
+                md.results.TransientSolution[3].SedimentHeadResidual,
+                md.results.TransientSolution[4].SedimentHead,
+                md.results.TransientSolution[4].EplHead,
+                md.results.TransientSolution[4].SedimentHeadResidual,
+                md.results.TransientSolution[8].SedimentHead,
+                md.results.TransientSolution[8].EplHead,
+                md.results.TransientSolution[8].SedimentHeadResidual]
Index: issm/trunk-jpl/test/NightlyRun/test905.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test905.m	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test905.m	(revision 24569)
@@ -0,0 +1,84 @@
+%Test Name: SquareNoDynHydrologyDCSmbCoupled
+md=triangle(model(),'../Exp/Square.exp',100000.);
+md=setmask(md,'','');
+md=parameterize(md,'../Par/SquareNoDyn.par');
+md.cluster = generic('name',oshostname(),'np',1);
+
+md.transient.ishydrology=1;
+md.transient.issmb=1;
+md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
+md.smb=(SMBgradientscomponents);
+
+md.hydrology.isefficientlayer=1;
+
+md.hydrology.sedimentlimit_flag=1;
+md.hydrology.sedimentlimit=400.0;
+md.hydrology.sediment_transmitivity=3.0*ones(md.mesh.numberofvertices,1);
+md.hydrology.mask_thawed_node=ones(md.mesh.numberofvertices,1);
+
+md.hydrology.mask_eplactive_node=zeros(md.mesh.numberofvertices,1);
+md.hydrology.epl_conductivity=3.;
+md.hydrology.epl_initial_thickness=20;
+md.hydrology.epl_colapse_thickness=1.0e-3;
+md.hydrology.epl_thick_comp=0;
+md.hydrology.epl_max_thickness=1;
+
+md.hydrology.spcsediment_head=NaN*ones(md.mesh.numberofvertices,1);
+md.hydrology.spcepl_head=NaN*ones(md.mesh.numberofvertices,1);
+
+md.initialization.sediment_head=zeros(md.mesh.numberofvertices,1);
+md.initialization.epl_head=zeros(md.mesh.numberofvertices,1);
+md.initialization.epl_thickness=np.ones(md.mesh.numberofvertices,1);
+
+md.hydrology.steps_per_step=5;
+md.smb.steps_per_step=7;
+md.timestepping.time_step=2;
+md.timestepping.final_time=20.0;
+
+smb_step=md.timestepping.time_step/md.smb.steps_per_step;
+duration=[md.timestepping.start_time:smb_step:md.timestepping.final_time];
+
+ddf=10.0e-3;
+md.smb.accuref=[[0.5 0.5];[0. 2.0]];
+md.smb.accualti=0.0;
+md.smb.accugrad=0.0;
+
+md.smb.runoffref=0.9*duration*ddf;
+md.smb.runoffref=[md.smb.runoffref;duration];
+md.smb.runoffalti=0.0;
+md.smb.runoffgrad=-6.5e-3*ddf;  %lapse rate *ddf*day per year
+
+md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+
+md=solve(md,'Transient');
+
+field_names={'SedimentWaterHead1','EplWaterHead1','SedimentHeadResidual1',...
+	     'SedimentWaterHead4','EplWaterHead4','SedimentHeadResidual4',...
+	     'SedimentWaterHead5','EplWaterHead5','SedimentHeadResidual5',...
+	     'SedimentWaterHead9','EplWaterHead9','SedimentHeadResidual9',...
+	     'EplWaterHead10','EplWaterHeadSubstep10','SedimentWaterHead10',...
+	     'SedimentWaterHeadSubstep10'}
+field_tolerances={1e-13,1e-13,1e-13,...
+		  1e-13,1e-13,1e-13,...
+		  1e-13,5e-12,1e-11,...
+		  1e-13,5e-12,1e-11,...
+		  1e-13,1e-13,1e-13,...
+		  1e-13};
+field_values={md.results.TransientSolution(1).SedimentHead,...
+	      md.results.TransientSolution(1).EplHead,...
+	      md.results.TransientSolution(1).SedimentHeadResidual,...
+	      md.results.TransientSolution(4).SedimentHead,...
+	      md.results.TransientSolution(4).EplHead,...
+	      md.results.TransientSolution(4).SedimentHeadResidual,...
+	      md.results.TransientSolution(5).SedimentHead,...
+	      md.results.TransientSolution(5).EplHead,...
+	      md.results.TransientSolution(5).SedimentHeadResidual,...
+	      md.results.TransientSolution(9).SedimentHead,...
+	      md.results.TransientSolution(9).EplHead,...
+	      md.results.TransientSolution(9).SedimentHeadResidual,...
+	      md.results.TransientSolution(end).EplHead,...
+	      md.results.TransientSolution(end).EplHeadSubstep,...
+	      md.results.TransientSolution(end).SedimentHead,...
+	      md.results.TransientSolution(end).SedimentHeadSubstep};
Index: issm/trunk-jpl/test/NightlyRun/test905.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test905.py	(revision 24569)
+++ issm/trunk-jpl/test/NightlyRun/test905.py	(revision 24569)
@@ -0,0 +1,94 @@
+#Test Name: SquareNoDynHydrologyDCSmbCoupled
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import triangle
+from setmask import setmask
+from parameterize import parameterize
+from solve import solve
+from generic import generic
+from SMBgradientscomponents import SMBgradientscomponents
+
+md = triangle(model(), '../Exp/Square.exp', 100000.)
+md = setmask(md, '', '')
+md = parameterize(md, '../Par/SquareNoDyn.py')
+md.cluster = generic('name', gethostname(), 'np', 1)
+
+md.transient.ishydrology = True
+md.transient.issmb = True
+md.hydrology = hydrologydc()
+md.hydrology = md.hydrology.initialize(md)
+md.smb = SMBgradientscomponents()
+
+md.hydrology.isefficientlayer = 1
+
+md.hydrology.sedimentlimit_flag = 1
+md.hydrology.sedimentlimit = 400.0
+md.hydrology.sediment_transmitivity = 3.0 * np.ones((md.mesh.numberofvertices))
+md.hydrology.mask_thawed_node = np.ones((md.mesh.numberofvertices))
+
+md.hydrology.mask_eplactive_node = np.zeros((md.mesh.numberofvertices))
+md.hydrology.epl_conductivity = 3.
+md.hydrology.epl_initial_thickness = 20
+md.hydrology.epl_colapse_thickness = 1.0e-3
+md.hydrology.epl_thick_comp = 0
+md.hydrology.epl_max_thickness = 1
+
+md.hydrology.spcsediment_head = np.nan * np.ones((md.mesh.numberofvertices))
+md.hydrology.spcepl_head = np.nan * np.ones((md.mesh.numberofvertices))
+
+md.initialization.sediment_head = np.zeros((md.mesh.numberofvertices))
+md.initialization.epl_head = np.zeros((md.mesh.numberofvertices))
+md.initialization.epl_thickness = np.ones((md.mesh.numberofvertices))
+
+md.hydrology.steps_per_step = 5
+md.smb.steps_per_step = 7  #md.hydrology.steps_per_step
+md.timestepping.time_step = 2
+md.timestepping.final_time = 20.0
+
+smb_step = md.timestepping.time_step / md.smb.steps_per_step
+duration = np.arange(md.timestepping.start_time, md.timestepping.final_time + smb_step, smb_step)
+
+ddf = 10.0e-3
+md.smb.accuref = np.array([[0.5, 0.5], [0., 2.0]])
+md.smb.accualti = 0.0
+md.smb.accugrad = 0.0
+
+md.smb.runoffref = 0.9 * duration * ddf
+md.smb.runoffref = np.vstack((md.smb.runoffref, duration))
+md.smb.runoffalti = 0.0
+md.smb.runoffgrad = -6.5e-3 * ddf  #lapse rate *ddf*day per year
+
+md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices))
+md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices))
+
+md = solve(md, 'Transient')
+
+field_names = ['SedimentWaterHead1', 'EplWaterHead1', 'SedimentHeadResidual1',
+               'SedimentWaterHead4', 'EplWaterHead4', 'SedimentHeadResidual4',
+               'SedimentWaterHead5', 'EplWaterHead5', 'SedimentHeadResidual5',
+               'SedimentWaterHead9', 'EplWaterHead9', 'SedimentHeadResidual9',
+               'EplWaterHead10', 'EplWaterHeadSubstep10', 'SedimentWaterHead10',
+               'SedimentWaterHeadSubstep10']
+field_tolerances = [1e-13, 1e-13, 1e-13,
+                    1e-13, 1e-13, 1e-13,
+                    1e-13, 5e-12, 1e-11,
+                    1e-13, 5e-12, 1e-11,
+                    1e-13, 1e-13, 1e-13,
+                    1e-13]
+field_values = [md.results.TransientSolution[0].SedimentHead,
+                md.results.TransientSolution[0].EplHead,
+                md.results.TransientSolution[0].SedimentHeadResidual,
+                md.results.TransientSolution[3].SedimentHead,
+                md.results.TransientSolution[3].EplHead,
+                md.results.TransientSolution[3].SedimentHeadResidual,
+                md.results.TransientSolution[4].SedimentHead,
+                md.results.TransientSolution[4].EplHead,
+                md.results.TransientSolution[4].SedimentHeadResidual,
+                md.results.TransientSolution[8].SedimentHead,
+                md.results.TransientSolution[8].EplHead,
+                md.results.TransientSolution[8].SedimentHeadResidual,
+                md.results.TransientSolution[-1].EplHead,
+                md.results.TransientSolution[-1].EplHeadSubstep,
+                md.results.TransientSolution[-1].SedimentHead,
+                md.results.TransientSolution[-1].SedimentHeadSubstep]
