Index: /issm/trunk-jpl/scripts/ConvertAllNetcdfArchivesToArch.m
===================================================================
--- /issm/trunk-jpl/scripts/ConvertAllNetcdfArchivesToArch.m	(revision 21151)
+++ /issm/trunk-jpl/scripts/ConvertAllNetcdfArchivesToArch.m	(revision 21151)
@@ -0,0 +1,20 @@
+root=pwd;
+cd('../test/Archives/');
+ncdir=pwd;
+cd(root);
+
+file_regexp=strcat(ncdir,'/*.nc');
+flist=dir(file_regexp);
+fnames={flist.name}';
+
+pattern='.nc';
+replacement='.arch';
+for i=1:numel(fnames)
+	% Remove the extension '.nc' and replace with '.arch'
+	arch_file_name=regexprep(fnames{i},pattern,replacement);
+	disp(['Converting ' fnames{i} ' to .arch format']);
+	% Reconstruct file path
+	fpath=strcat(ncdir,'/',fnames{i});
+	archpath=strcat(ncdir,'/',arch_file_name);
+	netcdf2arch(fpath,archpath);
+end
Index: /issm/trunk-jpl/src/m/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/m/Makefile.am	(revision 21150)
+++ /issm/trunk-jpl/src/m/Makefile.am	(revision 21151)
@@ -8,4 +8,5 @@
 if !DEVELOPMENT
 bin_SCRIPTS += ${ISSM_DIR}/src/m/qmu/*.m \
+					${ISSM_DIR}/src/m/archive/*.m \
 					${ISSM_DIR}/src/m/qmu/setupdesign/*.m \
 					${ISSM_DIR}/src/m/qmu/plot/*.m \
@@ -50,5 +51,6 @@
 if PYTHON
 if !DEVELOPMENT
-bin_SCRIPTS += ${ISSM_DIR}/src/m/classes/*.py \
+bin_SCRIPTS += ${ISSM_DIR}/src/m/archive/*.py \
+					${ISSM_DIR}/src/m/classes/*.py \
 					${ISSM_DIR}/src/m/classes/clusters/*.py \
 					${ISSM_DIR}/src/m/consistency/*.py \
Index: /issm/trunk-jpl/src/m/archive/arch.py
===================================================================
--- /issm/trunk-jpl/src/m/archive/arch.py	(revision 21151)
+++ /issm/trunk-jpl/src/m/archive/arch.py	(revision 21151)
@@ -0,0 +1,263 @@
+import numpy
+import math
+import struct
+import sys
+import os
+from collections import OrderedDict
+
+def archwrite(filename,*args): # {{{
+	"""
+	ARCHWRITE - Write data to a field, given the file name, field name, and data.
+
+		Usage:
+			archwrite('archive101.arch','variable_name',data);
+	"""
+	nargs=len(args);
+	if nargs % 2 != 0 :
+		raise ValueError('Incorrect number of arguments.')
+	# open file
+	try:
+		if not os.path.isfile(filename):
+			fid=open(filename,'wb')
+		else:
+			fid=open(filename,'ab')
+	except IOError as e:
+		raise IOError("archwrite error: could not open '%s' to write to." % filename)
+
+	nfields=len(args)/2
+	# generate data to write
+	for i in range(nfields):
+		# write field name
+		name=args[2*i]
+		write_field_name(fid,name)
+		
+		# write data associated with field name
+		data=args[2*i+1]
+		code=format_archive_code(data)
+		if code==1:
+			raise ValueError("archwrite : error writing data, string should not be written as field data")
+		elif code==2:
+			write_scalar(fid,data)
+		elif code==3:
+			write_vector(fid,data)
+		else:
+			raise ValueError("archwrite : error writing data, invalid code entered '%d'" % code)
+	
+	fid.close()
+
+# }}}
+def archread(filename,fieldname): # {{{
+	"""
+	ARCHREAD - Given an arch file name, and a field name, find and return the data
+					associated with that field name.
+		Usage:
+			archread('archive101.arch','field_var_1')
+	"""
+	try:
+		if os.path.isfile(filename):
+			fid=open(filename,'rb')
+		else:
+			raise IOError("archread error : file '%s' does not exist" % filename)
+	except IOError as e:
+		raise IOError("archread error : could not open file '%s' to read from" % filename)
+	
+	archive_results=[]
+
+	# read first result
+	result=read_field(fid)
+	while result:
+		if fieldname == result['field_name']:
+			# found the data we wanted
+			archive_results=result['data']; # we only want the data
+			break
+		
+		# read next result
+		result=read_field(fid)
+	
+	# close file
+	fid.close()
+	
+	return archive_results
+# }}}
+def archdisp(filename): # {{{
+	"""
+	ARCHDISP - Given an arch filename, display the contents of that file
+
+		Usage:
+			archdisp('archive101.arch')
+	"""
+	try:
+		if os.path.isfile(filename):
+			fid=open(filename,'rb')
+		else:
+			raise IOError("archread error : file '%s' does not exist" % filename)
+	except IOError as e:
+		raise IOError("archread error : could not open file '%s' to read from" % filename)
+	
+	print 'Source file: '
+	print '\t{0}'.format(filename)
+	print 'Variables: '
+
+	result=read_field(fid)
+	while result:
+		print '\t{0}'.format(result['field_name'])
+		print '\t\tSize:\t\t{0}'.format(result['size'])
+		print '\t\tDatatype:\t{0}'.format(result['data_type'])
+		# go to next result
+		result=read_field(fid)
+	
+	# close file
+	fid.close()
+
+# }}}
+
+# Helper functions 
+def write_field_name(fid,data): # {{{
+	"""
+	Routine to write field name (variable name) to an archive file.
+	"""
+	# write the length of the record
+	# length to write + string size (len) + format code
+	reclen=len(data)+4+4
+	fid.write(struct.pack('>i',reclen))
+	
+	# write format code
+	code=format_archive_code(data);
+	if code != 1:
+		raise TypeError("archwrite : error writing field name, expected string, but got %s" % type(data))
+	fid.write(struct.pack('>i',1))
+
+	# write string length, and then the string
+	fid.write(struct.pack('>i',len(data)))
+	fid.write(struct.pack('>%ds' % len(data),data))
+# }}}
+def write_scalar(fid,data): # {{{
+	"""
+	Procedure to write a double to an arch file pointed to by fid
+	"""
+	# write length of record
+	# double (8 bytes) + format code (4 bytes)
+	reclen=8+4
+	fid.write(struct.pack('>i',reclen))
+
+	# write the format code (2 for scalar)
+	fid.write(struct.pack('>i',2))
+	
+	# write the double
+	fid.write(struct.pack('>d',data))
+
+# }}}
+def write_vector(fid,data): # {{{
+	"""
+	Procedure to write a numpy array to an arch file
+	"""
+	# Make sure our vector is the correct shape.
+	# Reshape it into a row vector if it is not correct.
+	if isinstance(data,(bool,int,long,float)):
+		data=numpy.array([data])
+	elif isinstance(data,(list,tuple)):
+		data=numpy.array(data).reshape(-1,1)
+	
+	if numpy.ndim(data) == 1:
+		if numpy.size(data):
+			data=data.reshape(numpy.size(data),1)
+		else:
+			data=data.reshape(0,0)
+	
+	# get size of data
+	sz=data.shape
+
+	# write length of record
+	# format code + row size + col size + (double size * row amt * col amt)
+	reclen=4+4+4+8*sz[0]*sz[1]
+	# make sure we can fit data into file
+	if reclen>2**31:
+		raise ValueError("archwrite error : can not write vector to binary file because it is too large")
+	fid.write(struct.pack('>i',reclen))
+	
+	# write format code
+	fid.write(struct.pack('>i',3))
+
+	# write vector
+	fid.write(struct.pack('>i',sz[0]))
+	fid.write(struct.pack('>i',sz[1]))
+	for i in xrange(sz[0]):
+		for j in xrange(sz[1]):
+			fid.write(struct.pack('>d',float(data[i][j])))
+
+# }}}
+
+def read_field(fid): # {{{
+	"""
+	Procedure to read a field and return a results list with the following attributes:
+	result['field_name']	-> the name of the variable that was just read
+	result['size']			-> size (dimensions) of the variable just read
+	result['data_type']	-> the type of data that was just read
+	result['data']			-> the actual data
+	"""
+
+	try:
+		# first, read the string
+		reclen=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0]
+		check_name=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0]
+		if check_name != 1:
+			raise ValueError('archread error : a string was not present at the start of the arch file')
+		namelen=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0]
+		fieldname=struct.unpack('>%ds' % namelen,fid.read(namelen))[0]
+		
+		# then, read the data
+		datalen=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0]
+		data_type=struct.unpack('>i',fid.read(struct.calcsize('>i')))[0]
+
+		if data_type==2:
+			# unpack scalar
+			data=struct.unpack('>d',fid.read(struct.calcsize('>d')))[0]
+		elif data_type==3:
+			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)
+			for i in xrange(rows):
+				data[i,:]=struct.unpack('>%dd' % cols,fid.read(cols*struct.calcsize('>d')))
+		else:
+			raise TypeError("Cannot read data type %d" % data_type)
+			
+		# give additional data to user
+		if data_type==2:
+			data_size='1x1'
+			data_type_str='double'
+		elif data_type==3:
+			data_size='{0}x{1}'.format(rows,cols)
+			data_type_str='vector'
+
+		result=OrderedDict()
+		result['field_name']=fieldname
+		result['size']=data_size
+		result['data_type']=data_type_str
+		result['data']=data
+
+	except struct.error as e:
+		result=None
+
+	return result
+# }}}
+
+def format_archive_code(format): # {{{
+	"""
+	Given a variable, determine it's type and return
+	an integer value:
+
+	1 : string
+	2 : double (scalar)
+	3 : vector (of type double)
+
+	"""
+	if isinstance(format,basestring):
+		code=1
+	elif format.shape[0] == 1 and format.shape[1] == 1:
+		code=2
+	elif isinstance(format,(list,tuple,numpy.ndarray)):
+		code=3
+	else:
+		raise TypeError("archwrite error: data type '%s' is not valid." % type(format))
+	return code
+# }}}
Index: /issm/trunk-jpl/src/m/archive/archdisp.m
===================================================================
--- /issm/trunk-jpl/src/m/archive/archdisp.m	(revision 21151)
+++ /issm/trunk-jpl/src/m/archive/archdisp.m	(revision 21151)
@@ -0,0 +1,65 @@
+function archdisp(filename) % {{{
+%ARCHDISP - Display the contents of a .arch file
+%
+%	Usage:
+%		archdisp('archive101.arch');
+
+	if nargin~=1
+		help archdisplay
+		error('Error: Invalid number of arguments, only one file can be displayed at a time');
+	end
+
+	if ~exist(filename,'file'),
+		error('Error: File does not exist!');
+	end
+	
+	fid=fopen(filename,'rb');
+	fprintf('Source file: \n');
+	fprintf('\t%s\n',filename);
+	fprintf('Variables: \n');
+
+	field_names={};
+	field_sizes={};
+	field_types={};
+	archive_data=[];
+
+	% Read file
+	while ~feof(fid)
+		[reclen,count]=fread(fid,1,'int','ieee-be');
+		if count == 0, % reached eof
+			break;
+		end
+		% Read variable name
+		fread(fid,1,'int','ieee-be'); % this will always return a string
+		field_name_length=fread(fid,1,'int','ieee-be');
+		field_names{end+1}=char(fread(fid,field_name_length,'char','ieee-be')');
+		% Read variable traits
+		data_reclen=fread(fid,1,'int','ieee-be');
+		field_type=fread(fid,1,'int','ieee-be');
+		% return related type
+		if field_type==2
+			field_types{end+1}='double';
+			field_sizes{end+1}='1x1'; % a scalar will always be a 1x1 cell array
+			archive_data{end+1}=fread(fid,1,'double','ieee-be');
+		elseif field_type==3
+			field_types{end+1}='vector';
+			rows=fread(fid,1,'int','ieee-be');
+			cols=fread(fid,1,'int','ieee-be');
+			data=fread(fid,[rows,cols],'double','ieee-be');
+			field_size_cat=strcat(num2str(rows),'x',num2str(cols)); 
+			field_sizes{end+1}=field_size_cat;
+			archive_data{end+1}=data;
+		else
+			fclose(fid);
+			error('Encountered invalid data type while reading archive.');
+		end
+	end
+	fclose(fid);
+
+	% Display contents
+	for i=1:numel(field_names)
+		fprintf('\t%s\n',field_names{i});
+		fprintf('\t\tSize:\t\t%s\n',field_sizes{i});
+		fprintf('\t\tDatatype:\t%s\n',field_types{i});	
+	end
+end % }}}
Index: /issm/trunk-jpl/src/m/archive/archread.m
===================================================================
--- /issm/trunk-jpl/src/m/archive/archread.m	(revision 21151)
+++ /issm/trunk-jpl/src/m/archive/archread.m	(revision 21151)
@@ -0,0 +1,100 @@
+function [archive_data]=archread(filename,varargin) % {{{
+%ARCHREAD - Given a variable name list of variable names (given as separate arguments),
+%	find and return the data associated with those variable names
+%
+%	Usage:
+%		archive_data = archread('archive101.arch','archive101_field1');
+
+	nvarargs=length(varargin);
+
+	if nargin<2
+		help archread
+		error('Error: Must provide correct number of arguments.');
+	end
+
+	if ~exist(filename,'file')
+		error(['Error: File: ' filename ' does not exist.']);
+	end
+
+	fid=fopen(filename,'rb');
+	% Get structure of file
+	field_names={};
+	field_data_positions=[]; % starting position (in bytes) of data
+	archive_data=[];
+
+	% pass 1: get field names and relate them to field data 
+	while ~feof(fid)
+		[reclen,count]=fread(fid,1,'int','ieee-be'); 
+		if count == 0, % reached eof 
+			break;
+		end
+		record_type=fread(fid,1,'int','ieee-be');
+		if record_type~=1
+			fclose(fid);
+			error(['Error: Reading failed. Expected variable of type string, type given was ' format_archive_code_to_type(record_type)]);
+		else
+			field_name_length=fread(fid,1,'int','ieee-be');
+			field_names{end+1}=char(fread(fid,field_name_length,'char','ieee-be')');
+			% once we have a field name at index i, we will associate
+			% the starting positiion of the data related to that field name
+			% by using field data positions at index i
+			field_data_positions(end+1)=ftell(fid);
+			data_reclen=fread(fid,1,'int','ieee-be');
+			fseek(fid,data_reclen,'cof'); % advance to next field name
+		end
+	end
+	fseek(fid,0,-1); % rewind
+
+	% Determine what fields we need (compare to varargin)
+	idx=find(ismember(field_names,varargin));
+	if size(idx) == 0,
+		fclose(fid);
+		error('Error: No matching variables found in archive.');
+	end
+	% Get data from field_data_positions
+	for i=1:length(idx)
+		fseek(fid,field_data_positions(idx(i)),0); % get starting position of data
+		reclen=fread(fid,1,'int','ieee-be');
+		field_type=fread(fid,1,'int','ieee-be');
+		if field_type==2
+			archive_data{i}=fread(fid,1,'double','ieee-be');
+		elseif field_type==3
+			archive_data{i}=read_vector(fid);
+		else
+			fclose(fid);
+			error('Error: Encountered invalid field type when reading data.');
+		end
+	end
+	fclose(fid);
+end%}}}
+
+%Helper functions{{{
+function [code]=format_archive_code(format) % {{{
+%Given a format, return corresponding code (for reading and writing)
+	if ischar(format)
+		code=1; % string
+	elseif isscalar(format)
+		code=2; % scalar
+	elseif isvector(format)
+		code=3; % vector 
+	else
+		error('Error! Please ensure arguments are strings, scalars, or vectors.');
+	end
+end%}}}
+function [str]=format_archive_code_to_type(code) % {{{
+	if code == 1
+		str='string';
+	elseif code == 2
+		str='scalar';
+	elseif code == 3
+		str='vector';
+	else
+		error(['Code ' num2str(code) ' is not associated with any known format.']);
+	end
+end%}}}
+function [data]=read_vector(fid) % {{{
+	rows=fread(fid,1,'int','ieee-be'); 
+	cols=fread(fid,1,'int','ieee-be');
+	data=fread(fid,[rows,cols],'double','ieee-be');
+end%}}}
+%}}}
Index: /issm/trunk-jpl/src/m/archive/archwrite.m
===================================================================
--- /issm/trunk-jpl/src/m/archive/archwrite.m	(revision 21151)
+++ /issm/trunk-jpl/src/m/archive/archwrite.m	(revision 21151)
@@ -0,0 +1,91 @@
+function archwrite(filename,varargin) % {{{
+%ARCHWRITE - Write data to a field, given the file name, field name, and data.
+%
+%	Usage:
+%		archwrite('archive101.arch','variable_name',data);
+%		archwrite('archive102.arch','archive102_field1',transpose(field_value{1}),...
+%					'archive_field2',transpose(field_value{2}),...
+%					'archive_field3',-6.2420512521312);
+
+	nvarargs=length(varargin);
+
+	if nargin==1 || mod(nvarargs,2)~=0
+		help archwrite
+		error('Error: Invalid number of arguments provided. Please check the usage above.');
+	end
+
+	if ~exist(filename,'file'),
+		fid=fopen(filename,'wb');	% create new file
+	else
+		fid=fopen(filename,'ab+');	% append to file
+	end
+
+	numfields=nvarargs/2;
+	% Generate data to write
+	for i=1:numfields,
+		% write field name
+		name=varargin{2*i-1};
+		write_field_name(fid,name);
+
+		% write data
+		data=varargin{2*i};
+		code=format_archive_code(data);
+		if code==1
+			fclose(fid);
+			error('Error: Field data should be a vector or scalar, string was detected');
+		elseif code==2 % scalar
+			write_scalar(fid,data);
+		elseif code==3 % vector
+			write_vector(fid,data);
+		end
+	end
+	% clean up
+	fclose(fid);
+end%}}}
+
+%Helper functions{{{
+function [code]=format_archive_code(format) % {{{
+	%Given a format, return corresponding code (for reading and writing)
+	if ischar(format)
+		code=1; % string
+	elseif isscalar(format)
+		code=2; % scalar
+	elseif isvector(format)
+		code=3; % vector
+	else
+		error('Error! Please ensure arguments are strings, scalars, or vectors.');
+	end
+end%}}}
+function write_field_name(fid,field_name) % {{{
+	% write record length for name
+	rclen=4+4+length(field_name); % format code + length to write + amt of chars to write
+	fwrite(fid,rclen,'int','ieee-be');
+
+	% write field name to file
+	fwrite(fid,1,'int','ieee-be');							% data will be a string
+	fwrite(fid,length(field_name),'int','ieee-be');	% show how many chars need to be read
+	fwrite(fid,field_name,'char','ieee-be');			% write field name
+end%}}}
+function write_scalar(fid,data) % {{{
+	% write record length for scalar
+	rclen=4+8; % format code + size of double (8 bytes)
+	fwrite(fid,rclen,'int','ieee-be');
+
+	% write scalar to file
+	fwrite(fid,2,'int','ieee-be');			% data will be scalar
+	fwrite(fid,data,'double','ieee-be');	% write scalar to file
+end%}}}
+function write_vector(fid,data) % {{{
+	% write record length for vector
+	row_size=size(data,1);
+	col_size=size(data,2);
+	rclen=4+4+4+8*row_size*col_size; % code+rowsz+colsz+(doublesz*rowamt*colamt)
+	fwrite(fid,rclen,'int','ieee-be');
+	
+	% write vector to file
+	fwrite(fid,3,'int','ieee-be');
+	fwrite(fid,row_size,'int','ieee-be');
+	fwrite(fid,col_size,'int','ieee-be');
+	fwrite(fid,data,'double','ieee-be');
+end%}}}
+%}}}
Index: /issm/trunk-jpl/src/m/archive/netcdf2arch.m
===================================================================
--- /issm/trunk-jpl/src/m/archive/netcdf2arch.m	(revision 21151)
+++ /issm/trunk-jpl/src/m/archive/netcdf2arch.m	(revision 21151)
@@ -0,0 +1,29 @@
+function netcdf2arch(netcdf_filename,arch_filename) % {{{
+%NETCDF2ARCH - Convert a netcdf file into an archive file.
+%
+%	Usage:
+%		netcdf2arch('Archive101.nc','Archive101.arch');
+	ncdata=readnetcdf(netcdf_filename);
+	ncstruct=struct();
+
+	variables = ncdata.VarArray;
+	for i=1:size(variables,2),
+		fieldname=deblank(variables(i).Str);
+		fieldvalue=double(squeeze(variables(i).Data));
+		ncstruct.(fieldname)=fieldvalue;
+	end
+
+	variables=ncdata.AttArray;
+	for i=1:size(variables,2),
+		fieldname=deblank(variables(i).Str);
+		fieldvalue=double(variables(i).Val);
+		ncstruct.(fieldname)=fieldvalue;
+	end
+	ncfields=fieldnames(ncstruct);
+	% first, remove old arch file
+	delete(arch_filename);
+	% write data to file 
+	for i=1:numel(ncfields)
+		archwrite(arch_filename,ncfields{i},ncstruct.(ncfields{i}));
+	end
+end%}}}
Index: /issm/trunk-jpl/src/m/miscellaneous/netcdf2struct.m
===================================================================
--- /issm/trunk-jpl/src/m/miscellaneous/netcdf2struct.m	(revision 21150)
+++ /issm/trunk-jpl/src/m/miscellaneous/netcdf2struct.m	(revision 21151)
@@ -5,25 +5,25 @@
 %      S=netcdf2struct(File);
 
-%Read netcdf file
-data=readnetcdf(File);
+	%Read netcdf file
+	data=readnetcdf(File);
 
-%initialize output
-S=struct();
+	%initialize output
+	S=struct();
 
-%All the variables are in VarArray field
-variables=data.VarArray;
-for i=1:size(variables,2),
-	fieldname=deblank(variables(i).Str);
-	fieldvalue=double(squeeze(variables(i).Data));
-	S.(fieldname)=fieldvalue;
-end
+	%All the variables are in VarArray field
+	variables=data.VarArray;
+	for i=1:size(variables,2),
+		fieldname=deblank(variables(i).Str);
+		fieldvalue=double(squeeze(variables(i).Data));
+		S.(fieldname)=fieldvalue;
+	end
 
-%All the variables are in AttArray field
-variables=data.AttArray;
-for i=1:size(variables,2),
-	fieldname=deblank(variables(i).Str);
-	fieldvalue=double(variables(i).Val);
-	S.(fieldname)=fieldvalue;
-end
+	%All the variables are in AttArray field
+	variables=data.AttArray;
+	for i=1:size(variables,2),
+		fieldname=deblank(variables(i).Str);
+		fieldvalue=double(variables(i).Val);
+		S.(fieldname)=fieldvalue;
+	end
 end
 
@@ -113,21 +113,26 @@
 % Size of NetCDF data type, ID, in bytes
    S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}}));
+end
 
 function T = Type(ID)
 % Matlab string for CDF data type, ID
    T = subsref({'int8','char','int16','int32','single','double'},...
-               struct('type','{}','subs',{{ID}}));
+		struct('type','{}','subs',{{ID}}));
+end
 
 function N = Pad(Num,ID)
 % Number of elements to read after padding to 4 bytes for type ID
    N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0);
+end
 
 function S = String(fp)
 % Read a CDF string; Size,[String,[Padding]]
    S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').';
+end
 
 function A = ReOrder(A,S)
 % Rearrange CDF array A to size S with matlab ordering
    A = permute(reshape(A,fliplr(S)),fliplr(1:length(S)));
+end
 
 function S = DimArray(fp)
@@ -139,4 +144,5 @@
       end
    else fread(fp,1,'uint32=>uint32'); S = []; end
+end
 
 function S = AttArray(fp)
@@ -150,4 +156,5 @@
       end
    else fread(fp,1,'uint32=>uint32'); S = []; end
+end
 
 function S = VarArray(fp)
Index: /issm/trunk-jpl/src/m/miscellaneous/readnetcdf.m
===================================================================
--- /issm/trunk-jpl/src/m/miscellaneous/readnetcdf.m	(revision 21151)
+++ /issm/trunk-jpl/src/m/miscellaneous/readnetcdf.m	(revision 21151)
@@ -0,0 +1,185 @@
+function S = readnetcdf(File,varargin)
+% Function to read NetCDF files
+%   S = netcdf(File)
+% Input Arguments
+%   File = NetCDF file to read
+% Optional Input Arguments:
+%   'Var',Var - Read data for VarArray(Var), default [1:length(S.VarArray)]
+%   'Rec',Rec - Read data for Record(Rec), default [1:S.NumRecs]
+% Output Arguments:
+%   S    = Structure of NetCDF data organised as per NetCDF definition
+% Notes:
+%   Only version 1, classic 32bit, NetCDF files are supported. By default
+% data are extracted into the S.VarArray().Data field for all variables.
+% To read the header only call S = netcdf(File,'Var',[]);
+%
+% SEE ALSO
+% ---------------------------------------------------------------------------
+S = [];
+
+try
+   if exist(File,'file') 
+		fp = fopen(File,'r','b');
+   else 
+		fp = []; 
+		error('File not found'); 
+	end
+   
+	if fp == -1   
+		error('Unable to open file'); 
+	end
+
+% Read header
+   Magic = fread(fp,4,'uint8=>char');
+   if strcmp(Magic(1:3),'CDF') 
+		error('Not a NetCDF file'); 
+	end
+   if uint8(Magic(4))~=1       
+		error('Version not supported'); 
+	end
+  
+	S.NumRecs  = fread(fp,1,'uint32=>uint32');
+   S.DimArray = DimArray(fp);
+   S.AttArray = AttArray(fp);
+   S.VarArray = VarArray(fp);
+
+% Setup indexing to arrays and records
+   Var = ones(1,length(S.VarArray));
+   Rec = ones(1,S.NumRecs);
+   for i = 1:2:length(varargin)
+      if     strcmp(upper(varargin{i}),'VAR') 
+			Var=Var*0; 
+			Var(varargin{i+1})=1;
+      elseif strcmp(upper(varargin{i}),'REC') 
+			Rec=Rec*0; 
+			Rec(varargin{i+1})=1;
+      else 
+			error('Optional input argument not recognised'); 
+		end
+   end
+   
+	if sum(Var)==0 
+		fclose(fp); 
+		return; 
+	end
+
+% Read non-record variables
+   Dim = double(cat(2,S.DimArray.Dim));
+   ID  = double(cat(2,S.VarArray.Type));
+
+   for i = 1:length(S.VarArray)
+      D = Dim(S.VarArray(i).DimID+1); N = prod(D); RecID{i}=find(D==0);
+      if isempty(RecID{i})
+         if length(D)==0 
+				D = [1,1]; 
+				N = 1; 
+			elseif length(D)==1 
+				D=[D,1]; 
+			end
+
+         if Var(i)
+            S.VarArray(i).Data = ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D);
+            fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8');
+         else 
+				fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); 
+			end
+      else 
+			S.VarArray(i).Data = []; 
+		end
+   end
+
+% Read record variables
+   for k = 1:S.NumRecs
+      for i = 1:length(S.VarArray)
+         if ~isempty(RecID{i})
+            D = Dim(S.VarArray(i).DimID+1); D(RecID{i}) = 1; N = prod(D);
+            if length(D)==1 
+					D=[D,1]; 
+				end
+            if Var(i) & Rec(k)
+               S.VarArray(i).Data = cat(RecID{i},S.VarArray(i).Data,...
+                  ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D));
+               if N > 1 
+						fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); 
+					end
+            else 
+					fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); 
+				end
+         end
+      end
+   end
+
+   fclose(fp);
+catch
+   Err = lasterror; fprintf('%s\n',Err.message);
+   if ~isempty(fp) && fp ~= -1 
+		fclose(fp); 
+	end
+end
+end
+
+% ---------------------------------------------------------------------------------------
+% Utility functions
+
+function S = Size(ID)
+% Size of NetCDF data type, ID, in bytes
+   S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}}));
+end
+
+function T = Type(ID)
+% Matlab string for CDF data type, ID
+   T = subsref({'int8','char','int16','int32','single','double'},...
+		struct('type','{}','subs',{{ID}}));
+end
+
+function N = Pad(Num,ID)
+% Number of elements to read after padding to 4 bytes for type ID
+   N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0);
+end
+
+function S = String(fp)
+% Read a CDF string; Size,[String,[Padding]]
+   S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').';
+end
+
+function A = ReOrder(A,S)
+% Rearrange CDF array A to size S with matlab ordering
+   A = permute(reshape(A,fliplr(S)),fliplr(1:length(S)));
+end
+
+function S = DimArray(fp)
+% Read DimArray into structure
+   if fread(fp,1,'uint32=>uint32') == 10 % NC_DIMENSION
+      for i = 1:fread(fp,1,'uint32=>uint32')
+         S(i).Str = String(fp);
+         S(i).Dim = fread(fp,1,'uint32=>uint32');
+      end
+   else fread(fp,1,'uint32=>uint32'); S = []; end
+end
+
+function S = AttArray(fp)
+% Read AttArray into structure
+   if fread(fp,1,'uint32=>uint32') == 12 % NC_ATTRIBUTE
+      for i = 1:fread(fp,1,'uint32=>uint32')
+         S(i).Str = String(fp);
+         ID       = fread(fp,1,'uint32=>uint32');
+         Num      = fread(fp,1,'uint32=>uint32');
+         S(i).Val = fread(fp,Pad(Num,ID),[Type(ID),'=>',Type(ID)]).';
+      end
+   else fread(fp,1,'uint32=>uint32'); S = []; end
+end
+
+function S = VarArray(fp)
+% Read VarArray into structure
+   if fread(fp,1,'uint32=>uint32') == 11 % NC_VARIABLE
+      for i = 1:fread(fp,1,'uint32=>uint32')
+         S(i).Str      = String(fp);
+         Num           = double(fread(fp,1,'uint32=>uint32'));
+         S(i).DimID    = double(fread(fp,Num,'uint32=>uint32'));
+         S(i).AttArray = AttArray(fp);
+         S(i).Type     = fread(fp,1,'uint32=>uint32');
+         S(i).VSize    = fread(fp,1,'uint32=>uint32');
+         S(i).Begin    = fread(fp,1,'uint32=>uint32'); % Classic 32 bit format only
+      end
+   else fread(fp,1,'uint32=>uint32'); S = []; end
+end
Index: /issm/trunk-jpl/test/NightlyRun/README
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/README	(revision 21150)
+++ /issm/trunk-jpl/test/NightlyRun/README	(revision 21151)
@@ -34,3 +34,3 @@
 To modify some characteristics, do it in the testxxx.m file.
 Specify the field_names the field_values and the tolerance at the end of the test file.
-Don't forget to commit the archive (Archivexxx.nc in ../Archives)
+Don't forget to commit the archive (Archivexxx.arch in ../Archives)
Index: /issm/trunk-jpl/test/NightlyRun/runme.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/runme.m	(revision 21150)
+++ /issm/trunk-jpl/test/NightlyRun/runme.m	(revision 21151)
@@ -140,14 +140,11 @@
 		archive_name=['Archive' num2str(id) ];
 		if strcmpi(procedure,'update'),
-			delete(['../Archives/' archive_name '.nc'])
+			delete(['../Archives/' archive_name '.arch'])
 			for k=1:length(field_names),
 				field=field_values{k};
 				% matlab writes the dimensions reversed and matrices transposed into netcdf, so compensate for that
-				nccreate(['../Archives/' archive_name '.nc'],[archive_name '_field' num2str(k)],...
-				         'Dimensions',{[archive_name '_field' num2str(k) '_2'] size(field,2) [archive_name '_field' num2str(k) '_1'] size(field,1)},...
-				         'Format','classic');
-				ncwrite(['../Archives/' archive_name '.nc'],[archive_name '_field' num2str(k)],transpose(field));
+				archwrite(['../Archives/' archive_name '.arch'],[archive_name '_field' num2str(k)], transpose(field));
 			end
-			disp(sprintf(['File ./../Archives/' archive_name '.nc saved\n']));
+			disp(sprintf(['File ./../Archives/' archive_name '.arch saved\n']));
 
 		%CHECK for memory leaks?
@@ -215,7 +212,8 @@
 
 					%compare to archive
-					% matlab reads the dimensions reversed and matrices transposed from netcdf, so compensate for that
-					archive=transpose(ncread(['../Archives/' archive_name '.nc'],[archive_name '_field' num2str(k)]));
-					error_diff=full(max(abs(archive(:)-field(:)))/(max(abs(archive))+eps));
+					%our output is in the correct order (n,1) or (1,1), so we do not need to transpose again
+					archive_cell=archread(['../Archives/' archive_name '.arch'],[archive_name '_field' num2str(k)]);
+					archive=archive_cell{1};
+					error_diff=full(max(abs(archive(:)-field(:)))/(max(abs(archive(:)))+eps));
 
 					%disp test result
Index: /issm/trunk-jpl/test/NightlyRun/runme.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/runme.py	(revision 21150)
+++ /issm/trunk-jpl/test/NightlyRun/runme.py	(revision 21151)
@@ -4,5 +4,4 @@
 import socket
 import numpy
-import netCDF4
 import sys
 import traceback
@@ -47,4 +46,7 @@
 	from MatlabFuncs import strcmpi
 	from MatlabFuncs import ismember
+	from arch import archread
+	from arch import archwrite
+	from arch import archdisp
 
 	#Get ISSM_DIR variable
@@ -140,8 +142,7 @@
 			archive_name='Archive'+str(id)
 			if strcmpi(procedure,'update'):
-				archive=os.path.join('..','Archives',archive_name+'.nc')
-				if os.path.isfile(archive):
-					os.remove(archive)
-				f = netCDF4.Dataset(archive,'w',format='NETCDF3_CLASSIC')
+				archive_file=os.path.join('..','Archives',archive_name+'.arch')
+				if os.path.isfile(archive_file):
+					os.remove(archive_file)
 				for k,fieldname in enumerate(field_names):
 					field=numpy.array(field_values[k],dtype=float)
@@ -154,10 +155,6 @@
 						field=field.reshape(1,1)
 					# Matlab uses base 1, so use base 1 in labels
-					f.createDimension(archive_name+'_field'+str(k+1)+'_1',numpy.size(field,0))
-					f.createDimension(archive_name+'_field'+str(k+1)+'_2',numpy.size(field,1))
-					v = f.createVariable(archive_name+'_field'+str(k+1),field.dtype,(archive_name+'_field'+str(k+1)+'_1',archive_name+'_field'+str(k+1)+'_2'))
-					v[:] = field
-				f.close()
-				print "File '%s' saved.\n" % os.path.join('..','Archives',archive_name+'.nc')
+					archwrite(archive_file,archive_name+'_field'+str(k+1),field)
+				print "File '%s' saved.\n" % os.path.join('..','Archives',archive_name+'.arch')
 
 			#ELSE: CHECK TEST
@@ -165,8 +162,8 @@
 
 				#load archive
-				if os.path.exists(os.path.join('..','Archives',archive_name+'.nc')):
-					f = netCDF4.Dataset(os.path.join('..','Archives',archive_name+'.nc'),'r')
+				if os.path.exists(os.path.join('..','Archives',archive_name+'.arch')):
+					archive_file=os.path.join('..','Archives',archive_name+'.arch')
 				else:
-					raise IOError("Archive file '"+os.path.join('..','Archives',archive_name+'.nc')+"' does not exist.")
+					raise IOError("Archive file '"+os.path.join('..','Archives',archive_name+'.arch')+"' does not exist.")
 
 				for k,fieldname in enumerate(field_names):
@@ -184,7 +181,6 @@
 						#compare to archive
 						# Matlab uses base 1, so use base 1 in labels
-						if archive_name+'_field'+str(k+1) in f.variables:
-							archive=f.variables[archive_name+'_field'+str(k+1)][:]
-						else:
+						archive=numpy.array(archread(archive_file,archive_name+'_field'+str(k+1)))
+						if archive == None:
 							raise NameError("Field name '"+archive_name+'_field'+str(k+1)+"' does not exist in archive file.")
 						error_diff=numpy.amax(numpy.abs(archive-field),axis=0)/ \
@@ -214,5 +210,4 @@
 							raise RuntimeError(message)
 
-				f.close()
 
 		except Exception as message: