diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 59a883f0e808653433216bb16f3378a14cb705ba..9060c26a4125adb98cf90f024c4274d2477e41db 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -20,6 +20,7 @@
 #include "NumLib/ODESolver/PETScNonlinearSolver.h"
 #include "NumLib/ODESolver/TimeDiscretizedODESystem.h"
 #include "NumLib/StaggeredCoupling/StaggeredCoupling.h"
+#include "NumLib/TimeStepping/Time.h"
 #include "ProcessData.h"
 
 namespace
@@ -50,7 +51,8 @@ bool isOutputStep(std::vector<ProcessLib::Output> const& outputs,
 }
 
 void preOutputForAllProcesses(
-    int const timestep, double const t, double const dt, const double end_time,
+    int const timestep, NumLib::Time const& t, double const dt,
+    const NumLib::Time& end_time,
     std::vector<std::unique_ptr<ProcessLib::ProcessData>> const&
         per_process_data,
     std::vector<GlobalVector*> const& process_solutions,
@@ -67,7 +69,7 @@ void preOutputForAllProcesses(
         auto const process_id = process_data->process_id;
         auto& pcs = process_data->process;
 
-        pcs.preOutput(t, dt, process_solutions, process_solutions_prev,
+        pcs.preOutput(t(), dt, process_solutions, process_solutions_prev,
                       process_id);
     }
 }
@@ -76,7 +78,7 @@ void preOutputForAllProcesses(
 namespace ProcessLib
 {
 void preTimestepForAllProcesses(
-    double const t, double const dt,
+    NumLib::Time const& t, double const dt,
     std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
     std::vector<GlobalVector*> const& _process_solutions)
 {
@@ -84,12 +86,12 @@ void preTimestepForAllProcesses(
     {
         auto const process_id = process_data->process_id;
         auto& pcs = process_data->process;
-        pcs.preTimestep(_process_solutions, t, dt, process_id);
+        pcs.preTimestep(_process_solutions, t(), dt, process_id);
     }
 }
 
 void postTimestepForAllProcesses(
-    double const t, double const dt,
+    NumLib::Time const& t, double const dt,
     std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
     std::vector<GlobalVector*> const& process_solutions,
     std::vector<GlobalVector*> const& process_solutions_prev)
@@ -99,10 +101,10 @@ void postTimestepForAllProcesses(
         auto const process_id = process_data->process_id;
         auto& pcs = process_data->process;
 
-        pcs.computeSecondaryVariable(t, dt, process_solutions,
+        pcs.computeSecondaryVariable(t(), dt, process_solutions,
                                      *process_solutions_prev[process_id],
                                      process_id);
-        pcs.postTimestep(process_solutions, process_solutions_prev, t, dt,
+        pcs.postTimestep(process_solutions, process_solutions_prev, t(), dt,
                          process_id);
     }
 }
@@ -168,7 +170,7 @@ void setTimeDiscretizedODESystem(ProcessData& process_data)
 
 std::pair<std::vector<GlobalVector*>, std::vector<GlobalVector*>>
 setInitialConditions(
-    double const t0,
+    NumLib::Time const& t0,
     std::vector<std::unique_ptr<ProcessData>> const& per_process_data)
 {
     std::vector<GlobalVector*> process_solutions;
@@ -192,11 +194,11 @@ setInitialConditions(
     {
         auto& pcs = process_data->process;
         auto const process_id = process_data->process_id;
-        pcs.setInitialConditions(process_solutions, process_solutions_prev, t0,
-                                 process_id);
+        pcs.setInitialConditions(process_solutions, process_solutions_prev,
+                                 t0(), process_id);
 
         auto& time_disc = *process_data->time_disc;
-        time_disc.setInitialState(t0);  // push IC
+        time_disc.setInitialState(t0());  // push IC
     }
 
     return {process_solutions, process_solutions_prev};
@@ -245,8 +247,8 @@ NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
         // lead to some inconsistencies in the data compared to regular output.
         for (auto const& output : outputs)
         {
-            output.doOutputNonlinearIteration(process, process_id, timestep, t,
-                                              iteration, x);
+            output.doOutputNonlinearIteration(process, process_id, timestep,
+                                              NumLib::Time(t), iteration, x);
         }
     };
 
@@ -267,7 +269,7 @@ TimeLoop::TimeLoop(
     std::vector<Output>&& outputs,
     std::vector<std::unique_ptr<ProcessData>>&& per_process_data,
     std::unique_ptr<NumLib::StaggeredCoupling>&& staggered_coupling,
-    const double start_time, const double end_time)
+    const NumLib::Time& start_time, const NumLib::Time& end_time)
     : _outputs{std::move(outputs)},
       _per_process_data(std::move(per_process_data)),
       _start_time(start_time),
@@ -277,7 +279,8 @@ TimeLoop::TimeLoop(
 }
 
 bool computationOfChangeNeeded(
-    NumLib::TimeStepAlgorithm const& timestep_algorithm, double const time)
+    NumLib::TimeStepAlgorithm const& timestep_algorithm,
+    NumLib::Time const& time)
 {
     // for the first time step we can't compute the changes to the previous
     // time step
@@ -289,9 +292,9 @@ bool computationOfChangeNeeded(
 }
 
 std::pair<double, bool> TimeLoop::computeTimeStepping(
-    const double prev_dt, double& t, std::size_t& accepted_steps,
+    const double prev_dt, NumLib::Time& t, std::size_t& accepted_steps,
     std::size_t& rejected_steps,
-    std::vector<std::function<double(double, double)>> const&
+    std::vector<std::function<double(NumLib::Time const&, double)>> const&
         time_step_constraints)
 {
     bool all_process_steps_accepted = true;
@@ -327,11 +330,7 @@ std::pair<double, bool> TimeLoop::computeTimeStepping(
             solution_error, ppd.nonlinear_solver_status.number_iterations,
             ppd.timestep_previous, ppd.timestep_current);
 
-        if (!previous_step_accepted &&
-            // In case of FixedTimeStepping, which makes
-            // timestep_algorithm.next(...) return false when the ending time
-            // is reached.
-            t + eps < timestep_algorithm.end())
+        if (!previous_step_accepted)
         {
             // Not all processes have accepted steps.
             all_process_steps_accepted = false;
@@ -345,8 +344,7 @@ std::pair<double, bool> TimeLoop::computeTimeStepping(
             all_process_steps_accepted = false;
         }
 
-        if (timestepper_dt > eps ||
-            std::abs(t - timestep_algorithm.end()) < eps)
+        if (timestepper_dt > eps || t < timestep_algorithm.end())
         {
             dt = std::min(timestepper_dt, dt);
         }
@@ -371,7 +369,7 @@ std::pair<double, bool> TimeLoop::computeTimeStepping(
         }
         else
         {
-            if (t < _end_time || std::abs(t - _end_time) < eps)
+            if (t <= _end_time)
             {
                 t -= prev_dt;
                 rejected_steps++;
@@ -429,7 +427,7 @@ std::pair<double, bool> TimeLoop::computeTimeStepping(
         }
         else
         {
-            if (t < _end_time || std::abs(t - _end_time) < eps)
+            if (t <= _end_time)
             {
                 WARN(
                     "Time step {:d} was rejected {:d} times and it will be "
@@ -443,21 +441,21 @@ std::pair<double, bool> TimeLoop::computeTimeStepping(
     return {dt, last_step_rejected};
 }
 
-std::vector<std::function<double(double, double)>>
+std::vector<std::function<double(NumLib::Time const&, double)>>
 TimeLoop::generateOutputTimeStepConstraints(
     std::vector<double>&& fixed_times) const
 {
-    std::vector<std::function<double(double, double)>> const
+    std::vector<std::function<double(NumLib::Time const&, double)>> const
         time_step_constraints{
-            [fixed_times = std::move(fixed_times)](double t, double dt) {
+            [fixed_times = std::move(fixed_times)](NumLib::Time t, double dt) {
                 return NumLib::possiblyClampDtToNextFixedTime(t, dt,
                                                               fixed_times);
             },
-            [this](double t, double dt)
+            [this](NumLib::Time t, double dt) -> double
             {
-                if (t < _end_time && t + dt > _end_time)
+                if (t < _end_time && _end_time < t + dt)
                 {
-                    return _end_time - t;
+                    return _end_time() - t();
                 }
                 return dt;
             }};
@@ -495,7 +493,7 @@ void TimeLoop::initialize()
         _staggered_coupling->initializeCoupledSolutions(_process_solutions);
     }
 
-    updateDeactivatedSubdomains(_per_process_data, _start_time);
+    updateDeactivatedSubdomains(_per_process_data, _start_time());
 
     auto const time_step_constraints = generateOutputTimeStepConstraints(
         calculateUniqueFixedTimesForAllOutputs(_outputs));
@@ -507,7 +505,7 @@ void TimeLoop::initialize()
     // Output initial conditions
     {
         preOutputInitialConditions(_start_time, _dt);
-        outputSolutions(0, _start_time, &Output::doOutput);
+        outputSolutions(0, _start_time(), &Output::doOutput);
     }
 
     calculateNonEquilibriumInitialResiduum(
@@ -526,9 +524,9 @@ bool TimeLoop::executeTimeStep()
     INFO(
         "=== Time stepping at step #{:d} and time {:.18g} with step size "
         "{:.18g}",
-        timesteps, _current_time, _dt);
+        timesteps, _current_time(), _dt);
 
-    updateDeactivatedSubdomains(_per_process_data, _current_time);
+    updateDeactivatedSubdomains(_per_process_data, _current_time());
 
     successful_time_step =
         preTsNonlinearSolvePostTs(_current_time, _dt, timesteps);
@@ -540,7 +538,7 @@ bool TimeLoop::executeTimeStep()
 bool TimeLoop::calculateNextTimeStep()
 {
     const double prev_dt = _dt;
-    double const current_time = _current_time;
+    auto const current_time = _current_time;
 
     const std::size_t timesteps = _accepted_steps + 1;
 
@@ -554,22 +552,19 @@ bool TimeLoop::calculateNextTimeStep()
 
     if (!_last_step_rejected)
     {
-        outputSolutions(timesteps, current_time, &Output::doOutput);
+        outputSolutions(timesteps, current_time(), &Output::doOutput);
     }
 
-    if (std::abs(_current_time - _end_time) <
-            std::numeric_limits<double>::epsilon() ||
-        _current_time + _dt > _end_time)
+    if (current_time == (_current_time + _dt))
     {
+        ERR("Time step size of {:.18g} is too small.\n"
+            "Time stepping stops at step {:d} and at time of {:.18g}.",
+            _dt, timesteps, _current_time());
         return false;
     }
 
-    if (_dt < std::numeric_limits<double>::epsilon())
+    if (_current_time >= _end_time || _current_time + _dt > _end_time)
     {
-        WARN(
-            "Time step size of {:.18g} is too small.\n"
-            "Time stepping stops at step {:d} and at time of {:.18g}.",
-            _dt, timesteps, _current_time);
         return false;
     }
 
@@ -586,12 +581,12 @@ void TimeLoop::outputLastTimeStep() const
     // output last time step
     if (successful_time_step)
     {
-        outputSolutions(_accepted_steps + _rejected_steps, _current_time,
+        outputSolutions(_accepted_steps + _rejected_steps, _current_time(),
                         &Output::doOutputLastTimestep);
     }
 }
 
-bool TimeLoop::preTsNonlinearSolvePostTs(double const t, double const dt,
+bool TimeLoop::preTsNonlinearSolvePostTs(NumLib::Time const& t, double const dt,
                                          std::size_t const timesteps)
 {
     preTimestepForAllProcesses(t, dt, _per_process_data, _process_solutions);
@@ -629,7 +624,7 @@ bool TimeLoop::preTsNonlinearSolvePostTs(double const t, double const dt,
 }
 
 static NumLib::NonlinearSolverStatus solveMonolithicProcess(
-    const double t, const double dt, const std::size_t timestep_id,
+    const NumLib::Time& t, const double dt, const std::size_t timestep_id,
     ProcessData const& process_data, std::vector<GlobalVector*>& x,
     std::vector<GlobalVector*> const& x_prev,
     std::vector<Output> const& outputs)
@@ -638,7 +633,7 @@ static NumLib::NonlinearSolverStatus solveMonolithicProcess(
     time_timestep_process.start();
 
     auto const nonlinear_solver_status = solveOneTimeStepOneProcess(
-        x, x_prev, timestep_id, t, dt, process_data, outputs);
+        x, x_prev, timestep_id, t(), dt, process_data, outputs);
 
     INFO("[time] Solving process #{:d} took {:g} s in time step #{:d}",
          process_data.process_id, time_timestep_process.elapsed(), timestep_id);
@@ -650,7 +645,7 @@ static constexpr std::string_view timestepper_cannot_reduce_dt =
     "Time stepper cannot reduce the time step size further.";
 
 NumLib::NonlinearSolverStatus TimeLoop::solveUncoupledEquationSystems(
-    const double t, const double dt, const std::size_t timestep_id)
+    const NumLib::Time& t, const double dt, const std::size_t timestep_id)
 {
     NumLib::NonlinearSolverStatus nonlinear_solver_status;
 
@@ -666,7 +661,7 @@ NumLib::NonlinearSolverStatus TimeLoop::solveUncoupledEquationSystems(
         {
             ERR("The nonlinear solver failed in time step #{:d} at t = {:.18g} "
                 "s for process #{:d}.",
-                timestep_id, t, process_id);
+                timestep_id, t(), process_id);
 
             if (!process_data->timestep_algorithm->canReduceTimestepSize(
                     process_data->timestep_current,
@@ -692,11 +687,11 @@ NumLib::NonlinearSolverStatus TimeLoop::solveUncoupledEquationSystems(
 
 NumLib::NonlinearSolverStatus
 TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
-    const double t, const double dt, const std::size_t timestep_id)
+    const NumLib::Time& t, const double dt, const std::size_t timestep_id)
 {
     auto const nonlinear_solver_status =
         _staggered_coupling->execute<ProcessData, Output>(
-            t, dt, timestep_id, _process_solutions, _process_solutions_prev,
+            t(), dt, timestep_id, _process_solutions, _process_solutions_prev,
             _per_process_data, _outputs, &solveOneTimeStepOneProcess);
 
     _last_step_rejected = nonlinear_solver_status.error_norms_met;
@@ -708,7 +703,7 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
             int const process_id = process_data->process_id;
             auto& ode_sys = *process_data->tdisc_ode_sys;
             pcs.solveReactionEquation(_process_solutions,
-                                      _process_solutions_prev, t, dt, ode_sys,
+                                      _process_solutions_prev, t(), dt, ode_sys,
                                       process_id);
         }
     }
@@ -735,7 +730,7 @@ void TimeLoop::outputSolutions(unsigned timestep, const double t,
         for (auto const& output_object : _outputs)
         {
             (output_object.*output_class_member)(
-                pcs, process_id, timestep, t,
+                pcs, process_id, timestep, NumLib::Time(t),
                 process_data->nonlinear_solver_status.number_iterations,
                 _process_solutions);
         }
@@ -754,7 +749,8 @@ TimeLoop::~TimeLoop()
     }
 }
 
-void TimeLoop::preOutputInitialConditions(const double t, const double dt) const
+void TimeLoop::preOutputInitialConditions(NumLib::Time const& t,
+                                          const double dt) const
 {
     for (auto const& process_data : _per_process_data)
     {
@@ -768,16 +764,16 @@ void TimeLoop::preOutputInitialConditions(const double t, const double dt) const
         auto const process_id = process_data->process_id;
         auto& pcs = process_data->process;
 
-        process_data->time_disc->nextTimestep(t, dt);
+        process_data->time_disc->nextTimestep(t(), dt);
 
-        pcs.preTimestep(_process_solutions, _start_time, dt, process_id);
+        pcs.preTimestep(_process_solutions, _start_time(), dt, process_id);
 
-        pcs.preOutput(_start_time, dt, _process_solutions,
+        pcs.preOutput(_start_time(), dt, _process_solutions,
                       _process_solutions_prev, process_id);
 
         // Update secondary variables, which might be uninitialized, before
         // output.
-        pcs.computeSecondaryVariable(_start_time, dt, _process_solutions,
+        pcs.computeSecondaryVariable(_start_time(), dt, _process_solutions,
                                      *_process_solutions_prev[process_id],
                                      process_id);
     }
diff --git a/ProcessLib/TimeLoop.h b/ProcessLib/TimeLoop.h
index e96047a97925c01f2417d4f660ebcd5d8f385ca4..6445d2ce01c71d040f74c4cfae8f088f7471e794 100644
--- a/ProcessLib/TimeLoop.h
+++ b/ProcessLib/TimeLoop.h
@@ -22,6 +22,7 @@ namespace NumLib
 {
 class ConvergenceCriterion;
 class StaggeredCoupling;
+struct Time;
 }
 
 namespace ChemistryLib
@@ -40,7 +41,7 @@ public:
     TimeLoop(std::vector<Output>&& outputs,
              std::vector<std::unique_ptr<ProcessData>>&& per_process_data,
              std::unique_ptr<NumLib::StaggeredCoupling>&& staggered_coupling,
-             const double start_time, const double end_time);
+             const NumLib::Time& start_time, const NumLib::Time& end_time);
 
     void initialize();
     void outputLastTimeStep() const;
@@ -58,12 +59,12 @@ public:
     /// otherwise.
     bool calculateNextTimeStep();
 
-    double endTime() const { return _end_time; }
-    double currentTime() const { return _current_time; }
+    double endTime() const { return _end_time(); }
+    double currentTime() const { return _current_time(); }
     bool successful_time_step = false;
 
 private:
-    bool preTsNonlinearSolvePostTs(double const t, double const dt,
+    bool preTsNonlinearSolvePostTs(NumLib::Time const& t, double const dt,
                                    std::size_t const timesteps);
 
     /**
@@ -78,7 +79,7 @@ private:
      *                    false: if any of nonlinear solvers divergences.
      */
     NumLib::NonlinearSolverStatus solveUncoupledEquationSystems(
-        const double t, const double dt, const std::size_t timestep_id);
+        const NumLib::Time& t, const double dt, const std::size_t timestep_id);
 
     /**
      * \brief Member to solver coupled systems of equations by the staggered
@@ -91,7 +92,7 @@ private:
      *                    false:  if any of nonlinear solvers divergences.
      */
     NumLib::NonlinearSolverStatus solveCoupledEquationSystemsByStaggeredScheme(
-        const double t, const double dt, const std::size_t timestep_id);
+        const NumLib::Time& t, const double dt, const std::size_t timestep_id);
 
     /**
      *  Find the minimum time step size among the predicted step sizes of
@@ -109,9 +110,9 @@ private:
      *  rejected
      */
     std::pair<double, bool> computeTimeStepping(
-        const double prev_dt, double& t, std::size_t& accepted_steps,
+        const double prev_dt, NumLib::Time& t, std::size_t& accepted_steps,
         std::size_t& rejected_steps,
-        std::vector<std::function<double(double, double)>> const&
+        std::vector<std::function<double(NumLib::Time const&, double)>> const&
             time_step_constraints);
 
     template <typename OutputClassMember>
@@ -120,17 +121,18 @@ private:
                          OutputClassMember output_class_member) const;
 
 private:
-    std::vector<std::function<double(double, double)>>
+    std::vector<std::function<double(NumLib::Time const&, double)>>
     generateOutputTimeStepConstraints(std::vector<double>&& fixed_times) const;
-    void preOutputInitialConditions(const double t, const double dt) const;
+    void preOutputInitialConditions(NumLib::Time const& t,
+                                    const double dt) const;
     std::vector<GlobalVector*> _process_solutions;
     std::vector<GlobalVector*> _process_solutions_prev;
     std::vector<Output> _outputs;
     std::vector<std::unique_ptr<ProcessData>> _per_process_data;
 
-    const double _start_time;
-    const double _end_time;
-    double _current_time = _start_time;
+    const NumLib::Time _start_time;
+    const NumLib::Time _end_time;
+    NumLib::Time _current_time = _start_time;
     std::size_t _accepted_steps = 0;
     std::size_t _rejected_steps = 0;
     double _dt = 0;