Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 3836)
+++ /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 3837)
@@ -59,5 +59,5 @@
 
 	displaystring(verbose,'\n%s',['computing horizontal velocities...']);
-	[u_g m_dh.loads]=diagnostic_core_nonlinear(m_dh,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+	[u_g m_dh.loads]=diagnostic_core_nonlinear(m_dh,m_dh.loads,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 
 	if dim==2,
Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m	(revision 3836)
+++ /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m	(revision 3837)
@@ -1,39 +1,35 @@
-function [u_g varargout]=diagnostic_core_nonlinear(m,analysis_type,sub_analysis_type)
+function [u_g varargout]=diagnostic_core_nonlinear(fem,input_loads,analysis_type,sub_analysis_type)
 %DIAGNOSTIC_CORE_NONLINEAR - core solution of diagnostic non-linear
 %
 %   Usage:
-%      [u_g varargout]=diagnostic_core_nonlinear(m,analysis_type,sub_analysis_type)
+%      [u_g varargout]=diagnostic_core_nonlinear(fem,analysis_type,sub_analysis_type)
 	
 %   first off! We are going to modify the loads dataset. We need to shield the loads from those changes once we return;
-	loads=m.loads;
+	loads=fem.loads;
 
 %	stiffness and load generation only:
-	m.parameters.kflag=1; m.parameters.pflag=1;
+	fem.parameters.Kflag=1; fem.parameters.Pflag=1;
+
+%	Were loads requested as outputs ?
+	if isempty(loads),
+		loads=fem.loads;
+	else
+		loads=input_loads;
+	end
 
 %   initialize solution vector
 	converged=0; count=1;
-	soln(count).u_g=get(inputs,'velocity',[1 1 (m.parameters.NumberOfDofsPerNode>=3) (m.parameters.NumberOfDofsPerNode>=3)]); %only keep vz if running with more than 3 dofs per node
 
-	soln(count).u_f= Reducevectorgtof(soln(count).u_g,m.nodesets);
-		
-	%add converged=0 flag first in inputs.
-	inputs=add(inputs,'converged',converged,'double');
+	ug=GetSolutionFromInputs(fem.elements, fem.nodes, fem.vertices, fem.loads, fem.materials, fem.parameters, analysis_type, sub_analysis_type);
+	uf=Reducevectorgtof( ug, fem.nodesets);
 
-	displaystring(m.parameters.Verbose,'\n%s',['   starting direct shooting method']);
+	displaystring(fem.parameters.Verbose,'\n%s',['   starting direct shooting method']);
 	while(~converged),
 		
-		%add velocity to inputs
-		inputs=add(inputs,'velocity',soln(count).u_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes);
-		if count>1,
-			inputs=add(inputs,'old_velocity',soln(count-1).u_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes);
-		else
-			inputs=add(inputs,'old_velocity',soln(count).u_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes);
-		end
-		
-		%Update inputs in datasets
-		[m.elements,m.nodes,m.vertices, loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes,m.vertices, loads,m.materials,m.parameters,inputs);
+		old_ug=ug;
+		old_uf=uf;
 		
 		%system matrices 
-		[K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+		[K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(fem.elements,fem.nodes,fem.vertices,loads,fem.materials,fem.parameters,analysis_type,sub_analysis_type);
 	
 		%penalties
@@ -88,5 +84,5 @@
 	if nout==2,
 		inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes);
-		m.parameters.kflag=1; m.parameters.pflag=0; 
+		m.parameters.Kflag=1; m.parameters.Pflag=0; 
 		[K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
 		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets);