Index: /issm/trunk/src/m/solutions/cielo/SpawnCore.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/SpawnCore.m	(revision 959)
+++ /issm/trunk/src/m/solutions/cielo/SpawnCore.m	(revision 959)
@@ -0,0 +1,61 @@
+function responses=SpawnCore(models,inputs,variables,variabledescriptors,analysis_type,sub_analysis_type,counter);
+%SPAWNCORE - for Qmu analysis, using Dakota. Spawn the core solution.
+%
+%   Usage:
+%      responses=SpawnCore(models,inputs,variables,variabledescriptors)
+%
+
+%recover parameters
+debug=models.dh.parameters.debug;
+parameters=models.dh.parameters;
+responsedescriptors=models.dh.parameters.responsedescriptors; 
+npart=models.dh.parameters.qmu_npart;
+part=models.dh.parameters.qmu_part+1; %C indexing
+numberofnodes=models.dh.parameters.numberofnodes;
+
+disp(['   qmu iteration: ',num2str(counter)]);
+
+%first update the inputs to the models using the variables provided to us by dakota.
+count=1;
+while count<=numel(variables)
+	
+	descriptor=variabledescriptors{count};
+	if ~qmuisdistributed(descriptor),
+		inputs=add(inputs,descriptor,variables(count),'double');
+
+		count=count+1;
+	else
+		root=qmuroot(descriptor);
+		param=parameters.(root);
+
+		%next npart values in variables are partition values for this param, collect them.
+		partition=variables(count:count+npart-1);
+
+		%update parameter:
+		param=param.*partition(part);
+
+		%add parameter to inputs
+		inputs=add(inputs,root,param,'doublevec',1,numberofnodes);
+
+		%skip next npart iterations, they all deal with the same parameter descriptor
+		count=count+npart;
+	end
+end
+
+
+%now run the core solution
+if analysis_type==diagnosticenum(),
+
+	results=diagnostic_core(models,inputs);
+
+else
+	error(['SpawnCore error message: could not find core solutoin for analysis type: ' analysis_type]);
+end
+
+%now process the results to get response function values
+responses=zeros(numel(responsedescriptors),1);
+for i=1:numel(responsedescriptors),
+	descriptor=responsedescriptors{i};
+	responses(i)=qmuresponse(models,results,descriptor);
+end
+
Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 958)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 959)
@@ -33,8 +33,13 @@
 
 	%compute solution
-	results=diagnostic_core(models,inputs);
+	if ~models.dh.parameters.qmu_analysis,
+		%launch core of diagnostic solution.
+		results=diagnostic_core(models,inputs);
+	else
+		%launch dakota driver for diagnostic core solution
+		Qmu(models,inputs,models.dh.parameters);
+	end
 
-	%load results onto model
-	md=loadresults(md,models,results);
+	md.results.diagnostic=processresults(models,results);
 
 	%stop timing
Index: sm/trunk/src/m/solutions/cielo/loadresults.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/loadresults.m	(revision 958)
+++ 	(revision )
@@ -1,53 +1,0 @@
-function md=loadresults(md,models,rawresults);
-%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,models,rawresults)
-
-%recover models
-m_dh=models.dh;
-m_ds=models.ds;
-m_dhu=models.dhu;
-
-%recover parameters
-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;
-
-%recover results
-u_g=rawresults.u_g;
-p_g=rawresults.p_g;
-
-if ~isstruct(md.results), md.results=struct(); end 
-if dim==2,
-	results.step=1;
-	results.time=1;
-	results.vx=u_g(1:2:end)*md.yts;
-	results.vy=u_g(2:2:end)*md.yts;
-	results.vel=sqrt(results.vx.^2+results.vy.^2);
-	results.pressure=p_g;
-	md.results.diagnostic=results;
-
-else
-	
-	results.step=1;
-	results.time=1;
-	if isstokes,
-		results.vx=u_g(1:4:end)*md.yts;
-		results.vy=u_g(2:4:end)*md.yts;
-		results.vz=u_g(3:4:end)*md.yts;
-		results.pressure=p_g;
-	else
-		results.vx=u_g(1:3:end)*md.yts;
-		results.vy=u_g(2:3:end)*md.yts;
-		results.vz=u_g(3:3:end)*md.yts;
-		results.pressure=p_g;
-	end
-	results.vel=sqrt(results.vx.^2+results.vy.^2+results.vz.^2);
-	md.results.diagnostic=results;
-end
Index: /issm/trunk/src/m/solutions/cielo/processresults.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 959)
+++ /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 959)
@@ -0,0 +1,52 @@
+function newresults=processresults(models,results);
+%PROCESSRESULTS - process results from core solution
+%
+%   The solution (vel,pressure,...) are on the g-set,
+%   use index to get the value on each node
+%
+%   Usage:
+%      results=processresults(md,results,m_dh,m_ds,m_dhu)
+
+%recover models
+m_dh=models.dh;
+m_ds=models.ds;
+m_dhu=models.dhu;
+
+%recover parameters
+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;
+yts=m_ds.parameters.yts;
+
+%recover results
+u_g=results.u_g;
+p_g=results.p_g;
+
+newresults=struct();
+if dim==2,
+	newresults.step=1;
+	newresults.time=1;
+	newresults.vx=u_g(1:2:end)*yts;
+	newresults.vy=u_g(2:2:end)*yts;
+	newresults.vel=sqrt(newresults.vx.^2+newresults.vy.^2);
+	newresults.pressure=p_g;
+
+else
+	
+	newresults.step=1;
+	newresults.time=1;
+	if isstokes,
+		newresults.vx=u_g(1:4:end)*yts;
+		newresults.vy=u_g(2:4:end)*yts;
+		newresults.vz=u_g(3:4:end)*yts;
+		newresults.pressure=p_g;
+	else
+		newresults.vx=u_g(1:3:end)*yts;
+		newresults.vy=u_g(2:3:end)*yts;
+		newresults.vz=u_g(3:3:end)*yts;
+		newresults.pressure=p_g;
+	end
+	newresults.vel=sqrt(newresults.vx.^2+newresults.vy.^2+newresults.vz.^2);
+end