Index: /issm/trunk-jpl/test/NightlyRun/test355.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test355.m	(revision 23976)
+++ /issm/trunk-jpl/test/NightlyRun/test355.m	(revision 23976)
@@ -0,0 +1,94 @@
+%Test Name: SquareSheetHydrologyGlaDS
+cluster=generic('name',oshostname(),'np',1);
+
+%create model: 
+md=triangle(model(),'../Exp/Square.exp',50000.);
+md.mesh.x = md.mesh.x/100;
+md.mesh.y = md.mesh.y/100;
+md.miscellaneous.name='testChannels';
+
+%miscellaneous
+md=setmask(md,'',''); %everywhere grounded
+md=setflowequation(md,'SSA','all');
+
+%Some constants
+md.constants.g=9.8;
+md.materials.rho_ice=910;
+md.materials.rho_water=1000;
+
+%Geometry
+md.geometry.surface= -.02*md.mesh.x + 320;
+md.geometry.bed = zeros(md.mesh.numberofvertices,1);
+md.geometry.base=md.geometry.bed; 
+md.geometry.thickness = md.geometry.surface-md.geometry.bed;
+
+%Define initial conditions
+md.initialization.vx = 10^-6*md.constants.yts*ones(md.mesh.numberofvertices,1);
+md.initialization.vy = zeros(md.mesh.numberofvertices,1);
+md.initialization.temperature=(273.-20.)*ones(md.mesh.numberofvertices,1);
+md.initialization.watercolumn=0.03*ones(md.mesh.numberofvertices,1);
+md.initialization.hydraulic_potential = md.materials.rho_ice*md.constants.g*md.geometry.thickness;
+
+%Materials
+md.materials.rheology_B    = paterson(md.initialization.temperature);
+md.materials.rheology_B(:) = (5e-25)^(-1/3);
+md.materials.rheology_n    = 3.*ones(md.mesh.numberofelements,1);
+
+%Friction
+md.friction.coefficient=zeros(md.mesh.numberofvertices,1);
+md.friction.p=ones(md.mesh.numberofelements,1);
+md.friction.q=ones(md.mesh.numberofelements,1);
+md.friction.coupling=0;
+
+%Boundary conditions:
+md=SetIceSheetBC(md);
+
+md.inversion.iscontrol=0;
+md.transient=deactivateall(md.transient);
+md.transient.ishydrology=1;
+
+% Set numerical conditions
+md.timestepping.time_step=.1/365;
+md.timestepping.final_time=.4/365;
+
+%Change hydrology class to Glads model
+md.hydrology=hydrologyglads();
+md.hydrology.ischannels=1;
+md.hydrology.englacial_void_ratio=1e-5;
+md.hydrology.moulin_input=zeros(md.mesh.numberofvertices,1);
+md.hydrology.neumannflux=zeros(md.mesh.numberofelements,1);
+md.hydrology.bump_height = 1e-1 * ones(md.mesh.numberofvertices,1);
+md.hydrology.sheet_conductivity= 1e-3 * ones(md.mesh.numberofvertices,1);
+
+% BCs for hydrology
+pos=find(md.mesh.x==100 & md.mesh.vertexonboundary);
+md.hydrology.spcphi=NaN(md.mesh.numberofvertices,1);
+md.hydrology.spcphi(pos) = md.materials.rho_ice * md.constants.g * md.geometry.thickness(pos);
+
+md.cluster=cluster;
+md=solve(md,'Transient'); %or 'tr'
+
+%Fields and tolerances to track changes
+field_names ={...
+	'HydrologySheetThickness1','HydraulicPotential1','EffectivePressure1',...
+	'HydrologySheetThickness2','HydraulicPotential2','EffectivePressure2',...
+	'HydrologySheetThickness3','HydraulicPotential3','EffectivePressure3',...
+	'HydrologySheetThickness4','HydraulicPotential4','EffectivePressure4',};
+field_tolerances={...
+	1e-13,1e-13,1e-13,...
+	1e-13,1e-13,1e-13,...
+	1e-13,1e-13,1e-13,...
+	1e-13,1e-13,1e-13};
+field_values={...
+	md.results.TransientSolution(1).HydrologySheetThickness, ...
+	md.results.TransientSolution(1).HydraulicPotential,...
+	md.results.TransientSolution(1).EffectivePressure,...
+	md.results.TransientSolution(2).HydrologySheetThickness, ...
+	md.results.TransientSolution(2).HydraulicPotential,...
+	md.results.TransientSolution(2).EffectivePressure,...
+	md.results.TransientSolution(3).HydrologySheetThickness, ...
+	md.results.TransientSolution(3).HydraulicPotential,...
+	md.results.TransientSolution(3).EffectivePressure,...
+	md.results.TransientSolution(4).HydrologySheetThickness, ...
+	md.results.TransientSolution(4).HydraulicPotential,...
+	md.results.TransientSolution(4).EffectivePressure};