function md=icethermal(md,analysis_type)
%ICETHERMAL - thermal matlab driven solution
%
%   Usage:
%      md=icethermal(md,analysis_type)

%define global variables
iceglobal

%determine if run is parallel
if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end;

%for now, only serial support is in
if cluster,
	error('icethermal error message: parallel support not implemented yet');
end

%First, build elements,grids,loads, etc ... for thermal and melting 
m_t=CreateFemModel(md,analysis_type);
if strcmpi(analysis_type,'thermalsteady'),
	m_m=CreateFemModel(md,'meltingsteady');
else
	m_m=CreateFemModel(md,'meltingtransient');
end
%initialize velocity,pressure and material properties
velocity=zeros(m_t.gridset.gsize,1);
velocity(1:6:m_t.gridset.gsize)=md.vx(m_t.part)/md.yts;
velocity(2:6:m_t.gridset.gsize)=md.vy(m_t.part)/md.yts;
velocity(3:6:m_t.gridset.gsize)=md.vz(m_t.part)/md.yts;

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

if strcmpi(analysis_type,'thermalsteady'),

	%Call core thermal computation
	disp('   computing temperature...');
	gridset=m_t.gridset;
	inputs=struct('pressure',pressure,'velocity',velocity,'dt',md.dt);
	[t_g loads_t melting_offset]=icethermal_core(m_t,analysis_type,inputs);
	
	%Call core melting computation
	disp('   computing melting...');
	gridset=m_m.gridset;
	inputs=struct('pressure',pressure,'temperature',t_g,'melting_offset',melting_offset);
	melting_g=icemelting_core(m_m,analysis_type,inputs);
	
	%Wrap up
	indx=1:6:gridset.gsize;
	indx=indx(m_t.tpart);
	md.temperature=t_g(indx);
	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(gridset.gsize,1);
	melting=zeros(gridset.gsize,1);
	temperature(1:6:gridset.gsize)=md.temperature;
	melting(1:6:gridset.gsize)=md.melting/md.yts; %in m/s
	soln.t_g=temperature;
	soln.melting_g=melting;
	soln.time=0;


	for n=1:nsteps, 

		soln(n+1).time=(n)*md.dt;
		disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt));

		%Call core thermal computation
		disp('   computing temperature...');
		gridset=m_t.gridset;
		inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'velocity',velocity,'dt',md.dt);
		[soln(n+1).t_g loads_t  melting_offset]=icethermal_core(m_t,analysis_type,inputs);

		%Call core melting computation
		disp('   computing melting...');
		gridset=m_m.gridset;
		inputs=struct('pressure',pressure,'temperature',soln(n).t_g,'melting_offset',melting_offset,'dt',md.dt);
		soln(n+1).melting_g=icemelting_core(m_m,analysis_type,inputs);
		
	end
	
	%Wrap up
	indx=1:6:m_t.gridset.gsize;
	indx=indx(m_t.tpart);
 	solution=struct('time',{},'temperature',{},'melting',{});

	for n=1:nsteps+1,
		solution(n).temperature=soln(n).t_g(indx);
		solution(n).melting=soln(n).melting_g(indx)*md.yts; %in m/a
		solution(n).time=soln(n).time/md.yts; %in year
	end
	md.thermaltransient_results=solution;
end
