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

global gridset

%recover inputs
inputs=fem.inputs;

%Compute surface slope if there are hutter elements
if fem.ishutter,
	m_ss=fem.m_ss;
	m_dhu=fem.m_dhu;

	%Compute slope of the surface
	if m_ss.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);

	%Project slopes computed at surface, onto the 3d volume
	slopex=project3d(md,project2d(md,slopex(1:6:end),md.numlayers),'node');
	slopey=project3d(md,project2d(md,slopey(1:6:end),md.numlayers),'node');

	%Feed hutter solution with the slope input
	slopesurface=zeros(m_ss.gridset.gsize,1);
	slopesurface(1:6:m_ss.gridset.gsize,1)=slopex;
	slopesurface(2:6:m_ss.gridset.gsize,1)=slopey;
	inputs.slopesurface=slopesurface;

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

if fem.ismacayealpattyn,
	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);
end

%Extrude velocities for collapsed penta elements
u_g=VelocityExtrude(md,u_g);

%Compute depth averaged velocity and add it to inputs
gridset=fem.m_dv.gridset;
velocity_average=HorizontalVelocityDepthAverage(md,u_g);
inputs.velocity=u_g;
inputs.velocity_average=velocity_average;

%Compute basal vertical velocity, in 3d
disp(sprintf('\n%s',['computing basal vertical velocity...']));
u_g_basevert=icediagnostic_core_linear(fem.m_dbv,'diagnostic_basevert',inputs);

%update dirichlet boundary conditions with these new base vertical velocities
m_dv=fem.m_dv;
gridset=m_dv.gridset;
m_dv.y_g=u_g_basevert; m_dv.ys=Reducevector_g(m_dv.y_g);

%compute vertical velocities
disp(sprintf('\n%s',['computing vertical velocity...']));
u_g_vert=icediagnostic_core_linear(m_dv,'diagnostic_vert',inputs);

%add contribution to vertical velocity
u_g=u_g+u_g_vert;

%Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
u_g(4:6:m_dv.gridset.gsize)=md.rho_ice*md.g*(md.surface-md.z)/md.stokesreconditioning;

if fem.isstokes,
	m_bs=fem.m_bs;
	m_ds=fem.m_ds;

	%Compute slope of the bed
	if m_bs.params.debug
		disp(sprintf('\n%s',['computing bed slope (x and y derivatives)...']));
	end
	slopex=icediagnostic_core_linear(m_bs,'bed_slope_compute_x',inputs);
	slopey=icediagnostic_core_linear(m_bs,'bed_slope_compute_y',inputs);

	slopex=project3d(md,project2d(md,slopex(1:6:end),1),'node');
	slopey=project3d(md,project2d(md,slopey(1:6:end),1),'node');

	slopebed(1:6:m_bs.gridset.gsize,1)=slopex;
	slopebed(2:6:m_bs.gridset.gsize,1)=slopey;
	inputs.slopebed=slopebed;
	inputs.velocity=u_g;

	%update dirichlet boundary conditions with the velocities computed previously
	gridset=m_ds.gridset;
	m_ds.y_g=u_g; m_ds.ys=Reducevector_g(m_ds.y_g);

	%compute stokes velocities and pressure
	disp(sprintf('\n%s',['computing stokes velocities and pressure ...']));
	u_g=icediagnostic_core_nonlinear(m_ds,'diagnostic_stokes',inputs);
end

%Reconditioning of the pressure
u_g(4:6:m_dv.gridset.gsize)=u_g(4:6:m_dv.gridset.gsize)*md.stokesreconditioning;
