function md=slopecompute(md,field)
%SLOPECOMPUTE - compute the slope of the bed or the surface on the grids
%
%   Usage:
%      md=slopecompute(md,field)

% Define global variables
iceglobal
global gridset

if ~strcmpi(field, 'bed') & ~strcmpi(field, 'surface')
	error('slope compute error message: bed and surface computing only ');
end

%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('slopecompute error message: parallel support not implemented yet');
end

%First, build elements,grids,loads, etc ... for horizontal, base vertical and vertical model
if strcmpi(field, 'bed'),
	m_bs=CreateFemModel(md,'bed_slope_compute');
elseif strcmpi(field, 'surface'),
	m_ss=CreateFemModel(md,'surface_slope_compute');
end

if strcmpi(field, 'bed'),
	slopex=icediagnostic_core_linear(m_bs,'bed_slope_compute_x');
	slopey=icediagnostic_core_linear(m_bs,'bed_slope_compute_y');
	
	slopex=slopex(1:6:end);
	slopey=slopey(1:6:end);

	if strcmpi(md.type,'3d'),
		slopex=project3d(md,project2d(md,slopex,1),'node');
		slopey=project3d(md,project2d(md,slopey,1),'node');
	end

	md.bed_slopex=slopex;
	md.bed_slopey=slopey;
	
elseif strcmpi(field, 'surface'),
	slopex=icediagnostic_core_linear(m_ss,'surface_slope_compute_x');
	slopey=icediagnostic_core_linear(m_ss,'surface_slope_compute_y');
		
	slopex=slopex(1:6:end);
	slopey=slopey(1:6:end);

	if strcmpi(md.type,'3d'),
		slopex=project3d(md,project2d(md,slopex,md.numlayers),'node');
		slopey=project3d(md,project2d(md,slopey,md.numlayers),'node');
	end

	md.surface_slopex=slopex;
	md.surface_slopey=slopey;
end
	
