Index: /issm/trunk/src/m/solutions/cielo/CreateFemModel.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/CreateFemModel.m	(revision 301)
+++ /issm/trunk/src/m/solutions/cielo/CreateFemModel.m	(revision 302)
@@ -8,32 +8,32 @@
 function  m=CreateFEMModel(md)
 
-	disp(sprintf('\n   reading data from model %s...',md.name));
+	displaystring(md.debug,'\n   reading data from model %s...',md.name);
 	[m.elements,m.nodes,m.constraints,m.loads,m.materials,parameters]=ModelProcessor(md);
 	
-	disp('   generating degrees of freedom...');
+	displaystring(md.debug,'%s','   generating degrees of freedom...');
 	[m.nodes,m.part,m.tpart]=Dof(m.elements,m.nodes,parameters);
 
-	disp('   generating single point constraints...');
+	displaystring(md.debug,'%s','   generating single point constraints...');
 	[m.nodes,m.yg]=SpcNodes(m.nodes,m.constraints);
 
-	disp('   generating rigid body constraints...');
+	displaystring(md.debug,'%s','   generating rigid body constraints...');
 	[m.Rmg,m.nodes]=MpcNodes(m.nodes,m.constraints);
 
-	disp('   generating node sets...');
+	displaystring(md.debug,'%s','   generating node sets...');
 	m.nodesets=BuildNodeSets(m.nodes);
 
-	disp('   reducing single point constraints vector...');
+	displaystring(md.debug,'%s','   reducing single point constraints vector...');
 	[m.ys m.ys0]=Reducevectorgtos(m.yg,m.nodesets);
 
-	disp('   normalizing rigid body constraints matrix...');
+	displaystring(md.debug,'%s','   normalizing rigid body constraints matrix...');
 	m.Gmn = NormalizeConstraints(m.Rmg,m.nodesets);
 	
-	disp('   configuring element and loads...');
+	displaystring(md.debug,'%s','   configuring element and loads...');
 	[m.elements,m.loads,m.nodes] = ConfigureObjects( m.elements, m.loads, m.nodes, m.materials);
 
-	disp('   processing parameters...');
+	displaystring(md.debug,'%s','   processing parameters...');
 	parameters= ProcessParams(parameters,m.part);
 
-	disp('   creating parameters in matlab workspace...');
+	displaystring(md.debug,'%s','   creating parameters in matlab workspace...');
 	m.parameters= ParamsToWorkspace(parameters);
 
Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 301)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 302)
@@ -1,7 +1,7 @@
-function md=diagnostic(md)
-%DIAGNOSTIC - diagnostic solution sequence
-%
+function md=diagnostic_proto(md);
+%DIAGNOSTIC - compute the velocity field of a model
 %   Usage:
 %      md=diagnostic(md)
+%
 
 	%timing
@@ -9,10 +9,20 @@
 
 	%Build all models requested for diagnostic simulation
+	displaystring(md.debug,'%s',['reading diagnostic horiz model data']);
 	md.analysis_type='diagnostic_horiz'; m_dh=CreateFemModel(md);
+	
+	displaystring(md.debug,'\n%s',['reading diagnostic vert model data']);
+	md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md);
+	
+	displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']);
+	md.analysis_type='diagnostic_stokes'; m_ds=CreateFemModel(md);
 
-	if strcmpi(md.type,'3d'),
-		md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md);
-	end
-
+	displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']);
+	md.analysis_type='diagnostic_hutter'; m_dhu=CreateFemModel(md);
+	
+	displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']);
+	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
@@ -22,41 +32,11 @@
 	inputs=add(inputs,'velocity',m_dh.parameters.u_g,'doublevec',3,m_dh.parameters.numberofnodes);
 
-	%Get horizontal solution. 
-	disp(sprintf('\n%s',['computing horizontal velocities...']));
+	%compute solution
+	[u_g,p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_ss,m_bs,inputs);
 
-	u_g=diagnostic_core_nonlinear(m_dh,inputs);
+	%Load results onto model
+	md=loadresults(md,u_g,p_g,m_dh,m_ds,m_dhu);
 
-	if strcmpi(md.type,'3d'),
-	
-		%extrude velocities for collapsed penta elements
-		disp(sprintf('\n%s',['extruding horizontal velocities...']));
-		u_g=VelocityExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g);
-
-		disp(sprintf('\n%s',['computing vertical velocities...']));
-	
-		inputs=add(inputs,'velocity',u_g,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
-
-		u_g_vert=diagnostic_core_linear(m_dv,inputs);
-
-		%load results onto model: 
-		md.vx=u_g(1:2:end)*md.yts;
-		md.vy=u_g(2:2:end)*md.yts;
-		md.vz=u_g_vert*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
-		%load results onto model
-		md.vx=u_g(1:2:end)*md.yts;
-		md.vy=u_g(2:2:end)*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
+	%stop timing
 	t2=clock;
-	disp(sprintf('\n%s\n',['   solution converged in ' num2str(etime(t2,t1)) ' seconds']));
+	displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);
Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 302)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 302)
@@ -0,0 +1,99 @@
+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;
+dim=m_dhu.parameters.dim;
+ishutter=m_dhu.parameters.ishutter;
+ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
+isstokes=m_ds.parameters.isstokes;
+
+if ishutter,
+
+	displaystring(debug,'\n%s',['computing surface slope (x and y derivatives)...']);
+	slopex=diagnostic_core_linear(m_ss,'surface_slope_compute_x',inputs);
+	slopey=diagnostic_core_linear(m_ss,'surface_slope_compute_y',inputs);
+
+	if dim==3,
+		displaystring(debug,'\n%s',['extruding slopes in 3d...']);
+		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
+
+	displaystring(debug,'\n%s',['computing slopes...']);
+	inputs=add(inputs,'surfaceslopex',slopex,'doublevec',m_ss.parameters.numberofdofspernode,m_ss.parameters.numberofnodes);
+	inputs=add(inputs,'surfaceslopey',slopey,'doublevec',m_ss.parameters.numberofdofspernode,m_ss.parameters.numberofnodes);
+
+	displaystring(debug,'\n%s',['computing hutter velocities...']);
+	u_g=diagnostic_core_linear(m_dhu,'diagnostic_hutter',inputs);
+
+	displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']);
+	p_g=ComputePressure(m_dhu.elements,m_dhu.nodes,m_dhu.loads,m_dhu.materials,m_dhu.parameters);
+
+	displaystring(debug,'\n%s',['update boundary conditions for macyeal pattyn using hutter results...']);
+	if ismacayealpattyn,
+		m_dh.y_g=u_g;
+		[m_dh.ys m_dh.ys0]=Reducevectorgtos(m_dh.y_g,m_dh.nodesets);
+	end
+
+end
+
+if ismacayealpattyn,
+
+	displaystring(debug,'\n%s',['computing horizontal velocities...']);
+	u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic_horiz');
+
+	displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']);
+	p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters);
+end
+	
+if dim==3,
+
+	displaystring(debug,'\n%s',['extruding horizontal velocities...']);
+	u_g_horiz=VelocityExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g);
+
+	displaystring(debug,'\n%s',['computing vertical velocities...']);
+	inputs=add(inputs,'velocity',u_g,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
+	u_g_vert=diagnostic_core_linear(m_dv,inputs);
+
+	displaystring(debug,'\n%s',['combining horizontal and vertical velocities...']);
+	u_g=zeros(m_dh.parameters.numberofnodes*3,1);
+	u_g(dofsetgen([1,2],3,m_dh.parameters.numberofnodes*3))=u_g_horiz;
+	u_g(dofsetgen([3],3,m_dh.parameters.numberofnodes*3))=u_g_vert;
+
+	displaystring(debug,'\n%s',['computing pressure according to Pattyn...']);
+	p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters);
+	
+	if isstokes,
+
+		%"recondition" pressure 
+		p_g=p_g/m_ds.parameters.stokesreconditioning;
+
+		displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']);
+		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);
+
+		displaystring(debug,'\n%s',['update boundary conditions for stokes using velocities previously computed...']);
+		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);
+
+		displaystring(debug,'\n%s',['computing stokes velocities and pressure ...']);
+		u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic_stokes');
+	
+		%"decondition" pressure
+		p_g=u_g(4:4:end)*m_dh.parameters.stokesreconditioning;
+	end
+end
Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m	(revision 301)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m	(revision 302)
@@ -1,3 +1,3 @@
-function [u_g varargout]=diagnostic_core_nonlinear(m,inputs)
+function [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type)
 %INPUT function [ru_g varargout]=cielodiagnostic_core_nonlinear(m,inputs)
 	
Index: /issm/trunk/src/m/solutions/cielo/loadresults.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/loadresults.m	(revision 302)
+++ /issm/trunk/src/m/solutions/cielo/loadresults.m	(revision 302)
@@ -0,0 +1,39 @@
+function md=loadresults(md,u_g,p_g,m_dh,m_ds,m_dhu);
+%LOADRESULTS - load results onto model
+%
+%   The solution (vel,pressure,...) are on the g-set,
+%   use index to get the value on each node
+%
+%   Usage:
+%      md=Loadresults(md,u_g,p_g,m_dh,m_ds,m_dhu)
+
+debug=m_dhu.parameters.debug;
+dim=m_dhu.parameters.dim;
+ishutter=m_dhu.parameters.ishutter;
+ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
+isstokes=m_ds.parameters.isstokes;
+
+if dim==2,
+	
+	md.vx=u_g(1:2:end)*md.yts;
+	md.vy=u_g(2:2:end)*md.yts;
+	md.vel=sqrt(md.vx.^2+md.vy.^2);
+	md.pressure=p_g;
+
+else
+
+	if ismacayealpattyn | ishutter,
+	
+		md.vx=u_g(1:3:end)*md.yts;
+		md.vy=u_g(2:3:end)*md.yts;
+		md.vz=u_g(3:3:end)*md.yts;
+		md.pressure=p_g;
+	
+	else
+		md.vx=u_g(1:4:end)*md.yts;
+		md.vy=u_g(2:4:end)*md.yts;
+		md.vz=u_g(3:4:end)*md.yts;
+		md.pressure=u_g(4:4:end)
+		md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
+	end
+end