function md=diagnostic(md);
%DIAGNOSTIC - compute the velocity field of a model
%
%   Usage:
%      md=diagnostic(md)
%
	%timing
	t1=clock;

	analyses=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopecomputeAnalysisEnum];

	displaystring(md.verbose,'%s',['create fem model']);
	femmodel=CreateFemModel(md,analyses);

	%compute solution
	if ~femmodel.parameters.QmuAnalysis,
		if femodel.parameters.control_analysis,
			
			%launch core of control solution.
			md.results.DiagnosticAnalysis=control_core(femmodel);

		else,
			
			%launch core of diagnostic solution.
			md.results.DiagnosticAnalysis=diagnostic_core(femmodel);

		end
	else
		%launch dakota driver for diagnostic core solution
		Qmu(femmodel);
	end

	%stop timing
	t2=clock;
	displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);