function femmodel=transient3d_core(femmodel)
%TRANSIENT3D_CORE - core of transient 2d solution
%
%   Usage:
%      femmodel=transient3d_core(femmodel)

	%recover parameters common to all solutions
	dim=femmodel.parameters.MeshDimension;
	ndt=femmodel.parameters.TimesteppingFinalTime;
	dt=femmodel.parameters.TimesteppingTimeStep;
	yts=femmodel.parameters.ConstantsYts;
	control_analysis=femmodel.parameters.InversionIscontrol;
	solution_type=femmodel.parameters.SolutionType;
	output_frequency=femmodel.parameters.SettingsOutputFrequency;
	time_adapt=femmodel.parameters.TimesteppingTimeAdapt;
	dakota_analysis=femmodel.parameters.QmuIsdakota;
	isdiagnostic=femmodel.parameters.TransientIsdiagnostic;
	isprognostic=femmodel.parameters.TransientIsprognostic;
	isthermal=femmodel.parameters.TransientIsthermal;
	isgroundingline=femmodel.parameters.TransientIsgroundingline;
	isenthalpy=femmodel.parameters.ThermalIsenthalpy;
	groundinglinemigration=femmodel.parameters.GroundinglineMigration;

	%Initialize
	time=0;
	step=0;

	%for qmu analysis, be sure the velocity input we are starting from  is the one in the parameters: 
	if dakota_analysis,
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVxEnum,VxEnum);
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVyEnum,VyEnum);
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVzEnum,VzEnum);
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuPressureEnum,PressureEnum);
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuBedEnum,BedEnum);
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuThicknessEnum,ThicknessEnum);      
		femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuSurfaceEnum,SurfaceEnum);
		if (isthermal & dim==3)
			femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuTemperatureEnum,TemperatureEnum);
			femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuMeltingEnum,BasalforcingsMeltingRateEnum);
		end
	end

	%Loop through time
	while (time < ndt - (yts*eps)),

		%Increment
		if(time_adapt),
			dt=TimeAdapt(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
			if(time+dt>ndt), dt=ndt-time; end
			femmodel.parameters.TimesteppingTimeStep=dt;
		end
		step=step+1;
		time=time+dt;
		femmodel.parameters.Time=time;

		issmprintf(VerboseSolution,'\n%s%g%s%i%s%g\n','time [yr] ',time/yts,' iteration number: ',step,'/',floor(ndt/dt));

		if (isthermal & dim==3)
			issmprintf(VerboseSolution,'\n%s',['   computing temperature']);
			if (isenthalpy==0),
				femmodel=thermal_core_step(femmodel); 
			else
				femmodel=enthalpy_core_step(femmodel); 
			end
		end

		if (isdiagnostic)
			issmprintf(VerboseSolution,'\n%s',['   computing new velocities']);
			femmodel=diagnostic_core(femmodel); 
		end

		if (isprognostic)
			issmprintf(VerboseSolution,'\n%s',['   computing new thickness']);
			femmodel=prognostic_core(femmodel); 
			issmprintf(VerboseSolution,'\n%s',['   updating vertices position']);
			[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=...
				UpdateVertexPositions(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
		end

		if (isgroundingline)
			issmprintf(VerboseSolution,'\n%s',['   computing new grounding line']);
			[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=...
				GroundinglineMigration(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
		end

		if (solution_type==TransientSolutionEnum & ~control_analysis & (mod(step,output_frequency)==0 | time==ndt)),
			issmprintf(VerboseSolution,'\n%s',['   saving results']);
			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,step,time);
			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,step,time);
			if (dim==3), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,step,time); end
			if (isdiagnostic), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VelEnum,step,time); end
			if (isdiagnostic), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,PressureEnum,step,time); end
			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ThicknessEnum,step,time);
			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceEnum,step,time);
			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BedEnum,step,time);
			if (dim==3 & isthermal), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum,step,time);end
			if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,WaterfractionEnum,step,time);end
			if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,EnthalpyEnum,step,time);end
			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BasalforcingsMeltingRateEnum,step,time);
			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceforcingsMassBalanceEnum,step,time);
			femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,MaskElementonfloatingiceEnum,step,time);
		end
	end
end %end of function
