Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.par
===================================================================
--- /issm/trunk-jpl/test/Par/SquareSheetConstrained.par	(revision 13635)
+++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.par	(revision 13636)
@@ -6,6 +6,7 @@
 ymin=min(md.mesh.y);
 ymax=max(md.mesh.y);
+
 md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin);
-md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+20;
+md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+20.;
 md.geometry.surface=md.geometry.bed+md.geometry.thickness;
 
@@ -16,4 +17,5 @@
 vy    = transpose(ncread('../Data/SquareSheetConstrained.nc','vy'));
 index = transpose(ncread('../Data/SquareSheetConstrained.nc','index'));
+
 md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y);
 md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y);
@@ -23,11 +25,11 @@
 
 %Materials
-md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
+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);
+md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
 
 %Friction
 pos=find(md.mask.elementonfloatingice);
-md.friction.coefficient=20*ones(md.mesh.numberofvertices,1);
+md.friction.coefficient=20.*ones(md.mesh.numberofvertices,1);
 md.friction.coefficient(md.mesh.elements(pos,:))=0;
 md.friction.p=ones(md.mesh.numberofelements,1);
@@ -36,6 +38,6 @@
 %Numerical parameters
 md.diagnostic.viscosity_overshoot=0.0;
-md.prognostic.stabilization=1;
-md.thermal.stabilization=1;
+md.prognostic.stabilization=1.;
+md.thermal.stabilization=1.;
 md.verbose=verbose(0);
 md.settings.waitonlock=30;
@@ -44,6 +46,6 @@
 md.diagnostic.reltol=0.05;
 md.diagnostic.abstol=NaN;
-md.timestepping.time_step=1;
-md.timestepping.final_time=3;
+md.timestepping.time_step=1.;
+md.timestepping.final_time=3.;
 
 %Boundary conditions:
Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 13636)
+++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 13636)
@@ -0,0 +1,67 @@
+import os.path
+import netCDF4
+import numpy
+import inspect
+from verbose import *
+from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
+from paterson import *
+from SetIceSheetBC import *
+
+#Start defining model parameters here
+
+#Geometry
+hmin=300
+hmax=1000
+ymin=numpy.min(md.mesh.y)
+ymax=numpy.max(md.mesh.y)
+
+md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)
+md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+20.
+md.geometry.surface=md.geometry.bed+md.geometry.thickness
+
+#Initial velocity 
+f = netCDF4.Dataset('../Data/SquareSheetConstrained.nc','r')
+x     = numpy.reshape(f.variables['x'][:],(-1))
+y     = numpy.reshape(f.variables['y'][:],(-1))
+vx    = f.variables['vx'][:]
+vy    = f.variables['vy'][:]
+index = f.variables['index'][:]
+f.close()
+
+[md.initialization.vx]=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
+[md.initialization.vy]=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y)
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))
+
+#Materials
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.materials.rheology_B=paterson(md.initialization.temperature)
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+
+#Friction
+pos=numpy.nonzero(md.mask.elementonfloatingice)
+md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient[md.mesh.elements[pos,:].astype(int)-1]=0.
+md.friction.p=numpy.ones((md.mesh.numberofelements,1))
+md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+
+#Numerical parameters
+md.diagnostic.viscosity_overshoot=0.0
+md.prognostic.stabilization=1.
+md.thermal.stabilization=1.
+md.verbose=verbose(0)
+md.settings.waitonlock=30.
+md.diagnostic.restol=0.05
+md.steadystate.reltol=0.05
+md.diagnostic.reltol=0.05
+md.diagnostic.abstol=float('NaN')
+md.timestepping.time_step=1.
+md.timestepping.final_time=3.
+
+#Boundary conditions:
+md=SetIceSheetBC(md)
+
+#Change name so that no test have the same name
+if len(inspect.stack()) > 2:
+	md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]
+
Index: /issm/trunk-jpl/test/Par/SquareShelfConstrained.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 13635)
+++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 13636)
@@ -1,6 +1,6 @@
 import os.path
+import netCDF4
+import numpy
 import inspect
-import netCDF4
-from numpy import *
 from verbose import *
 from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
@@ -10,12 +10,12 @@
 #Start defining model parameters here
 #Geometry
-hmin = 300.
-hmax = 1000.
-ymin = min(md.mesh.y)
-ymax = max(md.mesh.y)
+hmin=300.
+hmax=1000.
+ymin=numpy.min(md.mesh.y)
+ymax=numpy.max(md.mesh.y)
 
-md.geometry.thickness = hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)
+md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)
 md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
-md.geometry.surface = md.geometry.bed+md.geometry.thickness
+md.geometry.surface=md.geometry.bed+md.geometry.thickness
 
 #Initial velocity 
@@ -29,35 +29,41 @@
 f.close()
 
-[md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)
-[md.initialization.vy] = InterpFromMeshToMesh2d(index, x, y, vy, md.mesh.x, md.mesh.y)
-md.initialization.vz = zeros((md.mesh.numberofvertices,1))
-md.initialization.pressure = zeros((md.mesh.numberofvertices,1))
+[md.initialization.vx]=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
+[md.initialization.vy]=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y)
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.pressure=numpy.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))
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.materials.rheology_B=paterson(md.initialization.temperature)
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+
 #Surface mass balance and basal melting
-md.surfaceforcings.mass_balance = 10.*ones((md.mesh.numberofvertices,1))
-md.basalforcings.melting_rate = 5.*ones((md.mesh.numberofvertices,1))
+md.surfaceforcings.mass_balance=10.*numpy.ones((md.mesh.numberofvertices,1))
+md.basalforcings.melting_rate=5.*numpy.ones((md.mesh.numberofvertices,1))
+
 #Friction
-pos = nonzero(md.mask.elementonfloatingice)
-md.friction.coefficient = 20.*ones((md.mesh.numberofvertices,1))
-md.friction.coefficient[md.mesh.elements[pos,:].astype(int)-1] =0.
-md.friction.p = ones((md.mesh.numberofelements,1))
-md.friction.q = ones((md.mesh.numberofelements,1))
+pos=numpy.nonzero(md.mask.elementonfloatingice)
+md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient[md.mesh.elements[pos,:].astype(int)-1]=0.
+md.friction.p=numpy.ones((md.mesh.numberofelements,1))
+md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+
 #Numerical parameters
-md.diagnostic.viscosity_overshoot = 0.0
-md.prognostic.stabilization = 1.
-md.thermal.stabilization = 1.
-md.verbose = verbose()
-md.settings.waitonlock = 30.
-md.diagnostic.restol = 0.05
-md.diagnostic.reltol = 0.05
-md.steadystate.reltol = 0.05
-md.diagnostic.abstol = float('nan')
-md.timestepping.time_step = 1.
-md.timestepping.final_time = 3.
+md.diagnostic.viscosity_overshoot=0.0
+md.prognostic.stabilization=1.
+md.thermal.stabilization=1.
+md.verbose = verbose(0)
+md.settings.waitonlock=30.
+md.diagnostic.restol=0.05
+md.diagnostic.reltol=0.05
+md.steadystate.reltol=0.05
+md.diagnostic.abstol=float('nan')
+md.timestepping.time_step=1.
+md.timestepping.final_time=3.
+
 #Deal with boundary conditions:
 md = SetIceShelfBC(md)
+
 #Change name so that no tests have the same name
 if len(inspect.stack()) > 2: