function md=iceprognostic(md)
%ICEPROGNOSTIC - prognostic solution
%
%   This routine is used to compute the evolution of the geometry
%   of a model
%
%   Usage:
%      md=iceprognostic(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('iceprognostic error message: parallel support not implemented yet');
end

%setup some parameters to be passed to the core solution
params.sparsity=md.sparsity;
params.solver_type=md.solver_type;
params.eps_rel=md.eps_rel;
params.eps_abs=md.eps_abs;
params.debug=md.debug;

%First, build elements,grids,loads, etc ... for prognostic
m=CreateFemModel(md,'prognostic');

%compute depth averaged horizontal velocity
u_g=zeros(gridset.gsize,1);
u_g(1:6:end)=md.vx/md.yts; %from m/yr to m/s
u_g(2:6:end)=md.vy/md.yts; %from m/yr to m/s
u_g(3:6:end)=md.vz/md.yts; %from m/yr to m/s

if strcmpi(md.type,'3d'),
	velocity_average=HorizontalVelocityDepthAverage(md,u_g);
else
	%the velocity is already depthaverdged
	velocity_average=u_g;
end

%move surface and bed velocities to first lower layer. This will be needed by the prognostic model
melting=zeros(gridset.gsize,1);
accumulation=zeros(gridset.gsize,1);
if strcmpi(md.type,'3d'),
	melting(1:6:end)=ShiftLayers(md,md.melting,1,1);
	accumulation(1:6:end)=ShiftLayers(md,md.accumulation,md.numlayers,1);
else
	melting(1:6:end)=md.melting/md.yts;           %from m/yr to m/s
	accumulation(1:6:end)=md.accumulation/md.yts; %from m/yr to m/s
end

%prepare inputs
thickness=zeros(gridset.gsize,1);thickness(1:6:end)=md.thickness;
%inputs=struct('thickness',thickness,'melting',melting,'accumulation',accumulation,...
%'surface_vertical_velocity',ws,'basal_vertical_velocity',wb,'dt',md.dt,'velocity_average',velocity_average);
inputs=struct('thickness',thickness,'melting',melting,'accumulation',accumulation,'dt',md.dt,'velocity_average',velocity_average);

%Run core solution
h_g=iceprognostic_core(m,'prognostic',inputs);

%project collapsed thickness onto 3d mesh
indx=1:6:m.gridset.gsize; indx=indx(m.tpart);
if strcmpi(md.type,'3d'),
	md.new_thickness=project3d(md,project2d(md,h_g(indx),1),'node');
else
	md.new_thickness=h_g(indx);
end
