Index: /issm/trunk-jpl/test/Par/SquareShelfConstrained.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 13237)
+++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 13238)
@@ -1,10 +1,10 @@
+import os.path
+import inspect
+import netCDF4
 from numpy import *
 from verbose import *
-import scipy.io as matio
 from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
-from paterson import  *
+from paterson import *
 from SetIceShelfBC import *
-import inspect
-import os.path
 
 #Start defining model parameters here
@@ -16,17 +16,18 @@
 
 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;
+md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
 md.geometry.surface = md.geometry.bed+md.geometry.thickness
 
 #Initial velocity 
-mat=matio.loadmat('../Data/SquareShelfConstrained.data')
+f = netCDF4.Dataset('../Data/SquareShelfConstrained.nc','r')
 #Reshape as Rank-1 arrays
-x=reshape(mat['x'],(-1))
-y=reshape(mat['y'],(-1))
-vx=mat['vx']
-vy=mat['vy']
+x=reshape(f.variables['x'][:],(-1))
+y=reshape(f.variables['y'][:],(-1))
+vx=f.variables['vx'][:]
+vy=f.variables['vy'][:]
 #deal with 'F' oriented matlab matrices!
-index=mat['index'].astype(float)
+index=f.variables['index'][:].astype(float)
 index=reshape(index.T,(len(index),3),order='F')
+f.close()
 
 [md.initialization.vx] = InterpFromMeshToMesh2d(index, x, y, vx, md.mesh.x, md.mesh.y)