Index: /issm/trunk/src/m/solutions/jpl/bedslope.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/bedslope.m	(revision 4112)
+++ /issm/trunk/src/m/solutions/jpl/bedslope.m	(revision 4113)
@@ -1,7 +1,7 @@
-function md=slopecompute(md);
-%SLOPECOMPUTE - compute the slope of a model
+function md=bedslopecompute(md);
+%SLOPECOMPUTE - compute the bed slope of a model
 %
 %   Usage:
-%      md=slopecompute(md)
+%      md=bedslope(md)
 %
 	%timing
@@ -21,5 +21,9 @@
 	if ~qmu_analysis,
 		displaystring(verbose,'%s',['call computational core']);
-		bedslope_core(femmodel);
+		femmodel=bedslope_core(femmodel);
+	
+		displaystring(verbose,'%s',['write results']);
+		OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters);
+
 	else
 		%launch dakota driver for diagnostic core solution
Index: /issm/trunk/src/m/solutions/jpl/bedslope_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/bedslope_core.m	(revision 4113)
+++ /issm/trunk/src/m/solutions/jpl/bedslope_core.m	(revision 4113)
@@ -0,0 +1,30 @@
+function femmodel=bedslope_core(femmodel)
+%SURFACESLOPE_CORE - core of the bed slope computation solution
+%
+%   Usage:
+%      femmodel=bedslope_core(femmodel)
+%
+
+
+	%Recover some parameters:
+	verbose=femmodel.parameters.Verbose;
+	dim=femmodel.parameters.Dim;
+
+	displaystring(verbose,'\n%s',['computing bed slope...']);
+
+	%Call on core computations: 
+	femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum);
+	femmodel=solver_linear(femmodel);
+	femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum);
+	femmodel=solver_linear(femmodel);
+	
+	%extrude inputs if we are in 3D: */
+	if dim==3,
+		displaystring(verbose,'\n%s',['extruding bed slopein 3d...']);
+		InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
+		InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
+	}
+	
+	displaystring(verbose,'\n%s',['saving results...']);
+	InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
+	InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
Index: /issm/trunk/src/m/solutions/jpl/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/diagnostic.m	(revision 4112)
+++ /issm/trunk/src/m/solutions/jpl/diagnostic.m	(revision 4113)
@@ -24,12 +24,16 @@
 			
 			displaystring(verbose,'%s',['call computational core']);
-			diagnostic_core(femmodel);
+			femmodel=diagnostic_core(femmodel);
 
 		else,
 			
 			displaystring(verbose,'%s',['call computational core']);
-			control_core(femmodel);
+			femmodel=control_core(femmodel);
 
 		end
+		
+		displaystring(verbose,'%s',['write results']);
+		OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters);
+
 	else
 		%launch dakota driver for diagnostic core solution
Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 4112)
+++ /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 4113)
@@ -1,3 +1,3 @@
-function results=diagnostic_core(model);
+function femmodel=diagnostic_core(femmodel);
 %DIAGNOSTIC_CORE - compute the core velocity field 
 %
@@ -28,6 +28,6 @@
 
 	%Compute slopes: 
-	if(ishutter)surfaceslope_core(femmodel);end
-	if(isstokes)bedslope_core(femmodel);end
+	if(ishutter)femmodel=surfaceslope_core(femmodel);end
+	if(isstokes)femmodel=bedslope_core(femmodel);end
 
 	if ishutter,
@@ -59,5 +59,5 @@
 			%"recondition" pressure computed previously:
 			InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum);
-			ScaleInput(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning);
+			InputScale(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning);
 
 			displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']);
Index: sm/trunk/src/m/solutions/jpl/slope_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/slope_core.m	(revision 4112)
+++ 	(revision )
@@ -1,56 +1,0 @@
-function [slopex, slopey]=slope_core(fem,AnalysisEnum)
-%SLOPE_CORE - Core solution to compute the slope
-%
-%    Usage: 
-%       [slopex, slopey]=slope_core(fem,AnalysisEnum)
-%
-%    Examples: 
-%       [bedslopex, bedslopey]=slope_core(fem,BedEnum);
-%       [surfaceslopex, surfaceslopey]=slope_core(fem,SurfaceEnum);
-%
-
-verbose=fem.parameters.Verbose;
-dim=fem.parameters.Dim;
-isstokes=fem.parameters.IsStokes;
-ishutter=fem.parameters.IsHutter;
-
-displaystring(verbose,'\n%s',['computing slope (x and y derivatives...)']);
-
-%Specify type of computations: 
-if AnalysisEnum==SurfaceAnalysisEnum,
-	xanalysis=SurfaceXAnalysisEnum;
-	yanalysis=SurfaceYAnalysisEnum;
-elseif AnalysisEnum==BedAnalysisEnum,
-	xanalysis=BedXAnalysisEnum;
-	yanalysis=BedYAnalysisEnum;
-else
-	displaystring(verbose,'%s%s%s','analysis ',EnumAsString(AnalysisEnum),' not supported yet!');
-end
-
-%Early return possible?
-if ~ishutter && AnalysisEnum==SurfaceAnalysisEnum,
-	%no need to compute Surface Slope except for Hutter: 
-	slopex=[];
-	slopey=[];
-	return;
-end
-	
-if ~isstokes && AnalysisEnum==BedAnalysisEnum,
-	%no need to compute Bed Slope except for full Stokes:
-
-	slopex=[];
-	slopey=[];
-	return;
-end
-	
-%Call on core computations: 
-slopex=diagnostic_core_linear(fem,SlopecomputeAnalysisEnum,xanalysis);
-slopey=diagnostic_core_linear(fem,SlopecomputeAnalysisEnum,yanalysis);
-
-%extrude if we are in 3D: 
-if dim==3,
-	displaystring(verbose,'%s\n','extruding slopes in 3d...');
-	slopex=FieldExtrude(fem.elements,fem.nodes,fem.vertices,fem.loads,fem.materials,fem.parameters,slopex,'slopex',1);
-	slopey=FieldExtrude(fem.elements,fem.nodes,fem.vertices,fem.loads,fem.materials,fem.parameters,slopey,'slopey',1);
-end
-
Index: sm/trunk/src/m/solutions/jpl/slopecompute.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/slopecompute.m	(revision 4112)
+++ 	(revision )
@@ -1,29 +1,0 @@
-function md=slopecompute(md);
-%SLOPECOMPUTE - compute the slope of a model
-%
-%   Usage:
-%      md=slopecompute(md)
-%
-	%timing
-	t1=clock;
-
-	%Build all models requested for diagnostic simulation
-	models.analysis_type=SlopecomputeAnalysisEnum; %needed for processresults
-
-	displaystring(md.verbose,'\n%s',['reading surface and bed slope computation model data']);
-	md.analysis_type=SlopecomputeAnalysisEnum; md.sub_analysis_type=NoneAnalysisEnum; models.sl=CreateFemModel(md);
-	
-	% figure out number of dof: just for information purposes.
-	md.dof=modelsize(models);
-
-	%compute solution
-	displaystring(md.verbose,'\n%s',['call computational core:']);
-	results=slopecompute_core(models,SlopecomputeAnalysisEnum(),NoneAnalysisEnum());
-
-	displaystring(md.verbose,'\n%s',['load results...']);
-	if ~isstruct(md.results), md.results=struct(); end
-	md.results.SlopecomputeAnalysis=processresults(models,results);
-
-	%stop timing
-	t2=clock;
-	displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);
Index: sm/trunk/src/m/solutions/jpl/slopecompute_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/slopecompute_core.m	(revision 4112)
+++ 	(revision )
@@ -1,20 +1,0 @@
-function results=slopecompute_core(models,analysis_type,sub_analysis_type)
-%SLOPECOMPUTE_CORE - linear solution sequence
-%
-%   Usage:
-%      [sx_g sy_g]=slopecompute_core(m,analysis_type,sub_analysis_type)
-
-	%get FE model
-	m=models.sl;
-	results.time=0;
-	results.step=1;
-
-	displaystring(m.parameters.Verbose,'\n%s',['call computational core:']);
-	results.sx_g=diagnostic_core_linear(m,analysis_type,SurfaceXAnalysisEnum());
-	results.sy_g=diagnostic_core_linear(m,analysis_type,SurfaceYAnalysisEnum());
-
-	displaystring(m.parameters.Verbose,'\n%s',['extrude computed slope on all layers:']);
-	results.sx_g=FieldExtrude(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.sx_g,'slopex',0);
-	results.sy_g=FieldExtrude(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.sy_g,'slopey',0);
-
-end %end function
Index: /issm/trunk/src/m/solutions/jpl/surfaceslope.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/surfaceslope.m	(revision 4112)
+++ /issm/trunk/src/m/solutions/jpl/surfaceslope.m	(revision 4113)
@@ -1,7 +1,7 @@
-function md=slopecompute(md);
-%SLOPECOMPUTE - compute the slope of a model
+function md=surfaceslopecompute(md);
+%SLOPECOMPUTE - compute the surface slope of a model
 %
 %   Usage:
-%      md=slopecompute(md)
+%      md=surfaceslope(md)
 %
 	%timing
@@ -21,5 +21,9 @@
 	if ~qmu_analysis,
 		displaystring(verbose,'%s',['call computational core']);
-		surfaceslope_core(femmodel);
+		femmodel=surfaceslope_core(femmodel);
+	
+		displaystring(verbose,'%s',['write results']);
+		OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters);
+
 	else
 		%launch dakota driver for diagnostic core solution
Index: /issm/trunk/src/m/solutions/jpl/surfaceslope_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/surfaceslope_core.m	(revision 4113)
+++ /issm/trunk/src/m/solutions/jpl/surfaceslope_core.m	(revision 4113)
@@ -0,0 +1,30 @@
+function femmodel=surfaceslope_core(femmodel)
+%SURFACESLOPE_CORE - core of the surface slope computation solution
+%
+%   Usage:
+%      femmodel=surfaceslope_core(femmodel)
+%
+
+
+	%Recover some parameters:
+	verbose=femmodel.parameters.Verbose;
+	dim=femmodel.parameters.Dim;
+
+	displaystring(verbose,'\n%s',['computing surface slope...']);
+
+	%Call on core computations: 
+	femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum);
+	femmodel=solver_linear(femmodel);
+	femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum);
+	femmodel=solver_linear(femmodel);
+	
+	%extrude inputs if we are in 3D: */
+	if dim==3,
+		displaystring(verbose,'\n%s',['extruding surface slopein 3d...']);
+		InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
+		InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
+	}
+	
+	displaystring(verbose,'\n%s',['saving results...']);
+	InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
+	InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);