Index: /issm/trunk-jpl/test/NightlyRun/test102.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test102.js	(revision 21174)
+++ /issm/trunk-jpl/test/NightlyRun/test102.js	(revision 21175)
@@ -7,5 +7,5 @@
 setflowequation(md,'SSA','all');
 //md.cluster=generic('name',oshostname(),'np',3);
-md=solve(md,'Stressbalance');
+md=solve(md,'Stressbalance', 'checkconsistency', 'no');
 
 //Fields and tolerances to track changes
Index: /issm/trunk-jpl/test/NightlyRun/test103.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test103.js	(revision 21174)
+++ /issm/trunk-jpl/test/NightlyRun/test103.js	(revision 21175)
@@ -8,5 +8,5 @@
 //md.cluster=generic('name',oshostname(),'np',3);
 md.stressbalance.requested_outputs=['default','StressTensorxx','StressTensoryy','StressTensorzz','StressTensorxy','StressTensorxz','StressTensoryz'];
-md=solve(md,'Stressbalance');
+md=solve(md,'Stressbalance', 'checkconsistency', 'no');
 
 //Fields and tolerances to track changes
Index: /issm/trunk-jpl/test/NightlyRun/test110.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test110.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test110.js	(revision 21175)
@@ -0,0 +1,42 @@
+//Test Name: SquareShelfConstrainedTranSSA2d
+var md = new model();
+triangle(md,square[0],150000.);
+setmask(md,'all','');
+parameterize(md);
+setflowequation(md,'SSA','all');
+//md.cluster=generic('name',oshostname(),'np',3);
+md.trans.requested_outputs=['IceVolume'];
+
+md=solve(md,'transient');
+
+//Fields and tolerances to track changes
+field_names     =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Volume1','Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Volume2','Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Volume3'];
+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,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[0](1).Vx),
+	(md.results.TransientSolution[0](1).Vy),
+	(md.results.TransientSolution[0](1).Vel),
+	(md.results.TransientSolution[0](1).Pressure),
+	(md.results.TransientSolution[0](1).Base),
+	(md.results.TransientSolution[0](1).Surface),
+	(md.results.TransientSolution[0](1).Thickness),
+	(md.results.TransientSolution[0](1).IceVolume),
+	(md.results.TransientSolution[0](2).Vx),
+	(md.results.TransientSolution[0](2).Vy),
+	(md.results.TransientSolution[0](2).Vel),
+	(md.results.TransientSolution[0](2).Pressure),
+	(md.results.TransientSolution[0](2).Base),
+	(md.results.TransientSolution[0](2).Surface),
+	(md.results.TransientSolution[0](2).Thickness),
+	(md.results.TransientSolution[0](2).IceVolume),
+	(md.results.TransientSolution[0](3).Vx),
+	(md.results.TransientSolution[0](3).Vy),
+	(md.results.TransientSolution[0](3).Vel),
+	(md.results.TransientSolution[0](3).Pressure),
+	(md.results.TransientSolution[0](3).Base),
+	(md.results.TransientSolution[0](3).Surface),
+	(md.results.TransientSolution[0](3).Thickness),
+	(md.results.TransientSolution[0](3).IceVolume),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test116.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test116.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test116.js	(revision 21175)
@@ -0,0 +1,26 @@
+//Test Name: SquareShelfConstrainedBalThic2d
+var md = new model();
+triangle(md,square[0],150000.);
+setmask(md,'all','');
+parameterize(md);
+//Add boundary conditions on thickness on the border
+pos=[];
+for (var i = 0; i < md.mesh.vertexonboundary.length; ++i) {
+    if ((md.mesh.vertexonboundary[i] !== 0)) {
+            pos.push(i);
+    };
+}
+
+for (var i = 0; i < pos.length; ++i) {
+    md.balancethickness.spcthickness[pos[i]] = md.geometry.thickness[pos[i]];
+}
+setflowequation(md,'SSA','all');
+//md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Balancethickness');
+
+//Fields and tolerances to track changes
+field_names     =['Thickness'];
+field_tolerances=[1e-13];
+field_values=[
+	(md.results.BalancethicknessSolution[0].Thickness),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test120.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test120.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test120.js	(revision 21175)
@@ -0,0 +1,31 @@
+//Test Name: SquareShelfConstrainedEnthalpyStea
+var md = new model();
+function zeros(...args) {
+	var array = [];
+	for (var i = 0; i < args[0]; ++i) {
+		array.push(args.length == 1 ? 0 : zeros(args.slice(1)));
+	}
+	return array;
+}
+triangle(md,square[0],180000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,1.);
+setflowequation(md,'SSA','all');
+md.timestepping.time_step=0;
+md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1);
+md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1);
+md.thermal.isenthalpy = 1;
+md.thermal.isdynamicbasalspc = 1;
+
+//md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Thermal');
+
+//Fields and tolerances to track changes
+field_names     =['Enthalpy','Waterfraction','Temperature'];
+field_tolerances=[1e-13,3e-10,1e-13];
+field_values=[
+	(md.results.ThermalSolution[0].Enthalpy),
+	(md.results.ThermalSolution[0].Waterfraction),
+	(md.results.ThermalSolution[0].Temperature),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test121.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test121.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test121.js	(revision 21175)
@@ -0,0 +1,42 @@
+//Test Name: SquareShelfConstrainedEnthalpyTran
+function zeros(...args) {
+	var array = [];
+	for (var i = 0; i < args[0]; ++i) {
+		array.push(args.length == 1 ? 0 : zeros(args.slice(1)));
+	}
+	return array;
+}
+var md = new model();
+triangle(md,square[0],180000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,1.);
+setflowequation(md,'SSA','all');
+//md.cluster=generic('name',oshostname(),'np',3);
+md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1);
+md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1);
+md.transient.isstressbalance=0;
+md.transient.ismasstransport=0;
+md.transient.issmb=1;
+md.transient.isthermal=1;
+md.transient.isgroundingline=0;
+md.thermal.isenthalpy=1;
+md.thermal.isdynamicbasalspc=1;
+md=solve(md,'Transient');
+
+//Fields and tolerances to track changes
+field_names     =['Enthalpy1','Waterfraction1','Temperature1',
+	'Enthalpy2','Waterfraction2','Temperature2',
+	'Enthalpy3','Waterfraction3','Temperature3'];
+field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-9,1e-13];
+field_values=[
+	(md.results.TransientSolution[0](1).Enthalpy),
+	(md.results.TransientSolution[0](1).Waterfraction),
+	(md.results.TransientSolution[0](1).Temperature),
+	(md.results.TransientSolution[0](2).Enthalpy),
+	(md.results.TransientSolution[0](2).Waterfraction),
+	(md.results.TransientSolution[0](2).Temperature),
+	(md.results.TransientSolution[0](3).Enthalpy),
+	(md.results.TransientSolution[0](3).Waterfraction),
+	(md.results.TransientSolution[0](3).Temperature),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test122.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test122.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test122.js	(revision 21175)
@@ -0,0 +1,65 @@
+//Test Name: SquareShelfConstrainedTransHOEnth
+function zeros(...args) {
+	var array = [];
+	for (var i = 0; i < args[0]; ++i) {
+		array.push(args.length == 1 ? 0 : zeros(args.slice(1)));
+	}
+	return array;
+}
+var md = new model();
+triangle(md,square[0],200000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,1.);
+setflowequation(md,'HO','all');
+md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1);
+md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1);
+md.thermal.isenthalpy=1;
+md.thermal.isdynamicbasalspc=1;
+md.thermal.stabilization=2;
+//md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Transient');
+
+//Fields and tolerances to track changes
+field_names     =['Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','Enthalpy1','Waterfraction1', 
+	'Vx2','Vy2','Vz2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Temperature2','Enthalpy2','Waterfraction2', 
+	'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Temperature3','Enthalpy3','Waterfraction3'];
+field_tolerances=[
+	1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,
+	1e-09,1e-09,1e-08,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,
+	1e-09,1e-09,1e-08,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09];
+field_values=[
+	(md.results.TransientSolution[0](1).Vx),
+	(md.results.TransientSolution[0](1).Vy),
+	(md.results.TransientSolution[0](1).Vz),
+	(md.results.TransientSolution[0](1).Vel),
+	(md.results.TransientSolution[0](1).Pressure),
+	(md.results.TransientSolution[0](1).Base),
+	(md.results.TransientSolution[0](1).Surface),
+	(md.results.TransientSolution[0](1).Thickness),
+	(md.results.TransientSolution[0](1).Temperature),
+	(md.results.TransientSolution[0](1).Enthalpy),
+	(md.results.TransientSolution[0](1).Waterfraction),
+	(md.results.TransientSolution[0](2).Vx),
+	(md.results.TransientSolution[0](2).Vy),
+	(md.results.TransientSolution[0](2).Vz),
+	(md.results.TransientSolution[0](2).Vel),
+	(md.results.TransientSolution[0](2).Pressure),
+	(md.results.TransientSolution[0](2).Base),
+	(md.results.TransientSolution[0](2).Surface),
+	(md.results.TransientSolution[0](2).Thickness),
+	(md.results.TransientSolution[0](2).Temperature),
+	(md.results.TransientSolution[0](2).Enthalpy),
+	(md.results.TransientSolution[0](2).Waterfraction),
+	(md.results.TransientSolution[0](3).Vx),
+	(md.results.TransientSolution[0](3).Vy),
+	(md.results.TransientSolution[0](3).Vz),
+	(md.results.TransientSolution[0](3).Vel),
+	(md.results.TransientSolution[0](3).Pressure),
+	(md.results.TransientSolution[0](3).Base),
+	(md.results.TransientSolution[0](3).Surface),
+	(md.results.TransientSolution[0](3).Thickness),
+	(md.results.TransientSolution[0](3).Temperature),
+	(md.results.TransientSolution[0](3).Enthalpy),
+	(md.results.TransientSolution[0](3).Waterfraction),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test123.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test123.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test123.js	(revision 21175)
@@ -0,0 +1,23 @@
+//Test Name: SquareShelfConstrainedTranMisfitSurface
+var md = new model();
+triangle(md,square[0],180000.);
+setmask(md,'all','');
+parameterize(md);
+setflowequation(md,'SSA','all');
+//md.cluster=generic('name',oshostname(),'np',3);
+
+fake_surface=[[md.geometry.surface+100,1.1],[md.geometry.surface+200,2.1],[md.geometry.surface+300,2.5]];
+
+md.transient.requested_outputs=['default','SurfaceMisfit'];
+md.outputdefinition.definitions=[misfit('name','SurfaceMisfit', 'definitionenum',Outputdefinition1Enum,'model_enum',SurfaceEnum,'observation_enum',SurfaceObservationEnum,'observation',fake_surface,'timeinterpolation','nearestneighbor','weights',ones(md.mesh.numberofvertices,1),'weights_enum',WeightsSurfaceObservationEnum)];
+
+md=solve(md,'Transient');
+
+//Fields and tolerances to track changes
+field_names     =['SurfaceMisfitFirstStep','SurfaceMisfitSecondStep','SurfaceMisfitThirdStep'];
+field_tolerances=[1e-13,1e-13,1e-13];
+field_values=[
+	(md.results.TransientSolution[0](1).SurfaceMisfit),
+	(md.results.TransientSolution[0](2).SurfaceMisfit),
+	(md.results.TransientSolution[0](3).SurfaceMisfit)
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test124.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test124.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test124.js	(revision 21175)
@@ -0,0 +1,49 @@
+//Test Name: SquareShelfConstrainedTranFSFreeSurface
+var md = new model();
+triangle(md,square[0],150000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,1.);
+setflowequation(md,'FS','all');
+
+//Free surface
+md.masstransport.isfreesurface=1;
+md.timestepping.time_step=0.00001;
+md.timestepping.final_time=0.00005;
+
+//Go solve
+//md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'TransientSolution');
+
+//Fields and tolerances to track changes
+field_names     =[
+	'Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1',
+	'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2',
+	'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3'];
+field_tolerances=[
+	2e-09,3e-9,3e-9,3e-9,1e-13,1e-12,1e-12,
+	2e-09,3e-9,3e-9,3e-9,1e-10,1e-10,1e-10,
+	2e-09,3e-9,3e-9,3e-9,1e-10,1e-10,1e-10];
+field_values=[
+	(md.results.TransientSolution[0](1).Vx),
+	(md.results.TransientSolution[0](1).Vy),
+	(md.results.TransientSolution[0](1).Vel),
+	(md.results.TransientSolution[0](1).Pressure),
+	(md.results.TransientSolution[0](1).Base),
+	(md.results.TransientSolution[0](1).Surface),
+	(md.results.TransientSolution[0](1).Thickness),
+	(md.results.TransientSolution[0](2).Vx),
+	(md.results.TransientSolution[0](2).Vy),
+	(md.results.TransientSolution[0](2).Vel),
+	(md.results.TransientSolution[0](2).Pressure),
+	(md.results.TransientSolution[0](2).Base),
+	(md.results.TransientSolution[0](2).Surface),
+	(md.results.TransientSolution[0](2).Thickness),
+	(md.results.TransientSolution[0](3).Vx),
+	(md.results.TransientSolution[0](3).Vy),
+	(md.results.TransientSolution[0](3).Vel),
+	(md.results.TransientSolution[0](3).Pressure),
+	(md.results.TransientSolution[0](3).Base),
+	(md.results.TransientSolution[0](3).Surface),
+	(md.results.TransientSolution[0](3).Thickness),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test125.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test125.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test125.js	(revision 21175)
@@ -0,0 +1,49 @@
+//Test Name: SquareShelfConstrainedRestartTranSSA2d
+var md = new model();
+triangle(md,square[0],150000.);
+setmask(md,'all','');
+parameterize(md);
+setflowequation(md,'SSA','all');
+//md.cluster=generic('name',oshostname(),'np',1);
+md.transient.requested_outputs=['IceVolume','TotalSmb'];
+
+//md.verbose=verbose('solution',true);
+md.settings.recording_frequency=5;
+
+// time steps and resolution
+md.timestepping.final_time=8;
+
+md=solve(md,'TransientSolution');
+md2=solve(md,'TransientSolution','restart',1);
+
+//Fields and tolerances to track changes
+field_names     =['Vx1','Vy1','Vel1','TotalSmb1','Bed1','Surface1','Thickness1','Volume1','Vx2','Vy2','Vel2','TotalSmb2','Bed2','Surface2','Thickness2','Volume2','Vx3','Vy3','Vel3','TotalSmb3','Bed3','Surface3','Thickness3','Volume3'];
+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,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(6).Vx)-(md2.results.TransientSolution[0](1).Vx),
+	(md.results.TransientSolution(6).Vy)-(md2.results.TransientSolution[0](1).Vy),
+	(md.results.TransientSolution(6).Vel)-(md2.results.TransientSolution[0](1).Vel),
+	(md.results.TransientSolution(6).TotalSmb)-(md2.results.TransientSolution[0](1).TotalSmb),
+	(md.results.TransientSolution(6).Base)-(md2.results.TransientSolution[0](1).Base),
+	(md.results.TransientSolution(6).Surface)-(md2.results.TransientSolution[0](1).Surface),
+	(md.results.TransientSolution(6).Thickness)-(md2.results.TransientSolution[0](1).Thickness),
+	(md.results.TransientSolution(6).IceVolume)-(md2.results.TransientSolution[0](1).IceVolume),
+	(md.results.TransientSolution(7).Vx)-(md2.results.TransientSolution[0](2).Vx),
+	(md.results.TransientSolution(7).Vy)-(md2.results.TransientSolution[0](2).Vy),
+	(md.results.TransientSolution(7).Vel)-(md2.results.TransientSolution[0](2).Vel),
+	(md.results.TransientSolution(7).TotalSmb)-(md2.results.TransientSolution[0](2).TotalSmb),
+	(md.results.TransientSolution(7).Base)-(md2.results.TransientSolution[0](2).Base),
+	(md.results.TransientSolution(7).Surface)-(md2.results.TransientSolution[0](2).Surface),
+	(md.results.TransientSolution(7).Thickness)-(md2.results.TransientSolution[0](2).Thickness),
+	(md.results.TransientSolution(7).IceVolume)-(md2.results.TransientSolution[0](2).IceVolume),
+	(md.results.TransientSolution(8).Vx)-(md2.results.TransientSolution[0](3).Vx),
+	(md.results.TransientSolution(8).Vy)-(md2.results.TransientSolution[0](3).Vy),
+	(md.results.TransientSolution(8).Vel)-(md2.results.TransientSolution[0](3).Vel),
+	(md.results.TransientSolution(8).TotalSmb)-(md2.results.TransientSolution[0](3).TotalSmb),
+	(md.results.TransientSolution(8).Base)-(md2.results.TransientSolution[0](3).Base),
+	(md.results.TransientSolution(8).Surface)-(md2.results.TransientSolution[0](3).Surface),
+	(md.results.TransientSolution(8).Thickness)-(md2.results.TransientSolution[0](3).Thickness),
+	(md.results.TransientSolution(8).IceVolume)-(md2.results.TransientSolution[0](3).IceVolume),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test126.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test126.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test126.js	(revision 21175)
@@ -0,0 +1,57 @@
+//Test Name: SquareShelfConstrainedRestartTranHO3d
+var md = new model();
+triangle(md,square[0],150000.);
+setmask(md,'all','');
+parameterize(md);
+//md.cluster=generic('name',oshostname(),'np',1);
+md.extrude(md,3,1.);
+setflowequation(md,'HO','all');
+md.transient.requested_outputs=['IceVolume','TotalSmb'];
+
+//md.verbose=verbose('solution',true);
+md.settings.recording_frequency=5;
+
+// time steps and resolution
+md.timestepping.final_time=8;
+
+md=solve(md,'TransientSolution');
+md2=solve(md,'TransientSolution','restart',1);
+
+//Fields and tolerances to track changes
+field_names     =['Vx1','Vy1','Vel1','TotalSmb1','Bed1','Surface1','Thickness1','Volume1','Temperature1','Pressure1','Vx2','Vy2','Vel2','TotalSmb2','Bed2','Surface2','Thickness2','Volume2','Temperature2','Pressure2','Vx3','Vy3','Vel3','TotalSmb3','Bed3','Surface3','Thickness3','Volume3','Temperature3','Pressure3'];
+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,1e-13,1e-13,1e-13,1e-13,
+						1e-13,1e-13,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(6).Vx)-(md2.results.TransientSolution[0](1).Vx),
+	(md.results.TransientSolution(6).Vy)-(md2.results.TransientSolution[0](1).Vy),
+	(md.results.TransientSolution(6).Vel)-(md2.results.TransientSolution[0](1).Vel),
+	(md.results.TransientSolution(6).TotalSmb)-(md2.results.TransientSolution[0](1).TotalSmb),
+	(md.results.TransientSolution(6).Base)-(md2.results.TransientSolution[0](1).Base),
+	(md.results.TransientSolution(6).Surface)-(md2.results.TransientSolution[0](1).Surface),
+	(md.results.TransientSolution(6).Thickness)-(md2.results.TransientSolution[0](1).Thickness),
+	(md.results.TransientSolution(6).IceVolume)-(md2.results.TransientSolution[0](1).IceVolume),
+	(md.results.TransientSolution(6).Temperature)-(md2.results.TransientSolution[0](1).Temperature),
+	(md.results.TransientSolution(6).Pressure)-(md2.results.TransientSolution[0](1).Pressure),
+	(md.results.TransientSolution(7).Vx)-(md2.results.TransientSolution[0](2).Vx),
+	(md.results.TransientSolution(7).Vy)-(md2.results.TransientSolution[0](2).Vy),
+	(md.results.TransientSolution(7).Vel)-(md2.results.TransientSolution[0](2).Vel),
+	(md.results.TransientSolution(7).TotalSmb)-(md2.results.TransientSolution[0](2).TotalSmb),
+	(md.results.TransientSolution(7).Base)-(md2.results.TransientSolution[0](2).Base),
+	(md.results.TransientSolution(7).Surface)-(md2.results.TransientSolution[0](2).Surface),
+	(md.results.TransientSolution(7).Thickness)-(md2.results.TransientSolution[0](2).Thickness),
+	(md.results.TransientSolution(7).IceVolume)-(md2.results.TransientSolution[0](2).IceVolume),
+	(md.results.TransientSolution(7).Temperature)-(md2.results.TransientSolution[0](2).Temperature),
+	(md.results.TransientSolution(7).Pressure)-(md2.results.TransientSolution[0](2).Pressure),
+	(md.results.TransientSolution(8).Vx)-(md2.results.TransientSolution[0](3).Vx),
+	(md.results.TransientSolution(8).Vy)-(md2.results.TransientSolution[0](3).Vy),
+	(md.results.TransientSolution(8).Vel)-(md2.results.TransientSolution[0](3).Vel),
+	(md.results.TransientSolution(8).TotalSmb)-(md2.results.TransientSolution[0](3).TotalSmb),
+	(md.results.TransientSolution(8).Base)-(md2.results.TransientSolution[0](3).Base),
+	(md.results.TransientSolution(8).Surface)-(md2.results.TransientSolution[0](3).Surface),
+	(md.results.TransientSolution(8).Thickness)-(md2.results.TransientSolution[0](3).Thickness),
+	(md.results.TransientSolution(8).IceVolume)-(md2.results.TransientSolution[0](3).IceVolume),
+	(md.results.TransientSolution(8).Temperature)-(md2.results.TransientSolution[0](3).Temperature),
+	(md.results.TransientSolution(8).Pressure)-(md2.results.TransientSolution[0](3).Pressure),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test210.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test210.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test210.js	(revision 21175)
@@ -0,0 +1,51 @@
+//Test Name: SquareShelfTranHO
+var md = new model();
+triangle(md,square[0],200000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,1.);
+setflowequation(md,'HO','all');
+//md.cluster=generic('name',oshostname(),'np',3);
+md.transient.requested_outputs=['default','GroundedArea'];
+md=solve(md,'Transient');
+
+//Fields and tolerances to track changes
+field_names     =['Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','BasalforcingsGroundediceMeltingRate1', 
+				      'Vx2','Vy2','Vz2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Temperature2','BasalforcingsGroundediceMeltingRate2', 
+					   'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Temperature3','BasalforcingsGroundediceMeltingRate3'];
+field_tolerances=[
+	1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,
+	1e-09,1e-09,1e-08,1e-09,1e-09,1e-09,1e-09,1e-09,1e-09,1e-06,
+	1e-09,1e-09,1e-08,1e-08,1e-09,1e-09,1e-09,1e-09,1e-09,1e-05];
+field_values=[
+	(md.results.TransientSolution[0](1).Vx),
+	(md.results.TransientSolution[0](1).Vy),
+	(md.results.TransientSolution[0](1).Vz),
+	(md.results.TransientSolution[0](1).Vel),
+	(md.results.TransientSolution[0](1).Pressure),
+	(md.results.TransientSolution[0](1).Base),
+	(md.results.TransientSolution[0](1).Surface),
+	(md.results.TransientSolution[0](1).Thickness),
+	(md.results.TransientSolution[0](1).Temperature),
+	(md.results.TransientSolution[0](1).BasalforcingsGroundediceMeltingRate),
+	(md.results.TransientSolution[0](2).Vx),
+	(md.results.TransientSolution[0](2).Vy),
+	(md.results.TransientSolution[0](2).Vz),
+	(md.results.TransientSolution[0](2).Vel),
+	(md.results.TransientSolution[0](2).Pressure),
+	(md.results.TransientSolution[0](2).Base),
+	(md.results.TransientSolution[0](2).Surface),
+	(md.results.TransientSolution[0](2).Thickness),
+	(md.results.TransientSolution[0](2).Temperature),
+	(md.results.TransientSolution[0](2).BasalforcingsGroundediceMeltingRate),
+	(md.results.TransientSolution[0](3).Vx),
+	(md.results.TransientSolution[0](3).Vy),
+	(md.results.TransientSolution[0](3).Vz),
+	(md.results.TransientSolution[0](3).Vel),
+	(md.results.TransientSolution[0](3).Pressure),
+	(md.results.TransientSolution[0](3).Base),
+	(md.results.TransientSolution[0](3).Surface),
+	(md.results.TransientSolution[0](3).Thickness),
+	(md.results.TransientSolution[0](3).Temperature),
+	(md.results.TransientSolution[0](3).BasalforcingsGroundediceMeltingRate),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test211.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test211.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test211.js	(revision 21175)
@@ -0,0 +1,52 @@
+//Test Name: SquareShelfTranFS
+var md = new model();
+triangle(md,square[0],200000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,1.);
+setflowequation(md,'FS','all');
+md.stressbalance.reltol=null;
+//md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Transient');
+
+//Fields and tolerances to track changes
+field_names=[
+	'Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','BasalforcingsGroundediceMeltingRate1', 
+	'Vx2','Vy2','Vz2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Temperature2','BasalforcingsGroundediceMeltingRate2', 
+	'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Temperature3','BasalforcingsGroundediceMeltingRate3'];
+field_tolerances=[
+	2e-08,2e-08,5e-05,2e-08,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08,
+	5e-06,5e-06,8e-05,5e-06,1e-07,5e-07,5e-07,5e-07,3e-06,5e-05,
+	8e-06,8e-06,8e-05,8e-06,5e-07,8e-07,8e-07,8e-07,5e-06,8e-05];
+field_values=[
+	(md.results.TransientSolution[0](1).Vx),
+	(md.results.TransientSolution[0](1).Vy),
+	(md.results.TransientSolution[0](1).Vz),
+	(md.results.TransientSolution[0](1).Vel),
+	(md.results.TransientSolution[0](1).Pressure),
+	(md.results.TransientSolution[0](1).Base),
+	(md.results.TransientSolution[0](1).Surface),
+	(md.results.TransientSolution[0](1).Thickness),
+	(md.results.TransientSolution[0](1).Temperature),
+	(md.results.TransientSolution[0](1).BasalforcingsGroundediceMeltingRate),
+	(md.results.TransientSolution[0](2).Vx),
+	(md.results.TransientSolution[0](2).Vy),
+	(md.results.TransientSolution[0](2).Vz),
+	(md.results.TransientSolution[0](2).Vel),
+	(md.results.TransientSolution[0](2).Pressure),
+	(md.results.TransientSolution[0](2).Base),
+	(md.results.TransientSolution[0](2).Surface),
+	(md.results.TransientSolution[0](2).Thickness),
+	(md.results.TransientSolution[0](2).Temperature),
+	(md.results.TransientSolution[0](2).BasalforcingsGroundediceMeltingRate),
+	(md.results.TransientSolution[0](3).Vx),
+	(md.results.TransientSolution[0](3).Vy),
+	(md.results.TransientSolution[0](3).Vz),
+	(md.results.TransientSolution[0](3).Vel),
+	(md.results.TransientSolution[0](3).Pressure),
+	(md.results.TransientSolution[0](3).Base),
+	(md.results.TransientSolution[0](3).Surface),
+	(md.results.TransientSolution[0](3).Thickness),
+	(md.results.TransientSolution[0](3).Temperature),
+	(md.results.TransientSolution[0](3).BasalforcingsGroundediceMeltingRate),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test212.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test212.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test212.js	(revision 21175)
@@ -0,0 +1,36 @@
+//Test Name: SquareShelfCMBSSA2d
+var md = new model();
+triangle(md,square[0],200000.);
+setmask(md,'all','');
+parameterize(md);
+setflowequation(md,'SSA','all');
+
+//control parameters
+md.inversion.iscontrol=1;
+md.inversion.control_parameters=['MaterialsRheologyBbar'];
+md.inversion.min_parameters=Math.pow(10,6)*ones(md.mesh.numberofvertices,1);
+md.inversion.max_parameters=2*Math.pow(10,9)*ones(md.mesh.numberofvertices,1);
+md.inversion.nsteps=2;
+md.inversion.cost_functions=101;
+md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1);
+md.inversion.gradient_scaling=Math.pow(10,7)*ones(md.inversion.nsteps,1);
+md.inversion.maxiter_per_step=2*ones(md.inversion.nsteps,1);
+md.inversion.step_threshold=0.3*ones(md.inversion.nsteps,1);
+md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy;
+md.verbose.control=true;
+
+//md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Stressbalance');
+
+//Fields and tolerances to track changes
+field_names     =['Gradient','Misfits','MaterialsRheologyBbar','Pressure','Vel','Vx','Vy'];
+field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13];
+field_values=[
+	(md.results.StressbalanceSolution[0].Gradient1),
+	(md.results.StressbalanceSolution[0].J),
+	(md.results.StressbalanceSolution[0].MaterialsRheologyBbar),
+	(md.results.StressbalanceSolution[0].Pressure),
+	(md.results.StressbalanceSolution[0].Vel),
+	(md.results.StressbalanceSolution[0].Vx),
+	(md.results.StressbalanceSolution[0].Vy)
+];
Index: /issm/trunk-jpl/test/NightlyRun/test213.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test213.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test213.js	(revision 21175)
@@ -0,0 +1,39 @@
+//Test Name: SquareShelfCMBSSA3d
+var md = new model();
+triangle(md,square[0],200000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,1.);
+setflowequation(md,'SSA','all');
+
+//control parameters
+md.inversion.iscontrol=1;
+md.inversion.control_parameters=['MaterialsRheologyBbar'];
+md.inversion.min_parameters=Math.pow(10,6)*ones(md.mesh.numberofvertices,1);
+md.inversion.max_parameters=2*Math.pow(10,9)*ones(md.mesh.numberofvertices,1);
+md.inversion.nsteps=2;
+md.inversion.cost_functions=101;
+md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1);
+md.inversion.gradient_scaling=Math.pow(10,7)*ones(md.inversion.nsteps,1);
+md.inversion.maxiter_per_step=2*ones(md.inversion.nsteps,1);
+md.inversion.step_threshold=0.3*ones(md.inversion.nsteps,1);
+md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy;
+
+md.verbose.control=true;
+md.verbose.solution=true;
+
+//md.cluster=generic('name',oshostname(),'np',1);
+md=solve(md,'Stressbalance');
+
+//Fields and tolerances to track changes
+field_names     =['Gradient','Misfits','MaterialsRheologyBbar','Pressure','Vel','Vx','Vy'];
+field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13];
+field_values=[
+	(md.results.StressbalanceSolution[0].Gradient1),
+	md.results.StressbalanceSolution[0].J,
+	(md.results.StressbalanceSolution[0].MaterialsRheologyBbar),
+	(md.results.StressbalanceSolution[0].Pressure),
+	(md.results.StressbalanceSolution[0].Vel),
+	(md.results.StressbalanceSolution[0].Vx),
+	(md.results.StressbalanceSolution[0].Vy)
+];
Index: /issm/trunk-jpl/test/NightlyRun/test214.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test214.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test214.js	(revision 21175)
@@ -0,0 +1,36 @@
+//Test Name: SquareShelfCMBHO
+var md = new model();
+triangle(md,square[0],200000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,1.);
+setflowequation(md,'HO','all');
+
+//control parameters
+md.inversion.iscontrol=1;
+md.inversion.control_parameters=['MaterialsRheologyBbar'];
+md.inversion.min_parameters=Math.pow(10,6)*ones(md.mesh.numberofvertices,1);
+md.inversion.max_parameters=2*Math.pow(10,9)*ones(md.mesh.numberofvertices,1);
+md.inversion.nsteps=2;
+md.inversion.cost_functions=101;
+md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1);
+md.inversion.gradient_scaling=Math.pow(10,7)*ones(md.inversion.nsteps,1);
+md.inversion.maxiter_per_step=2*ones(md.inversion.nsteps,1);
+md.inversion.step_threshold=0.3*ones(md.inversion.nsteps,1);
+md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy;
+
+//md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Stressbalance');
+
+//Fields and tolerances to track changes
+field_names     =['Gradient','Misfits','MaterialsRheologyBbar','Pressure','Vel','Vx','Vy'];
+field_tolerances=[1e-07,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08];
+field_values=[
+	(md.results.StressbalanceSolution[0].Gradient1),
+	(md.results.StressbalanceSolution[0].J),
+	(md.results.StressbalanceSolution[0].MaterialsRheologyBbar),
+	(md.results.StressbalanceSolution[0].Pressure),
+	(md.results.StressbalanceSolution[0].Vel),
+	(md.results.StressbalanceSolution[0].Vx),
+	(md.results.StressbalanceSolution[0].Vy)
+];
Index: /issm/trunk-jpl/test/NightlyRun/test215.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test215.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test215.js	(revision 21175)
@@ -0,0 +1,36 @@
+//Test Name: SquareShelfCMBFS
+var md = new model();
+triangle(md,square[0],200000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,1.);
+setflowequation(md,'FS','all');
+
+//control parameters
+md.inversion.iscontrol=1;
+md.inversion.control_parameters=['MaterialsRheologyBbar'];
+md.inversion.min_parameters=Math.pow(10,6)*ones(md.mesh.numberofvertices,1);
+md.inversion.max_parameters=2*Math.pow(10,9)*ones(md.mesh.numberofvertices,1);
+md.inversion.nsteps=2;
+md.inversion.cost_functions=101;
+md.inversion.cost_functions_coeffecients=ones(md.mesh.numberofvertices,1);
+md.inversion.gradient_scaling=Math.pow(10,7)*ones(md.inversion.nsteps,1);
+md.inversion.maxiter_per_step=2*ones(md.inversion.nsteps,1);
+md.inversion.step_threshold=0.3*ones(md.inversion.nsteps,1);
+md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy;
+
+//md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Stressbalance');
+
+//Fields and tolerances to track changes
+field_names     =['Gradient','Misfits','MaterialsRheologyBbar','Pressure','Vel','Vx','Vy'];
+field_tolerances=[4.6e-08,1e-08,2e-08,1e-09,2e-09,5e-09,2e-09];
+field_values=[
+	(md.results.StressbalanceSolution[0].Gradient1),
+	(md.results.StressbalanceSolution[0].J),
+	(md.results.StressbalanceSolution[0].MaterialsRheologyBbar),
+	(md.results.StressbalanceSolution[0].Pressure),
+	(md.results.StressbalanceSolution[0].Vel),
+	(md.results.StressbalanceSolution[0].Vx),
+	(md.results.StressbalanceSolution[0].Vy)
+];
Index: /issm/trunk-jpl/test/NightlyRun/test216.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test216.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test216.js	(revision 21175)
@@ -0,0 +1,26 @@
+//Test Name: SquareShelfStressSSA2dRift
+var md = new model();
+triangle(md,rifts[0],50000.);
+meshprocessrifts(md,'../Exp/Square.exp');
+setmask(md,'all','');
+parameterize(md);
+setflowequation(md,'SSA','all');
+//md.cluster=generic('name',oshostname(),'np',3);
+
+//rift settings
+md.rifts.riftstruct.fill='Melange';
+md.rifts.riftstruct.fraction=0;
+md.stressbalance.rift_penalty_lock=2;
+md.stressbalance.rift_penalty_threshold=0;
+md.rifts.riftstruct.fractionincrement=.1;
+md=solve(md,'Stressbalance');
+
+//Fields and tolerances to track changes
+field_names     =['Vx','Vy','Vel','Pressure'];
+field_tolerances=[8e-7,5e-8,7e-8,2e-11];
+field_values=[
+	(md.results.StressbalanceSolution[0].Vx),
+	(md.results.StressbalanceSolution[0].Vy),
+	(md.results.StressbalanceSolution[0].Vel),
+	(md.results.StressbalanceSolution[0].Pressure),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test217.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test217.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test217.js	(revision 21175)
@@ -0,0 +1,96 @@
+//Test Name: SquareShelfConstrained
+var md = new model();
+function zeros(...args) {
+	var array = [];
+	for (var i = 0; i < args[0]; ++i) {
+		array.push(args.length == 1 ? 0 : zeros(args.slice(1)));
+	}
+	return array;
+}
+var md = new model();
+triangle(md,square[0],150000.);
+setmask(md,'all','');
+parameterize(md);
+setflowequation(md,'SSA','all');
+//md.cluster=generic('name',oshostname(),'np',3);
+
+//redo the parameter file for this special shelf. 
+//constant thickness, constrained (vy=0) flow into an icefront, 
+//from 0 m/yr at the grounding line.
+
+//tighten
+md.stressbalance.restol=Math.pow(10,-4);
+
+//needed later
+ymin=Math.min.apply(null, md.mesh.y);
+ymax=Math.max.apply(null, md.mesh.y);
+xmin=Math.min.apply(null, md.mesh.x);
+xmax=Math.max.apply(null, md.mesh.x);
+
+di=md.materials.rho_ice/md.materials.rho_water;
+
+h=1000;
+md.geometry.thickness=h*ones(md.mesh.numberofvertices,1);
+md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness;
+md.geometry.surface=md.geometry.base+md.geometry.thickness;
+
+//Initial velocity and pressure
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
+
+//Materials
+md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
+md.materials.rheology_B=paterson(md.initialization.temperature);
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+
+//Boundary conditions:
+md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
+md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
+
+//constrain flanks to 0 normal velocity
+pos=[];
+for (var i = 0; i < md.mesh.x.length; ++i) {
+    if (md.mesh.x[i]==xmin || md.mesh.x[i]==xmax) {
+            pos.push(i);
+    };
+}
+
+for (var i = 0; i < pos.length; ++i) {
+    md.stressbalance.spcvx[pos[i]] = 0;
+    md.stressbalance.spcvz[pos[i]] = NaN;
+}
+
+//constrain grounding line to 0 velocity
+pos=[];
+for (var i = 0; i < md.mesh.y.length; ++i) {
+    if (md.mesh.y[i]==ymin) {
+            pos.push(i);
+    };
+}
+for (var i = 0; i < pos.length; ++i) {
+    md.stressbalance.spcvx[pos[i]] = 0;
+    md.stressbalance.spcvy[pos[i]] = 0;
+}
+
+//icefront
+nodeonicefront=zeros(md.mesh.numberofvertices,1);
+for (var i = 0; i < md.mesh.y.length; ++i) {
+    if (md.mesh.y[i]==ymax) {
+        nodeonicefront[i] = 1;
+    }
+}
+md.mask.ice_levelset = nodeonicefront.map(function(x) { return -1 + x; });
+
+md=solve(md,'Stressbalance');
+
+//create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)Math.pow(,3) (Paterson, 4th Edition, page 292.
+//ey_c=(md.materials.rho_ice*md.constants.g*(1-di)*md.geometry.thickness./(4*md.materials.rheology_B)).Math.pow(,3);
+//vy_c=ey_c.*md.mesh.y*md.constants.yts;
+
+//Fields and tolerances to track changes
+field_names     =['Vy'];
+field_tolerances=[1e-13];
+field_values=[(md.results.StressbalanceSolution[0].Vy)];
Index: /issm/trunk-jpl/test/NightlyRun/test218.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test218.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test218.js	(revision 21175)
@@ -0,0 +1,114 @@
+//Test Name: SquareShelfConstrainedDakotaB
+var md = new model();
+function zeros(...args) {
+	var array = [];
+	for (var i = 0; i < args[0]; ++i) {
+		array.push(args.length == 1 ? 0 : zeros(args.slice(1)));
+	}
+	return array;
+}
+var md = new model();
+squaremesh(model(),1000000,1000000,5,5);
+setmask(md,'all','');
+parameterize(md);
+setflowequation(md,'SSA','all');
+//md.cluster=generic('name',oshostname(),'np',3);
+
+//redo the parameter file for this special shelf. 
+//constant thickness, constrained (vy=0) flow into an icefront, 
+//from 0 m/yr at the grounding line.
+
+//needed later
+ymin=Math.min.apply(null, md.mesh.y);
+ymax=Math.max.apply(null, md.mesh.y);
+xmin=Math.min.apply(null, md.mesh.x);
+xmax=Math.max.apply(null, md.mesh.x);
+
+di=md.materials.rho_ice/md.materials.rho_water;
+
+h=1000;
+md.geometry.thickness=h*ones(md.mesh.numberofvertices,1);
+md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness;
+md.geometry.surface=md.geometry.base+md.geometry.thickness;
+
+//Initial velocity and pressure
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
+
+//Materials
+md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
+md.materials.rheology_B=paterson(md.initialization.temperature);
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+
+//Boundary conditions:
+md.stressbalance.spcvx=null*ones(md.mesh.numberofvertices,1);
+md.stressbalance.spcvy=null*ones(md.mesh.numberofvertices,1);
+md.stressbalance.spcvz=null*ones(md.mesh.numberofvertices,1);
+
+//constrain flanks to 0 normal velocity
+pos=[];
+for (var i = 0; i < md.mesh.x==xmin | md.mesh.x==xmax.length; ++i) {
+if ((md.mesh.x==xmin | md.mesh.x==xmax[i] !== 0)) {
+		pos.push(i);
+};
+md.stressbalance.spcvx(pos)=0;
+md.stressbalance.spcvz(pos)=null;
+
+//constrain grounding line to 0 velocity
+pos=[];
+for (var i = 0; i < md.mesh.y==ymin.length; ++i) {
+if ((md.mesh.y==ymin[i] !== 0)) {
+		pos.push(i);
+};
+md.stressbalance.spcvx(pos)=0;
+md.stressbalance.spcvy(pos)=0;
+
+//partitioning
+md.qmu.numberofpartitions=md.mesh.numberofvertices;
+partitioner(md,'package','linear','npart',md.qmu.numberofpartitions);
+md.qmu.partition=md.qmu.partition-1;
+
+//Dakota options
+
+//dakota version
+version=IssmConfig('_DAKOTA_VERSION_'); version=version.slice(0,2); version=str2num(version);
+
+//variables
+md.qmu.variables.rheology_B=normal_uncertain('scaled_MaterialsRheologyB',1,.05);
+
+//responses
+md.qmu.responses.MaxVel=response_function('MaxVel',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+
+//method
+md.qmu.method     =dakota_method('nond_l');
+
+//parameters
+md.qmu.params.direct=true;
+md.qmu.params.interval_type='forward';
+
+if (version>=6,) {
+	md.qmu.params.analysis_driver='matlab';
+	md.qmu.params.evaluation_scheduling='master';
+	md.qmu.params.processors_per_evaluation=2;
+} else {
+	md.qmu.params.analysis_driver='stressbalance';
+	md.qmu.params.evaluation_concurrency=1;
+}
+
+//imperative! 
+md.stressbalance.reltol=10^-10; //tighten for qmu analysese
+md.qmu.isdakota=1;
+
+//solve
+md=solve(md,'Stressbalance','overwrite','y');
+
+//Fields and tolerances to track changes
+md.qmu.results=md.results.dakota;
+md.results.dakota.importancefactors=importancefactors(md,'scaled_MaterialsRheologyB','MaxVel')';
+field_names     =['importancefactors'];
+field_tolerances=[1e-10];
+field_values=[
+         md.results.dakota.importancefactors,
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test219.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test219.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test219.js	(revision 21175)
@@ -0,0 +1,20 @@
+//Test Name: SquareShelfStressSSAHOTiling
+var md = new model();
+triangle(md,square[0],150000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,3,2.);
+setflowequation(md,'HO','../Exp/SquareHalfRight.exp','fill','SSA');
+//md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Stressbalance');
+
+//Fields and tolerances to track changes
+field_names     =['Vx','Vy','Vz','Vel','Pressure'];
+field_tolerances=[1e-09,1e-09,1e-09,1e-09,1e-13];
+field_values=[
+	(md.results.StressbalanceSolution[0].Vx),
+	(md.results.StressbalanceSolution[0].Vy),
+	(md.results.StressbalanceSolution[0].Vz),
+	(md.results.StressbalanceSolution[0].Vel),
+	(md.results.StressbalanceSolution[0].Pressure),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test220.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test220.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test220.js	(revision 21175)
@@ -0,0 +1,21 @@
+//Test Name: SquareShelfStressHOFS3dTiling
+var md = new model();
+triangle(md,square[0],120000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,2,1.);
+setflowequation(md,'FS','../Exp/SquareHalfRight.exp','fill','HO');
+//md.cluster=generic('name',oshostname(),'np',3);
+md.verbose=verbose('all');
+md=solve(md,'Stressbalance');
+
+//Fields and tolerances to track changes
+field_names     =['Vx','Vy','Vz','Vel','Pressure'];
+field_tolerances=[1e-09,1e-09,5e-06,1e-09,1e-09];
+field_values=[
+	(md.results.StressbalanceSolution[0].Vx),
+	(md.results.StressbalanceSolution[0].Vy),
+	(md.results.StressbalanceSolution[0].Vz),
+	(md.results.StressbalanceSolution[0].Vel),
+	(md.results.StressbalanceSolution[0].Pressure),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test221.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test221.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test221.js	(revision 21175)
@@ -0,0 +1,21 @@
+//Test Name: SquareShelfStressSSAFS3dTiling
+var md = new model();
+triangle(md,square[0],120000.);
+setmask(md,'all','');
+parameterize(md);
+md.extrude(md,2,1.);
+setflowequation(md,'FS','../Exp/SquareHalfRight.exp','fill','SSA');
+//md.cluster=generic('name',oshostname(),'np',3);
+md.stressbalance.viscosity_overshoot=0;
+md=solve(md,'Stressbalance');
+
+//Fields and tolerances to track changes
+field_names     =['Vx','Vy','Vz','Vel','Pressure'];
+field_tolerances=[1e-09,1e-09,5e-06,1e-09,1e-09];
+field_values=[
+	(md.results.StressbalanceSolution[0].Vx),
+	(md.results.StressbalanceSolution[0].Vy),
+	(md.results.StressbalanceSolution[0].Vz),
+	(md.results.StressbalanceSolution[0].Vel),
+	(md.results.StressbalanceSolution[0].Pressure),
+	];
Index: /issm/trunk-jpl/test/NightlyRun/test222.js
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test222.js	(revision 21175)
+++ /issm/trunk-jpl/test/NightlyRun/test222.js	(revision 21175)
@@ -0,0 +1,36 @@
+//Test Name: SquareShelfStressSSA2dTransientIncrHydro
+var md = new model();
+triangle(md,square[0],150000.);
+setmask(md,'all','');
+parameterize(md);
+setflowequation(md,'SSA','all');
+//md.cluster=generic('name',oshostname(),'np',3);
+md.masstransport.hydrostatic_adjustment='Incremental';
+md=solve(md,'Transient');
+
+//Fields and tolerances to track changes
+field_names     =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3'];
+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,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13];
+field_values=[
+	(md.results.TransientSolution[0](1).Vx),
+	(md.results.TransientSolution[0](1).Vy),
+	(md.results.TransientSolution[0](1).Vel),
+	(md.results.TransientSolution[0](1).Pressure),
+	(md.results.TransientSolution[0](1).Base),
+	(md.results.TransientSolution[0](1).Surface),
+	(md.results.TransientSolution[0](1).Thickness),
+	(md.results.TransientSolution[0](2).Vx),
+	(md.results.TransientSolution[0](2).Vy),
+	(md.results.TransientSolution[0](2).Vel),
+	(md.results.TransientSolution[0](2).Pressure),
+	(md.results.TransientSolution[0](2).Base),
+	(md.results.TransientSolution[0](2).Surface),
+	(md.results.TransientSolution[0](2).Thickness),
+	(md.results.TransientSolution[0](3).Vx),
+	(md.results.TransientSolution[0](3).Vy),
+	(md.results.TransientSolution[0](3).Vel),
+	(md.results.TransientSolution[0](3).Pressure),
+	(md.results.TransientSolution[0](3).Base),
+	(md.results.TransientSolution[0](3).Surface),
+	(md.results.TransientSolution[0](3).Thickness),
+	];