Index: /issm/trunk/src/m/solutions/cielo/proto/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/proto/diagnostic.m	(revision 275)
+++ /issm/trunk/src/m/solutions/cielo/proto/diagnostic.m	(revision 275)
@@ -0,0 +1,33 @@
+function md=diagnostic_proto(md);
+%DIAGNOSTIC - compute the velocity field of a model
+%   Usage:
+%      md=diagnostic(md)
+%
+
+	%timing
+	t1=clock;
+
+	%Build all models requested for diagnostic simulation
+	md.analysis_type='diagnostic_horiz'; m_dh=CreateFemModel(md);
+	md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md);
+	md.analysis_type='diagnostic_stokes'; m_ds=CreateFemModel(md);
+	md.analysis_type='diagnostic_hutter'; m_dhu=CreateFemModel(md);
+	md.analysis_type='surface_slope_compute'; m_ss=CreateFemModel(md);
+	md.analysis_type='bed_slope_compute'; m_bs=CreateFemModel(md);
+	
+	% figure out number of dof: just for information purposes.
+	md.dof=m_dh.nodesets.fsize; %biggest dof number
+
+	%initialize inputs
+	inputs=inputlist;
+	inputs=add(inputs,'velocity',m_dh.parameters.u_g,'doublevec',3,m_dh.parameters.numberofnodes);
+
+	%compute solution
+	[u_g,pg]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_ss,m_bs,inputs);
+
+	%Load results onto model
+	md=loadresults(md,u_g,p_g,m_dh,m_ds,m_dhu);
+
+	%stop timing
+	t2=clock;
+	disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']));
Index: /issm/trunk/src/m/solutions/cielo/proto/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/proto/diagnostic_core.m	(revision 275)
+++ /issm/trunk/src/m/solutions/cielo/proto/diagnostic_core.m	(revision 275)
@@ -0,0 +1,96 @@
+function [u_g p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_ss,m_bs,inputs);
+%DIAGNOSTIC_CORE - compute the core velocity field 
+%
+%   Usage:
+%      u_g=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_ss,m_bs,inputs);
+%
+
+%recover parameters common to all solutions
+debug=m_dhu.parameters.debug;
+
+%Compute surface slope if there are hutter elements
+if m_dhu.parameters.ishutter,
+
+	display('\n%s',['computing surface slope (x and y derivatives)...'],debug);
+	
+	slopex=diagnostic_core_linear(m_ss,'surface_slope_compute_x',inputs);
+	slopey=diagnostic_core_linear(m_ss,'surface_slope_compute_y',inputs);
+
+	%Extrude slopes in 3d 
+	if m_dhu.parameters.dim==3,
+		slopex=SlopeExtrude(m_ss.elements,m_ss.nodes,m_ss.loads,m_ss.materials,slopex);
+		slopey=SlopeExtrude(m_ss.elements,m_ss.nodes,m_ss.loads,m_ss.materials,slopey);
+	end
+
+	%initialize slope surface input for hutter solution sequence
+	inputs=add(inputs,'surfaceslopex',slopex,'doublevec',1,m_ss.parameters.numberofnodes);
+	inputs=add(inputs,'surfaceslopey',slopey,'doublevec',1,m_ss.parameters.numberofnodes);
+
+	display('\n%s',['computing hutter velocities...']);
+	u_g=diagnostic_core_linear(m_dhu,'diagnostic_hutter',inputs);
+
+	%Computation of depth averaged pressure with MacAyeal's assumptions (P_bar=1/2*rho_ice*g*H in Pa)
+	p_g=Pressure(m_dhu.elements,m_dhu.nodes,m_dhu.materials,m_dhu.parameters);
+
+	%update dirichlet boundary conditions for macyeal pattyn, with the velocities computed by hutter
+	if m_dh.ismacayealpattyn,
+		m_dh.y_g=u_g;
+		[m_dh.ys m_dh.ys0]=Reducevectorgtos(m_dh.y_g,m_dh.nodesets);
+	end
+
+end
+
+if m_dh.ismacayealpattyn,
+
+	if debug, disp(sprintf('\n%s',['computing horizontal velocities...']));end;
+	u_g=diagnostic_core_nonlinear(m_dh,'diagnostic_horiz',inputs);
+
+	%Computation of depth averaged pressure with MacAyeal's assumptions (P_bar=1/2*rho_ice*g*H in Pa)
+	p_g=Pressure(m_dh.elements,m_dh.nodes,m_dh.materials,m_dh.parameters);
+end
+	
+if m_dh.parameters.dim==3,,
+
+	if debug, disp(sprintf('\n%s',['extruding horizontal velocities...']));end;
+	u_g=VelocityExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g);
+
+	if debug, disp(sprintf('\n%s',['computing vertical velocities...']));end;
+	inputs=add(inputs,'velocity',u_g,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
+	u_g_vert=diagnostic_core_linear(m_dv,inputs);
+
+	%Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
+	p_g=Pressure(m_dh.elements,m_dh.nodes,m_dh.materials,m_dh.parameters);
+	
+	if m_ds.parameters.isstokes,
+
+		%recondition pressure 
+		p_g=p_g/m_ds.parameters.stokesreconditioning;
+
+		%Compute slope of the bed
+		if debug, disp(sprintf('\n%s',['computing bed slope (x and y derivatives)...'])); end;
+		
+		slopex=diagnostic_core_linear(m_bs,'bed_slope_compute_x',inputs);
+		slopey=diagnostic_core_linear(m_bs,'bed_slope_compute_y',inputs);
+
+		slopex=SlopeExtrude(m_bs.elements,m_bs.nodes,m_bs.loads,m_bs.materials,slopex);
+		slopey=SlopeExtrude(m_bs.elements,m_bs.nodes,m_bs.loads,m_bs.materials,slopey);
+
+		inputs=add(inputs,'bedslopex',slopex,'doublevec',m_ss.parameters.numberofdofspernode,m_ss.parameters.numberofnodes);
+		inputs=add(inputs,'bedslopey',slopey,'doublevec',m_ss.parameters.numberofdofspernode,m_ss.parameters.numberofnodes);
+		inputs=add(inputs,'velocity',u_g,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
+		inputs=add(inputs,'pressure',p_g,'doublevec',1,m_dh.parameters.numberofnodes);
+
+		%update dirichlet boundary conditions with the velocities computed previously
+		m_ds.y_g(dofsetgen([1,2],4,m_ds.parameters.gsize))=u_g;
+		m_ds.y_g(dofsetgen([3],4,m_ds.parameters.gsize))=u_g_vert;
+		
+		[m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets);
+
+		%compute stokes velocities and pressure
+		if debug, disp(sprintf('\n%s',['computing stokes velocities and pressure ...']));end;
+		u_g=diagnostic_core_nonlinear(m_ds,'diagnostic_stokes',inputs);
+	
+		%decondition pressure
+		p_g=u_g(4:4:end)*m_dh.parameters.stokesreconditioning;
+	end
+end