function femmodel=diagnostic_core(femmodel);
%DIAGNOSTIC_CORE - compute the core velocity field 
%
%   Usage:
%      results=diagnostic_core(model);
%

	%some parameters
	modify_loads=true;
	conserve_loads=true;

	%recover parameters common to all solutions
	dim=femmodel.parameters.MeshDimension;
	ishutter=femmodel.parameters.FlowequationIshutter;
	ismacayealpattyn=femmodel.parameters.FlowequationIsmacayealpattyn;
	isstokes=femmodel.parameters.FlowequationIsstokes;
	isnewton=femmodel.parameters.DiagnosticIsnewton;
	dakota_analysis=femmodel.parameters.QmuIsdakota;
	control_analysis=femmodel.parameters.InversionIscontrol;
	solution_type=femmodel.parameters.SolutionType;

	%for qmu analysis, be sure the velocity input we are starting from  is the one in the parameters: 
	if dakota_analysis & solution_type==DiagnosticSolutionEnum,
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVxEnum,VxEnum);
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVyEnum,VyEnum);
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVzEnum,VzEnum);
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuPressureEnum,PressureEnum);
	end

	%Compute slopes: 
	if(ishutter), femmodel=surfaceslope_core(femmodel); end
	if(isstokes),
		femmodel=bedslope_core(femmodel); 
		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
		femmodel.nodes=ResetCoordinateSystem(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
	end

	if ishutter,

		issmprintf(VerboseSolution,'\n%s',['   computing hutter velocities']);

		%Take the last velocity into account so that the velocity on the MacAyeal domain is not zero
		if(ismacayealpattyn)femmodel=ResetBoundaryConditions(femmodel,DiagnosticHutterAnalysisEnum); end 

		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHutterAnalysisEnum);
		femmodel=solver_linear(femmodel);

		if(ismacayealpattyn)femmodel=ResetBoundaryConditions(femmodel,DiagnosticHorizAnalysisEnum); end

	end
			
	if xor(ismacayealpattyn,isstokes), %if macayealpattyn and stokes, we must use stokescoupling

		issmprintf(VerboseSolution,'\n%s',['   computing horizontal velocities']);
		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
		if isnewton,
			femmodel=solver_newton(femmodel); 
		else
			femmodel=solver_nonlinear(femmodel,modify_loads); 
		end
	end
	
	if ismacayealpattyn & isstokes,

		issmprintf(VerboseSolution,'\n%s',['   computing coupling velocities for macayealpattyn and stokes']);
		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
		femmodel=solver_stokescoupling_nonlinear(femmodel,conserve_loads); 
	end

	if dim==3 & (ismacayealpattyn | ishutter) & ~isstokes,
	
		issmprintf(VerboseSolution,'\n%s',['   computing vertical velocities']);
		femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
		femmodel=solver_linear(femmodel);
	end

	if (solution_type==DiagnosticSolutionEnum & ~control_analysis)
		issmprintf(VerboseSolution,'\n%s',['   saving results']);
		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum);
		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum);
		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VelEnum);
		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,PressureEnum);
		if(dim==3) femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum); end
	end
