Index: /issm/trunk-jpl/src/m/classes/frictionpism.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/frictionpism.py	(revision 25023)
+++ /issm/trunk-jpl/src/m/classes/frictionpism.py	(revision 25023)
@@ -0,0 +1,77 @@
+import numpy as np
+from fielddisplay import fielddisplay
+from project3d import project3d
+from checkfield import checkfield
+from WriteData import WriteData
+
+
+class frictionpism(object):
+    """
+    frictionpism class definition
+
+    Usage:
+      frictionpism = frictionpism()
+    """
+    def __init__(self):
+        self.pseudoplasticity_exponent = 0.
+        self.threshold_speed = 0.
+        self.delta = 0.
+        self.void_ratio = 0.
+        self.till_friction_angle = np.nan
+        self.sediment_compressibility_coefficient = np.nan
+
+        self.setdefaultparameters()
+        self.requested_outputs = []
+    #}}}
+
+    def extrude(self, md):  # {{{
+        self.till_friction_angle = project3d(md, 'vector', self.till_friction_angle, 'type', 'node', 'layer', 1)
+        self.sediment_compressibility_coefficient = project3d(md, 'vector', self.sediment_compressibility_coefficient, 'type', 'node', 'layer', 1)
+        return self
+        # }}}
+
+    def setdefaultparameters(self):  # {{{
+        self.pseudoplasticity_exponent = 0.6
+        self.threshold_speed = 100.
+        self.delta = 0.02
+        self.void_ratio = 0.69
+        return self
+    # }}}
+
+    def checkconsistency(self, md, solution, analyses):  # {{{
+
+        #Early return
+        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
+            return md
+        if solution == 'TransientSolution' and md.transient.isstressbalance == 0 and md.transient.isthermal == 0:
+            return md
+
+        md = checkfield(md, 'fieldname', 'friction.pseudoplasticity_exponent', 'numel', [1], '>', 0, 'NaN', 1, 'Inf', 1)
+        md = checkfield(md, 'fieldname', 'friction.threshold_speed', 'numel', [1], '>', 0, 'NaN', 1, 'Inf', 1)
+        md = checkfield(md, 'fieldname', 'friction.delta', 'numel', [1], '>', 0, '<', 1, 'NaN', 1, 'Inf', 1)
+        md = checkfield(md, 'fieldname', 'friction.void_ratio', 'numel', [1], '>', 0, '<', 1, 'NaN', 1, 'Inf', 1)
+        md = checkfield(md, 'fieldname', 'friction.till_friction_angle', 'NaN', 1, 'Inf', 1, '<', 360., '>', 0., 'size', [md.mesh.numberofvertices])  #User should give angle in degrees, Matlab calculates in rad
+        md = checkfield(md, 'fieldname', 'friction.sediment_compressibility_coefficient', 'NaN', 1, 'Inf', 1, '<', 1., '>', 0., 'size', [md.mesh.numberofvertices])
+    # }}}
+
+    def __repr__(self):  # {{{
+        string = 'Basal shear stress parameters for the PISM friction law (See Aschwanden et al. 2016 for more details)'
+        string = "%s\n%s" % (string, fielddisplay(self, 'pseudoplasticity_exponent', 'pseudoplasticity exponent [dimensionless]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'threshold_speed', 'threshold speed [m / yr]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'delta', 'lower limit of the effective pressure, expressed as a fraction of overburden pressure [dimensionless]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'void_ratio', 'void ratio at a reference effective pressure [dimensionless]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'till_friction_angle', 'till friction angle [deg], recommended default: 30 deg'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'sediment_compressibility_coefficient', 'coefficient of compressibility of the sediment [dimensionless], recommended default: 0.12'))
+        # }}}
+
+    def marshall(self, prefix, md, fid):  # {{{
+        yts = md.constants.yts
+
+        WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 10, 'format', 'Integer')
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'pseudoplasticity_exponent', 'format', 'Double')
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'threshold_speed', 'format', 'Double', 'scale', 1. / yts)
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'delta', 'format', 'Double')
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'void_ratio', 'format', 'Double')
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'till_friction_angle', 'format', 'DoubleMat', 'mattype', 1)
+        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'sediment_compressibility_coefficient', 'format', 'DoubleMat', 'mattype', 1)
+        # }}}
Index: /issm/trunk-jpl/src/m/classes/hydrologyglads.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/hydrologyglads.py	(revision 25023)
+++ /issm/trunk-jpl/src/m/classes/hydrologyglads.py	(revision 25023)
@@ -0,0 +1,141 @@
+import numpy as np
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
+
+
+class hydrologyglads(object):
+    """
+    HydrologyGlaDSclass definition
+
+    Usage:
+      hydrologyglads = hydrologyglads()
+    """
+
+    def __init__(self):  # {{{
+        #Sheet
+        self.pressure_melt_coefficient = 0.
+        self.sheet_conductivity = np.nan
+        self.cavity_spacing = 0.
+        self.bump_height = np.nan
+
+        #Channels
+        self.ischannels = 0
+        self.channel_conductivity = 0.
+        self.channel_sheet_width = 0.
+
+        #Other
+        self.spcphi = np.nan
+        self.moulin_input = np.nan
+        self.neumannflux = np.nan
+        self.englacial_void_ratio = 0.
+        self.requested_outputs = []
+        self.melt_flag = 0
+
+        # set defaults
+        self.setdefaultparameters()
+        # }}}
+
+    def __repr__(self):  # {{{
+        string = '   GlaDS (hydrologyglads) solution parameters:'
+        string = "%s\n\t%s" % (string, '--SHEET')
+        string = "%s\n%s" % (string, fielddisplay(self, 'pressure_melt_coefficient', 'Pressure melt coefficient (c_t) [K Pa^ - 1]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'sheet_conductivity', 'sheet conductivity (k) [m^(7 / 4) kg^(- 1 / 2)]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'cavity_spacing', 'cavity spacing (l_r) [m]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'bump_height', 'typical bump height (h_r) [m]'))
+        string = "%s\n\t%s" % (string, '--CHANNELS')
+        string = "%s\n%s" % (string, fielddisplay(self, 'ischannels', 'Do we allow for channels? 1: yes, 0: no'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'channel_conductivity', 'channel conductivity (k_c) [m^(3 / 2) kg^(- 1 / 2)]'))
+        string = "%s\n\t%s" % (string, '--OTHER')
+        string = "%s\n%s" % (string, fielddisplay(self, 'spcphi', 'Hydraulic potential Dirichlet constraints [Pa]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'neumannflux', 'water flux applied along the model boundary (m^2 / s)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'moulin_input', 'moulin input (Q_s) [m^3 / s]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'englacial_void_ratio', 'englacial void ratio (e_v)'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'melt_flag', 'User specified basal melt? 0: no (default), 1: use md.basalforcings.groundedice_melting_rate'))
+        return string
+    # }}}
+
+    def defaultoutputs(self, md):  # {{{
+        list = ['EffectivePressure', 'HydraulicPotential', 'HydrologySheetThickness', 'ChannelArea']
+        return list
+    # }}}
+
+    def setdefaultparameters(self):  # {{{
+
+        #Sheet parameters
+        self.pressure_melt_coefficient = 7.5e-8  #K / Pa (See table 1 in Erder et al. 2013)
+        self.cavity_spacing = 2.  #m
+
+        #Channel parameters
+        self.ischannels = False
+        self.channel_conductivity = 5.e-2  #Dow's default, Table uses 0.1
+        self.channel_sheet_width = 2.  #m
+
+        #Other
+        self.englacial_void_ratio = 1.e-5  #Dow's default, Table from Werder et al. uses 1e-3
+        self.requested_outputs = ['default']
+        self.melt_flag = 0
+
+        return self
+    # }}}
+
+    def checkconsistency(self, md, solution, analyses):  # {{{
+
+        #Early return
+        if 'HydrologyGladsAnalysis' not in analyses:
+            return md
+
+        #Sheet
+        md = checkfield(md, 'fieldname', 'hydrology.pressure_melt_coefficient', 'numel', [1], '>=', 0)
+        md = checkfield(md, 'fieldname', 'hydrology.sheet_conductivity', 'size', [md.mesh.numberofvertices], '>', 0, 'np.nan', 1, 'Inf', 1)
+        md = checkfield(md, 'fieldname', 'hydrology.cavity_spacing', 'numel', [1], '>', 0)
+        md = checkfield(md, 'fieldname', 'hydrology.bump_height', 'size', [md.mesh.numberofvertices], '>=', 0, 'np.nan', 1, 'Inf', 1)
+
+        #Channels
+        md = checkfield(md, 'fieldname', 'hydrology.ischannels', 'numel', [1], 'values', [0, 1])
+        md = checkfield(md, 'fieldname', 'hydrology.channel_conductivity', 'numel', [1], '>', 0)
+        md = checkfield(md, 'fieldname', 'hydrology.channel_sheet_width', 'numel', [1], '>=', 0)
+
+        #Other
+        md = checkfield(md, 'fieldname', 'hydrology.spcphi', 'Inf', 1, 'timeseries', 1)
+        md = checkfield(md, 'fieldname', 'hydrology.englacial_void_ratio', 'numel', [1], '>=', 0)
+        md = checkfield(md, 'fieldname', 'hydrology.moulin_input', 'size', [md.mesh.numberofvertices], '>=', 0, 'NaN', 1, 'Inf', 1)
+        md = checkfield(md, 'fieldname', 'hydrology.neumannflux', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
+        md = checkfield(md, 'fieldname', 'hydrology.requested_outputs', 'stringrow', 1)
+        md = checkfield(md, 'fieldname', 'hydrology.melt_flag', 'numel', [1], 'values', [0, 1])
+        if self.melt_flag:
+            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
+    # }}}
+
+    def marshall(self, prefix, md, fid):  # {{{
+        yts = md.constants.yts
+        #Marshall model code first
+        WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 5, 'format', 'Integer')
+
+        #Sheet
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'pressure_melt_coefficient', 'format', 'Double')
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'sheet_conductivity', 'format', 'DoubleMat', 'mattype', 1)
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'cavity_spacing', 'format', 'Double')
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'bump_height', 'format', 'DoubleMat', 'mattype', 1)
+
+        #Channels
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'ischannels', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'channel_conductivity', 'format', 'Double')
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'channel_sheet_width', 'format', 'Double')
+
+        #Others
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'spcphi', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'neumannflux', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'moulin_input', 'format', 'DoubleMat', 'mattype', 1)
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'englacial_void_ratio', 'format', 'Double')
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'melt_flag', 'format', 'Boolean')
+
+        outputs = self.requested_outputs
+        indices = [i for i, x in enumerate(outputs) if x == 'default']
+        if len(indices) > 0:
+            outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
+            outputs = outputscopy
+        WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray')
+
+    # }}}
Index: /issm/trunk-jpl/src/m/classes/hydrologypism.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/hydrologypism.py	(revision 25023)
+++ /issm/trunk-jpl/src/m/classes/hydrologypism.py	(revision 25023)
@@ -0,0 +1,57 @@
+import numpy as np
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
+
+
+class hydrologypism(object):
+    """
+    HydrologyPism class definition
+
+    Usage:
+      hydrologypism = hydrologypism()
+    """
+
+    def __init__(self):  # {{{
+        self.drainage_rate = np.nan
+        self.watercolumn_max = np.nan
+
+        # set defaults
+        self.setdefaultparameters()
+    # }}}
+
+    def extrude(self, md):  # {{{
+        return self
+    # }}}
+
+    def defaultoutputs(self, md):  # {{{
+        list = ['Watercolumn']
+        return list
+
+    def setdefaultparameters(self):  # {{{
+        return self
+
+    def checkconsistency(self, md, solution, analyses):  # {{{
+
+        #Early return
+        if 'HydrologyPismAnalysis' not in analyses:
+            return
+
+        md = checkfield(md, 'fieldname', 'hydrology.drainage_rate', 'Inf', 1, 'NaN', 1, '>=', 0, 'size', [md.mesh.numberofvertices])
+        md = checkfield(md, 'fieldname', 'friction.watercolumn_max', 'NaN', 1, 'Inf', 1, '>', 0., 'size', [md.mesh.numberofvertices])
+
+    def __repr__(self):  # {{{
+        string = '   hydrologypism solution parameters:'
+        string = "%s\n%s" % (string, fielddisplay(self, 'drainage_rate', 'fixed drainage rate [mm / yr]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'watercolumn_max', 'maximum water column height [m], recommended default: 2 m'))
+        return string
+        # }}}
+
+    def marshall(self, prefix, md, fid):  # {{{
+        yts = md.constants.yts
+
+        WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 4, 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'drainage_rate', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / (1000. * yts))  #from mm / yr to m / s
+        WriteData(fid, prefix, 'class', 'hydrology', 'object', self, 'fieldname', 'watercolumn_max', 'format', 'DoubleMat', 'mattype', 1)
+        WriteData(fid, prefix, 'data', {'Watercolumn'}, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray')
+    # }}}
Index: /issm/trunk-jpl/src/m/classes/initialization.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/initialization.py	(revision 25022)
+++ /issm/trunk-jpl/src/m/classes/initialization.py	(revision 25023)
@@ -28,4 +28,6 @@
         self.epl_head = float('NaN')
         self.epl_thickness = float('NaN')
+        self.hydraulic_potential = float('NaN')
+        self.channelarea = float('NaN')
 
     #set defaults
@@ -48,4 +50,6 @@
         string = "%s\n%s" % (string, fielddisplay(self, 'epl_head', 'epl water head of subglacial system [m]'))
         string = "%s\n%s" % (string, fielddisplay(self, 'epl_thickness', 'thickness of the epl [m]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'hydraulic_potential', 'Hydraulic potential (for GlaDS) [Pa]'))
+        string = "%s\n%s" % (string, fielddisplay(self, 'channelarea', 'subglaciale water channel area (for GlaDS) [m2]'))
 
         return string
@@ -97,5 +101,5 @@
                                      np.sum(np.abs(md.initialization.vy[md.mesh.elements - 1]), axis=1) == 0)):
                 md.checkmessage("at least one triangle has all its vertices with a zero velocity")
-        if 'ThermalAnalysis' in analyses:
+        if 'ThermalAnalysis' in analyses and solution != 'TransientSolution' and md.transient.isthermal == 0:
             md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
             md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
@@ -121,5 +125,9 @@
                     md = checkfield(md, 'fieldname', 'initialization.epl_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
                     md = checkfield(md, 'fieldname', 'initialization.epl_thickness', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
-
+        if 'HydrologyGlaDSAnalysis' in analyses:
+            if hasattr(md.hydrology, '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])
+                md = checkfield(md, 'fieldname', 'initialization.channelarea', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [md.mesh.numberofelements])
         return md
     # }}}
@@ -139,5 +147,6 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'epl_thickness', 'format', 'DoubleMat', 'mattype', 1)
         WriteData(fid, prefix, 'object', self, 'fieldname', 'watercolumn', 'format', 'DoubleMat', 'mattype', 1)
-
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'channelarea', 'format', 'DoubleMat', 'mattype', 1)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'hydraulic_potential', 'format', 'DoubleMat', 'mattype', 1)
         if md.thermal.isenthalpy:
             if (np.size(self.enthalpy) <= 1):
Index: /issm/trunk-jpl/src/m/classes/model.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.py	(revision 25022)
+++ /issm/trunk-jpl/src/m/classes/model.py	(revision 25023)
@@ -46,4 +46,7 @@
 from hydrologyshreve import hydrologyshreve
 from hydrologydc import hydrologydc
+from hydrologyglads import hydrologyglads
+from hydrologypism import hydrologypism
+from hydrologyshakti import hydrologyshakti
 from masstransport import masstransport
 from thermal import thermal
Index: /issm/trunk-jpl/src/m/classes/stressbalance.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/stressbalance.py	(revision 25022)
+++ /issm/trunk-jpl/src/m/classes/stressbalance.py	(revision 25023)
@@ -88,5 +88,5 @@
         self.maxiter = 100
         #Convergence criterion: absolute, relative and residual
-        self.restol = 10**-4
+        self.restol = 1.0e-4
         self.reltol = 0.01
         self.abstol = 10
Index: /issm/trunk-jpl/test/NightlyRun/test350.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test350.py	(revision 25022)
+++ /issm/trunk-jpl/test/NightlyRun/test350.py	(revision 25023)
@@ -4,12 +4,11 @@
 from model import *
 from socket import gethostname
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-from frictionshakti import *
-from hydrologyshakti import *
-from transient import *
+from triangle import triangle
+from setmask import setmask
+from parameterize import parameterize
+from setflowequation import setflowequation
+from solve import solve
+from frictionshakti import frictionshakti
+from transient import transient
 
 md = triangle(model(), '../Exp/Square.exp', 50000.)
@@ -18,5 +17,5 @@
 md = setmask(md, '', '')
 md = parameterize(md, '../Par/SquareSheetConstrained.py')
-md.transient = transient().deactivateall()
+md.transient = transient.deactivateall(md.transient)
 md.transient.ishydrology = 1
 md = setflowequation(md, 'SSA', 'all')
Index: /issm/trunk-jpl/test/NightlyRun/test355.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test355.py	(revision 25023)
+++ /issm/trunk-jpl/test/NightlyRun/test355.py	(revision 25023)
@@ -0,0 +1,97 @@
+#Test Name: SquareSheetHydrologyGlaDS
+import numpy as np
+from model import *
+from triangle import triangle
+from setmask import setmask
+from setflowequation import setflowequation
+from parameterize import parameterize
+from solve import solve
+from SetIceSheetBC import SetIceSheetBC
+from generic import generic
+
+#create model:
+md = triangle(model(), '../Exp/Square.exp', 50000.)
+md.mesh.x = md.mesh.x / 100
+md.mesh.y = md.mesh.y / 100
+#miscellaneous
+md = setmask(md, '', '')   #everywhere grounded
+md = setflowequation(md, 'SSA', 'all')
+md.miscellaneous.name = 'testChannels'
+md.stressbalance.maxiter = 2  #Make sure it runs quickly...
+
+#Some constants
+md.constants.g = 9.8
+md.materials.rho_ice = 910
+
+#Geometry
+md.geometry.surface = -0.02 * md.mesh.x + 320
+md.geometry.bed = np.zeros((md.mesh.numberofvertices))
+md.geometry.base = md.geometry.bed
+md.geometry.thickness = md.geometry.surface - md.geometry.bed
+
+#Define initial conditions
+md.initialization.vx = 1.0e-6 * md.constants.yts * np.ones((md.mesh.numberofvertices))
+md.initialization.vy = np.zeros((md.mesh.numberofvertices))
+md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices))
+md.initialization.watercolumn = 0.03 * np.ones((md.mesh.numberofvertices))
+md.initialization.hydraulic_potential = md.materials.rho_ice * md.constants.g * md.geometry.thickness
+
+#Materials
+md.materials.rheology_B = (5e-25)**(-1 / 3) * np.ones((md.mesh.numberofvertices))
+md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements))
+
+#Friction
+md.friction.coefficient = np.zeros((md.mesh.numberofvertices))
+md.friction.p = np.ones((md.mesh.numberofelements))
+md.friction.q = np.ones((md.mesh.numberofelements))
+#md.friction.coupling = 0
+
+#Boundary conditions:
+md = SetIceSheetBC(md)
+
+md.inversion.iscontrol = 0
+md.transient = transient.deactivateall(md.transient)
+md.transient.ishydrology = 1
+
+# Set numerical conditions
+md.timestepping.time_step = 0.1 / 365
+md.timestepping.final_time = 0.4 / 365
+
+#Change hydrology class to Glads model
+md.hydrology = hydrologyglads()
+md.hydrology.ischannels = 1
+md.hydrology.englacial_void_ratio = 1.e-5
+md.hydrology.moulin_input = np.zeros((md.mesh.numberofvertices))
+md.hydrology.neumannflux = np.zeros((md.mesh.numberofelements))
+md.hydrology.bump_height = 1.e-1 * np.ones((md.mesh.numberofvertices))
+md.hydrology.sheet_conductivity = 1.e-3 * np.ones((md.mesh.numberofvertices))
+
+# BCs for hydrology
+pos = np.where(np.logical_and(md.mesh.x == 100, md.mesh.vertexonboundary))
+md.hydrology.spcphi = np.nan * np.ones((md.mesh.numberofvertices))
+md.hydrology.spcphi[pos] = md.materials.rho_ice * md.constants.g * md.geometry.thickness[pos]
+
+md.cluster = generic('np', 2)
+md = solve(md, 'Transient')   #or 'tr'
+
+#Fields and tolerances to track changes
+field_names = ['HydrologySheetThickness1', 'HydraulicPotential1', 'ChannelArea1',
+               'HydrologySheetThickness2', 'HydraulicPotential2', 'ChannelArea2',
+               'HydrologySheetThickness3', 'HydraulicPotential3', 'ChannelArea3',
+               'HydrologySheetThickness4', 'HydraulicPotential4', 'ChannelArea4']
+field_tolerances = [1e-13, 1e-13, 1e-13,
+                    1e-13, 1e-13, 2e-13,
+                    1e-13, 1e-13, 7e-13,
+                    1e-13, 1e-13, 2e-12]
+field_values = [md.results.TransientSolution[0].HydrologySheetThickness,
+                md.results.TransientSolution[0].HydraulicPotential,
+                md.results.TransientSolution[0].ChannelArea,
+                md.results.TransientSolution[1].HydrologySheetThickness,
+                md.results.TransientSolution[1].HydraulicPotential,
+                md.results.TransientSolution[1].ChannelArea,
+                md.results.TransientSolution[2].HydrologySheetThickness,
+                md.results.TransientSolution[2].HydraulicPotential,
+                md.results.TransientSolution[2].ChannelArea,
+                md.results.TransientSolution[3].HydrologySheetThickness,
+                md.results.TransientSolution[3].HydraulicPotential,
+                md.results.TransientSolution[3].ChannelArea]
Index: /issm/trunk-jpl/test/NightlyRun/test610.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test610.py	(revision 25023)
+++ /issm/trunk-jpl/test/NightlyRun/test610.py	(revision 25023)
@@ -0,0 +1,43 @@
+#Test Name: 79NorthPISMhydro
+import numpy as np
+from model import *
+from setmask import setmask
+from triangle import triangle
+from parameterize import parameterize
+from setflowequation import setflowequation
+from solve import solve
+from socket import gethostname
+from generic import generic
+
+
+md = triangle(model(), '../Exp/79North.exp', 10000.)
+md = setmask(md, '../Exp/79NorthShelf.exp', '')
+md = parameterize(md, '../Par/79North.py')
+md = setflowequation(md, 'SSA', 'all')
+
+#Hydrology
+md.hydrology = hydrologypism()
+md.hydrology.drainage_rate = 0.001 * np.ones((md.mesh.numberofvertices))
+md.hydrology.watercolumn_max = 20 * np.ones((md.mesh.numberofvertices))
+md.initialization.pressure = np.zeros((md.mesh.numberofvertices))
+md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices))
+md.basalforcings.groundedice_melting_rate = np.arange(0, md.mesh.numberofvertices) + 1
+md.transient.ishydrology = 1
+md.transient.issmb = 0
+md.transient.ismasstransport = 0
+md.transient.isstressbalance = 0
+md.transient.isthermal = 0
+md.transient.isgroundingline = 0
+
+md.cluster = generic('name', gethostname(), 'np', 3)
+md = solve(md, 'Transient')
+
+#Plot to check result
+# plotmodel(md, 'data', md.results.TransientSolution(3).Watercolumn - 3 * (md.materials.rho_freshwater / md.materials.rho_ice * [1:md.mesh.numberofvertices]' - 1))
+
+#Fields and tolerances to track changes
+field_names = ['WaterColumn1', 'WaterColumn2', 'WaterColumn3']
+field_tolerances = [1e-12, 1e-12, 1e-12]
+field_values = [md.results.TransientSolution[0].Watercolumn,
+                md.results.TransientSolution[1].Watercolumn,
+                md.results.TransientSolution[2].Watercolumn]
Index: /issm/trunk-jpl/test/NightlyRun/test612.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test612.py	(revision 25023)
+++ /issm/trunk-jpl/test/NightlyRun/test612.py	(revision 25023)
@@ -0,0 +1,56 @@
+#Test Name: 79NorthPISMhydrofriction
+import numpy as np
+from model import *
+from setmask import *
+from triangle import triangle
+from parameterize import parameterize
+from solve import solve
+from frictionpism import frictionpism
+from socket import gethostname
+from generic import generic
+
+md = triangle(model(), '../Exp/79North.exp', 10000.)
+md = setmask(md, '../Exp/79NorthShelf.exp', '')
+md = parameterize(md, '../Par/79North.py')
+md = setflowequation(md, 'SSA', 'all')
+
+#Hydrology
+md.hydrology = hydrologypism()
+md.hydrology.drainage_rate = 0.001 * np.ones((md.mesh.numberofvertices))
+md.hydrology.watercolumn_max = 2 * np.ones((md.mesh.numberofvertices))
+md.initialization.pressure = np.zeros((md.mesh.numberofvertices))
+md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices))
+md.basalforcings.groundedice_melting_rate = np.arange(0, md.mesh.numberofvertices) + 1
+
+#Friction
+md.friction = frictionpism()
+md.friction.till_friction_angle = 30 * np.ones((md.mesh.numberofvertices))
+md.friction.sediment_compressibility_coefficient = 0.12 * np.ones((md.mesh.numberofvertices))
+
+md.transient.ishydrology = 1
+md.transient.issmb = 0
+md.transient.ismasstransport = 0
+md.transient.isstressbalance = 1
+md.transient.isthermal = 0
+md.transient.isgroundingline = 0
+
+md.cluster = generic('name', gethostname(), 'np', 3)
+md = solve(md, 'Transient')
+
+#Plot to check result
+# plotmodel(md, 'data', md.results.TransientSolution(3).Vx - md.results.TransientSolution(2).Vx, 'caxis', [ - 1 1])
+# plotmodel(md, 'data', md.results.TransientSolution(3).Vy - md.results.TransientSolution(2).Vy, 'caxis', [ - 1 1])
+# plotmodel(md, 'data', md.results.TransientSolution(3).Vel - md.results.TransientSolution(2).Vel)
+
+#Fields and tolerances to track changes
+field_names = ['Vx1', 'Vx2', 'Vx3', 'Vy1', 'Vy2', 'Vy3', 'Vel1', 'Vel2', 'Vel3']
+field_tolerances = [2e-11, 2e-11, 2e-11, 2e-11, 2e-11, 2e-11, 2e-11, 2e-11, 2e-11]
+field_values = [md.results.TransientSolution[0].Vx,
+                md.results.TransientSolution[1].Vx,
+                md.results.TransientSolution[2].Vx,
+                md.results.TransientSolution[0].Vy,
+                md.results.TransientSolution[1].Vy,
+                md.results.TransientSolution[2].Vy,
+                md.results.TransientSolution[0].Vel,
+                md.results.TransientSolution[1].Vel,
+                md.results.TransientSolution[2].Vel]