Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 1786)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 1787)
@@ -33,4 +33,5 @@
 int SteadyAnalysisEnum(void){           return          231; }
 int TransientAnalysisEnum(void){        return          232; }
+int ThermalstaticAnalysisEnum(void){    return          233; }
 //slope
 int SlopeComputeAnalysisEnum(void){     return          240; }
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 1786)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 1787)
@@ -34,4 +34,5 @@
 int SteadyAnalysisEnum(void);
 int TransientAnalysisEnum(void);
+int ThermalstaticAnalysisEnum(void);
 //slope
 int SlopeComputeAnalysisEnum(void);
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 1786)
+++ /issm/trunk/src/c/Makefile.am	(revision 1787)
@@ -604,5 +604,5 @@
 bin_PROGRAMS = 
 else 
-dnl bin_PROGRAMS = diagnostic.exe  control.exe thermal.exe prognostic.exe transient.exe
+bin_PROGRAMS = diagnostic.exe  control.exe thermal.exe prognostic.exe transient.exe
 endif
 
Index: /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp
===================================================================
--- /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp	(revision 1786)
+++ /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp	(revision 1787)
@@ -40,4 +40,5 @@
 		/*Do nothing, no constraints management!:*/
 		num_unstable_constraints=0;
+		converged=1;
 	}
 
Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 1786)
+++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 1787)
@@ -217,4 +217,28 @@
 		bool=0; return;
 	end
+end
+
+%THERMALSTATIC
+if md.analysis_type==ThermalstaticAnalysisEnum,
+	%PRESSURE
+	if isnans(md.pressure),
+		disp(['For a thermalstatic computation, the model must have an initial pressure, even lithostatic will do.']);
+		bool=0;return;
+	end
+
+	%eps: 
+	if isnan(md.eps_rel),
+		disp(['For a thermalstatic computation, eps_rel (relative convergence criterion) must be defined!']);
+		bool=0;return;
+	end
+
+	%dim: 
+	if strcmpi(md.type,'2d'),
+		disp(['For a thermalstatic computation, model needs to be 3d']);
+		bool=0;return;
+	end
+
+
+
 end
 
Index: /issm/trunk/src/m/classes/public/process_solve_options.m
===================================================================
--- /issm/trunk/src/m/classes/public/process_solve_options.m	(revision 1786)
+++ /issm/trunk/src/m/classes/public/process_solve_options.m	(revision 1787)
@@ -73,5 +73,5 @@
 
 %check solution type is supported
-if ~ismemberi(analysis_type,{'control','diagnostic','prognostic','thermal','parameters','mesh2grid','transient'}),
+if ~ismemberi(analysis_type,{'control','diagnostic','prognostic','thermal','thermalstatic','parameters','mesh2grid','transient'}),
 	error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']);
 else
Index: /issm/trunk/src/m/classes/public/solve.m
===================================================================
--- /issm/trunk/src/m/classes/public/solve.m	(revision 1786)
+++ /issm/trunk/src/m/classes/public/solve.m	(revision 1787)
@@ -76,4 +76,7 @@
 	md=parameters(md);
 
+elseif md.analysis_type==ThermalStaticAnalysisEnum,
+	md=thermalstatic(md);
+
 else
 	error('solution type not supported for this model configuration.');
Index: sm/trunk/src/m/classes/public/transfervel.m
===================================================================
--- /issm/trunk/src/m/classes/public/transfervel.m	(revision 1786)
+++ 	(revision )
@@ -1,16 +1,0 @@
-function md=transfervel(md,string)
-%TRANSFERVEL transfer results for velocity from results to md.vx and md.vy fields.
-%
-%    Usage: md=transfervel(md,string)
-%
-%
-%    Example: md=transfervel(md,'diagnostic');
-%             md=transfervel(md,'control');
-
-if strcmpi(string,'diagnostic'),
-	md.vx=md.results.diagnostic.vx;
-	md.vy=md.results.diagnostic.vy;
-	md.vel=md.results.diagnostic.vel;
-else 
-	error(['transfervel error message: analysis ' string ' not supported yet!']);
-end
Index: /issm/trunk/src/m/classes/public/tres.m
===================================================================
--- /issm/trunk/src/m/classes/public/tres.m	(revision 1786)
+++ /issm/trunk/src/m/classes/public/tres.m	(revision 1787)
@@ -13,5 +13,5 @@
 	md.vel=md.results.diagnostic.vel;
 	md.pressure=md.results.diagnostic.pressure;
-else if strcmpi(string,'thermalstatic'),
+elseif strcmpi(string,'thermalstatic'),
 	md.vx=md.results.thermalstatic.vx;
 	md.vy=md.results.thermalstatic.vy;
@@ -20,5 +20,5 @@
 	md.temperature=md.results.thermalstatic.temperature;
 	md.melting=md.results.thermalstatic.melting;
-else if strcmpi(string,'thermal'),
+elseif strcmpi(string,'thermal'),
 	md.temperature=md.results.thermalstatic.temperature;
 	md.melting=md.results.thermalstatic.melting;
Index: /issm/trunk/src/m/enum/AnalysisTypeAsEnum.m
===================================================================
--- /issm/trunk/src/m/enum/AnalysisTypeAsEnum.m	(revision 1786)
+++ /issm/trunk/src/m/enum/AnalysisTypeAsEnum.m	(revision 1787)
@@ -65,4 +65,8 @@
 end
 
+if enum==ThermalstaticAnalysisEnum(),
+	string='thermalstatic';
+end
+
 if enum==SlopeComputeAnalysisEnum(),
 	string='slopecompute';
Index: /issm/trunk/src/m/enum/ThermalstaticAnalysisEnum.m
===================================================================
--- /issm/trunk/src/m/enum/ThermalstaticAnalysisEnum.m	(revision 1787)
+++ /issm/trunk/src/m/enum/ThermalstaticAnalysisEnum.m	(revision 1787)
@@ -0,0 +1,9 @@
+function macro=ThermalstaticAnalysisEnum()
+%THERMALSTATICANALYSISENUM - Enum of ThermalstaticAnalysis
+%
+%   file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum
+%
+%   Usage:
+%      macro=ThermalstaticAnalysisEnum()
+
+macro=233;
Index: /issm/trunk/src/m/solutions/cielo/processresults.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 1786)
+++ /issm/trunk/src/m/solutions/cielo/processresults.m	(revision 1787)
@@ -15,5 +15,5 @@
 
 %recover models first
-if (analysis_type==DiagnosticAnalysisEnum | analysis_type==TransientAnalysisEnum),
+if (analysis_type==DiagnosticAnalysisEnum | analysis_type==TransientAnalysisEnum | analysis_type==ThermalstaticAnalysisEnum()),
 	m_dh=models.dh;
 	m_ds=models.ds;
@@ -35,5 +35,5 @@
 	isstokes=m_dh.parameters.isstokes;
 end
-if (analysis_type==ThermalAnalysisEnum() | analysis_type==TransientAnalysisEnum()),
+if (analysis_type==ThermalAnalysisEnum() | analysis_type==TransientAnalysisEnum() | analysis_type==ThermalstaticAnalysisEnum()),
 	m_m=models.m;
 end
Index: /issm/trunk/src/m/solutions/cielo/thermalstatic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermalstatic.m	(revision 1787)
+++ /issm/trunk/src/m/solutions/cielo/thermalstatic.m	(revision 1787)
@@ -0,0 +1,61 @@
+function md=thermalstatic(md);
+%THERMALSTATIC - compute the velocity and temperature field of a model in steady state.
+%
+%   Usage:
+%      md=thermalstatic(md)
+%
+
+	%timing
+	t1=clock;
+
+	models.analysis_type=ThermalstaticAnalysisEnum; %needed for processresults
+	
+	%Build all models requested for diagnostic simulation
+	displaystring(md.debug,'%s',['reading diagnostic horiz model data']);
+	md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HorizAnalysisEnum; models.dh=CreateFemModel(md);
+	
+	displaystring(md.debug,'\n%s',['reading diagnostic vert model data']);
+	md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=VertAnalysisEnum; models.dv=CreateFemModel(md);
+	
+	displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']);
+	md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=StokesAnalysisEnum; models.ds=CreateFemModel(md);
+	
+	displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']);
+	md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HutterAnalysisEnum; models.dhu=CreateFemModel(md);
+	
+	displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']);
+	md.analysis_type=SlopeComputeAnalysisEnum; md.sub_analysis_type=NoneAnalysisEnum; models.sl=CreateFemModel(md);
+
+	%Build all models requested for thermal simulation
+	displaystring(md.debug,'%s',['reading thermal model data']);
+	md.analysis_type=ThermalAnalysisEnum(); md.sub_analysis_type=SteadyAnalysisEnum(); models.t=CreateFemModel(md);
+
+	displaystring(md.debug,'%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);
+
+	%initialize inputs
+	displaystring(md.debug,'\n%s',['setup inputs...']);
+	inputs=inputlist;
+	inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,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
+	if ~models.dh.parameters.qmu_analysis,
+		%launch core of thermalstatic solution.
+		results=thermalstatic_core(models,inputs);
+	
+		%process results
+		if ~isstruct(md.results), md.results=struct(); end
+		md.results.thermalstatic=processresults(models,results);
+	else
+		%launch dakota driver for thermalstatic core solution
+		Qmu(models,inputs,models.dh.parameters);
+	end
+
+	%stop timing
+	t2=clock;
+	displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);
Index: /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m	(revision 1787)
+++ /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m	(revision 1787)
@@ -0,0 +1,90 @@
+function results=thermalstatic_core(models,inputs);
+%THERMALSTATIC_CORE - compute the core temperature and velocity field  at thermal steady state.
+%
+%   Usage:
+%      results=thermalstatic_core(models,inputs);
+%
+
+%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
+debug=m_dhu.parameters.debug;debug=1;
+dim=m_dhu.parameters.dim;
+eps_rel=m_dhu.parameters.eps_rel;
+ishutter=m_dhu.parameters.ishutter;
+ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
+isstokes=m_ds.parameters.isstokes;
+
+%convergence
+converged=0;
+
+step=1;
+
+%initialize: 
+t_g_old=0;
+u_g_old=0;
+
+if isstokes, ndof=4; else ndof=3; end
+
+while(~converged),
+	
+	displaystring(debug,'%s%i','computing temperature and velocity for step ',step);
+	
+	%first compute temperature at steady state.
+	if (step>1),
+		inputs=add(inputs,'velocity',results.u_g,'doublevec',ndof,m_t.parameters.numberofnodes);
+	end
+	results_thermal=thermal_core(models,inputs);
+	
+	%save results from thermal
+	results.t_g=results_thermal.t_g;
+	results.m_g=results_thermal.m_g;
+
+	%add temperature to inputs.
+	%Compute depth averaged temperature
+	temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,results.t_g,'temperature');
+	inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,m_t.parameters.numberofnodes);
+	inputs=add(inputs,'temperature',results.t_g,'doublevec',1,m_t.parameters.numberofnodes);
+	
+	%now compute diagnostic velocity using the steady state temperature.
+	results_diagnostic=diagnostic_core(models,inputs);
+	
+	%save results from thermal
+	results.u_g=results_diagnostic.u_g;
+	results.p_g=results_diagnostic.p_g;
+
+	%convergence? 
+	du_g=results.u_g-u_g_old;
+	ndu=norm(du_g,2); 
+	nu=norm(u_g_old,2);
+	
+	dt_g=results.t_g-t_g_old;
+	ndt=norm(dt_g,2); 
+	nt=norm(t_g_old,2); 
+
+	u_g_old=results.u_g;
+	t_g_old=results.t_g;
+	
+	displaystring(debug,'%-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;
+	end
+
+	step=step+1;
+	if converged,
+		break;
+	end
+end
+
+results.step=step;
+results.time=0;