Index: /issm/trunk-jpl/test/Par/79North.py
===================================================================
--- /issm/trunk-jpl/test/Par/79North.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/79North.py	(revision 21409)
@@ -26,20 +26,20 @@
 
 #Materials
-md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices))
 md.materials.rheology_B=paterson(md.initialization.temperature)
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 md.initialization.temperature=md.initialization.temperature
 
 #Friction
-md.friction.coefficient=50.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=50.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 #Ice shelf melting and surface mass balance 
-md.basalforcings.floatingice_melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
+md.basalforcings.floatingice_melting_rate=numpy.zeros((md.mesh.numberofvertices))
 md.basalforcings.floatingice_melting_rate[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.basalforcings.groundedice_melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
-md.smb.mass_balance=15*numpy.ones((md.mesh.numberofvertices,1))
+md.basalforcings.groundedice_melting_rate=numpy.zeros((md.mesh.numberofvertices))
+md.smb.mass_balance=15*numpy.ones((md.mesh.numberofvertices))
 
 #Numerical parameters
Index: /issm/trunk-jpl/test/Par/ISMIPA.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPA.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/ISMIPA.py	(revision 21409)
@@ -5,17 +5,17 @@
 
 print "      creating thickness"
-md.geometry.surface=-md.mesh.x.reshape(-1,1)*numpy.tan(0.5*numpy.pi/180.)
-md.geometry.base=md.geometry.surface-1000.+500.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))*numpy.sin(md.mesh.y.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))
+md.geometry.surface=-md.mesh.x*numpy.tan(0.5*numpy.pi/180.)
+md.geometry.base=md.geometry.surface-1000.+500.*numpy.sin(md.mesh.x*2.*numpy.pi/numpy.max(md.mesh.x))*numpy.sin(md.mesh.y*2.*numpy.pi/numpy.max(md.mesh.x))
 md.geometry.thickness=md.geometry.surface-md.geometry.base
 
 print "      creating drag"
-md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 print "      creating flow law parameter"
-md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 print "      boundary conditions for stressbalance model"
Index: /issm/trunk-jpl/test/Par/ISMIPB.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPB.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/ISMIPB.py	(revision 21409)
@@ -5,17 +5,17 @@
 
 print "      creating thickness"
-md.geometry.surface=-md.mesh.x.reshape(-1,1)*numpy.tan(0.5*numpy.pi/180.)
-md.geometry.base=md.geometry.surface-1000.+500.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))
+md.geometry.surface=-md.mesh.x*numpy.tan(0.5*numpy.pi/180.)
+md.geometry.base=md.geometry.surface-1000.+500.*numpy.sin(md.mesh.x*2.*numpy.pi/numpy.max(md.mesh.x))
 md.geometry.thickness=md.geometry.surface-md.geometry.base
 
 print "      creating drag"
-md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 print "      creating flow law parameter"
-md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 print "      boundary conditions for stressbalance model"
Index: /issm/trunk-jpl/test/Par/ISMIPC.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPC.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/ISMIPC.py	(revision 21409)
@@ -5,5 +5,5 @@
 
 print "      creating thickness"
-md.geometry.surface=2000.-md.mesh.x.reshape(-1,1)*numpy.tan(0.1*numpy.pi/180.)    #to have z>0
+md.geometry.surface=2000.-md.mesh.x*numpy.tan(0.1*numpy.pi/180.)    #to have z>0
 md.geometry.base=md.geometry.surface-1000.
 md.geometry.thickness=md.geometry.surface-md.geometry.base
@@ -11,12 +11,12 @@
 print "      creating drag"
 #md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x/2.)).*sin(md.mesh.y*2.*pi/max(md.mesh.x/2.)))./(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.base)));
-md.friction.coefficient=numpy.sqrt(md.constants.yts*(1000.+1000.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))*numpy.sin(md.mesh.y.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))))
+md.friction.coefficient=numpy.sqrt(md.constants.yts*(1000.+1000.*numpy.sin(md.mesh.x*2.*numpy.pi/numpy.max(md.mesh.x))*numpy.sin(md.mesh.y*2.*numpy.pi/numpy.max(md.mesh.x))))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.zeros((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.zeros((md.mesh.numberofelements))
 
 print "      creating flow law parameter"
-md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 print "      boundary conditions for stressbalance model:"
Index: /issm/trunk-jpl/test/Par/ISMIPD.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPD.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/ISMIPD.py	(revision 21409)
@@ -5,17 +5,17 @@
 
 print "      creating thickness"
-md.geometry.surface=2000.-md.mesh.x.reshape(-1,1)*numpy.tan(0.1*numpy.pi/180.)    #to have z>0
+md.geometry.surface=2000.-md.mesh.x*numpy.tan(0.1*numpy.pi/180.)    #to have z>0
 md.geometry.base=md.geometry.surface-1000.
 md.geometry.thickness=md.geometry.surface-md.geometry.base
 
 print "      creating drag"
-md.friction.coefficient=numpy.sqrt(md.constants.yts*(1000.+1000.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))))
+md.friction.coefficient=numpy.sqrt(md.constants.yts*(1000.+1000.*numpy.sin(md.mesh.x*2.*numpy.pi/numpy.max(md.mesh.x))))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.zeros((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.zeros((md.mesh.numberofelements))
 
 print "      creating flow law parameter"
-md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 print "      boundary conditions for stressbalance model:"
Index: /issm/trunk-jpl/test/Par/ISMIPE.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPE.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/ISMIPE.py	(revision 21409)
@@ -7,6 +7,6 @@
 print "      creating thickness"
 data = numpy.array(archread('../Data/ISMIPE.arch','data'));
-md.geometry.surface=numpy.zeros((md.mesh.numberofvertices,1))
-md.geometry.base=numpy.zeros((md.mesh.numberofvertices,1))
+md.geometry.surface=numpy.zeros((md.mesh.numberofvertices))
+md.geometry.base=numpy.zeros((md.mesh.numberofvertices))
 for i in xrange(0,md.mesh.numberofvertices):
 	y=md.mesh.y[i]
@@ -21,11 +21,11 @@
 
 print "      creating drag"
-md.friction.coefficient=numpy.zeros((md.mesh.numberofvertices,1))
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.coefficient=numpy.zeros((md.mesh.numberofvertices))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 print "      creating flow law parameter"
-md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 print "      boundary conditions for stressbalance model:"
Index: /issm/trunk-jpl/test/Par/ISMIPF.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPF.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/ISMIPF.py	(revision 21409)
@@ -6,17 +6,17 @@
 
 print "      creating thickness"
-md.geometry.surface=-md.mesh.x.reshape(-1,1)*numpy.tan(3.*numpy.pi/180.)
+md.geometry.surface=-md.mesh.x*numpy.tan(3.*numpy.pi/180.)
 #md.geometry.base=md.geometry.surface-1000.
-md.geometry.base=md.geometry.surface-1000.+100.*numpy.exp(-((md.mesh.x.reshape(-1,1)-numpy.max(md.mesh.x)/2.)**2+(md.mesh.y.reshape(-1,1)-numpy.max(md.mesh.y)/2.)**2)/(10000.**2))
+md.geometry.base=md.geometry.surface-1000.+100.*numpy.exp(-((md.mesh.x-numpy.max(md.mesh.x)/2.)**2+(md.mesh.y-numpy.max(md.mesh.y)/2.)**2)/(10000.**2))
 md.geometry.thickness=md.geometry.surface-md.geometry.base
 
 print "      creating drag"
-md.friction.coefficient=numpy.sqrt(md.constants.yts/(2.140373*10**-7*1000.))*numpy.ones((md.mesh.numberofvertices,1))
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.zeros((md.mesh.numberofelements,1))
+md.friction.coefficient=numpy.sqrt(md.constants.yts/(2.140373*10**-7*1000.))*numpy.ones((md.mesh.numberofvertices))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.zeros((md.mesh.numberofelements))
 
 print "      creating flow law parameter"
-md.materials.rheology_B=1.4734*10**14*numpy.ones((md.mesh.numberofvertices,1))
-md.materials.rheology_n=1.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_B=1.4734*10**14*numpy.ones((md.mesh.numberofvertices))
+md.materials.rheology_n=1.*numpy.ones((md.mesh.numberofelements))
 md.materials.rheology_law='None'
 
@@ -24,18 +24,18 @@
 #Create node on boundary first (because we cannot use mesh)
 md=SetIceSheetBC(md)
-md.stressbalance.spcvx=100.*numpy.ones((md.mesh.numberofvertices,1))
-md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vel=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.temperature=255.*numpy.ones((md.mesh.numberofvertices,1))
+md.stressbalance.spcvx=100.*numpy.ones((md.mesh.numberofvertices))
+md.initialization.vx=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vy=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vel=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.temperature=255.*numpy.ones((md.mesh.numberofvertices))
 pos=numpy.nonzero(numpy.logical_or(numpy.logical_or(md.mesh.x==numpy.min(md.mesh.x),md.mesh.x==numpy.max(md.mesh.x)),numpy.logical_or(md.mesh.y==numpy.min(md.mesh.y),md.mesh.y==numpy.max(md.mesh.y))))
-md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices))
 md.balancethickness.spcthickness[pos]=md.geometry.thickness[pos]
-md.masstransport.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+md.masstransport.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices))
 md.masstransport.spcthickness[pos]=md.geometry.thickness[pos]
-md.thermal.spctemperature=255.*numpy.ones((md.mesh.numberofvertices,1))
-md.basalforcings.geothermalflux=0.4*numpy.ones((md.mesh.numberofvertices,1))
+md.thermal.spctemperature=255.*numpy.ones((md.mesh.numberofvertices))
+md.basalforcings.geothermalflux=0.4*numpy.ones((md.mesh.numberofvertices))
 
 #Parallel options
Index: /issm/trunk-jpl/test/Par/Pig.py
===================================================================
--- /issm/trunk-jpl/test/Par/Pig.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/Pig.py	(revision 21409)
@@ -26,18 +26,18 @@
 md.initialization.vx=md.inversion.vx_obs
 md.initialization.vy=md.inversion.vy_obs
-md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 #Materials
-md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices))
 md.materials.rheology_B=paterson(md.initialization.temperature)
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 md.initialization.temperature=md.initialization.temperature
 
 #Friction
-md.friction.coefficient=50.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=50.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 #Numerical parameters
Index: /issm/trunk-jpl/test/Par/RoundSheetEISMINT.py
===================================================================
--- /issm/trunk-jpl/test/Par/RoundSheetEISMINT.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/RoundSheetEISMINT.py	(revision 21409)
@@ -4,23 +4,23 @@
 #Ok, start defining model parameters here
 print "      creating thickness"
-md.geometry.thickness=10.*numpy.ones((md.mesh.numberofvertices,1))
-md.geometry.base=numpy.zeros((md.mesh.numberofvertices,1))
+md.geometry.thickness=10.*numpy.ones((md.mesh.numberofvertices))
+md.geometry.base=numpy.zeros((md.mesh.numberofvertices))
 md.geometry.surface=md.geometry.base+md.geometry.thickness
 
 print "      creating drag"
-md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1)) 
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices)) 
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 print "      creating temperatures"
 tmin=238.15    #K
 st=1.67*10**-2/1000.    #k/m
-radius=numpy.sqrt((md.mesh.x)**2+(md.mesh.y)**2).reshape(-1,1)
+radius=numpy.sqrt((md.mesh.x)**2+(md.mesh.y)**2)
 md.initialization.temperature=tmin+st*radius
-md.basalforcings.geothermalflux=4.2*10**-2*numpy.ones((md.mesh.numberofvertices,1))
+md.basalforcings.geothermalflux=4.2*10**-2*numpy.ones((md.mesh.numberofvertices))
 
 print "      creating flow law parameter"
-md.materials.rheology_B=6.81*10**7*numpy.ones((md.mesh.numberofvertices,1))    #to have the same B as the analytical solution 
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_B=6.81*10**7*numpy.ones((md.mesh.numberofvertices))    #to have the same B as the analytical solution 
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 print "      creating surface mass balance"
@@ -32,11 +32,11 @@
 print "      creating velocities"
 constant=0.3
-md.inversion.vx_obs=constant/2.*md.mesh.x.reshape(-1,1)*(md.geometry.thickness)**-1
-md.inversion.vy_obs=constant/2.*md.mesh.y.reshape(-1,1)*(md.geometry.thickness)**-1
+md.inversion.vx_obs=constant/2.*md.mesh.x*(md.geometry.thickness)**-1
+md.inversion.vy_obs=constant/2.*md.mesh.y*(md.geometry.thickness)**-1
 md.inversion.vel_obs=numpy.sqrt((md.inversion.vx_obs)**2+(md.inversion.vy_obs)**2)
-md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.vx=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vy=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 #Deal with boundary conditions:
Index: /issm/trunk-jpl/test/Par/RoundSheetShelf.py
===================================================================
--- /issm/trunk-jpl/test/Par/RoundSheetShelf.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/RoundSheetShelf.py	(revision 21409)
@@ -14,5 +14,5 @@
 hmin=300.
 hmax=1000.
-radius=numpy.sqrt(md.mesh.x.reshape(-1,1)*md.mesh.x.reshape(-1,1)+md.mesh.y.reshape(-1,1)*md.mesh.y.reshape(-1,1))
+radius=numpy.sqrt(md.mesh.x*md.mesh.x+md.mesh.y*md.mesh.y.reshape(-1))
 ymin=numpy.min(radius)
 ymax=numpy.max(radius)
@@ -58,26 +58,26 @@
 
 #Initial velocity 
-md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.vx=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vy=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 #Materials
-md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices))
 md.materials.rheology_B=paterson(md.initialization.temperature)
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 #Surface mass balance and basal melting
-md.smb.mass_balance=-10.*numpy.ones((md.mesh.numberofvertices,1))
-md.basalforcings.groundedice_melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
+md.smb.mass_balance=-10.*numpy.ones((md.mesh.numberofvertices))
+md.basalforcings.groundedice_melting_rate=numpy.zeros((md.mesh.numberofvertices))
 pos=numpy.nonzero(md.mask.groundedice_levelset>0.)[0]
 md.basalforcings.groundedice_melting_rate[pos]=10.
-md.basalforcings.floatingice_melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
-md.basalforcings.geothermalflux=numpy.ones((md.mesh.numberofvertices,1))
+md.basalforcings.floatingice_melting_rate=numpy.zeros((md.mesh.numberofvertices))
+md.basalforcings.geothermalflux=numpy.ones((md.mesh.numberofvertices))
 
 #Friction
 radius=1.e6
 shelfextent=2.e5
-md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices))
 xelem=numpy.mean(md.mesh.x[md.mesh.elements.astype(int)-1],axis=1)
 yelem=numpy.mean(md.mesh.y[md.mesh.elements.astype(int)-1],axis=1)
@@ -86,6 +86,6 @@
 pos=numpy.nonzero(rad>=(radius-shelfextent))
 md.friction.coefficient[md.mesh.elements[pos,:]-1]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 #Numerical parameters
@@ -109,7 +109,7 @@
 
 #Deal with boundary conditions:
-md.stressbalance.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-md.stressbalance.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-md.stressbalance.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+md.stressbalance.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices))
+md.stressbalance.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices))
+md.stressbalance.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices))
 
 pos=numpy.nonzero(numpy.logical_and(md.mesh.x==0,md.mesh.y==0))
@@ -119,9 +119,9 @@
 pos=numpy.nonzero(md.mesh.vertexonboundary)
 md.mask.ice_levelset[pos]=0
-md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices))
+md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices))
 md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
 md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
-md.thermal.spctemperature=737.*numpy.ones((md.mesh.numberofvertices,1))
+md.thermal.spctemperature=737.*numpy.ones((md.mesh.numberofvertices))
 
 #Change name so that no test have the same name
Index: /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.py
===================================================================
--- /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.py	(revision 21409)
@@ -5,16 +5,16 @@
 hmin=0.01
 hmax=2756.7
-radius=numpy.sqrt((md.mesh.x)**2+(md.mesh.y)**2).reshape(-1,1)
+radius=numpy.sqrt((md.mesh.x)**2+(md.mesh.y)**2)
 radiusmax=numpy.max(radius)
 radius[numpy.nonzero(radius>(1.-10**-9)*radiusmax)]=radiusmax    #eliminate roundoff issues in next statement
-md.geometry.thickness=hmin*numpy.ones((numpy.size(md.mesh.x),1))+hmax*(4.*((1./2.)**(4./3.)*numpy.ones((numpy.size(md.mesh.x),1))-((radius)/(2.*radiusmax))**(4./3.)))**(3./8.)
+md.geometry.thickness=hmin*numpy.ones((numpy.size(md.mesh.x)))+hmax*(4.*((1./2.)**(4./3.)*numpy.ones((numpy.size(md.mesh.x)))-((radius)/(2.*radiusmax))**(4./3.)))**(3./8.)
 md.geometry.base=0.*md.geometry.thickness
 md.geometry.surface=md.geometry.base+md.geometry.thickness
 
 print "      creating drag"
-md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 print "      creating temperatures"
@@ -22,9 +22,9 @@
 st=1.67*10**-2/1000.    #k/m
 md.initialization.temperature=tmin+st*radius
-md.basalforcings.geothermalflux=4.2*10**-2*numpy.ones((md.mesh.numberofvertices,1))
+md.basalforcings.geothermalflux=4.2*10**-2*numpy.ones((md.mesh.numberofvertices))
 
 print "      creating flow law parameter"
-md.materials.rheology_B=6.81*10**7*numpy.ones((md.mesh.numberofvertices,1))    #to have the same B as the analytical solution 
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_B=6.81*10**7*numpy.ones((md.mesh.numberofvertices))    #to have the same B as the analytical solution 
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 print "      creating surface mass balance"
@@ -36,11 +36,11 @@
 print "      creating velocities"
 constant=0.3
-md.inversion.vx_obs=constant/2.*md.mesh.x.reshape(-1,1)*(md.geometry.thickness)**-1
-md.inversion.vy_obs=constant/2.*md.mesh.y.reshape(-1,1)*(md.geometry.thickness)**-1
+md.inversion.vx_obs=constant/2.*md.mesh.x*(md.geometry.thickness)**-1
+md.inversion.vy_obs=constant/2.*md.mesh.y*(md.geometry.thickness)**-1
 md.inversion.vel_obs=numpy.sqrt((md.inversion.vx_obs)**2+(md.inversion.vy_obs)**2)
-md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.vx=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vy=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 #Deal with boundary conditions:
Index: /issm/trunk-jpl/test/Par/SquareEISMINT.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareEISMINT.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/SquareEISMINT.py	(revision 21409)
@@ -7,30 +7,30 @@
 ymin=numpy.min(md.mesh.y)
 ymax=numpy.max(md.mesh.y)
-md.geometry.thickness=500.*numpy.ones((md.mesh.numberofvertices,1))
+md.geometry.thickness=500.*numpy.ones((md.mesh.numberofvertices))
 md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
 md.geometry.surface=md.geometry.base+md.geometry.thickness
 
 print "      creating drag"
-md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 print "      creating initial values"
-md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
-md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vel=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices))
+md.initialization.vx=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vy=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vel=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 print "      creating flow law parameter"
-md.materials.rheology_B=1.7687*10**8*numpy.ones((md.mesh.numberofvertices,1))
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_B=1.7687*10**8*numpy.ones((md.mesh.numberofvertices))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 print "      creating surface mass balance"
-md.smb.mass_balance=0.2*numpy.ones((md.mesh.numberofvertices,1))    #0m/a
-md.basalforcings.floatingice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices,1))    #0m/a
-md.basalforcings.groundedice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices,1))    #0m/a
+md.smb.mass_balance=0.2*numpy.ones((md.mesh.numberofvertices))    #0m/a
+md.basalforcings.floatingice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices))    #0m/a
+md.basalforcings.groundedice_melting_rate=0.*numpy.ones((md.mesh.numberofvertices))    #0m/a
 
 print "      boundary conditions"
@@ -39,7 +39,7 @@
 #Evolution of the ice shelf
 pos=numpy.nonzero(md.mesh.y==200000.)    #nodes on the upper boundary condition
-md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices))
 md.balancethickness.spcthickness[pos]=500.
-md.masstransport.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+md.masstransport.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices))
 md.masstransport.spcthickness[pos]=500.
 md.masstransport.stabilization=0    #Better result with no artificial diffusivity
Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 21409)
@@ -17,5 +17,5 @@
 xmin=min(md.mesh.x)
 xmax=max(md.mesh.x)
-md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x.reshape(-1,1)-xmin)/(xmax-xmin)
+md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin)
 md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness+20.
 md.geometry.surface=md.geometry.base+md.geometry.thickness
@@ -30,21 +30,21 @@
 [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))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 #Materials
-md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices))
 md.materials.rheology_B=paterson(md.initialization.temperature)
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 #Calving
-md.calving.calvingrate=0.*numpy.ones((md.mesh.numberofvertices,1))
-md.levelset.spclevelset=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+md.calving.calvingrate=0.*numpy.ones((md.mesh.numberofvertices))
+md.levelset.spclevelset=float('NaN')*numpy.ones((md.mesh.numberofvertices))
 
 #Friction
-md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 #Numerical parameters
@@ -62,6 +62,6 @@
 
 #GIA:
-md.gia.lithosphere_thickness=100.*numpy.ones((md.mesh.numberofvertices,1)); # in km
-md.gia.mantle_viscosity=1.*10**21*numpy.ones((md.mesh.numberofvertices,1)); # in Pa.s
+md.gia.lithosphere_thickness=100.*numpy.ones((md.mesh.numberofvertices)); # in km
+md.gia.mantle_viscosity=1.*10**21*numpy.ones((md.mesh.numberofvertices)); # in Pa.s
 md.materials.lithosphere_shear_modulus=6.7*10**10;                          # in Pa
 md.materials.lithosphere_density=3.32;                                      # in g/cm^-3
Index: /issm/trunk-jpl/test/Par/SquareSheetShelf.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareSheetShelf.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/SquareSheetShelf.py	(revision 21409)
@@ -17,5 +17,5 @@
 xmin=min(md.mesh.x)
 xmax=max(md.mesh.x)
-md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x.reshape(-1,1)-xmin)/(xmax-xmin)
+md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin)
 md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
 bed_sheet=-md.materials.rho_ice/md.materials.rho_water*(hmax+(hmin-hmax)*(ymax/2-ymin)/(ymax-ymin))
@@ -33,22 +33,22 @@
 [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))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 #Materials
-md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices))
 md.materials.rheology_B=paterson(md.initialization.temperature)
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 #Accumulation and melting
-md.smb.mass_balance=10.*numpy.ones((md.mesh.numberofvertices,1))
-md.basalforcings.groundedice_melting_rate=5.*numpy.ones((md.mesh.numberofvertices,1))
-md.basalforcings.floatingice_melting_rate=5.*numpy.ones((md.mesh.numberofvertices,1))
+md.smb.mass_balance=10.*numpy.ones((md.mesh.numberofvertices))
+md.basalforcings.groundedice_melting_rate=5.*numpy.ones((md.mesh.numberofvertices))
+md.basalforcings.floatingice_melting_rate=5.*numpy.ones((md.mesh.numberofvertices))
 
 #Friction
-md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 #Numerical parameters
Index: /issm/trunk-jpl/test/Par/SquareShelf.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareShelf.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/SquareShelf.py	(revision 21409)
@@ -16,5 +16,5 @@
 xmin=min(md.mesh.x)
 xmax=max(md.mesh.x)
-md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x.reshape(-1,1)-xmin)/(xmax-xmin)
+md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin)
 md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
 md.geometry.surface=md.geometry.base+md.geometry.thickness
@@ -35,6 +35,6 @@
 [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))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 #dbg - begin
@@ -51,13 +51,13 @@
 
 #Materials
-md.initialization.temperature = (273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.temperature = (273.-20.)*numpy.ones((md.mesh.numberofvertices))
 md.materials.rheology_B = paterson(md.initialization.temperature)
-md.materials.rheology_n = 3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n = 3.*numpy.ones((md.mesh.numberofelements))
 
 #Friction
-md.friction.coefficient = 20.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient = 20.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p = numpy.ones((md.mesh.numberofelements,1))
-md.friction.q = numpy.ones((md.mesh.numberofelements,1))
+md.friction.p = numpy.ones((md.mesh.numberofelements))
+md.friction.q = numpy.ones((md.mesh.numberofelements))
 
 #Numerical parameters
Index: /issm/trunk-jpl/test/Par/SquareShelf2.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareShelf2.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/SquareShelf2.py	(revision 21409)
@@ -16,5 +16,5 @@
 xmin=min(md.mesh.x)
 xmax=max(md.mesh.x)
-md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x.reshape(-1,1)-xmin)/(xmax-xmin)
+md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin)
 md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
 md.geometry.surface=md.geometry.base+md.geometry.thickness
@@ -34,6 +34,6 @@
 [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))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 #dbg - begin
@@ -50,13 +50,13 @@
 
 #Materials
-md.initialization.temperature = (273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.temperature = (273.-20.)*numpy.ones((md.mesh.numberofvertices))
 md.materials.rheology_B = paterson(md.initialization.temperature)
-md.materials.rheology_n = 3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n = 3.*numpy.ones((md.mesh.numberofelements))
 
 #Friction
-md.friction.coefficient = 20.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient = 20.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p = numpy.ones((md.mesh.numberofelements,1))
-md.friction.q = numpy.ones((md.mesh.numberofelements,1))
+md.friction.p = numpy.ones((md.mesh.numberofelements))
+md.friction.q = numpy.ones((md.mesh.numberofelements))
 
 #Numerical parameters
Index: /issm/trunk-jpl/test/Par/SquareShelfConstrained.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 21409)
@@ -1,5 +1,5 @@
 import os.path
 from arch import *
-import numpy
+import numpy as np
 import inspect
 from verbose import verbose
@@ -12,9 +12,9 @@
 hmin=300.
 hmax=1000.
-ymin=numpy.min(md.mesh.y)
-ymax=numpy.max(md.mesh.y)
+ymin=np.min(md.mesh.y)
+ymax=np.max(md.mesh.y)
 xmin=min(md.mesh.x)
 xmax=max(md.mesh.x)
-md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y.reshape(-1,1)-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x.reshape(-1,1)-xmin)/(xmax-xmin)
+md.geometry.thickness=hmax+(hmin-hmax)*(md.mesh.y-ymin)/(ymax-ymin)+0.1*(hmin-hmax)*(md.mesh.x-xmin)/(xmax-xmin)
 md.geometry.base=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
 md.geometry.surface=md.geometry.base+md.geometry.thickness
@@ -22,32 +22,32 @@
 
 #Initial velocity 
-#x         = numpy.reshape(numpy.array(archread('../Data/SquareShelfConstrained.arch','x')),(-1))
-#y         = numpy.reshape(numpy.array(archread('../Data/SquareShelfConstrained.arch','y')),(-1))
-x         = numpy.array(archread('../Data/SquareShelfConstrained.arch','x'))
-y         = numpy.array(archread('../Data/SquareShelfConstrained.arch','y'))
-vx        = numpy.array(archread('../Data/SquareShelfConstrained.arch','vx'))
-vy        = numpy.array(archread('../Data/SquareShelfConstrained.arch','vy'))
-index     = numpy.array(archread('../Data/SquareShelfConstrained.arch','index').astype(int))
+#x         = np.reshape(np.array(archread('../Data/SquareShelfConstrained.arch','x')),(-1))
+#y         = np.reshape(np.array(archread('../Data/SquareShelfConstrained.arch','y')),(-1))
+x         = np.array(archread('../Data/SquareShelfConstrained.arch','x'))
+y         = np.array(archread('../Data/SquareShelfConstrained.arch','y'))
+vx        = np.array(archread('../Data/SquareShelfConstrained.arch','vx'))
+vy        = np.array(archread('../Data/SquareShelfConstrained.arch','vy'))
+index     = np.array(archread('../Data/SquareShelfConstrained.arch','index').astype(int))
 
 [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))
+md.initialization.vz=np.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=np.zeros((md.mesh.numberofvertices))
 
 #Materials
-md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.temperature=(273.-20.)*np.ones((md.mesh.numberofvertices))
 md.materials.rheology_B=paterson(md.initialization.temperature)
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n=3.*np.ones((md.mesh.numberofelements))
 
 #Surface mass balance and basal melting
-md.smb.mass_balance=10.*numpy.ones((md.mesh.numberofvertices,1))
-md.basalforcings.groundedice_melting_rate=5.*numpy.ones((md.mesh.numberofvertices,1))
-md.basalforcings.floatingice_melting_rate=5.*numpy.ones((md.mesh.numberofvertices,1))
+md.smb.mass_balance=10.*np.ones((md.mesh.numberofvertices))
+md.basalforcings.groundedice_melting_rate=5.*np.ones((md.mesh.numberofvertices))
+md.basalforcings.floatingice_melting_rate=5.*np.ones((md.mesh.numberofvertices))
 
 #Friction
-md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
-md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.coefficient=20.*np.ones((md.mesh.numberofvertices))
+md.friction.coefficient[np.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
+md.friction.p=np.ones((md.mesh.numberofelements))
+md.friction.q=np.ones((md.mesh.numberofelements))
 
 #Numerical parameters
Index: /issm/trunk-jpl/test/Par/SquareThermal.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareThermal.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/SquareThermal.py	(revision 21409)
@@ -9,29 +9,29 @@
 print "      creating thickness"
 h=1000.
-md.geometry.thickness=h*numpy.ones((md.mesh.numberofvertices,1))
-md.geometry.base=-1000.*numpy.ones((md.mesh.numberofvertices,1))
+md.geometry.thickness=h*numpy.ones((md.mesh.numberofvertices))
+md.geometry.base=-1000.*numpy.ones((md.mesh.numberofvertices))
 md.geometry.surface=md.geometry.base+md.geometry.thickness;
 
 print "      creating velocities"
-md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.vx=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vy=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
 
 print "      creating drag"
-md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 print "      creating temperatures"
-md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices))
 
 print "      creating flow law parameter"
 md.materials.rheology_B=paterson(md.initialization.temperature)
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 print "      creating surface mass balance"
-md.smb.mass_balance=numpy.ones((md.mesh.numberofvertices,1))/md.constants.yts    #1m/a
-md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices,1))/md.constants.yts    #1m/a
+md.smb.mass_balance=numpy.ones((md.mesh.numberofvertices))/md.constants.yts    #1m/a
+md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices))/md.constants.yts    #1m/a
 
 #Deal with boundary conditions:
@@ -41,5 +41,5 @@
 
 print "      boundary conditions for thermal model"
-md.thermal.spctemperature[:]=md.initialization.temperature
-md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) 
+md.thermal.spctemperature=md.initialization.temperature
+md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices)) 
 md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)[0]]=1.*10**-3    #1 mW/m^2
Index: /issm/trunk-jpl/test/Par/ValleyGlacierShelf.py
===================================================================
--- /issm/trunk-jpl/test/Par/ValleyGlacierShelf.py	(revision 21408)
+++ /issm/trunk-jpl/test/Par/ValleyGlacierShelf.py	(revision 21409)
@@ -10,6 +10,6 @@
 
 #Start defining model parameters here
-x=md.mesh.x.reshape(-1,1)
-y=md.mesh.y.reshape(-1,1)
+x=md.mesh.x
+y=md.mesh.y
 xmin, xmax = min(x), max(x)
 ymin, ymax = min(y), max(y)
@@ -33,22 +33,22 @@
 #Mask
 md.mask.ice_levelset=x - alpha*Lx
-md.mask.groundedice_levelset= numpy.ones((md.mesh.numberofvertices,1))
+md.mask.groundedice_levelset= numpy.ones((md.mesh.numberofvertices))
 
 #Initial velocity 
-md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.vx=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vy=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices))
 
 #Materials
-md.initialization.temperature=(273.15-5.)*numpy.ones((md.mesh.numberofvertices,1))
-md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1))
-md.initialization.watercolumn=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.temperature=(273.15-5.)*numpy.ones((md.mesh.numberofvertices))
+md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices))
+md.initialization.watercolumn=numpy.zeros((md.mesh.numberofvertices))
 md.materials.rheology_B=paterson(md.initialization.temperature)
-md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements))
 
 #Thermal
 md.thermal.isenthalpy=False
-md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices))
 
 #Groundingline
@@ -56,13 +56,13 @@
 
 #Surface mass balance and basal melting
-md.smb.mass_balance=0.3*numpy.ones((md.mesh.numberofvertices,1))
+md.smb.mass_balance=0.3*numpy.ones((md.mesh.numberofvertices))
 md.basalforcings.groundedice_melting_rate=md.smb.mass_balance
 md.basalforcings.floatingice_melting_rate=md.smb.mass_balance
 
 #Friction
-md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.coefficient=20.*numpy.ones((md.mesh.numberofvertices))
 md.friction.coefficient[numpy.nonzero(md.mask.groundedice_levelset<0.)[0]]=0.
-md.friction.p=numpy.ones((md.mesh.numberofelements,1))
-md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements))
+md.friction.q=numpy.ones((md.mesh.numberofelements))
 
 #Transient
@@ -82,7 +82,7 @@
 
 #Masstransport
-md.calving.calvingrate=0.*numpy.ones((md.mesh.numberofvertices,1))
-md.calving.meltingrate=0.*numpy.ones((md.mesh.numberofvertices,1))
-md.levelset.spclevelset=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+md.calving.calvingrate=0.*numpy.ones((md.mesh.numberofvertices))
+md.calving.meltingrate=0.*numpy.ones((md.mesh.numberofvertices))
+md.levelset.spclevelset=float('NaN')*numpy.ones((md.mesh.numberofvertices))
 md.masstransport.stabilization=1.