function md=cieloprognostic(md)
%diagnostic solution for horizontal velocity

	%timing
	t1=clock;

	analysis_p='prognostic';
	
	%configure, and initialise data model and parameter defaults:
	SetParameterDefaults;

	%Build model requested for prognostic simulation
	m=CreateFEMModel(md,analysis_p);
		
	% figure out number of dof: just for information purposes.
	md.dof=m.uset.fsize; 

	%Build partitioning vectors. We only need one partitioning for all our models, we made sure of it in Imp!
	indx=1:6:m.uset.gsize; indx=indx(tpart);

	%build some parameters needed in core solution.
	m.rifts=md.rifts;

	%compute depth averaged horizontal velocity in Matlab workspace partitioning
	u_g=zeros(m.uset.gsize,1); u_g(1:6:end)=md.vx/md.yts; u_g(2:6:end)=md.vy/md.yts;
	velocity_average=CieloHorizontalVelocityDepthAverage(md,u_g);

	%move surface and bed velocities to first lower layer, in workspace partitioning
	ws=ShiftLayers(md,md.vz/md.yts,md.numlayers,1);
	wb=ShiftLayers(md,md.vz/md.yts,1,1);
	melting=ShiftLayers(md,md.melting,1,1);
	accumulation=ShiftLayers(md,md.accumulation,md.numlayers,1);

	%plug into g set, with cluster partitioning
	ws_g=zeros(m.uset.gsize,1);ws_g(1:6:m.uset.gsize)=ws(part);
	wb_g=zeros(m.uset.gsize,1);wb_g(1:6:m.uset.gsize)=wb(part);
	melting_g=zeros(m.uset.gsize,1);melting_g(1:6:m.uset.gsize)=melting(part);
	accumulation_g=zeros(m.uset.gsize,1);accumulation_g(1:6:m.uset.gsize)=accumulation(part);
	
	%initialize thickness
	thickness=zeros(m.uset.gsize,1);thickness(1:6:end)=md.thickness(part);

	%prepare inputs strcuture
	inputs=struct('thickness',thickness,'melting',melting_g,'accumulation',accumulation_g,...
	'surface_vertical_velocity',ws_g,'basal_vertical_velocity',wb_g,'dt',md.dt,'velocity_average',velocity_average);

	%Run core solution
	h_g=cieloprognostic_core(m,m.params,inputs,analysis_p);

	%recover thickness and project onto every grid of the mesh
	md.new_thickness=project3d(md,project2d(md,h_g(indx),1),'node');
