function md=cielothermal(md,solutiontype)
%thermal  solution: steady state and transient

	analysis_t=solutiontype;
	analysis_m='melting';
	
	%configure, and initialise data model and parameter defaults:
	SetParameterDefaults;

	%Build thermal and melting models
	m_t=CreateFEMModel(md,analysis_t);
	m_m=CreateFEMModel(md,analysis_m);


	%Build partitioning vectors. 
	indx=1:6:mt.uset.gsize; indx=indx(tpart);

	% figure out number of dof: just for information purposes.
	md.dof=m_t.uset.fsize; %biggest dof number

	%initialize velocity,pressure and material properties
	velocity=zeros(m_t.uset.gsize,1);
	velocity(1:6:m_t.uset.gsize)=md.vx(part)/md.yts;
	velocity(2:6:m_t.uset.gsize)=md.vy(part)/md.yts;
	velocity(3:6:m_t.uset.gsize)=md.vz(part)/md.yts;

	pressure=zeros(m_t.uset.gsize,1);
	pressure(1:6:m_t.uset.gsize)=md.pressure(part);

	thermalparams.beta=md.beta;
	thermalparams.meltingpoint=md.meltingpoint;
	thermalparams.latentheat=md.latentheat;
	thermalparams.heatcapacity=md.heatcapacity;
	thermalparams.penalty=md.penalty_thermal;
	thermalparams.min_thermal_constraints=md.min_thermal_constraints;
	thermalparams.mintemp=min(md.dirichletvalues_thermal(find(md.gridondirichlet_thermal)));

	if strcmpi(analysis_t,'thermalsteady'),

		inputs=struct('pressure',pressure,'velocity',velocity,'dt',md.dt);

		%Call core thermal computation
		disp(sprintf('\n%s\n','computing temperature...'));
		[t_g,m_t]=cielothermal_core(m_t,m_t.params,thermalparams,inputs,analysis_t);

		%Call core melting computation: 
		%disp(sprintf('\n%s\n','computing melting...'));
		%SetUset(m_m);
		%inputs=struct('pressure',pressure,'temperature',rt_g,'dt',md.dt);
		%rmelting_g=cielomelting_core(m_m,m_m.params,thermalparams,inputs,analysis_m);

		%load results onto model
		md.temperature=t_g(indx);
		%melting_g=IMdb('select matrix from rmelting_g');
		%md.melting=melting_g(indx)*md.yts; %from m/s to m/a
	else

		nsteps=md.ndt/md.dt;

		%initialize temperature and melting
		temperature=zeros(m_t.uset.gsize,1);
		temperature(1:6:m_t.uset.gsize)=md.temperature(part);
		soln.t_g=temperature;
		%melting=zeros(uset.gsize,1);
		%melting(1:6:uset.gsize)=md.melting(part)/md.yts; %in m/s
		%soln.melting_g=melting;

		for n=1:nsteps, 

			disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt));
			%Call core thermal computation
			disp('   computing temperature...');
			inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'velocity',velocity,'dt',md.dt);
			[soln(n+1).t_g,m_t]=cielothermal_core(m_t,m_t.params,thermalparams,inputs,analysis_t);
			
			%Call core melting computation
			%disp('   computing melting...');
			%SetUset(m_m);
			%inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'dt',md.dt);
			%rmelting_g=cielomelting_core(m_m,m_m.params,thermalparams,inputs,analysis_m);
			%soln(n+1).melting_g=IMdb('select matrix from rmelting_g');
			
		end
		
		%Wrap up
		solution_temperature=struct('temperature',{});
		for n=1:nsteps+1,
			solution_temperature(n).temperature=soln(n).t_g(indx);
		end
		md.temperature=solution_temperature;
		
		%solution_melting=struct('melting',{});
		%for n=1:nsteps+1,
		%	solution_melting(n).melting=soln(n).melting_g(indx)*md.yts; %in m/a
		%end
		%md.melting=solution_melting;
	end
	
