function md=iceparameters(md);
%ICEPARAMETERS - compute parameters
%
%   This routine is caled by the solver to compute the parameters
%   that are given in md.outputparameters (stress, strain rate,...)
%
%   Usage:
%      md=iceparameters(md)

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

if strcmpi(md.type,'2d'),
	%First, build elements,grids,loads, etc ... for thermal and melting 
	m_dh=CreateFemModel(md,'diagnostic_horiz');

	%Are any more computations requested? 
	if ~isempty(md.parameteroutput),
		disp(sprintf('\n%s',['computing outputs...']));
		%velocity field
		u_g=zeros(gridset.gsize,1); 
		u_g(1:6:gridset.gsize)=md.vx/md.yts; %from m/a to m/s
		u_g(2:6:gridset.gsize)=md.vy/md.yts;

		%Go through fields of parameteroutput and compute
		for i=1:length(md.parameteroutput),
			if strcmpi(md.parameteroutput{i},'viscousheating'),
				disp(sprintf('%s',['   viscousheating']));
				md.viscousheating=ViscousHeatingCompute(m_dh,struct('velocity',u_g),md.type);
			elseif strcmpi(md.parameteroutput{i},'pressure_elem'),
				disp(sprintf('%s',['   pressure_elem']));
				md.pressure_elem=PressureElemCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
			elseif strcmpi(md.parameteroutput{i},'strainrate'),
				disp(sprintf('%s',['   strainrate']));
				md.strainrate=StrainRateCompute(m_dh,struct('velocity',u_g),md.type);
			elseif strcmpi(md.parameteroutput{i},'stress'),
				disp(sprintf('%s',['   stress']));
				md.stress=StressCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
			elseif strcmpi(md.parameteroutput{i},'stress_bed'),
				disp(sprintf('%s',['   stress_bed']));
				md.stress_bed=StressBedCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
			elseif strcmpi(md.parameteroutput{i},'stress_surface'),
				disp(sprintf('%s',['   stress_surface']));
				md.stress_surface=StressSurfaceCompute(m_dh,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
			elseif strcmpi(md.parameteroutput{i},'deviatoricstress'),
				disp(sprintf('%s',['   deviatoricstress']));
				md.deviatoricstress=DeviatoricStressCompute(m_dh,struct('velocity',u_g),md.type);
			end
		end
	end

else

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

	%Are any more computations requested? 
	if ~isempty(md.parameteroutput),
		disp(sprintf('\n%s',['computing outputs...']));
		%velocity field
		u_g=zeros(gridset.gsize,1); 
		u_g(1:6:gridset.gsize)=md.vx/md.yts; %from m/a to m/s
		u_g(2:6:gridset.gsize)=md.vy/md.yts;
		u_g(3:6:gridset.gsize)=md.vz/md.yts;
		if ~isnan(md.pressure),
			u_g(4:6:gridset.gsize)=md.pressure;
		end

		%Go through fields of parameteroutput and compute
		for i=1:length(md.parameteroutput),
			if strcmpi(md.parameteroutput{i},'viscousheating'),
				disp(sprintf('%s',['   viscousheating']));
				md.viscousheating=ViscousHeatingCompute(m_dv,struct('velocity',u_g),md.type);
			elseif strcmpi(md.parameteroutput{i},'pressure_elem'),
				disp(sprintf('%s',['   pressure_elem']));
				md.pressure_elem=PressureElemCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
			elseif strcmpi(md.parameteroutput{i},'strainrate'),
				disp(sprintf('%s',['   strainrate']));
				md.strainrate=StrainRateCompute(m_dv,struct('velocity',u_g),md.type);
			elseif strcmpi(md.parameteroutput{i},'stress'),
				disp(sprintf('%s',['   stress']));
				md.stress=StressCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
			elseif strcmpi(md.parameteroutput{i},'stress_bed'),
				disp(sprintf('%s',['   stress_bed']));
				md.stress_bed=StressBedCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
			elseif strcmpi(md.parameteroutput{i},'stress_surface'),
				disp(sprintf('%s',['   stress_surface']));
				md.stress_surface=StressSurfaceCompute(m_dv,struct('velocity',u_g,'elementonstokes',md.elements_type(:,2)),md.type);
			elseif strcmpi(md.parameteroutput{i},'deviatoricstress'),
				disp(sprintf('%s',['   deviatoricstress']));
				md.deviatoricstress=DeviatoricStressCompute(m_dv,struct('velocity',u_g),md.type);
			end
		end
	end

end
