function md=cielotransient(md)
%transient solution

	analysis_t='thermaltransient';
	analysis_m='melting';
	analysis_dh='diagnostic_horiz';
	analysis_dbv='diagnostic_basevert';
	analysis_dv='diagnostic_vert';
	analysis_p='prognostic';

	%configure, and initialise data model and parameter defaults:
	SetParameterDefaults;

	%Initialize model velocity at 0 before model processing
	md.vx=zeros(md.numberofgrids,1); md.vy=zeros(md.numberofgrids,1); md.vz=zeros(md.numberofgrids,1);

	%read input data for all analyses
	disp(sprintf('\n   reading data from model %s...',md.name));

	%Build a model for the diagnostic simulation
	m_t=CreateFEMModel(md,analysis_t);
	m_m=CreateFEMModel(md,analysis_m);
	m_dh=CreateFEMModel(md,analysis_dh);
	m_dbv=CreateFEMModel(md,analysis_dbv);
	m_dv=CreateFEMModel(md,analysis_dv);
	m_p=CreateFEMModel(md,analysis_p);

	% buid partitioning vectors. Same partitioning for all analyses. 
	tpart=m_t.tpart;
	part=m_t.part;
	indx=1:6:m_t.uset.gsize; indx=indx(tpart);
	indy=2:6:m_t.uset.gsize; indy=indy(tpart);
	indz=3:6:m_t.uset.gsize; indz=indz(tpart);
	indp=4:6:m_t.uset.gsize; indp=indp(tpart);

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

	nsteps=md.ndt/md.dt;

	%initialize (velocity,pressure) and (melting,temperature)
	solution.u_g=zeros(uset.gsize,1);
	if ~isempty(md.vx), solution.u_g(1:6:uset.gsize)=md.vx(part)/md.yts; end
	if ~isempty(md.vy), solution.u_g(2:6:uset.gsize)=md.vy(part)/md.yts; end
	if ~isempty(md.vz), solution.u_g(3:6:uset.gsize)=md.vz(part)/md.yts; end
	solution.pressure=zeros(uset.gsize,1); solution.pressure(1:6:uset.gsize)=md.rho_ice*md.g*(md.surface(part)-md.z(part)); %lithostatic pressure to spin up

	solution.t_g=zeros(uset.gsize,1);
	if ~isempty(md.temperature), 
		solution.t_g(1:6:uset.gsize)=md.temperature(part);
	else
		solution.t_g(1:6:uset.gsize)=md.meltingpoint;
	end
	if ~isempty(md.melting), solution.melting_g(1:6:uset.gsize)=md.melting(part);end;

	%initialize geometry and time
	solution.thickness=zeros(uset.gsize,1);solution.thickness(1:6:uset.gsize)=md.thickness(part);
	solution.surface=zeros(uset.gsize,1);solution.surface(1:6:uset.gsize)=md.surface(part);
	solution.bed=zeros(uset.gsize,1);solution.bed(1:6:uset.gsize)=md.bed(part);
	solution.time=0;

	%some thermal parameters
	thermalparams.beta=md.beta;
	thermalparams.meltingpoint=md.meltingpoint;
	thermalparams.latentheat=md.latentheat;
	thermalparams.heatcapacity=md.heatcapacity;
	thermalparams.penalty=md.penalty_thermal;

	for n=1:nsteps,

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

		%Deal with temperature first 
		disp(sprintf('%s','   computing temperature...'));

		%Call core thermal computation
		inputs=struct('velocity',solution(n).u_g,'pressure',solution(n).pressure,'temperature',solution(n).t_g,'dt',md.dt,...
					'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);

		[solution(n+1).t_g,m_t]=cielothermal_core(m_t,m_t.params,thermalparams,inputs,analysis_t); 
		
		%Call core melting computation
		disp(sprintf('%s','   computing melting...'));
		inputs=struct('pressure',solution(n).pressure,'temperature',solution(n+1).t_g,'dt',md.dt,...
		'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
	
		solution(n+1).melting_g=cielomelting_core(m_m,m_m.params,thermalparams,inputs,analysis_m); 
		
		%Compute depth averaged temperature for 2d type elements.
		temperature_workspace=solution(n+1).t_g(indx);
		temperature_average2d_workspace=DepthAverage(md,temperature_workspace);
		temperature_average_workspace=project3d(md,temperature_average2d_workspace,'node');
		temperature_average=zeros(uset.gsize,1);temperature_average(1:6:uset.gsize)=temperature_average_workspace(part);
		
		if md.debug,
			t_g_workspace=SwitchPartitioning(solution(n+1).t_g,'workspace',part,tpart,[1]); 
			melting_g_workspace=SwitchPartitioning(solution(n+1).melting_g,'workspace',part,tpart,[1]); 
			plot(md,'data',t_g_workspace(1:6:uset.gsize),'title','Temperature (K)','view',3,'data',melting_g_workspace(1:6:uset.gsize)*md.yts,'title','Melting (m/a)','view',3);pause(2);
		end

		%Deal with velocities.
		disp(sprintf('%s','   computing horizontal velocities...'));

		m_dh.rifts=md.rifts;
		inputs=struct('thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface,...
					'temperature',solution(n+1).t_g,'temperature_average',temperature_average);
		solution(n+1).u_g=cielodiagnostic_core_nonlinear(m_dh,m_dh.params,inputs,analysis_dh); 

		%Extrude velocities for collapsed penta elements
		u_g_workspace=SwitchPartitioning(solution(n+1).u_g,'workspace',part,tpart,[1 2]); u_g_workspace=CieloVelocityExtrude(md,u_g_workspace);
		
		%Plug velocities back into solution, with cluster partitioning: 
		solution(n+1).u_g=SwitchPartitioning(u_g_workspace,'cluster',part,tpart,[1 2]);
		
		if md.debug,
			plot(md,'data',sqrt(u_g_workspace(1:6:uset.gsize).^2+u_g_workspace(2:6:uset.gsize))*md.yts,'view',3,'title','Horizontal velocity (m/a)');pause(2);
		end

		%compute base vertical velocities
		disp(sprintf('%s','   computing vertical base velocity...'));
		inputs=struct('velocity',solution(n+1).u_g,'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
		u_g_basevert=cielodiagnostic_core_linear(m_dbv,m_dbv.params,inputs,analysis_dbv); 

		if md.debug,
			u_g_basevert_workspace=SwitchPartitioning(u_g_basevert,'workspace',part,tpart,[3]); 
			plot(md,'data',u_g_basevert_workspace(3:6:uset.gsize)*md.yts,'layer',1,'title','Basal Vertical velocity (m/a)');pause(2);
		end
	
		%use solution(n+1).u_g_basevert to constrain vertical velocity
		disp(sprintf('%s','   computing vertical velocity...'));
		m_dv.y_g=u_g_basevert;
	
		%reduce constraining vector to s-set (the one we solve on)
		m_dv.y_s = Reducevectorg( m_dv.y_g,m_dv.uset);

		%compute depth averaged horizontal velocity
		velocity_average_workspace=CieloHorizontalVelocityDepthAverage(md,u_g_workspace); velocity_average=SwitchPartitioning(velocity_average_workspace,'cluster',part,tpart,[1 2]);

		%compute vertical velocities
		inputs=struct('velocity',solution(n+1).u_g,'velocity_average',velocity_average,'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
		u_g_vert=cielodiagnostic_core_linear(m_dv,m_dv.params,inputs,analysis_dv); 
		
		
		%add contribution to vertical velocity
		solution(n+1).u_g(3:6:uset.gsize)=u_g_vert(3:6:uset.gsize);
	
		if md.debug,
			u_g_workspace=SwitchPartitioning(solution(n+1).u_g,'workspace',part,tpart,[3]);
			plot(md,'data',u_g_workspace(3:6:uset.gsize)*md.yts,'view',3,'title','Vertical velocity (m/a)');pause(2);
		end

		%compute pressure (for now, lithostatic)
		disp(sprintf('%s','   computing pressure...'));
		solution(n+1).pressure=zeros(uset.gsize,1); solution(n+1).pressure(1:6:uset.gsize)=md.rho_ice*md.g*(md.surface(part)-md.z(part));
		

		%move surface and bed velocities to first lower layer. This will be needed by the prognostic model
		u_g_vert_workspace=SwitchPartitioning(u_g_vert,'workspace',part,tpart,[3]);
		ws_workspace=ShiftLayers(md,u_g_vert_workspace(3:6:uset.gsize),md.numlayers,1);
		wb_workspace=ShiftLayers(md,u_g_vert_workspace(3:6:uset.gsize),1,1);
		
		melting_g_workspace=SwitchPartitioning(solution(n+1).melting_g,'workspace',part,tpart,[1]);
		melting_workspace=ShiftLayers(md,melting_g_workspace(1:6:uset.gsize),1,1);
		
		accumulation_workspace=ShiftLayers(md,md.accumulation,md.numlayers,1);
		

		%plug into g set
		ws_g=zeros(uset.gsize,1);ws_g(1:6:uset.gsize)=ws_workspace(part);
		wb_g=zeros(uset.gsize,1);wb_g(1:6:uset.gsize)=wb_workspace(part);
		melting_g=zeros(uset.gsize,1);melting_g(1:6:uset.gsize)=melting_workspace(part);
		accumulation_g=zeros(uset.gsize,1);accumulation_g(1:6:uset.gsize)=accumulation_workspace(part)/md.yts;
	
		
		if md.debug,
			plot(md,'data',ws_workspace*md.yts,'title','Surface velocity projected on layer 1','view',3,'data',wb_workspace*md.yts,'title','Basal velocity projected on layer 1','view',3,...
			'data',melting_workspace*md.yts,'title','Bed melting projected on layer 1 (m/yr)','view',3,'data',accumulation_workspace,'title','Surface accumulation projected on layer 1 (m/yr)','view',3);pause(2);
		end
		
		%compute new thickness
		disp(sprintf('%s','   computing new thickness...'));
		inputs=struct('thickness',solution(n).thickness,'melting',melting_g,'accumulation',accumulation_g,...
					   'surface_vertical_velocity',ws_g,'basal_vertical_velocity',wb_g,'dt',md.dt,'velocity_average',velocity_average);
		new_thickness=cieloprognostic_core(m_p,m_p.params,inputs,analysis_p);
	
		%project collapsed thickness onto 3d mesh
		new_thickness_workspace=SwitchPartitioning(new_thickness,'workspace',part,tpart,[1]); 
		new_thickness_workspace=project3d(md,project2d(md,new_thickness_workspace(1:6:end),1),'node');

		
		%update surface and bed using the new thickness
		[new_bed_workspace,new_surface_workspace,]=UpdateGeometry(md,new_thickness_workspace,solution(n).bed(indx),solution(n).surface(indx));

		if md.debug,
			plot(md,'data',new_thickness_workspace,'view',3,'title','New thickness (m)','data',new_surface_workspace,'view',3,'title','New surface (m)','data',new_bed_workspace,'view',3,'title','New bed (m)');pause(2);
		end

		%project onto g set
		solution(n+1).thickness=zeros(uset.gsize,1);solution(n+1).thickness(1:6:uset.gsize)=new_thickness_workspace(part);
		solution(n+1).bed=zeros(uset.gsize,1);solution(n+1).bed(1:6:uset.gsize)=new_bed_workspace(part);
		solution(n+1).surface=zeros(uset.gsize,1);solution(n+1).surface(1:6:uset.gsize)=new_surface_workspace(part);

		%update grid positions
		[m_t.bgpdt,m_t.bgpdtb]=CieloUpdateGridPosition(m_t.bgpdt,m_t.bgpdtb,solution(n+1).bed,solution(n+1).thickness);

		[m_dh.bgpdt,m_dh.bgpdtb]=CieloUpdateGridPosition(m_dh.bgpdt,m_dh.bgpdtb,solution(n+1).bed,solution(n+1).thickness);

		[m_dbv.bgpdt,m_dbv.bgpdtb]=CieloUpdateGridPosition(m_dbv.bgpdt,m_dbv.bgpdtb,solution(n+1).bed,solution(n+1).thickness);

		[m_dv.bgpdt,m_dv.bgpdtb]=CieloUpdateGridPosition(m_dv.bgpdt,m_dv.bgpdtb,solution(n+1).bed,solution(n+1).thickness);

		[m_p.bgpdt,m_p.bgpdtb]=CieloUpdateGridPosition(m_p.bgpdt,m_p.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
	end



%load results onto model
for i=1:length(solution),
	solution2(i).vx=solution(i).u_g(indx)*md.yts;
	solution2(i).vy=solution(i).u_g(indy)*md.yts;
	solution2(i).vz=solution(i).u_g(indz)*md.yts;
	solution2(i).vel=sqrt(solution2(i).vx.^2+solution2(i).vy.^2+solution2(i).vz.^2);
	solution2(i).pressure=solution(i).pressure(indx)/10^5;
	solution2(i).temperature=solution(i).t_g(indx);
	solution2(i).melting=solution(i).melting_g(indx);
	solution2(i).thickness=solution(i).thickness(indx);
	solution2(i).surface=solution(i).surface(indx);
	solution2(i).bed=solution(i).bed(indx);
	solution2(i).time=solution(i).time/md.yts;
end
md.transient_results=solution2;
