function bool=isresultconsistent(md,solutiontype),
%ISRESULTCONSISTENT: check that results are consistent
%
%   This routine is called after a computation and check
%   that the values t the penalties are consitent
%
%   Usage:
%      bool=IsModelSelfConsistent(model,solutiontype)

%tolerance we use in litmus tests for the consistency of the model
tolerance=10^-2;

if (nargin~=2  )
	IsResultConsistentUsage();
	error(' ');
end
%Check velocities i(no penalties for 2d meshes)
if (strcmp(solutiontype,'diagnostic') || strcmp(solutiontype,'diagnostic_horiz') ),

	if strcmpi(md.type,'3d'),
		if ~isnan(md.penalties),
			for i=1:size(md.penalties,1)
				for n=2:md.numlayers
					relativex=(md.vx(md.penalties(i,1))-md.vx(md.penalties(i,n)))./(md.vx(md.penalties(i,1)));
					if md.vx(md.penalties(i,1))==md.vx(md.penalties(i,n)), relativex=0; end
					relativey=(md.vy(md.penalties(i,1))-md.vy(md.penalties(i,n)))./(md.vy(md.penalties(i,1)));
					if md.vy(md.penalties(i,1))==md.vy(md.penalties(i,n)), relativey=0; end
					if abs(relativex)>tolerance | abs(relativey)>tolerance
						disp(['''diagnostic'' result not consistent: increase penalty_offset (discontinuous horizontal velocity)']);
						bool=0; return;
					end

				end
			end
		end
	end
end

if strcmp(solutiontype,'thermalsteady')

	%check melting 
	if any(abs(md.melting(md.numberofgrids2d+1:end))>tolerance)
		disp(['''thermalsteady'' result not consistent: increase penalty_melting (negative melting)']);
		bool=0; return; 
	end

end

if strcmp(solutiontype,'thermaltransient')

	for iter=1:length(md.thermaltransient_results)

	%check melting 
	if any(abs(md.thermaltransient_results(iter).melting(md.numberofgrids2d+1:end))>tolerance)
		disp(['''thermaltransient'' result not consistent: increase penalty_melting (negative melting)']);
		bool=0; return;
	end

	end
end


if (strcmp(solutiontype,'transient') & strcmpi(md.type,'3d')),
	for iter=1:length(md.transient_results)

		%check melting 
		if any(abs(md.transient_results(iter).melting(md.numberofgrids2d+1:end))>tolerance)
			disp(['''thermaltransient'' result not consistent: increase penalty_melting (negative melting)']);
			bool=0; return;
		end


		%check velocities
		if ~isnan(md.penalties),
			for i=1:size(md.penalties,1)
				for n=2:md.numlayers
					relativex=(md.transient_results(iter).vx(md.penalties(i,1))-md.transient_results(iter).vx(md.penalties(i,n)))./(md.transient_results(iter).vx(md.penalties(i,1)));
					if md.transient_results(iter).vx(md.penalties(i,1))==md.transient_results(iter).vx(md.penalties(i,n)), relativex=0; end
					relativey=(md.transient_results(iter).vy(md.penalties(i,1))-md.transient_results(iter).vy(md.penalties(i,n)))./(md.transient_results(iter).vy(md.penalties(i,1)));
					if md.transient_results(iter).vy(md.penalties(i,1))==md.transient_results(iter).vy(md.penalties(i,n)), relativey=0; end
					if abs(relativex)>tolerance | abs(relativey)>tolerance
						disp(['''transient'' result not consistent: increase penalty_offset (discontinuous horizontal velocity)']);
						bool=0; return;
					end

				end
			end
		end

	end

end


%No problems, just return 1;
bool=1;
return;
	
function IsResultConsistentUsage()
disp(' ');
disp('   IsResultConsistent usage: flag=IsResultConsistent(model,solutiontype)');
disp('         where model is an instance of the @model class, and flag a boolean');
