function md=diagnostic(md)

	%timing
	t1=clock;

	analysis_dh='diagnostic_horiz';
	analysis_dbv='diagnostic_basevert';
	analysis_dv='diagnostic_vert';

	%Build all models requested for diagnostic simulation
	%remark, partitions are all identical.
	md.analysis_type=analysis_dh; 
	m_dh=CreateFemModel(md);

	if strcmpi(md.type,'3d'),
		m_dbv=CreateFemModel(md,analysis_dbv);
		m_dv=CreateFemModel(md,analysis_dv);
	end

	% figure out number of dof: just for information purposes.
	md.dof=m_dh.nodesets.fsize; %biggest dof number

	%Get horizontal solution. 
	disp(sprintf('\n%s',['computing horizontal velocities...']));

	u_g=diagnostic_core_nonlinear(m_dh,{});

	if strcmpi(md.type,'3d'),
	
		%swith u_g to workspace partitioning.
		u_g=SwitchPartitioning(u_g,'workspace',part,tpart,[1 2]); 
		
		%extrude velocities for collapsed penta elements
		u_g=CieloVelocityExtrude(md,u_g);

		%Compute depth averaged velocity and add it to inputs
		velocity_average=CieloHorizontalVelocityDepthAverage(md,u_g);

		%swith u_g and velocity_average to cluster partitioning
		u_g=SwitchPartitioning(u_g,'cluster',part,tpart,[1 2]); 
		velocity_average=SwitchPartitioning(velocity_average,'cluster',part,tpart,[1 2]); 
		

		disp(sprintf('\n%s',['computing basal vertical velocities...']));
		SetUset(m_dbv);
		u_g_basevert=diagnostic_core_linear(m_dbv,m_dbv.parameters,struct('velocity',u_g,'velocity_average',velocity_average),analysis_dbv); 
		
		%use u_g_basevert to constrain vertical velocity
		SetUset(m_dv);
		m_dv.y_g=u_g_basevert;

		%reduce constraining vector to s-set (the one we solve on)
		m_dv.y_s = Reducevectorg( m_dv.y_g);

		%run core linear to solve for vertical velocity
		disp(sprintf('\n%s',['computing vertical velocities...']));
		u_g_vert=diagnostic_core_linear(m_dv,m_dv.parameters,struct('velocity',u_g),analysis_dv); 
		
		%load results onto model: carefull, u_g and u_g_vert are in cluster partitioning
		md.vx=u_g(indx)*md.yts;
		md.vy=u_g(indy)*md.yts;
		md.vz=u_g_vert(indz)*md.yts;
		md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
		
		%Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
		md.pressure=md.rho_ice*md.g*(md.surface-md.z);
	else
		
		%Build partitioning vectors to recover solution
		indx=m_dh.part(1:2:end);
		indy=m_dh.part(2:2:end);

		%load results onto model
		md.vx=u_g(indx)*md.yts;
		md.vy=u_g(indy)*md.yts;
		md.vz=zeros(md.numberofgrids,1);
		md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
	
		%Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
		md.pressure=md.rho_ice*md.g*(md.thickness);
	end		

	%timing
	t2=clock;
	disp(sprintf('\n%s\n',['   solution converged in ' num2str(etime(t2,t1)) ' seconds']));
