Index: /issm/trunk/src/m/solutions/jpl/balancedthickness2_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/balancedthickness2_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/balancedthickness2_core.m	(revision 4126)
@@ -1,33 +1,30 @@
-function results=balancedthickness2_core(models,analysis_type,sub_analysis_type)
+function femmodel=balancedthickness2_core(femmodel)
 %BALANCEDTHICKNESS_CORE - linear solution sequence
 %
 %   Usage:
-%      h_g=balancedthickness2_core(m,analysis_type,sub_analysis_type)
+%      femmodel=balancedthickness2_core(femmodel)
 
-	%get FE model
-	m=models.p;
-	results.time=0;
-	results.step=1;
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
+	dim=femmodel.parameters.Dim;
 
-	displaystring(m.parameters.Verbose,'\n%s',['depth averaging velocity...']);
-	%Take only the first two dofs of m.parameters.u_g
-	vx_g=get(inputs,'vx',1);
-	vy_g=get(inputs,'vy',1);
+	%Activate formulation
+	femmodel=SetCurrentAnalysis(femmodel,Balancedthickness2AnalysisEnum);
 
-	%NOT WORKING YET!!!!!
-	%vx_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vx_g,'vx');
-	%vy_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vy_g,'vy');
+	displaystring(verbose,'\n%s',['depth averaging velocities...']);
+	%[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum);
+	%[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum);
+	if dim==3,
+	%	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzAverageEnum);
+	end
 
-	inputs=add(inputs,'vx_average',vx_g,'doublevec',1,m.parameters.NumberOfVertices);
-	inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.NumberOfVertices);
+	displaystring(verbose,'\n%s',['call computational core...']);
+	femmodel=solver_linear(femmodel);
+	
+	displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
+	%femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
 
-	displaystring(m.parameters.Verbose,'\n%s',['call computational core:']);
-	results.h_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type);
-
-	displaystring(m.parameters.Verbose,'\n%s',['averaging over vertices']);
-	results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.h_g);
-
-	displaystring(m.parameters.Verbose,'\n%s',['extrude computed thickness on all layers:']);
-	%results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness');
-
+	displaystring(verbose,'\n%s',['saving results...']);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
+	
 end %end function
Index: /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m	(revision 4126)
@@ -1,27 +1,30 @@
-function results=balancedthickness_core(models,analysis_type,sub_analysis_type)
+function femmodel=balancedthickness_core(femmodel)
 %BALANCEDTHICKNESS_CORE - linear solution sequence
 %
 %   Usage:
-%      h_g=balancedthickness_core(m,analysis_type,sub_analysis_type)
+%      femmodel=balancedthickness_core(femmode)
 
-	%get FE model
-	verbose=models.bt.parameters.Verbose;
-	results.time=0;
-	results.step=1;
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
+	dim=femmodel.parameters.Dim;
 
-	displaystring(verbose,'\n%s',['depth averaging velocity...']);
-	[models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters]=DepthAverageInput(models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters,VxEnum);
-	[models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters]=DepthAverageInput(models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters,VyEnum);
+	%Activate formulation
+	femmodel=SetCurrentAnalysis(femmodel,BalancedthicknessAnalysisEnum);
 
-	displaystring(verbose,'\n%s',['call computational core:']);
-	h_g=diagnostic_core_linear(models.bt,analysis_type,sub_analysis_type);
+	displaystring(verbose,'\n%s',['depth averaging velocities...']);
+	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum);
+	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum);
+	if dim==3,
+		[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzAverageEnum);
+	end
 
-	%Update
-	models.bt.elements=UpdateInputsFromSolution(models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters,h_g,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
+	displaystring(verbose,'\n%s',['call computational core...']);
+	femmodel=solver_linear(femmodel);
+	
+	displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
+	femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
 
-	displaystring(verbose,'\n%s',['extrude computed thickness on all layers:']);
-	models.bt.elements=ExtrudeInput(models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters,ThicknessEnum);
-
-	%NEED TO BE CLEANED
-	results.h_g=h_g;
+	displaystring(verbose,'\n%s',['saving results...']);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
+	
 end %end function
Index: /issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m	(revision 4126)
@@ -1,24 +1,32 @@
-function results=balancedvelocities_core(models,analysis_type,sub_analysis_type)
+function femmodel=balancedvelocities_core(femmdoel)
 %BALANCEDVELOCITIES_CORE - linear solution sequence
 %
 %   Usage:
-%      v_g=balancedvelocities_core(m,analysis_type,sub_analysis_type)
+%      femmodel=balancedvelocities_core(femmodel)
 
-	%get FE model
-	m=models.p;
-	results.time=0;
-	results.step=1;
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
+	dim=femmodel.parameters.Dim;
 
-	displaystring(m.parameters.Verbose,'\n%s',['depth averaging velocity...']);
-	%Take only the first two dofs of m.parameters.u_g
-	u_g=get(inputs,'velocity',[1 1 0 0]);
-	u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,u_g,'velocity');
-	inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.NumberOfNodes);
+	%Activate formulation
+	femmodel=SetCurrentAnalysis(femmodel,BalancedvelocitiesAnalysisEnum);
 
-	displaystring(m.parameters.Verbose,'\n%s',['call computational core:']);
-	results.h_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type);
+	displaystring(verbose,'\n%s',['depth averaging velocities...']);
+	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum);
+	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum);
+	if dim==3,
+		[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzAverageEnum);
+	end
 
-	displaystring(m.parameters.Verbose,'\n%s',['extrude computed thickness on all layers:']);
-	results.v_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness',0);
+	displaystring(verbose,'\n%s',['call computational core...']);
+	femmodel=solver_linear(femmodel);
+	
+	displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
+	femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
+	femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
+
+	displaystring(verbose,'\n%s',['saving results...']);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
 
 end %end function
Index: /issm/trunk/src/m/solutions/jpl/bedslope_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/bedslope_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/bedslope_core.m	(revision 4126)
@@ -22,6 +22,6 @@
 	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);
+		femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
+		femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
 	}
 	
Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 4126)
@@ -21,8 +21,8 @@
 	%for qmu analysis, be sure the velocity input we are starting from  is the one in the parameters: 
 	if qmu_analysis,
-		InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVxEnum,VxEnum);
-		InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVyEnum,VyEnum);
-		InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVzEnum,VzEnum);
-		InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuPressureEnum,PressureEnum);
+		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVxEnum,VxEnum);
+		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVyEnum,VyEnum);
+		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVzEnum,VzEnum);
+		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuPressureEnum,PressureEnum);
 	end
 
@@ -58,5 +58,5 @@
 
 			%"recondition" pressure computed previously:
-			InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum);
+			femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum);
 			InputScale(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning);
 
Index: /issm/trunk/src/m/solutions/jpl/prognostic2_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/prognostic2_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/prognostic2_core.m	(revision 4126)
@@ -1,33 +1,26 @@
-function results=prognostic2_core(models,analysis_type,sub_analysis_type)
+function femmodel=prognostic2_core(femmodel)
 %PROGNOSTIC2_CORE - linear solution sequence
 %
 %   Usage:
-%      h_g=prognostic2_core(m,analysis_type,sub_analysis_type)
+%      femmodel=prognostic2_core(femmodel)
 
-	%get FE model
-	m=models.p;
-	results.time=0;
-	results.step=1;
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
 
-	displaystring(m.parameters.Verbose,'\n%s',['depth averaging velocity...']);
-	%Take only the first two dofs of m.parameters.u_g
-	vx_g=get(inputs,'vx',1);
-	vy_g=get(inputs,'vy',1);
+	%Activate formulation
+	femmodel=SetCurrentAnalysis(femmodel,Prognostic2AnalysisEnum);
 
-	%NOT WORKING YET!!!!!
-	%vx_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vx_g,'vx');
-	%vy_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vy_g,'vy');
+	displaystring(verbose,'\n%s',['depth averaging velocities...']);
+%	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum);
+%	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum);
 
-	inputs=add(inputs,'vx_average',vx_g,'doublevec',1,m.parameters.NumberOfVertices);
-	inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.NumberOfVertices);
+	displaystring(verbose,'\n%s',['call computational core...']);
+	femmodel=solver_linear(femmodel);
+	
+	displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
+%	femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
 
-	displaystring(m.parameters.Verbose,'\n%s',['call computational core:']);
-	results.h_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type);
-
-	displaystring(m.parameters.Verbose,'\n%s',['averaging over vertices']);
-	results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.h_g);
-
-	displaystring(m.parameters.Verbose,'\n%s',['extrude computed thickness on all layers:']);
-	%results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness');
-
+	displaystring(verbose,'\n%s',['saving results...']);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
+	
 end %end function
Index: /issm/trunk/src/m/solutions/jpl/prognostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/prognostic_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/prognostic_core.m	(revision 4126)
@@ -1,28 +1,26 @@
-function results=prognostic_core(models,analysis_type,sub_analysis_type)
+function femmodel=prognostic_core(femmodel)
 %PROGNOSTIC_CORE - linear solution sequence
 %
 %   Usage:
-%      results=prognostic_core(m,analysis_type,sub_analysis_type)
+%      femmodel=prognostic_core(femmodel)
 
-	%get FE model
-	verbose=models.p.parameters.Verbose;
-	results.time=0;
-	results.step=1;
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
 
-	displaystring(verbose,'\n%s',['depth averaging velocity...']);
-	[models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters]=DepthAverageInput(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,VxEnum);
-	[models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters]=DepthAverageInput(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,VyEnum);
+	%Activate formulation
+	femmodel=SetCurrentAnalysis(femmodel,PrognosticAnalysisEnum);
 
-	displaystring(verbose,'\n%s',['call computational core:']);
-	h_g=diagnostic_core_linear(models.p,analysis_type,sub_analysis_type);
+	displaystring(verbose,'\n%s',['depth averaging velocities...']);
+	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum);
+	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum);
 
-	%Update
-	models.p.elements=UpdateInputsFromSolution(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,h_g,PrognosticAnalysisEnum,NoneAnalysisEnum);
+	displaystring(verbose,'\n%s',['call computational core...']);
+	femmodel=solver_linear(femmodel);
+	
+	displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
+	femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
 
-	displaystring(verbose,'\n%s',['extrude computed thickness on all layers:']);
-	models.p.elements=ExtrudeInput(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,ThicknessEnum);
-
-	%NEED TO BE CLEANED
-	results.h_g=h_g;
+	displaystring(verbose,'\n%s',['saving results...']);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
 
 end %end function
Index: /issm/trunk/src/m/solutions/jpl/steadystate_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/steadystate_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/steadystate_core.m	(revision 4126)
@@ -1,87 +1,53 @@
-function results=steadystate_core(models);
+function femmodel=steadystate_core(femmodel);
 %STEADYSTATE_CORE - compute the core temperature and velocity field  at thermal steady state.
 %
 %   Usage:
-%      results=steadystate_core(models);
+%      femmodel=steadystate_core(femmodel);
 %
 
-%recover models
-m_dh=models.dh;
-m_dv=models.dv;
-m_ds=models.ds;
-m_dhu=models.dhu;
-m_sl=models.sl;
-m_t=models.t;
-m_m=models.m;
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
+	dim=femmodel.parameters.Dim;
 
-%recover parameters common to all solutions
-verbose=m_dhu.parameters.Verbose;
-dim=m_dhu.parameters.Dim;
-eps_rel=m_dhu.parameters.EpsRel;
-ishutter=m_dhu.parameters.IsHutter;
-ismacayealpattyn=m_dh.parameters.IsMacAyealPattyn;
-isstokes=m_ds.parameters.IsStokes;
+	%Initialize counter
+	step=1;
 
-%convergence
-converged=0;
+	while true,
 
-step=1;
+		displaystring(verbose,'\n%s%i/\n','computing velocities and temperatures for step: ',step);
+		femmodel=thermal_core(femmodel); 
 
-%initialize: 
-t_g_old=0;
-u_g_old=0;
+		displaystring(verbose,'\n%s',['computing depth average temperature...']);
+		[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum,TemperatureAverageEnum);
 
-if isstokes, ndof=4; else ndof=3; end
+		displaystring(verbose,'\n%s',['computing new velocity...']);
+		femmodel=diagnostic_core(femmodel); 
 
-while(~converged),
-	
-	displaystring(verbose,'%s%i','computing temperature and velocity for step ',step);
-	
-	%first compute temperature at steady state.
-	if (step>1),
-		inputs=add(inputs,'velocity',results_diagnostic.u_g,'doublevec',ndof,m_t.parameters.NumberOfNodes);
-	end
-	results_thermal=thermal_core(models);
-	
-	%add temperature to inputs.
-	%Compute depth averaged temperature
-	temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.vertices,m_t.loads,m_t.materials,m_t.parameters,results_thermal.t_g,'temperature');
-	inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,m_t.parameters.NumberOfNodes);
-	inputs=add(inputs,'temperature',results_thermal.t_g,'doublevec',1,m_t.parameters.NumberOfNodes);
-	
-	%now compute diagnostic velocity using the steady state temperature.
-	results_diagnostic=diagnostic_core(models);
-	
-	%convergence? 
-	du_g=results_diagnostic.u_g-u_g_old;
-	ndu=norm(du_g,2); 
-	nu=norm(u_g_old,2);
-	
-	dt_g=results_thermal.t_g-t_g_old;
-	ndt=norm(dt_g,2); 
-	nt=norm(t_g_old,2); 
+		displaystring(verbose,'\n%s',['checking temperature, velocity and pressure convergence...']);
+		if step>1,
+			if steadystateconvergence(femmodel),
+				break;
+			end
+		end 
 
-	u_g_old=results_diagnostic.u_g;
-	t_g_old=results_thermal.t_g;
-	
-	displaystring(verbose,'%-60s%g\n                                     %s%g\n                                     %s%g%s',...
-	              '      relative convergence criterion: velocity -> norm(du)/norm(u)=   ',ndu/nu*100,' temperature -> norm(dt)/norm(t)=',ndt/nt*100,' eps_rel:                        ',eps_rel*100,' %');
-	if (ndu/nu<=eps_rel)  & (ndt/nt<=eps_rel),
-		converged=1;
-	else
-		converged=0;
+		displaystring(verbose,'\n%s',['saving velocities, temperature and pressure for convergence...']);
+		femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxOldEnum);
+		femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyOldEnum);
+		femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzOldEnum);
+		femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureOlddnum);
+		femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum);
+
+		%Increase counter
+		step=step+1;
+
 	end
 
-	step=step+1;
-	if converged,
-		break;
+	displaystring(verbose,'\n%s',['saving results...']);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum);
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum);
 	end
-end
 
-%save results from thermal and diagnostic
-results.u_g=results_diagnostic.u_g;
-results.p_g=results_diagnostic.p_g;
-results.t_g=results_thermal.t_g;
-results.m_g=results_thermal.m_g;
-results.step=step;
-results.time=0;
+end %end of function
Index: /issm/trunk/src/m/solutions/jpl/surfaceslope_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/surfaceslope_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/surfaceslope_core.m	(revision 4126)
@@ -5,5 +5,4 @@
 %      femmodel=surfaceslope_core(femmodel)
 %
-
 
 	%Recover some parameters:
@@ -22,6 +21,6 @@
 	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);
+		femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
+		femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
 	}
 	
@@ -29,2 +28,4 @@
 	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);
+
+end %end function
Index: /issm/trunk/src/m/solutions/jpl/thermal_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/thermal_core.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/thermal_core.m	(revision 4126)
@@ -1,56 +1,34 @@
-function results=thermal_core(models)
+function femmodel=thermal_core(femmodel)
 %THERMAL_CORE - core of thermal solution
 %
 %   Usage:
-%      solution=thermal_core(models)
+%      femmodel=thermal_core(femmodel)
 
-%recover parameters common to all solutions
-verbose=models.t.parameters.Verbose;
 
-if models.t.parameters.Dt==0,
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
+	ndt=femmodel.parameters.Ndt;
+	dt=femmodel.parameters.Dt;
 
-	results.time=0;
-	results.step=1;
-
-	displaystring(verbose,'\n%s',['computing temperatures...']);
-	[t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,ThermalAnalysisEnum(),NoneAnalysisEnum());
-
-	displaystring(verbose,'\n%s',['computing melting...']);
-	models=ModelUpdateInputsFromVector(models,t_g,TemperatureEnum,VertexEnum);
-	[models.m.elements models.m.loads]=UpdateInputsFromConstant(models.m.elements,models.m.nodes,models.m.vertices,models.m.loads,models.m.materials,models.m.parameters,melting_offset,MeltingOffsetEnum);
-	m_g=diagnostic_core_linear(models.m,MeltingAnalysisEnum(),NoneAnalysisEnum());
-
-	%NEED TO BE CLEANED
-	results.t_g=t_g;
-	results.m_g=m_g;
-
-else
-
-	%initialize temperature and melting
-	nsteps=models.t.parameters.Ndt/models.t.parameters.Dt;
-	results.step=1;
-	results.time=0;
-	results.t_g=[]; %FIRST VALUE MISSING
-	results.m_g=[];
-
-	for n=1:nsteps, 
-
-		displaystring(verbose,'\n%s%i/%i\n','time step: ',n,nsteps);
-		results(n+1).step=n+1;
-		results(n+1).time=n*models.t.parameters.Dt;
-
-		displaystring(verbose,'\n%s',['    computing temperatures...']);
-		[t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,ThermalAnalysisEnum(),NoneAnalysisEnum());
-
-		displaystring(verbose,'\n%s',['    computing melting...']);
-		models=ModelUpdateInputsFromVector(models,t_g,TemperatureEnum,VertexEnum);
-		[models.m.elements models.m.loads]=UpdateInputsFromConstant(models.m.elements,models.m.nodes,models.m.vertices,models.m.loads,models.m.materials,models.m.parameters,melting_offset,MeltingOffsetEnum);
-		m_g=diagnostic_core_linear(models.m,MeltingAnalysisEnum(),NoneAnalysisEnum());
-
-		%NEED TO BE CLEANED
-		results(n+1).t_g=t_g;
-		results(n+1).m_g=m_g;
-
+	%Compute number of timesteps
+	if (dt==0 | ndt==0),
+		dt=0;
+		nsteps=1;
+	else
+		nsteps=floor(ndt/dt);
 	end
 
-end
+	%Loop through time
+	for i=1:nsteps+1,
+		displaystring(verbose,'\n%s%i/%i\n','time step: ',i,nsteps);
+		time=(i+1)*dt;
+
+		displaystring(verbose,'\n%s',['computing temperature...']);
+		femmodel=thermal_core_step(femmodel,i,time); 
+
+		displaystring(verbose,'\n%s',['saving results...']);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,i,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,i,time);
+	end
+
+end %end of function
Index: /issm/trunk/src/m/solutions/jpl/thermal_core_step.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/thermal_core_step.m	(revision 4126)
+++ /issm/trunk/src/m/solutions/jpl/thermal_core_step.m	(revision 4126)
@@ -0,0 +1,18 @@
+function femmodel=thermal_core_step(femmodel)
+%THERMAL_CORE_STEP - core of the thermal solution for one step 
+%
+%   Usage:
+%      femmodel=thermal_core_step(femmodel)
+
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
+
+	displaystring(verbose,'\n%s',['computing temperature...']);
+	femmodel=SetCurrentAnalysis(femmodel,ThermalAnalysisEnum);
+	femmodel=solver_thermal_nonlineat(femmodel);
+
+	displaystring(verbose,'\n%s',['compute melting...']);
+	femmodel=SetCurrentAnalysis(femmodel,MeltingAnalysisEnum);
+	femmodel=solver_linear(femmodel);
+	
+end %end function
Index: sm/trunk/src/m/solutions/jpl/transient.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/transient.m	(revision 4125)
+++ 	(revision )
@@ -1,18 +1,0 @@
-function md=transient(md);
-%TRANSIENT - transient solution
-%
-%   Usage:
-%      md=transient(md)
-
-%start timing
-t1=clock;
-
-if (md.dim==2),
-	md=transient2d(md);
-else
-	md=transient3d(md);
-end
-
-%stop timing
-t2=clock;
-disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']));
Index: /issm/trunk/src/m/solutions/jpl/transient2d.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/transient2d.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/transient2d.m	(revision 4126)
@@ -7,114 +7,28 @@
 %   See also: TRANSIENT3D, TRANSIENT
 
-%Build all models requested
-models.analysis_type='transient'; %needed for processresults
+	analysis_types=[DiagnosticHorizAnalysisEnum,PrognosticAnalysisEnum];
+	solution_type=Transient2DAnalysisEnum;
 
-displaystring(md.verbose,'%s',['reading diagnostic horiz model data']);
-md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HorizAnalysisEnum(); models.dh=CreateFemModel(md);
+	displaystring(md.verbose,'%s',['create finite element model']);
+	femmodel=NewFemModel(md,solution_type,analysis_types,2);
 
-displaystring(md.verbose,'\n%s',['reading diagnostic vert model data']);
-md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=VertAnalysisEnum(); models.dv=CreateFemModel(md);
+	%retrieve parameters
+	verbose=femmodel.parameters.Verbose;
+	qmu_analysis=femmodel.parameters.QmuAnalysis;
 
-displaystring(md.verbose,'\n%s',['reading diagnostic stokes model data']);
-md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=StokesAnalysisEnum(); models.ds=CreateFemModel(md);
+	%compute solution
+	if ~qmu_analysis,
+		displaystring(verbose,'%s',['call computational core']);
+		femmodel=transient2d_core(femmodel);
+		
+		displaystring(verbose,'%s',['write results']);
+		OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters);
 
-displaystring(md.verbose,'\n%s',['reading diagnostic hutter model data']);
-md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HutterAnalysisEnum(); models.dhu=CreateFemModel(md);
+	else
+		%launch dakota driver for diagnostic core solution
+		Qmu(femmodel);
+	end
 
-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);
-
-displaystring(md.verbose,'%s',['reading prognostic model data']);
-md.analysis_type=PrognosticAnalysisEnum(); models.p=CreateFemModel(md);
-
-%initialize solution
-solution=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]);
-solution.step=1;
-solution.time=0;
-%solution.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g)));
-solution.u_g=[];
-solution.p_g=[];
-%solution.h_g=models.p.parameters.h_g;
-solution.h_g=md.thickness;
-solution.s_g=md.surface;
-%solution.s_g=models.p.parameters.s_g;
-solution.b_g=md.bed;
-%solution.b_g=models.p.parameters.b_g;
-
-%initialize inputs
-displaystring(md.verbose,'\n%s',['setup inputs...']);
-inputs=inputlist;
-inputs=add(inputs,'velocity',solution.u_g,'doublevec',2,models.p.parameters.NumberOfNodes);
-inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-inputs=add(inputs,'dt',models.p.parameters.Dt*models.p.parameters.Yts,'double'); %change time from yrs to sec
-
-% figure out number of dof: just for information purposes.
-md.dof=modelsize(models);
-
-%first time step is given by model. 
-dt=models.p.parameters.Dt;
-finaltime=models.p.parameters.Ndt;
-time=dt;
-yts=models.p.parameters.Yts;
-n=1; %counter
-
-while  time<finaltime+dt, %make sure we run up to finaltime.
-
-	disp(sprintf('\n%s%g%s%g%s%g\n','time [yr]: ',time,'    iteration number: ',n,'/',floor(finaltime/dt)));
-
-	solution(n+1).step=n+1;
-	solution(n+1).time=time;
-
-	%update inputs
-	inputs=add(inputs,'thickness',solution(n).h_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-	inputs=add(inputs,'surface',solution(n).s_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-	inputs=add(inputs,'bed',solution(n).b_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-
-	%Deal with velocities.
-
-	%Get horizontal solution. 
-	rawresults=diagnostic_core(models);
-	solution(n+1).u_g=rawresults.u_g; solution(n+1).p_g=rawresults.p_g;
-
-	%compute new thickness
-	disp(sprintf('%s','   computing new thickness...'));
-	inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',2,models.p.parameters.NumberOfNodes);
-	rawresults=prognostic_core(models,PrognosticAnalysisEnum(),NoneAnalysisEnum());
-	new_thickness=rawresults.h_g;
-
-	%update surface and bed using the new thickness
-	disp(sprintf('%s','   updating geometry...'));
-	[new_thickness,new_bed,new_surface]=UpdateGeometry(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,new_thickness,solution(n).b_g,solution(n).s_g);
-
-	%Record bed surface and thickness in the solution
-	solution(n+1).h_g=new_thickness;
-	solution(n+1).b_g=new_bed;
-	solution(n+1).s_g=new_surface;
-
-	%figure out if time stepping is good
-	%disp(sprintf('%s','   checking time stepping...'));
-	%[back,dt,time]=TimeStepping(md,solution,dt,time);
-	%if back,
-	%	continue;
-	%end
-
-	%update time and counter
-	time=time+dt;
-	n=n+1;
-end
-
-%Load results onto model
-results=struct('step',[],'time',[],'vx',[],'vy',[],'vel',[],'pressure',[],'thickness',[],'surface',[],'bed',[]);
-for i=1:length(solution),
-	results(i).step=solution(i).step;
-	results(i).time=solution(i).time/yts;
-	results(i).vx=solution(i).u_g(1:2:end)*yts;
-	results(i).vy=solution(i).u_g(2:2:end)*yts;
-	results(i).vel=sqrt(solution(i).u_g(1:2:end).^2+solution(i).u_g(2:2:end).^2)*yts;
-	results(i).bed=solution(i).b_g;
-	results(i).surface=solution(i).s_g;
-	results(i).thickness=solution(i).h_g;
-end
-if ~isstruct(md.results), md.results=struct(); end
-md.results.TransientAnalysis=results;
+	%stop timing
+	t2=clock;
+	displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);
Index: /issm/trunk/src/m/solutions/jpl/transient2d_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/transient2d_core.m	(revision 4126)
+++ /issm/trunk/src/m/solutions/jpl/transient2d_core.m	(revision 4126)
@@ -0,0 +1,41 @@
+function femmodel=transient2d_core(femmodel)
+%TRANSIENT2D_CORE - core of transient 2d solution
+%
+%   Usage:
+%      femmodel=transient2d_core(femmodel)
+
+
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
+	ndt=femmodel.parameters.Ndt;
+	dt=femmodel.parameters.Dt;
+
+	%Initialize
+	time=0;
+	step=0;
+
+	%Loop through time
+	while time<ndt,
+		displaystring(verbose,'\n%s%g%s%i%s%g\n','time [yr]',time,'iteration number: ',step,'/',floor(ndt/dt));
+		step=step+1;
+		time=time+dt;
+
+		displaystring(verbose,'\n%s',['computing new velocities...']);
+		femmodel=diagnostic_core(femmodel); 
+
+		displaystring(verbose,'\n%s',['computing new thickness...']);
+		femmodel=prognostic_core(femmodel); 
+
+		displaystring(verbose,'\n%s',['updating geometry...']);
+		[femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=UpdateGeometry(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+		displaystring(verbose,'\n%s',['saving results...']);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
+	end
+
+end %end of function
Index: /issm/trunk/src/m/solutions/jpl/transient3d.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/transient3d.m	(revision 4125)
+++ /issm/trunk/src/m/solutions/jpl/transient3d.m	(revision 4126)
@@ -6,135 +6,28 @@
 %   See also: TRANSIENT2D, TRANSIENT
 
-%Build all models related to diagnostic
-models.analysis_type=TransientAnalysisEnum(); %needed for processresults
+	analysis_types=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,PrognosticAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum];
+	solution_type=Transient3DAnalysisEnum;
 
-displaystring(md.verbose,'%s',['reading diagnostic horiz model data']);
-md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HorizAnalysisEnum(); models.dh=CreateFemModel(md);
+	displaystring(md.verbose,'%s',['create finite element model']);
+	femmodel=NewFemModel(md,solution_type,analysis_types,8);
 
-displaystring(md.verbose,'\n%s',['reading diagnostic vert model data']);
-md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=VertAnalysisEnum(); models.dv=CreateFemModel(md);
+	%retrieve parameters
+	verbose=femmodel.parameters.Verbose;
+	qmu_analysis=femmodel.parameters.QmuAnalysis;
 
-displaystring(md.verbose,'\n%s',['reading diagnostic stokes model data']);
-md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=StokesAnalysisEnum(); models.ds=CreateFemModel(md);
+	%compute solution
+	if ~qmu_analysis,
+		displaystring(verbose,'%s',['call computational core']);
+		femmodel=transient3d_core(femmodel);
+		
+		displaystring(verbose,'%s',['write results']);
+		OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters);
 
-displaystring(md.verbose,'\n%s',['reading diagnostic hutter model data']);
-md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HutterAnalysisEnum(); models.dhu=CreateFemModel(md);
+	else
+		%launch dakota driver for diagnostic core solution
+		Qmu(femmodel);
+	end
 
-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);
-
-displaystring(md.verbose,'%s',['reading prognostic model data']);
-md.analysis_type=PrognosticAnalysisEnum(); models.p=CreateFemModel(md);
-
-%Build all models related to thermal
-displaystring(md.verbose,'%s',['reading thermal model data']);
-md.analysis_type=ThermalAnalysisEnum(); models.t=CreateFemModel(md);
-
-displaystring(md.verbose,'%s',['reading melting model data']);
-md.analysis_type=MeltingAnalysisEnum(); models.m=CreateFemModel(md);
-
-
-%initialize solution
-results=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]);
-results.step=1;
-results.time=0;
-results.u_g=models.dh.parameters.u_g;
-%results.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g)));
-results.p_g=models.p.parameters.p_g;
-results.h_g=models.p.parameters.h_g;
-results.s_g=models.p.parameters.s_g;
-results.b_g=models.p.parameters.b_g;
-results.t_g=models.p.parameters.t_g;
-results.m_g=models.p.parameters.m_g;
-results.a_g=models.p.parameters.a_g;
-
-%initialize inputs
-displaystring(md.verbose,'\n%s',['setup inputs...']);
-inputs=inputlist;
-inputs=add(inputs,'velocity',results.u_g,'doublevec',3,models.p.parameters.NumberOfNodes);
-inputs=add(inputs,'melting',results.m_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-inputs=add(inputs,'accumulation',results.a_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-inputs=add(inputs,'dt',models.p.parameters.Dt*models.p.parameters.Yts,'double');
-
-% figure out number of dof: just for information purposes.
-md.dof=modelsize(models);
-
-%first time step is given by model. 
-dt=models.p.parameters.Dt;
-finaltime=models.p.parameters.Ndt;
-time=dt;
-yts=models.p.parameters.Yts;
-n=1; %counter
-
-while  time<finaltime+dt, %make sure we run up to finaltime.
-
-	displaystring(md.verbose,'\n%s%g%s%g%s%g\n','time [yr]: ',time,'    iteration number: ',n,'/',floor(finaltime/dt));
-
-	results(n+1).step=n+1;
-	results(n+1).time=time;
-
-	%update inputs
-	inputs=add(inputs,'thickness',results(n).h_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-	inputs=add(inputs,'surface',results(n).s_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-	inputs=add(inputs,'bed',results(n).b_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-	inputs=add(inputs,'velocity',results(n).u_g,'doublevec',3,models.p.parameters.NumberOfNodes);
-	inputs=add(inputs,'pressure',results(n).p_g,'doublevec',1,models.p.parameters.NumberOfNodes);
-	inputs=add(inputs,'temperature',results(n).t_g,'doublevec',1,models.t.parameters.NumberOfNodes);
-
-	%Deal with temperature first 
-	displaystring(md.verbose,'\n%s',['    computing temperatures...']);
-	[results(n+1).t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,ThermalAnalysisEnum(),TransientAnalysisEnum());
-	inputs=add(inputs,'temperature',results(n+1).t_g,'doublevec',1,models.t.parameters.NumberOfNodes);
-	
-	displaystring(md.verbose,'\n%s',['    computing melting...']);
-	inputs=add(inputs,'melting_offset',melting_offset,'double');
-	results(n+1).m_g=diagnostic_core_linear(models.m,MeltingAnalysisEnum(),TransientAnalysisEnum());
-
-	%Compute depth averaged temperature
-	temperature_average=FieldDepthAverage(models.t.elements,models.t.nodes,models.t.vertices,models.t.loads,models.t.materials,models.t.parameters,results(n+1).t_g,'temperature');
-	inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,models.t.parameters.NumberOfNodes);
-
-	%Deal with velocities.
-	rawresults=diagnostic_core(models);
-	results(n+1).u_g=rawresults.u_g; results(n+1).p_g=rawresults.p_g;
-
-	%compute new thickness
-	displaystring(md.verbose,'\n%s',['    computing new thickness...']);
-	inputs=add(inputs,'velocity',results(n+1).u_g,'doublevec',3,models.p.parameters.NumberOfNodes);
-	rawresults=prognostic_core(models,PrognosticAnalysisEnum(),NoneAnalysisEnum());
-	new_thickness=rawresults.h_g;
-
-	%update surface and bed using the new thickness
-	displaystring(md.verbose,'\n%s',['    updating geometry...']);
-	[new_thickness,new_bed,new_surface]=UpdateGeometry(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,new_thickness,results(n).b_g,results(n).s_g);
-
-	%Record bed surface and thickness in the results
-	results(n+1).h_g=new_thickness;
-	results(n+1).b_g=new_bed;
-	results(n+1).s_g=new_surface;
-
-	%figure out if time stepping is good
-	%displaystring(md.verbose,'   checking time stepping...'));
-	%[back,dt,time]=TimeStepping(md,results,dt,time);
-	%if back,
-	%	continue;
-	%end
-
-	%update node positions
-	displaystring(md.verbose,'\n%s',['    updating node positions...']);
-	models.dh.vertices=UpdateVertexPositions(models.dh.vertices,new_bed,new_thickness);
-	models.dv.vertices=UpdateVertexPositions(models.dv.vertices,new_bed,new_thickness);
-	models.ds.vertices=UpdateVertexPositions(models.ds.vertices,new_bed,new_thickness);
-	models.sl.vertices=UpdateVertexPositions(models.sl.vertices,new_bed,new_thickness);
-	models.p.vertices =UpdateVertexPositions(models.p.vertices,new_bed,new_thickness);
-	models.t.vertices =UpdateVertexPositions(models.t.vertices,new_bed,new_thickness);
-	models.m.vertices =UpdateVertexPositions(models.m.vertices,new_bed,new_thickness);
-	
-	%update time and counter
-	time=time+dt;
-	n=n+1;
-end
-
-%process results
-if ~isstruct(md.results), md.results=struct(); end
-md.results.transient=processresults(models, results);
+	%stop timing
+	t2=clock;
+	displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);
Index: /issm/trunk/src/m/solutions/jpl/transient3d_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/transient3d_core.m	(revision 4126)
+++ /issm/trunk/src/m/solutions/jpl/transient3d_core.m	(revision 4126)
@@ -0,0 +1,53 @@
+function femmodel=transient3d_core(femmodel)
+%TRANSIENT3D_CORE - core of transient 2d solution
+%
+%   Usage:
+%      femmodel=transient3d_core(femmodel)
+
+
+	%recover parameters common to all solutions
+	verbose=femmodel.parameters.Verbose;
+	ndt=femmodel.parameters.Ndt;
+	dt=femmodel.parameters.Dt;
+
+	%Initialize
+	time=0;
+	step=0;
+
+	%Loop through time
+	while time<ndt,
+		displaystring(verbose,'\n%s%g%s%i%s%g\n','time [yr]',time,'iteration number: ',step,'/',floor(ndt/dt));
+		step=step+1;
+		time=time+dt;
+
+		displaystring(verbose,'\n%s',['computing temperature...']);
+		femmodel=thermal_core_step(femmodel); 
+
+		displaystring(verbose,'\n%s',['computing depth average temperature...']);
+		[femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=InputDepthAverage(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureAverageEnum);
+
+		displaystring(verbose,'\n%s',['computing new velocities...']);
+		femmodel=diagnostic_core(femmodel); 
+
+		displaystring(verbose,'\n%s',['computing new thickness...']);
+		femmodel=prognostic_core(femmodel); 
+
+		displaystring(verbose,'\n%s',['updating geometry...']);
+		[femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=UpdateGeometry(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+		displaystring(verbose,'\n%s',['updating vertices position...']);
+		[femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=UpdateVertexPostions(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+		displaystring(verbose,'\n%s',['saving results...']);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time);
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingEnum,step,time);
+	end
+
+end %end of function