function md=icetransient3d(md);
%ICETRANSIENT3D - transient solution of 3d models
%
%   Usage:
%      md=icetransient3d(md)

global gridset

%Create fem structure
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, base vertical and vertical model
if fem.ismacayealpattyn,
	fem.m_dh=CreateFemModel(md,'diagnostic_horiz');
end
if fem.ishutter,
	fem.m_ss=CreateFemModel(md,'surface_slope_compute');
	fem.m_dhu=CreateFemModel(md,'diagnostic_hutter');
end
if fem.isstokes,
	fem.m_bs=CreateFemModel(md,'bed_slope_compute');
	fem.m_ds=CreateFemModel(md,'diagnostic_stokes');
end
fem.m_dbv=CreateFemModel(md,'diagnostic_basevert');
fem.m_dv=CreateFemModel(md,'diagnostic_vert');

%Second, build elements,grids,loads, etc ... for thermal and prognostic model
fem.m_t=CreateFemModel(md,'thermaltransient');
fem.m_m=CreateFemModel(md,'meltingtransient');
fem.m_p=CreateFemModel(md,'prognostic');

%initialize (velocity,pressure,...)
solution.u_g=zeros(fem.m_dhu.gridset.gsize,1);
if ~isnan(md.vx), solution.u_g(1:6:fem.m_dhu.gridset.gsize)=md.vx/md.yts; end
if ~isnan(md.vy), solution.u_g(2:6:fem.m_dhu.gridset.gsize)=md.vy/md.yts; end
if ~isnan(md.vz), solution.u_g(3:6:fem.m_dhu.gridset.gsize)=md.vz/md.yts; end
solution.pressure=zeros(fem.m_dhu.gridset.gsize,1); solution.pressure(1:6:fem.m_dhu.gridset.gsize)=md.rho_ice*md.g*(md.surface-md.z); %lithostatic pressure to spin up
solution.t_g=zeros(fem.m_t.gridset.gsize,1);
if ~isempty(md.temperature), 
	solution.t_g(1:6:fem.m_t.gridset.gsize)=md.temperature;
else
	solution.t_g(1:6:fem.m_t.gridset.gsize)=md.meltingpoint;
end
pos=find(isnan(solution.t_g));
solution.t_g(pos)=md.meltingpoint;
if ~isempty(md.melting), solution.melting_g(1:6:fem.m_m.gridset.gsize)=md.melting;end;

%initialize geometry and time
solution.mesh=struct('elements',md.elements,'grids',[md.x md.y md.z]);
solution.thickness=zeros(fem.m_dhu.gridset.gsize,1);solution.thickness(1:6:fem.m_dhu.gridset.gsize)=md.thickness;
solution.surface=zeros(fem.m_dhu.gridset.gsize,1);solution.surface(1:6:fem.m_dhu.gridset.gsize)=md.surface;
solution.bed=zeros(fem.m_dhu.gridset.gsize,1);solution.bed(1:6:fem.m_dhu.gridset.gsize)=md.bed;
solution.time=0;

%build constant inputs
fem.inputs=struct();

fem.inputs.melting=zeros(fem.m_m.gridset.gsize,1); melting=ShiftLayers(md,md.melting,1,1); fem.inputs.melting(1:6:fem.m_m.gridset.gsize)=melting;
fem.inputs.accumulation=zeros(fem.m_m.gridset.gsize,1); accumulation=ShiftLayers(md,md.accumulation,md.numlayers,1); fem.inputs.accumulation(1:6:fem.m_m.gridset.gsize)=accumulation;
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.pressure=solution(n).pressure;
	fem.inputs.temperature=full(solution(n).t_g);
	fem.inputs.thickness=solution(n).thickness;
	fem.inputs.bed=solution(n).bed;
	fem.inputs.surface=solution(n).surface;

	%%Deal with temperature first 
	%disp('   computing temperature...');
	%[solution(n+1).t_g loads_t melting_offset]=icethermal_core(fem.m_t,'thermaltransient',fem.inputs);
	solution(n+1).t_g=solution(1).t_g;
	%disp('   computing melting...');
	%fem.inputs.temperature=full(solution(n+1).t_g);
	%fem.inputs.melting_offset=melting_offset;
	%solution(n+1).melting_g=icemelting_core(fem.m_m,'thermaltransient',fem.inputs);
	solution(n+1).melting_g=solution(1).melting_g;

	%Compute depth averaged temperature for 2d type elements.
	%temperature_average2d=DepthAverage(md,solution(n+1).t_g(1:6:fem.m_t.gridset.gsize));
	%temperature_average=project3d(md,temperature_average2d,'node');
	%temperature_average_g=zeros(fem.m_t.gridset.gsize,1);temperature_average_g(1:6:fem.m_t.gridset.gsize)=temperature_average;
	%fem.inputs.temperature_average=temperature_average_g;
	
	%Deal with velocities.
	disp(sprintf('%s','   computing velocities...'));
	solution(n+1).u_g=icediagnostic3d(md,fem);
	solution(n+1).pressure=zeros(fem.m_dv.gridset.gsize,1); solution(n+1).pressure(1:6:end)=solution(n+1).u_g(4:6:end);

	%compute new thickness
	disp(sprintf('\n%s','   computing new thickness...'));
	%compute depth averaged horizontal velocity
	velocity_average=HorizontalVelocityDepthAverage(md,solution(n+1).u_g);
	fem.inputs.velocity_average=velocity_average;
	new_thickness=iceprognostic_core(fem.m_p,'prognostic',fem.inputs);

	%project collapsed thickness onto 3d mesh
	indx=1:6:fem.m_p.gridset.gsize; indx=indx(fem.m_p.tpart);
	new_thickness=project3d(md,project2d(md,new_thickness(indx),1),'node');

	%update surface and bed using the new thickness
	disp(sprintf('%s','   updating geometry...'));
	[new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,solution(n).bed(indx),solution(n).surface(indx));

	%project onto 3d mesh
	[solution(n+1).mesh solution(n+1).u_g solution(n+1).t_g solution(n+1).pressure]=UpdateMesh(md,fem.m_dv.grids,solution(n).mesh,solution(n+1).u_g,solution(n+1).t_g,solution(n+1).pressure,new_bed,new_thickness);
	solution(n+1).thickness=zeros(fem.m_p.gridset.gsize,1);solution(n+1).thickness(1:6:fem.m_p.gridset.gsize)=new_thickness;
	solution(n+1).bed=zeros(fem.m_p.gridset.gsize,1);solution(n+1).bed(1:6:fem.m_p.gridset.gsize)=new_bed;
	solution(n+1).surface=zeros(fem.m_p.gridset.gsize,1);solution(n+1).surface(1:6:fem.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 grids
	disp(sprintf('%s','   updating grid positions...'));
	grids_t=UpdateGridPosition(md,fem.m_t.grids,new_bed,new_thickness);
	grids_dhu=UpdateGridPosition(md,fem.m_dhu.grids,new_bed,new_thickness);
	grids_dbv=UpdateGridPosition(md,fem.m_dbv.grids,new_bed,new_thickness);
	grids_dv=UpdateGridPosition(md,fem.m_dv.grids,new_bed,new_thickness);
	grids_p=UpdateGridPosition(md,fem.m_p.grids,new_bed,new_thickness);

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

end

%load results onto model
indx=1:6:fem.m_dv.gridset.gsize; indx=indx(fem.m_dv.tpart);
indy=2:6:fem.m_dv.gridset.gsize; indy=indy(fem.m_dv.tpart);
indz=3:6:fem.m_dv.gridset.gsize; indz=indz(fem.m_dv.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).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;
