Index: /issm/trunk-jpl/src/m/archive/arch.py
===================================================================
--- /issm/trunk-jpl/src/m/archive/arch.py	(revision 21169)
+++ /issm/trunk-jpl/src/m/archive/arch.py	(revision 21170)
@@ -216,7 +216,10 @@
 			rows=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0]
 			cols=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0]
-			data=numpy.zeros(shape=(rows,cols),dtype=float)
+			raw_data=numpy.zeros(shape=(rows,cols),dtype=float)
 			for i in xrange(rows):
-				data[i,:]=struct.unpack('>%dd' % cols,fid.read(cols*struct.calcsize('>d')))
+				raw_data[i,:]=struct.unpack('>%dd' % cols,fid.read(cols*struct.calcsize('>d')))
+			# The matrix will be unpacked in order and will be filled left -> right by column
+			# We need to reshape and transpose the matrix so it can be read correctly
+			data=raw_data.reshape(raw_data.shape[::-1]).T
 		else:
 			raise TypeError("Cannot read data type %d" % data_type)
Index: /issm/trunk-jpl/test/Par/79North.py
===================================================================
--- /issm/trunk-jpl/test/Par/79North.py	(revision 21169)
+++ /issm/trunk-jpl/test/Par/79North.py	(revision 21170)
@@ -1,5 +1,5 @@
 import os.path
 import inspect
-import netCDF4
+from arch import *
 import numpy
 from verbose import verbose 
@@ -11,13 +11,11 @@
 
 #Geometry and observation
-f         = netCDF4.Dataset('../Data/79North.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'][:]
-surface   = f.variables['surface'][:]
-thickness = f.variables['thickness'][:]
-f.close()
+x         = numpy.array(archread('../Data/79North.arch','x'))
+y         = numpy.array(archread('../Data/79North.arch','y'))
+vx        = numpy.array(archread('../Data/79North.arch','vx'));
+vy        = numpy.array(archread('../Data/79North.arch','vy'));
+index     = numpy.array(archread('../Data/79North.arch','index'));
+surface   = numpy.array(archread('../Data/79North.arch','surface')); 
+thickness = numpy.array(archread('../Data/79North.arch','thickness')); 
 
 [md.initialization.vx]  = InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
Index: /issm/trunk-jpl/test/Par/ISMIPE.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPE.py	(revision 21169)
+++ /issm/trunk-jpl/test/Par/ISMIPE.py	(revision 21170)
@@ -1,4 +1,4 @@
 import numpy
-import netCDF4
+from arch import *
 from SetIceSheetBC import SetIceSheetBC
 
@@ -6,7 +6,5 @@
 
 print "      creating thickness"
-f = netCDF4.Dataset('../Data/ISMIPE.nc','r')
-data = f.variables['data'][:]
-f.close()
+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))
Index: /issm/trunk-jpl/test/Par/Pig.py
===================================================================
--- /issm/trunk-jpl/test/Par/Pig.py	(revision 21169)
+++ /issm/trunk-jpl/test/Par/Pig.py	(revision 21170)
@@ -1,5 +1,5 @@
 import os.path
 import inspect
-import netCDF4
+from arch import *
 import numpy
 from verbose import verbose
@@ -11,13 +11,11 @@
 
 #Geometry and observation
-f = netCDF4.Dataset('../Data/Pig.nc','r')
-x         = numpy.reshape(f.variables['x'][:],(-1))
-y         = numpy.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()
+x         = numpy.array(archread('../Data/Pig.arch','x'))
+y         = numpy.array(archread('../Data/Pig.arch','y'))
+vx_obs    = numpy.array(archread('../Data/Pig.arch','vx_obs')) 
+vy_obs    = numpy.array(archread('../Data/Pig.arch','vy_obs')) 
+index     = numpy.array(archread('../Data/Pig.arch','index'))
+surface   = numpy.array(archread('../Data/Pig.arch','surface'))
+thickness = numpy.array(archread('../Data/Pig.arch','thickness')) 
 
 [md.inversion.vx_obs]   =InterpFromMeshToMesh2d(index,x,y,vx_obs,md.mesh.x,md.mesh.y)
Index: /issm/trunk-jpl/test/Par/SquareSheetConstrained.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 21169)
+++ /issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 21170)
@@ -1,4 +1,3 @@
 import os.path
-import netCDF4
 import numpy
 import inspect
@@ -7,4 +6,5 @@
 from paterson import paterson
 from SetIceSheetBC import SetIceSheetBC
+from arch import *
 
 #Start defining model parameters here
@@ -22,11 +22,9 @@
 
 #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()
+x         = numpy.array(archread('../Data/SquareSheetConstrained.arch','x'))
+y         = numpy.array(archread('../Data/SquareSheetConstrained.arch','y'))
+vx        = numpy.array(archread('../Data/SquareSheetConstrained.arch','vx'));
+vy        = numpy.array(archread('../Data/SquareSheetConstrained.arch','vy'));
+index     = archread('../Data/SquareSheetConstrained.arch','index').astype(int);
 
 [md.initialization.vx]=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
Index: /issm/trunk-jpl/test/Par/SquareSheetShelf.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareSheetShelf.py	(revision 21169)
+++ /issm/trunk-jpl/test/Par/SquareSheetShelf.py	(revision 21170)
@@ -1,5 +1,5 @@
 import os.path
 import inspect
-import netCDF4
+from arch import *
 import numpy
 from verbose import verbose
@@ -25,11 +25,9 @@
 
 #Initial velocity 
-f = netCDF4.Dataset('../Data/SquareSheetShelf.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()
+x         = numpy.array(archread('../Data/SquareSheetShelf.arch','x'))
+y         = numpy.array(archread('../Data/SquareSheetShelf.arch','y'))
+vx        = numpy.array(archread('../Data/SquareSheetShelf.arch','vx'));
+vy        = numpy.array(archread('../Data/SquareSheetShelf.arch','vy'));
+index     = numpy.array(archread('../Data/SquareSheetShelf.arch','index'));
 
 [md.initialization.vx]  = InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
Index: /issm/trunk-jpl/test/Par/SquareShelf.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareShelf.py	(revision 21169)
+++ /issm/trunk-jpl/test/Par/SquareShelf.py	(revision 21170)
@@ -1,5 +1,5 @@
 import os.path
 import inspect
-import netCDF4
+from arch import *
 import numpy
 from verbose import verbose
@@ -21,11 +21,9 @@
 
 #Initial velocity and pressure
-iVelF = netCDF4.Dataset('../Data/SquareShelf.nc','r')
-
-x=numpy.reshape(iVelF.variables['x'][:],(-1))
-y=numpy.reshape(iVelF.variables['y'][:],(-1))
-vx=iVelF.variables['vx'][:]
-vy=iVelF.variables['vy'][:]
-index=iVelF.variables['index'][:].astype(int)
+x         = numpy.array(archread('../Data/SquareShelf.arch','x'))
+y         = numpy.array(archread('../Data/SquareShelf.arch','y'))
+vx        = numpy.array(archread('../Data/SquareShelf.arch','vx'));
+vy        = numpy.array(archread('../Data/SquareShelf.arch','vy'));
+index     = archread('../Data/SquareShelf.arch','index').astype(int);
 
 #dbg - begin
@@ -34,6 +32,4 @@
 # #	print v
 #dbg - end 
-
-iVelF.close()
 
 [md.initialization.vx]=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
Index: /issm/trunk-jpl/test/Par/SquareShelf2.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareShelf2.py	(revision 21169)
+++ /issm/trunk-jpl/test/Par/SquareShelf2.py	(revision 21170)
@@ -1,5 +1,5 @@
 import os.path
 import inspect
-import netCDF4
+from arch import *
 import numpy
 from verbose import verbose
@@ -21,12 +21,9 @@
 
 #Initial velocity and pressure
-iVelF = netCDF4.Dataset('../Data/SquareShelf.nc','r')
-
-x=numpy.reshape(iVelF.variables['x'][:],(-1))
-y=numpy.reshape(iVelF.variables['y'][:],(-1))
-vx=iVelF.variables['vx'][:]
-vy=iVelF.variables['vy'][:]
-index=iVelF.variables['index'][:].astype(int)
-
+x         = numpy.array(archread('../Data/SquareShelf.arch','x'))
+y         = numpy.array(archread('../Data/SquareShelf.arch','y'))
+vx        = numpy.array(archread('../Data/SquareShelf.arch','vx'));
+vy        = numpy.array(archread('../Data/SquareShelf.arch','vy'));
+index     = archread('../Data/SquareShelf.arch','index').astype(int);
 #dbg - begin
 # #print 'vars in SquareShelf.nc:'
@@ -34,6 +31,4 @@
 # #	print v
 #dbg - end 
-
-iVelF.close()
 
 [md.initialization.vx]=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
Index: /issm/trunk-jpl/test/Par/SquareShelfConstrained.py
===================================================================
--- /issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 21169)
+++ /issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 21170)
@@ -1,4 +1,4 @@
 import os.path
-import netCDF4
+from arch import *
 import numpy
 import inspect
@@ -22,12 +22,11 @@
 
 #Initial velocity 
-f = netCDF4.Dataset('../Data/SquareShelfConstrained.nc','r')
-#Reshape as Rank-1 arrays
-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()
+#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))
 
 [md.initialization.vx]=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y)
Index: /issm/trunk-jpl/test/Par/ValleyGlacierShelf.py
===================================================================
--- /issm/trunk-jpl/test/Par/ValleyGlacierShelf.py	(revision 21169)
+++ /issm/trunk-jpl/test/Par/ValleyGlacierShelf.py	(revision 21170)
@@ -1,4 +1,4 @@
 import os.path
-import netCDF4
+from arch import *
 import numpy
 import inspect