Index: /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/DomainOutline.exp
===================================================================
--- /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/DomainOutline.exp	(revision 116)
+++ /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/DomainOutline.exp	(revision 116)
@@ -0,0 +1,10 @@
+## Name:domainoutline
+## Icon:0
+# Points Count  Value
+5 1.
+# X pos Y pos
+0 0
+1000000 0
+1000000 1000000
+0 1000000
+0 0
Index: /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/Square.par
===================================================================
--- /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/Square.par	(revision 116)
+++ /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/Square.par	(revision 116)
@@ -0,0 +1,36 @@
+%Start defining model parameters here
+
+md.viscosity_overshoot=0.3;
+
+disp('      creating thickness');
+hmin=300;
+hmax=1000;
+ymin=min(md.y);
+ymax=max(md.y);
+md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
+md.bed=-md.rho_ice/md.rho_water*md.thickness;
+md.surface=md.bed+md.thickness;
+
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
+
+disp('      creating temperature');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
+
+disp('      creating flow law paramter');
+md.B=paterson(md.observed_temperature);
+md.n=3*ones(md.numberofelements,1);
+
+%Deal with boundary conditions:
+md=SetIceSheetBC(md);
+
+%Parallel options
+md.np=8;
+md.time=50;
+md.waitonlock=1;
Index: /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/runme.m
===================================================================
--- /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/runme.m	(revision 116)
+++ /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/runme.m	(revision 116)
@@ -0,0 +1,87 @@
+% This file can be run to check that the current version of macayeal and issm models are giving 
+% coherent results. This test deals with an icesheet without icefront for a 2d model. The geometry 
+% is square. Just run this file in Matlab, with a properly setup ISSM code. 
+% The results of this test will indicate if there is a difference between current computations 
+% and archived results.
+
+% Errors  between archived results and the current version will get flagged if they are not within 
+% a certain tolerance. The current tolerance is 10^-14. If you have good reasons to believe this 
+% tolerance should be lowered (for example, if you are running single precision compilers?), feel 
+% free to tweak the tolerance variable.
+
+%packages and solutions to be tested
+packages={'macayeal','ice','cielo_serial','cielo_parallel'};
+
+%Initialize log message for nightly runs.
+testname='IceSheetNoIceFront2d_1';
+tolerance=10^-12;
+
+%go through the solutions requested
+for i=1:length(packages),
+	package=packages{i};
+
+	%initialize model
+	md=model;
+	md=mesh(md,'DomainOutline.exp',50000);
+	md=geography(md,'','');
+	md=parameterize(md,'Square.par');
+	md=setelementstype(md,'macayeal','all');
+
+	%Compute solution using requested package
+	if strcmpi(package,'macayeal'),
+
+		%DIAGNOSTIC
+		md=solve(md,'diagnostic','macayeal');
+		%check result
+		load ArchiveMacAyealDiag
+		pos=find(ArchiveMacAyealDiag); error_vel=abs(norm((ArchiveMacAyealDiag(pos)-md.vel(pos))./ArchiveMacAyealDiag(pos),2));
+		if (error_vel>tolerance);
+			disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> ERROR (difference=' num2str(error_vel) ')'])
+		else
+			disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> SUCCESS (difference=' num2str(error_vel) ')'])
+		end
+
+	elseif strcmpi(package,'ice'),
+
+		%DIAGNOSTIC
+		md=solve(md,'diagnostic','ice');
+		%check result
+		load ArchiveIceDiag
+		pos=find(ArchiveIceDiag); error_vel=abs(norm((ArchiveIceDiag(pos)-md.vel(pos))./ArchiveIceDiag(pos),2));
+		if (error_vel>tolerance);
+			disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> ERROR (difference=' num2str(error_vel) ')'])
+		else
+			disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> SUCCESS (difference=' num2str(error_vel) ')'])
+		end
+
+	elseif strcmpi(package,'cielo_serial'),
+
+		%DIAGNOSTIC
+		md=solve(md,'diagnostic','cielo');
+		%check result
+		load ArchiveCieloSerialDiag
+		pos=find(ArchiveCieloSerialDiag); error_vel=abs(norm((ArchiveCieloSerialDiag(pos)-md.vel(pos))./ArchiveCieloSerialDiag(pos),2));
+		if (error_vel>tolerance);
+			disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> ERROR (difference=' num2str(error_vel) ')'])
+		else
+			disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> SUCCESS (difference=' num2str(error_vel) ')'])
+		end
+
+	elseif strcmpi(package,'cielo_parallel'),
+
+		%DIAGNOSTIC
+		md.cluster='wilkes';
+		md=solve(md,'diagnostic_horiz','cielo');
+		%check result
+		load ArchiveCieloParallelDiag
+		pos=find(ArchiveCieloParallelDiag); error_vel=abs(norm((ArchiveCieloParallelDiag(pos)-md.vel(pos))./ArchiveCieloParallelDiag(pos),2));
+		if (error_vel>tolerance);
+			disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> ERROR (difference=' num2str(error_vel) ')'])
+		else
+			disp(['test: ' testname  ', solution: diagnostic, package: ' package ' -> SUCCESS (difference=' num2str(error_vel) ')'])
+		end
+
+	else
+		disp(sprintf(['\npackage ' package  ' not supported yet.\n']));
+	end
+end
Index: /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/updatearchive.m
===================================================================
--- /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/updatearchive.m	(revision 116)
+++ /issm/trunk/test/Verification/IceSheetNoIceFront2d_1/updatearchive.m	(revision 116)
@@ -0,0 +1,59 @@
+% This file can be run to update the velocity archives  of the test1.
+% This test deals with an icesheet without icefront for a 2d model. The geometry 
+% is square. Just run this file in Matlab, with a properly setup Ice code. 
+
+% The archive files will be saved in this directory but will not commited to ice1. 
+% Just commit the result if you want to.
+
+%packages and solutions to be tested
+packages={'macayeal','ice','cielo_serial','cielo_parallel'};
+
+%go through the solutions requested
+for i=1:length(packages),
+	package=packages{i};
+
+	%initialize model
+	md=model;
+	md=mesh(md,'DomainOutline.exp',50000);
+	md=geography(md,'','');
+	md=parameterize(md,'Square.par');
+	md=setelementstype(md,'macayeal','all');
+
+	%Compute solution using requested package
+	if strcmpi(package,'macayeal'),
+
+		%DIAGNOSTIC
+		md=solve(md,'diagnostic','macayeal');
+		%save archive
+		ArchiveMacAyealDiag=md.vel;
+		save ArchiveMacAyealDiag ArchiveMacAyealDiag
+
+	elseif strcmpi(package,'ice'),
+
+		%DIAGNOSTIC
+		md=solve(md,'diagnostic','ice');
+		%save archive
+		ArchiveIceDiag=md.vel;
+		save ArchiveIceDiag ArchiveIceDiag
+
+	elseif strcmpi(package,'cielo_serial'),
+
+		%DIAGNOSTIC
+		md=solve(md,'diagnostic','cielo');
+		%save archive
+		ArchiveCieloSerialDiag=md.vel;
+		save ArchiveCieloSerialDiag ArchiveCieloSerialDiag
+
+	elseif strcmpi(package,'cielo_parallel'),
+
+		%DIAGNOSTIC
+		md.cluster='wilkes';
+		md=solve(md,'diagnostic_horiz','cielo');
+		%save archive
+		ArchiveCieloParallelDiag=md.vel;
+		save ArchiveCieloParallelDiag ArchiveCieloParallelDiag
+
+	else
+		disp(sprintf(['\npackage ' package  ' not supported yet.\n']));
+	end
+end