function md=icediagnostic(md);
%ICEDIAGNOSTIC - compute the velocity field of a model
%
%   this routine is a Matlab drivn solution
%   
%   Usage:
%      u_g=icediagnostic_core_linear(m,analysis_type,varargin);
%
%   See also: ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR

%define global variables
iceglobal

%determine if run is parallel
if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end;

%for now, only serial support is in
if cluster,
	error('icediagnostic error message: parallel support not implemented yet');
end

%Create fem structure (input of icediagnostic3d)
fem=struct();
%Figure out which type of elements are present
[fem.ishutter,fem.ismacayealpattyn,fem.isstokes]=DiagnosticSolutionType(md.elements_type);

if strcmpi(md.type,'2d'),

	%First, build elements,grids,loads, etc ... for horizontal model
	if fem.ishutter,
		fem.m_ss=CreateFemModel(md,'surface_slope_compute');
		fem.m_dhu=CreateFemModel(md,'diagnostic_hutter');
	end
	if fem.ismacayealpattyn,
		fem.m_dh=CreateFemModel(md,'diagnostic_horiz');
	end

	%plug inputs in fem (existing velocity)
	if ~isnan(md.vx) & ~isnan(md.vy) & fem.ismacayealpattyn,
		%an input velocity is present in the model, use it to bootstrap the diagnostic core nonlinear solution.
		m_dh=fem.m_dh;
		fem.inputs.velocity=zeros(m_dh.gridset.gsize,1);
		fem.inputs.velocity(1:6:m_dh.gridset.gsize)=md.vx/md.yts;
		fem.inputs.velocity(2:6:m_dh.gridset.gsize)=md.vy/md.yts;
	else
		fem.inputs=struct();
	end

	%compute solution
	u_g=icediagnostic2d(md,fem);

	%Load results onto model
	md=Loadresults(md,fem,u_g);

else

	%First, build elements,grids,loads, etc ... for horizontal, base vertical and vertical model
	fem.m_dbv=CreateFemModel(md,'diagnostic_basevert');
	fem.m_dv=CreateFemModel(md,'diagnostic_vert');

	if fem.ismacayealpattyn,
		fem.m_dh=CreateFemModel(md,'diagnostic_horiz');
	end
	if fem.ishutter,
		fem.m_ss=CreateFemModel(md,'surface_slope_compute');
		fem.m_dhu=CreateFemModel(md,'diagnostic_hutter');
	end
	if fem.isstokes,
		fem.m_bs=CreateFemModel(md,'bed_slope_compute');
		fem.m_ds=CreateFemModel(md,'diagnostic_stokes');
	end

	%plug inputs in fem (existing velocity)
	if ~isnan(md.vx) & ~isnan(md.vy) & fem.ismacayealpattyn,
		%an input velocity is present in the model, use it to bootstrap the diagnostic core nonlinear solution.
		m_dh=fem.m_dh;
		fem.inputs.velocity=zeros(m_dh.gridset.gsize,1);
		fem.inputs.velocity(1:6:m_dh.gridset.gsize)=md.vx/md.yts;
		fem.inputs.velocity(2:6:m_dh.gridset.gsize)=md.vy/md.yts;
	else
		fem.inputs=struct();
	end

	%compute solution
	u_g=icediagnostic3d(md,fem);

	%Load results onto model
	md=Loadresults(md,fem,u_g);
end
