Index: /issm/trunk/src/m/solutions/cielo/control.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/control.m	(revision 1849)
+++ /issm/trunk/src/m/solutions/cielo/control.m	(revision 1850)
@@ -27,4 +27,11 @@
 	md.analysis_type=SlopeComputeAnalysisEnum(); md.sub_analysis_type=NoneAnalysisEnum(); models.sl=CreateFemModel(md);
 
+	displaystring(md.debug,'\n%s',['reading thermal model data']);
+	md.analysis_type=ThermalAnalysisEnum(); md.sub_analysis_type=SteadyAnalysisEnum(); models.t=CreateFemModel(md);
+
+	displaystring(md.debug,'\n%s',['reading melting model data']);
+	md.analysis_type=MeltingAnalysisEnum(); md.sub_analysis_type=SteadyAnalysisEnum(); models.m=CreateFemModel(md);
+
+
 	% figure out number of dof: just for information purposes.
 	md.dof=modelsize(models);
@@ -34,4 +41,6 @@
 	inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes);
 	inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes);
+	inputs=add(inputs,'pressure',models.t.parameters.p_g,'doublevec',1,models.t.parameters.numberofnodes);
+	inputs=add(inputs,'dt',models.t.parameters.dt*models.t.parameters.yts,'double');
 
 	%compute solution
Index: /issm/trunk/src/m/solutions/cielo/gradjcompute_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/gradjcompute_core.m	(revision 1850)
+++ /issm/trunk/src/m/solutions/cielo/gradjcompute_core.m	(revision 1850)
@@ -0,0 +1,51 @@
+function results=gradjcompute_core(models,inputs);
+
+%recover active models
+m=models.active;
+
+%recover parameters
+debug=m.parameters.debug;
+dim=m.parameters.dim;
+extrude_param=m.parameters.extrude_param;
+ishutter=m.parameters.ishutter;
+ismacayealpattyn=m.parameters.ismacayealpattyn;
+isstokes=m.parameters.isstokes;
+analysis_type=m.parameters.analysis_type;
+sub_analysis_type=m.parameters.sub_analysis_type;
+
+%Recover solution for this stiffness and right hand side: 
+displaystring(debug,'%s','         computing velocities...');
+[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
+inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
+
+
+%Buid Du, difference between observed velocity and model velocity.
+displaystring(debug,'%s','          computing Du...');
+[Du_g]=Du(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+
+%Reduce adjoint load from g-set to f-set
+[Du_f] = Reduceloadfromgtof( Du_g, m.Gmn, K_fs0, m.ys0, m.nodesets);
+
+%Solve for adjoint vector: 
+displaystring(debug,'%s','          computing adjoint state...');
+lambda_f=Solver(K_ff0,Du_f,[],m.parameters);
+
+%Merge back to g set
+lambda_g= Mergesolutionfromftog( lambda_f, m.Gmn, m.ys0, m.nodesets ); 
+inputs=add(inputs,'adjoint',lambda_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
+
+%Compute gradJ 
+grad_g=Gradj(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+if (dim==3 & extrude_param),
+	displaystring(debug,'%s','          extruding gradient...');
+	grad_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,grad_g,'gradj',0);
+end
+
+%compute initial velocity from diagnostic_core (horiz+vertical)
+displaystring(debug,'%s','          compute initial velocity...');
+results_diag=diagnostic_core(models,inputs);
+u_g=results_diag.u_g;
+
+%Assign output
+results.u_g=u_g;
+results.grad_g=grad_g;
Index: /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m	(revision 1849)
+++ /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m	(revision 1850)
@@ -1,8 +1,12 @@
-function J =objectivefunctionC(search_scalar,m,inputs,p_g,grad_g,n,analysis_type,sub_analysis_type);
-        
+function J =objectivefunctionC(search_scalar,models,inputs,p_g,grad_g,n,analysis_type,sub_analysis_type);
+
+%recover active model.
+m=models.active;
+
 %recover some parameters
 optscal=m.parameters.optscal(n);
 fit=m.parameters.fit(n);
 control_type=m.parameters.control_type;
+thermalstatic=m.parameters.thermalstatic;
 
 %Update along gradient using scalar supplied by fmincon optimization routine
@@ -13,5 +17,18 @@
 
 %Run diagnostic with updated parameters.
-u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
+if ~thermalstatic, 
+	%do a simple diagnostic, with the current temperature profile, do not look for steady state.
+	u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
+else
+	%do a full thermal mechanical steady state converged computation, much slower!
+	results=thermalstatic_core(models,inputs); u_g=results.u_g; 
+
+	%u_g ships with 3 or 4 dofs, we only want the horizontal ones!
+	if ~models.dh.parameters.isstokes, 
+		u_g=u_g(dofsetgen([1,2],3,m.parameters.numberofnodes*3));
+	end
+end
+
+%add velocity to inputs.
 inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
 
Index: /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m	(revision 1849)
+++ /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m	(revision 1850)
@@ -16,5 +16,5 @@
 
 %recover parameters common to all solutions
-debug=m_dhu.parameters.debug;debug=1;
+debug=m_dhu.parameters.debug;
 dim=m_dhu.parameters.dim;
 eps_rel=m_dhu.parameters.eps_rel;
@@ -40,5 +40,5 @@
 	%first compute temperature at steady state.
 	if (step>1),
-		inputs=add(inputs,'velocity',results.u_g,'doublevec',ndof,m_t.parameters.numberofnodes);
+		inputs=add(inputs,'velocity',results_diagnostic.u_g,'doublevec',ndof,m_t.parameters.numberofnodes);
 	end
 	results_thermal=thermal_core(models,inputs);