Index: /issm/trunk-jpl/test/NightlyRun/test121.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test121.py	(revision 13462)
+++ /issm/trunk-jpl/test/NightlyRun/test121.py	(revision 13463)
@@ -14,5 +14,5 @@
 md=setflowequation(md,'macayeal','all')
 md.cluster=generic('name',oshostname(),'np',3);
-md.initialization.waterfraction=numpy.zeros(md.mesh.numberofvertices)
+md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1))
 md.transient.isdiagnostic=0
 md.transient.isprognostic=0
Index: /issm/trunk-jpl/test/NightlyRun/test122.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test122.py	(revision 13462)
+++ /issm/trunk-jpl/test/NightlyRun/test122.py	(revision 13463)
@@ -13,5 +13,5 @@
 md.extrude(3,1)
 md=setflowequation(md,'pattyn','all')
-md.initialization.waterfraction=numpy.zeros(md.mesh.numberofvertices)
+md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1))
 md.thermal.isenthalpy=1
 md.thermal.stabilization=2
Index: /issm/trunk-jpl/test/Par/Pig.py
===================================================================
--- /issm/trunk-jpl/test/Par/Pig.py	(revision 13463)
+++ /issm/trunk-jpl/test/Par/Pig.py	(revision 13463)
@@ -0,0 +1,63 @@
+import os.path
+import inspect
+import netCDF4
+from numpy import *
+from verbose import *
+from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
+from paterson import *
+from SetIceShelfBC import *
+
+#Start defining model parameters here
+
+#Geometry and observation
+f = netCDF4.Dataset('../Data/Pig.nc','r')
+x         = reshape(f.variables['x'][:],(-1))
+y         = reshape(f.variables['y'][:],(-1))
+vx_obs    = f.variables['vx_obs'][:]
+vy_obs    = f.variables['vy_obs'][:]
+index     = f.variables['index'][:]
+surface   = f.variables['surface'][:]
+thickness = f.variables['thickness'][:]
+f.close()
+
+[md.inversion.vx_obs]   =InterpFromMeshToMesh2d(index,x,y,vx_obs,md.mesh.x,md.mesh.y)
+[md.inversion.vy_obs]   =InterpFromMeshToMesh2d(index,x,y,vy_obs,md.mesh.x,md.mesh.y)
+[md.geometry.surface]  =InterpFromMeshToMesh2d(index,x,y,surface,md.mesh.x,md.mesh.y)
+[md.geometry.thickness]=InterpFromMeshToMesh2d(index,x,y,thickness,md.mesh.x,md.mesh.y)
+md.geometry.bed=md.geometry.surface-md.geometry.thickness
+md.initialization.vx=md.inversion.vx_obs
+md.initialization.vy=md.inversion.vy_obs
+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))
+md.initialization.temperature=md.initialization.temperature
+
+#Friction
+pos=numpy.nonzero(md.mask.elementonfloatingice)
+md.friction.coefficient=50.*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))
+
+#Numerical parameters
+md.diagnostic.viscosity_overshoot=0.3
+md.prognostic.stabilization=1.
+md.verbose=verbose(0)
+md.settings.waitonlock=30.
+md.timestepping.time_step=1.
+md.timestepping.final_time=2.
+md.diagnostic.restol=0.05
+md.diagnostic.reltol=1.
+md.steadystate.reltol=1.
+md.diagnostic.abstol=float('nan')
+
+#Boundary conditions:
+#md=SetMarineIceSheetBC(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 13462)
+++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 13463)
@@ -15,5 +15,5 @@
 ymax = max(md.mesh.y)
 
-md.geometry.thickness = hmax+(hmin-hmax)*(md.mesh.y-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
@@ -26,26 +26,24 @@
 vx=f.variables['vx'][:]
 vy=f.variables['vy'][:]
-#deal with 'F' oriented matlab matrices!
-index=f.variables['index'][:].astype(float)
-index=reshape(index.T,(len(index),3),order='F')
+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 = zeros(md.mesh.numberofvertices)
-md.initialization.pressure = zeros(md.mesh.numberofvertices)
+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)
+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)
+md.materials.rheology_n = 3.*ones((md.mesh.numberofelements,1))
 #Surface mass balance and basal melting
-md.surfaceforcings.mass_balance = 10.*ones(md.mesh.numberofvertices)
-md.basalforcings.melting_rate = 5.*ones(md.mesh.numberofvertices)
+md.surfaceforcings.mass_balance = 10.*ones((md.mesh.numberofvertices,1))
+md.basalforcings.melting_rate = 5.*ones((md.mesh.numberofvertices,1))
 #Friction
 pos = nonzero(md.mask.elementonfloatingice)
-md.friction.coefficient = 20.*ones(md.mesh.numberofvertices)
+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)
-md.friction.q = ones(md.mesh.numberofelements)
+md.friction.p = ones((md.mesh.numberofelements,1))
+md.friction.q = ones((md.mesh.numberofelements,1))
 #Numerical parameters
 md.diagnostic.viscosity_overshoot = 0.0
@@ -57,5 +55,5 @@
 md.diagnostic.reltol = 0.05
 md.steadystate.reltol = 0.05
-md.diagnostic.abstol = nan
+md.diagnostic.abstol = float('nan')
 md.timestepping.time_step = 1.
 md.timestepping.final_time = 3.