Index: /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/DomainOutline.exp
===================================================================
--- /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/DomainOutline.exp	(revision 178)
+++ /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/DomainOutline.exp	(revision 178)
@@ -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/IceShelfIceFrontP3d_5/Front.exp
===================================================================
--- /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/Front.exp	(revision 178)
+++ /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/Front.exp	(revision 178)
@@ -0,0 +1,10 @@
+## Name:icefront
+## Icon:0
+# Points Count  Value
+5 1.
+# X pos Y pos
+-1000 900000
+-1000 1100000
+1100000 1100000
+1100000 900000
+-1000 900000
Index: /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/Square.par
===================================================================
--- /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/Square.par	(revision 178)
+++ /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/Square.par	(revision 178)
@@ -0,0 +1,39 @@
+%Start defining model parameters here
+
+%dynamics
+md.dt=1*md.yts; %1 year
+md.ndt=md.dt*10; 
+md.artificial_diffusivity=1;
+
+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=SetMarineIceSheetBC(md,'Front.exp');
+
+%Parallel options
+md.np=8;
+md.time=50;
+md.waitonlock=1;
Index: /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/runme.m
===================================================================
--- /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/runme.m	(revision 178)
+++ /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/runme.m	(revision 178)
@@ -0,0 +1,94 @@
+% This file can be run to check that the current version of issm is giving 
+% coherent results. This test deals with an icesheet with icefront for a 3d 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^-12. 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'};
+solutions={'diagnostic','thermalsteady','thermaltransient','prognostic','transient'};
+
+%Initialize log message for nightly runs.
+testname='IceShelfIceFrontM3d_2';
+tolerance=10^-12;
+
+%go through the solutions requested
+for i=1:length(packages),
+	package=packages{i};
+
+	for j=1:length(solutions),
+		solution=solutions{j};
+
+		%check package
+		if ~(strcmpi(package,'macayeal') | strcmpi(package,'ice') | strcmpi(package,'cielo_serial') | strcmpi(package,'cielo_parallel'));
+			error(['package: ' package  ' in test: ' testname  ' not supported yet']);
+		%check solution
+		elseif ~(strcmpi(solution,'diagnostic') | strcmpi(solution,'thermalsteady') | strcmpi(solution,'thermaltransient') |...
+				strcmpi(solution,'prognostic') | strcmpi(solution,'transient'));
+			error(['solution: ' solution  ' in test: ' testname  ' not supported yet']);
+		end
+		%check solution requested
+		if (~strcmpi(package,'ice') | strcmpi(solution,'transient')),
+			disp(sprintf(['\nsolution: ' solution  ', with package: ' package  ', in test: ' testname  ', not supported yet.\n']));
+			continue
+		end
+
+		%initialize model
+		md=model;
+		md=mesh(md,'DomainOutline.exp',100000);
+		md=geography(md,'all','');
+		md=parameterize(md,'Square.par');
+		md=extrude(md,8,4);
+		md=setelementstype(md,'pattyn','all');
+		if strcmpi(package,'cielo_parallel'), md.cluster='wilkes'; end
+		if md.numberofgrids==832
+			load Velocities; md.vx=0.8*vx; md.vy=0.8*vy; md.vz=0.8*vz; md.temperature=temperature-1; md.pressure=pressure;
+		end
+
+		%compute solution
+		if strcmpi(package,'cielo_parallel') & strcmpi(solution,'diagnostic'),
+			md=solve(md,'diagnostic_horiz','cielo');
+		elseif strcmpi(package,'cielo_serial'),
+			eval(['md=solve(md,''' solution ''',''cielo'');']);
+		else
+			eval(['md=solve(md,''' solution ''',''' package ''');']);
+		end
+
+		%compute field to be checked
+		if strcmpi(solution,'diagnostic'),
+			fields={'vy','vz'};
+		elseif strcmpi(solution,'thermalsteady'),
+			fields={'temperature','melting'};
+		elseif strcmpi(solution,'thermaltransient'),
+			fields={'thermaltransient_results(end).temperature','thermaltransient_results(end).melting'};
+		elseif strcmpi(solution,'prognostic'),
+			fields={'new_thickness'};
+		elseif strcmpi(solution,'transient'),
+			fields={'transient_results(end).vy','transient_results(end).vz','transient_results(end).temperature','transient_results(end).melting','transient_results(end).thickness'};
+		end
+
+		%load archive
+		eval(['load Archive' package solution ]);
+
+		for k=1:length(fields),
+			field=fields{k};
+
+			%compare to archive
+			eval(['Archive=Archive' package solution '_field' num2str(k) ';']);
+			eval(['error_diff=abs(norm((Archive(find(Archive))-md.' field  '(find(Archive)))./Archive(find(Archive)),2));']);
+
+			%disp test result
+			if (error_diff>tolerance);
+				disp(sprintf(['\n\nERROR (difference = %8.3g > %g)   --->   test: ' testname ', solution: ' solution  ', package: ' package ', field tested: ' field  '.\n\n'],error_diff,tolerance));
+			else
+				disp(sprintf(['\n\nSUCCESS (difference = %8.3g < %g)   --->   test: ' testname ', solution: ' solution  ', package: ' package ', field tested: ' field  '.\n\n'],error_diff,tolerance));
+			end
+
+		end
+	end
+end
Index: /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/updatearchive.m
===================================================================
--- /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/updatearchive.m	(revision 178)
+++ /issm/trunk/test/Verification/IceShelfIceFrontP3d_5/updatearchive.m	(revision 178)
@@ -0,0 +1,75 @@
+% This file can be run to update the velocity archives  of the test1.
+% This test deals with an icesheet with icefront for a 3d 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'};
+solutions={'diagnostic','thermalsteady','thermaltransient','prognostic','transient'};
+
+%go through the solutions requested
+testname='IceShelfIceFrontM3d_2';
+for i=1:length(packages),
+	package=packages{i};
+
+	for j=1:length(solutions),
+		solution=solutions{j};
+
+		%check package
+		if ~(strcmpi(package,'macayeal') | strcmpi(package,'ice') | strcmpi(package,'cielo_serial') | strcmpi(package,'cielo_parallel'));
+			error(['package: ' package  ' in test: ' testname  ' not supported yet']);
+			%check solution
+		elseif ~(strcmpi(solution,'diagnostic') | strcmpi(solution,'thermalsteady') | strcmpi(solution,'thermaltransient') |...
+				strcmpi(solution,'prognostic') | strcmpi(solution,'transient'));
+			error(['solution: ' solution  ' in test: ' testname  ' not supported yet']);
+		end
+		%check solution requested
+		if (~strcmpi(package,'ice') | strcmpi(solution,'transient')),
+			disp(sprintf(['\nsolution: ' solution  ', with package: ' package  ', in test: ' testname  ', not supported yet.\n']));
+			continue
+		end
+
+		%initialize model
+		md=model;
+		md=mesh(md,'DomainOutline.exp',100000);
+		md=geography(md,'all','');
+		md=parameterize(md,'Square.par');
+		md=extrude(md,8,4);
+		md=setelementstype(md,'pattyn','all');
+		if strcmpi(package,'cielo_parallel'), md.cluster='wilkes'; end
+		if md.numberofgrids==832
+			load Velocities; md.vx=0.8*vx; md.vy=0.8*vy; md.vz=0.8*vz; md.temperature=temperature-1; md.pressure=pressure;
+		end
+
+		%compute solution
+		if strcmpi(package,'cielo_parallel') & strcmpi(solution,'diagnostic'),
+			md=solve(md,'diagnostic_horiz','cielo');
+		elseif strcmpi(package,'cielo_serial'),
+			eval(['md=solve(md,''' solution ''',''cielo'');']);
+		else
+			eval(['md=solve(md,''' solution ''',''' package ''');']);
+		end
+
+		%compute field to be checked
+		if strcmpi(solution,'diagnostic'),
+			fields={'vy','vz'};
+		elseif strcmpi(solution,'thermalsteady'),
+			fields={'temperature','melting'};
+		elseif strcmpi(solution,'thermaltransient'),
+			fields={'thermaltransient_results(end).temperature','thermaltransient_results(end).melting'};
+		elseif strcmpi(solution,'prognostic'),
+			fields={'new_thickness'};
+		elseif strcmpi(solution,'transient'),
+			fields={'transient_results(end).vy','transient_results(end).vz','transient_results(end).temperature','transient_results(end).melting','transient_results(end).thickness'};
+		end
+
+		%save new archive
+		for k=1:length(fields),
+			field=fields{k};
+			eval(['Archive' package solution '_field' num2str(k) '=md.' field  ';']);
+		end
+		eval(['save Archive' package solution ' Archive' package solution '_field*']);
+	end
+end