function varargout=runme(varargin)
%RUNME - test deck for ISSM nightly runs
%
%   Usage:
%      varargout=runme(varargin);
%
%   Examples:
%      runme;
%      runme({'ice'});
%      runme({'ice','cielo_serial'},{'diagnostic'});
%      md=runme({'cielo_parallel'});
%
%   See Also: UPDATEARCHIVE

% 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.

%check arguments
if (nargin>2 | nargout>1)
	help runme
	error('runme error message: bad usage');
end

%packages and analysis_types to be tested
if nargin==2,
	solutions=varargin{2};
else
	solutions={'diagnostic','thermal_steady','thermal_transient','prognostic','transient'};
end
if nargin,
	packages=varargin{1};
else
	packages={'macayeal','ice','cielo_serial'};
end

%Initialize log message for nightly runs.
testname='IceShelfIceFrontMPS3d_4';

%go through the solutions requested
for i=1:length(packages),
	package=packages{i};

	for j=1:length(solutions),
		solution=solutions{j};

		%check solution requested
		if ((~strcmpi(package,'ice') & (strcmpi(solution,'transient') | strcmpi(solution,'diagnostic')))...
				| strcmpi(package,'macayeal'));
			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','Pattyn.exp','fill','macayeal','stokes','Stokes.exp');
		if md.numberofgrids==728
			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
		[md analysis_type sub_analysis_type packagefinal]=testsparameterization(md,solution,package);
		if ~testschecksolution(md.type,solution), continue; end
		md=solve(md,'analysis_type',analysis_type,'sub_analysis_type',sub_analysis_type,'package',packagefinal);

		%compute fields to be checked
		fields=testsgetfields(md,solution);

		%load archive
		eval(['load Archive' package solution ]);

		for k=1:length(fields),
			
			%Get field and tolerance
			field=fields{k};
			tolerance=testsgettolerance(md,package,solution,field);

			%compare to archive
			eval(['Archive=Archive' package solution '_field' num2str(k) ';']);
			eval(['error_diff=full(max(abs(Archive-md.results.' field '))/(max(abs(Archive))+eps));']);

			%disp test result
			if (error_diff>tolerance);
				disp(sprintf(['\n\nERROR   difference: %-7.2g > %7.2g test: %-25s solution: %-16s package: %-14s field: ' field  '.\n\n'],error_diff,tolerance,testname,solution,package));
			else
				disp(sprintf(['\n\nSUCCESS difference: %-7.2g < %7.2g test: %-25s solution: %-16s package: %-14s field: ' field  '.\n\n'],error_diff,tolerance,testname,solution,package));
			end

		end
	end
end
if nargout==1, 
	varargout{1}=md; 
end