function md=solve(md,solutiontype,varargin)
%SOLVE - apply solution sequence for this model
%
%   solutiontype is 'diagnostic','prognostic','transient','thermalsteady','thermaltransient','parameters' or 'control'
%   and varargin is an optional package name ('cielo', 'ice' or 'macayeal')
%
%   Usage:
%       md=solve(md,solutiontype,varargin)
%
%   Examples:
%      md=solve(md,'diagnostic','macayeal');
%      md=solve(md,'control','cielo');

%some checks on list of arguments

solutions={'mesh2grid','mesh','thermaltransient','thermalsteady','qmu','diagnostic','diagnostic_horiz','prognostic','transient','paramters','control'};
found=0;
for i=1:length(solutions),
	if strcmpi(solutiontype,solutions{i}),
		found=1;
		break;
	end
end

if found==0,
	error(['solve error messae: solution type ' solutiontype ' not supported yet!']);
end

%we are good with this solutiontype, put in the analysis_type field of md: 
md.analysis_type=solutiontype;

%Recover type of package being used: 
if nargin==2,
	package='Ice';
else
	package=varargin{1};
end

if ~ischar(package), 
	error('Package specified in varargin can only be ''ice'', or ''cielo''');
end

if ~(strcmpi(package,'ice') || strcmpi(package,'cielo') || strcmpi(package,'macayeal'))
	error('Package specified in varargin can only be ''ice'', ''macayeal'', or ''cielo''');
end

%Use package to set solution namespace
usenamespace(package);


%Check model is self-consistent
if ~strcmpi(solutiontype,'mesh2grid'),
	disp(sprintf('\n%s','checking model consistency'));
	if ~ismodelselfconsistent(md,solutiontype,package),
		error(' '); %previous error messages should explain what is going on.
	end
end
disp('launching solution sequence');

%If running in parallel, we have a different way of launching the solution
%sequences.
if ~strcmpi(solutiontype,'qmu'),
	if ~strcmpi(md.cluster,'none'),
		md=solveparallel(md,solutiontype,package);
		return;
	end
end

%Lauch correct solution sequence
if strcmpi(solutiontype,'diagnostic'),
	md=diagnostic(md);

elseif strcmpi(solutiontype,'mesh'),
	md=mesh(md);

elseif strcmpi(solutiontype,'transient'),
	md=transient(md);

elseif strcmpi(solutiontype,'qmu'),
	md=qmu(md,package);

elseif strcmpi(solutiontype,'parameters'),
	md=parameters(md);;

elseif strcmpi(solutiontype,'mesh2grid'),
	md=mesh2grid(md);;

elseif strcmpi(solutiontype,'prognostic'),
	md=prognostic(md);;

elseif strcmpi(solutiontype,'control'),
	md=control(md);

elseif strcmpi(solutiontype,'thermalsteady') | strcmpi(solutiontype,'thermaltransient'),
	md=thermal(md,solutiontype);
else
	error('solution type not supported for this model configuration.');
end

%Check result is consistent
disp(sprintf('%s\n','checking result consistency'));
if ~isresultconsistent(md,solutiontype),
	disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
end

function solveusage();
disp(' ');
disp('   Solve usage: md=solve(md,solutiontype,package)');
disp('         solutiontype can be ''diagnostic'',''transient'', ''thermalsteady'',''thermaltransient'',''parameters'',''mesh2grid''  or ''control''');
disp('         package is either ''cielo'' or ''ice''');
