function u_g=icediagnostic2d(md,fem),
%ICEDIAGNOSTIC2D - diagnostic of a 2d model
%
%   this routine computes the velocity field of a glacier in 2d 
%
%   Usage:
%      u_g=icediagnostic2d(md,fem)
%
%   See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR

%recover inputs
inputs=fem.inputs;

%Compute surface slope if there are hutter elements
if fem.ishutter,

	%Get field of fem
	m_ss=fem.m_ss;
	m_dhu=fem.m_dhu;

	%Compute slope of the surface
	if fem.m_dhu.params.debug
		disp(sprintf('\n%s',['computing surface slope (x and y derivatives)...']));
	end
	slopex=icediagnostic_core_linear(m_ss,'surface_slope_compute_x',inputs);
	slopey=icediagnostic_core_linear(m_ss,'surface_slope_compute_y',inputs);
	
	%initialize slope surface input for hutter solution sequence
	slopesurface=zeros(m_ss.gridset.gsize,1);
	slopesurface(1:6:m_ss.gridset.gsize,1)=slopex(1:6:end);
	slopesurface(2:6:m_ss.gridset.gsize,1)=slopey(1:6:end);

	%Compute Hutter solution
	disp(sprintf('\n%s',['computing hutter velocities...']));
	inputs.slopesurface=slopesurface;
	u_g=icediagnostic_core_linear(m_dhu,'diagnostic_hutter',inputs);

	%Computation of depth averaged pressure with MacAyeal's assumptions (P_bar=1/2*rho_ice*g*H in Pa)
	u_g(4:6:m_dhu.gridset.gsize)=1/2*md.rho_ice*md.g*md.thickness;
end

if fem.ismacayealpattyn,

	%Get field of fem
	m_dh=fem.m_dh;

	%Compute horizontal solution. 
	disp(sprintf('\n%s',['computing horizontal velocities...']));

	%Run core solution
	u_g=icediagnostic_core_nonlinear(m_dh,'diagnostic_horiz',inputs);

	%Computation of depth averaged pressure with MacAyeal's assumptions (P_bar=1/2*rho_ice*g*H in Pa)
	u_g(4:6:m_dh.gridset.gsize)=1/2*md.rho_ice*md.g*md.thickness;
end
