Index: /issm/trunk/src/m/solutions/jpl/ModelUpdateInputsFromVector.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/ModelUpdateInputsFromVector.m	(revision 3872)
+++ /issm/trunk/src/m/solutions/jpl/ModelUpdateInputsFromVector.m	(revision 3873)
@@ -1,3 +1,3 @@
-function ModelUpdateInputsFromVector(models, vector, enum, typeenum);
+function models=ModelUpdateInputsFromVector(models, vector, enum, typeenum);
 %MODELUPDATEINPUTSFROMVECTOR - Update inputs using a vector
 %
@@ -26,4 +26,5 @@
 	if isstruct(model), %there is an analysis_type model
 		[model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters] = UpdateInputsFromVector(model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters,vector,enum, typeenum);
+		models.(field)=model;
 	end
 
Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 3872)
+++ /issm/trunk/src/m/solutions/jpl/diagnostic_core.m	(revision 3873)
@@ -6,50 +6,43 @@
 %
 
-%recover fem models
-m_dh=model.dh;
-m_dv=model.dv;
-m_ds=model.ds;
-m_dhu=model.dhu;
-m_sl=model.sl;
-
 %recover parameters common to all solutions
-verbose=m_dhu.parameters.Verbose;
-dim=m_dh.parameters.Dim;
-ishutter=m_dhu.parameters.IsHutter;
-ismacayealpattyn=m_dh.parameters.IsMacAyealPattyn;
-isstokes=m_ds.parameters.IsStokes;
-numrifts=m_dhu.parameters.NumRifts;
-qmu_analysis=m_dh.parameters.QmuAnalysis;
+verbose=model.dhu.parameters.Verbose;
+dim=model.dh.parameters.Dim;
+ishutter=model.dhu.parameters.IsHutter;
+ismacayealpattyn=model.dh.parameters.IsMacAyealPattyn;
+isstokes=model.ds.parameters.IsStokes;
+numrifts=model.dhu.parameters.NumRifts;
+qmu_analysis=model.dh.parameters.QmuAnalysis;
 
 	
 %for qmu analysis, be sure the velocity input we are starting from  is the one in the parameters: 
 if qmu_analysis,
-	ModelUpdateInputsFromVector(model,m_dh.vx,VxEnum,VertexEnum);
-	ModelUpdateInputsFromVector(model,m_dh.vy,VyEnum,VertexEnum);
-	ModelUpdateInputsFromVector(model,m_dh.vz,VzEnum,VertexEnum);
+	model=ModelUpdateInputsFromVector(model,m_dh.vx,VxEnum,VertexEnum);
+	model=ModelUpdateInputsFromVector(model,m_dh.vy,VyEnum,VertexEnum);
+	model=ModelUpdateInputsFromVector(model,m_dh.vz,VzEnum,VertexEnum);
 end
 
 %Compute slopes: 
-[surfaceslopex,surfaceslopey]=slope_core(m_sl,SurfaceAnalysisEnum);
-[bedslopex,bedslopey]=slope_core(m_sl,BedAnalysisEnum);
+[surfaceslopex,surfaceslopey]=slope_core(model.sl,SurfaceAnalysisEnum);
+[bedslopex,bedslopey]=slope_core(model.sl,BedAnalysisEnum);
 
 %Update:
-ModelUpdateInputsFromVector(model,surfaceslopex,SurfaceSlopexEnum,VertexEnum);
-ModelUpdateInputsFromVector(model,surfaceslopey,SurfaceSlopeyEnum,VertexEnum);
-ModelUpdateInputsFromVector(model,bedslopex,BedSlopexEnum,VertexEnum);
-ModelUpdateInputsFromVector(model,bedslopey,BedSlopeyEnum,VertexEnum);
+model=ModelUpdateInputsFromVector(model,surfaceslopex,SurfaceSlopexEnum,VertexEnum);
+model=ModelUpdateInputsFromVector(model,surfaceslopey,SurfaceSlopeyEnum,VertexEnum);
+model=ModelUpdateInputsFromVector(model,bedslopex,BedSlopexEnum,VertexEnum);
+model=ModelUpdateInputsFromVector(model,bedslopey,BedSlopeyEnum,VertexEnum);
 
 if ishutter,
 
 	displaystring(verbose,'\n%s',['computing hutter velocities...']);
-	u_g=diagnostic_core_linear(m_dhu,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
+	u_g=diagnostic_core_linear(model.dhu,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
 
 	displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']);
-	p_g=ComputePressure(m_dhu.elements,m_dhu.nodes,m_dhu.vertices,m_dhu.loads,m_dhu.materials,m_dhu.parameters,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
+	p_g=ComputePressure(model.dhu.elements,model.dhu.nodes,model.dhu.vertices,model.dhu.loads,model.dhu.materials,model.dhu.parameters,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
 
 	displaystring(verbose,'\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);
+		model.dh.y_g=u_g;
+		[model.dh.ys model.dh.ys0]=Reducevectorgtos(model.dh.y_g,model.dh.nodesets);
 	end
 
@@ -59,9 +52,9 @@
 
 	displaystring(verbose,'\n%s',['computing horizontal velocities...']);
-	[u_g m_dh.loads]=diagnostic_core_nonlinear(m_dh,m_dh.loads,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+	[u_g model.dh.loads]=diagnostic_core_nonlinear(model.dh,model.dh.loads,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 
 	if dim==2,
 		displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']);
-		p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+		p_g=ComputePressure(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 	end
 end
@@ -70,17 +63,21 @@
 
 	displaystring(verbose,'\n%s',['extruding horizontal velocities...']);
-	u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,u_g,'velocity',1);
+	u_g_horiz=FieldExtrude(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,u_g,'velocity',1);
 
-	[vx,vy]=SplitSolutionVector(u_g_horiz,m_dh.parameters.NumberOfNodes,m_dh.parameters.NumberOfDofsPerNode);
-	ModelUpdateInputsFromVector(model,vx,VxEnum,VertexEnum);
-	ModelUpdateInputsFromVector(model,vy,VyEnum,VertexEnum);
+	[vx,vy]=SplitSolutionVector(u_g_horiz,model.dh.parameters.NumberOfNodes,model.dh.parameters.NumberOfDofsPerNode);
+	model=ModelUpdateInputsFromVector(model,vx,VxEnum,VertexEnum);
+	model=ModelUpdateInputsFromVector(model,vy,VyEnum,VertexEnum);
 		
 	displaystring(verbose,'\n%s',['computing vertical velocities...']);
-	u_g_vert=diagnostic_core_linear(m_dv,DiagnosticAnalysisEnum(),VertAnalysisEnum());
-	ModelUpdateInputsFromVector(model,u_g_vert,VzEnum,VertexEnum);
+	u_g_vert=diagnostic_core_linear(model.dv,DiagnosticAnalysisEnum(),VertAnalysisEnum());
+	model=ModelUpdateInputsFromVector(model,u_g_vert,VzEnum,VertexEnum);
 
 	displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']);
-	p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
-	ModelUpdateInputsFromVector(model,p_g,PressureEnum,VertexEnum);
+	p_g=ComputePressure(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+	model=ModelUpdateInputsFromVector(model,p_g,PressureEnum,VertexEnum);
+	u_g=zeros(model.dh.parameters.NumberOfNodes*3,1); %%%%%%%%%%%%%%%%%%%%%%% NEED TO BE CLEANED 
+	u_g(1:3:end)=vx;
+	u_g(2:3:end)=vy;
+	u_g(3:3:end)=u_g_vert;
 	
 	if isstokes,
@@ -90,14 +87,14 @@
 
 		displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']);
-		m_ds.y_g=zeros(m_ds.nodesets.gsize,1);
-		m_ds.y_g(dofsetgen([1,2],4,m_ds.nodesets.gsize))=u_g;
-		m_ds.y_g(dofsetgen([3],4,m_ds.nodesets.gsize))=u_g_vert;
-		[m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets);
+		model.ds.y_g=zeros(model.ds.nodesets.gsize,1);
+		model.ds.y_g(dofsetgen([1,2],4,model.ds.nodesets.gsize))=u_g;
+		model.ds.y_g(dofsetgen([3],4,model.ds.nodesets.gsize))=u_g_vert;
+		[model.ds.ys model.ds.ys0]=Reducevectorgtos(model.ds.y_g,model.ds.nodesets);
 
 		displaystring(verbose,'\n%s',['computing stokes velocities and pressure ...']);
-		u_g=diagnostic_core_nonlinear(m_ds,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
+		u_g=diagnostic_core_nonlinear(model.ds,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
 	
 		%"decondition" pressure
-		p_g=u_g(4:4:end)*m_dh.parameters.stokesreconditioning;
+		p_g=u_g(4:4:end)*model.dh.parameters.stokesreconditioning;
 	end
 end
@@ -109,4 +106,4 @@
 
 if numrifts,
-	results.riftproperties=OutputRifts(m_dh.loads,m_dh.parameters);
+	results.riftproperties=OutputRifts(model.dh.loads,model.dh.parameters);
 end
Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m	(revision 3872)
+++ /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m	(revision 3873)
@@ -8,7 +8,4 @@
 	m.parameters.Kflag=1; m.parameters.Pflag=1;
 
-	%Update inputs in datasets
-	[m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters);
-	
 	%system matrices
 	[K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
@@ -21,5 +18,5 @@
 	%Reduce load from g size to f size
 	[p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
-	
+
 	%Solve	
 	[u_f]=Solver(K_ff,p_f,[],m.parameters);