function md=transient2d(md);
%ICETRANSIENT2D - transient solution of 2d models
%
%   Usage:
%      md=transient2d(md)

global gridset

%Create fem structure (input of diagnostic2d)
fem=struct();

%Figure out which type of elements are present
[fem.ishutter,fem.ismacayealpattyn,fem.isstokes]=DiagnosticSolutionType(md.elements_type);

%First, build elements,grids,loads, etc ... for horizontal model
if fem.ishutter,
	fem.m_ss=CreateFemModel(md,'surface_slope_compute');
	fem.m_dhu=CreateFemModel(md,'diagnostic_hutter');
end
if fem.ismacayealpattyn,
	fem.m_dh=CreateFemModel(md,'diagnostic_horiz');
end
m_p=CreateFemModel(md,'prognostic');

%initialize (velocity,pressure)
solution.u_g=zeros(m_p.gridset.gsize,1);
if ~isnan(md.vx), solution.u_g(1:6:m_p.gridset.gsize)=md.vx/md.yts; end
if ~isnan(md.vy), solution.u_g(2:6:m_p.gridset.gsize)=md.vy/md.yts; end

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

%melting and accumulation are taken constant for all iterations
melting=zeros(m_p.gridset.gsize,1);
accumulation=zeros(m_p.gridset.gsize,1);
melting(1:6:end)=md.melting/md.yts;           %from m/yr to m/s
accumulation(1:6:end)=md.accumulation/md.yts; %from m/yr to m/s

%build constant inputs
fem.inputs=struct();
fem.inputs.accumulation=accumulation;
fem.inputs.melting=melting;
fem.inputs.dt=md.dt;

%first time step is given by model. 
dt=md.dt;
finaltime=md.ndt;
time=dt;
n=1; %counter

while  time<finaltime+dt, %make sure we run up to finaltime.

	disp(sprintf('\n%s%g%s%g%s%g\n','time [yr]: ',time/md.yts,'    iteration number: ',n,'/',floor(finaltime/dt)));
	
	solution(n+1).time=time;

	%update inputs
	fem.inputs.velocity=solution(n).u_g;
	fem.inputs.thickness=solution(n).thickness;
	fem.inputs.bed=solution(n).bed;
	fem.inputs.surface=solution(n).surface;

	%Deal with velocities.

	%Get horizontal solution. 
	solution(n+1).u_g=diagnostic2d(md,fem);

	%compute new thickness
	disp(sprintf('%s','   computing new thickness...'));
	fem.inputs.velocity_average=solution(n+1).u_g;
	new_thickness=prognostic_core(m_p,'prognostic',fem.inputs);

	%update surface and bed using the new thickness
	disp(sprintf('%s','   updating geometry...'));
	indx=1:6:m_p.gridset.gsize; indx=indx(m_p.tpart);
	[new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness(indx),solution(n).bed(indx),solution(n).surface(indx));
	
	%Record bed surface and thickness in the solution
	solution(n+1).thickness=zeros(m_p.gridset.gsize,1);solution(n+1).thickness(1:6:m_p.gridset.gsize)=new_thickness;
	solution(n+1).bed=zeros(m_p.gridset.gsize,1);solution(n+1).bed(1:6:m_p.gridset.gsize)=new_bed;
	solution(n+1).surface=zeros(m_p.gridset.gsize,1);solution(n+1).surface(1:6:m_p.gridset.gsize)=new_surface;

	%figure out if time stepping is good
	disp(sprintf('%s','   checking time stepping...'));
	[back,dt,time]=TimeStepping(md,solution,dt,time);
	if back,
		continue;
	end

	%update time and counter
	time=time+dt;
	n=n+1;

end

%load results onto model
indx=1:6:m_p.gridset.gsize; indx=indx(m_p.tpart);
indy=2:6:m_p.gridset.gsize; indy=indy(m_p.tpart);
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).vel=sqrt(solution2(i).vx.^2+solution2(i).vy.^2);
	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;