diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh
index c7e8d11f35..00b03445a0 100644
--- a/ebos/eclproblem.hh
+++ b/ebos/eclproblem.hh
@@ -582,15 +582,17 @@ public:
         const auto& vanguard = simulator.vanguard();
         SolventModule::initFromDeck(vanguard.deck(), vanguard.eclState());
         PolymerModule::initFromDeck(vanguard.deck(), vanguard.eclState());
-        if (EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput))
-            // create the ECL writer
-            eclWriter_.reset(new EclWriterType(simulator));
+
+        // create the ECL writer
+        eclWriter_.reset(new EclWriterType(simulator));
 
         if (enableExperiments)
             enableDriftCompensation_ = EWOMS_GET_PARAM(TypeTag, bool, EclEnableDriftCompensation);
         else
             enableDriftCompensation_ = false;
 
+        enableEclOutput_ = EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput);
+
         enableTuning_ = EWOMS_GET_PARAM(TypeTag, bool, EclEnableTuning);
         initialTimeStepSize_ = EWOMS_GET_PARAM(TypeTag, Scalar, InitialTimeStepSize);
         minTimeStepSize_ = EWOMS_GET_PARAM(TypeTag, Scalar, MinTimeStepSize);
@@ -608,6 +610,19 @@ public:
         ParentType::finishInit();
 
         auto& simulator = this->simulator();
+        const auto& eclState = simulator.vanguard().eclState();
+        const auto& schedule = simulator.vanguard().schedule();
+        const auto& timeMap = schedule.getTimeMap();
+
+        // Set the start time of the simulation
+        simulator.setStartTime(timeMap.getStartTime(/*reportStepIdx=*/0));
+
+        // We want the episode index to be the same as the report step index to make
+        // things simpler, so we have to set the episode index to -1 because it is
+        // incremented by endEpisode(). The size of the initial time step and
+        // length of the initial episode is set to zero for the same reason.
+        simulator.setEpisodeIndex(-1);
+        simulator.setEpisodeLength(0.0);
 
         // the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false
         // disables gravity, else the standard value of the gravity constant at sea level
@@ -632,7 +647,6 @@ public:
         }
 
         // deal with DRSDT
-        const auto& eclState = simulator.vanguard().eclState();
         unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables();
         maxDRs_.resize(ntpvt, 1e30);
         dRsDtOnlyFreeGas_.resize(ntpvt, false);
@@ -649,34 +663,17 @@ public:
         readThermalParameters_();
         transmissibilities_.finishInit();
 
-        const auto& timeMap = simulator.vanguard().schedule().getTimeMap();
-
         readInitialCondition_();
-        // Set the start time of the simulation
-        simulator.setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0));
-
-        // We want the episode index to be the same as the report step index to make
-        // things simpler, so we have to set the episode index to -1 because it is
-        // incremented inside beginEpisode(). The size of the initial time step and
-        // length of the initial episode is set to zero for the same reason.
-        simulator.setEpisodeIndex(-1);
-        simulator.setEpisodeLength(0.0);
-        simulator.setTimeStepSize(0.0);
 
         updatePffDofData_();
 
         if (GET_PROP_VALUE(TypeTag, EnablePolymer)) {
-            const auto& vanguard = this->simulator().vanguard();
+            const auto& vanguard = simulator.vanguard();
             const auto& gridView = vanguard.gridView();
             int numElements = gridView.size(/*codim=*/0);
             maxPolymerAdsorption_.resize(numElements, 0.0);
         }
 
-        if (eclWriter_) {
-            eclWriter_->writeInit();
-            this->simulator().vanguard().releaseGlobalTransmissibilities();
-        }
-
         tracerModel_.init();
 
         readBoundaryConditions_();
@@ -689,6 +686,16 @@ public:
         if (enableExperiments)
             checkDeckCompatibility_();
 
+        // write the static output files (EGRID, INIT, SMSPEC, etc.)
+        if (enableEclOutput_)
+            eclWriter_->writeInit();
+
+        simulator.vanguard().releaseGlobalTransmissibilities();
+
+        // after finishing the initialization and writing the initial solution, we move
+        // to the first "real" episode/report step
+        simulator.startNextEpisode(timeMap.getTimeStepLength(0));
+        simulator.setEpisodeIndex(0);
     }
 
     void prefetch(const Element& elem) const
@@ -737,37 +744,21 @@ public:
     void beginEpisode(bool isOnRestart = false)
     {
         // Proceed to the next report step
-        Simulator& simulator = this->simulator();
-        auto& eclState = this->simulator().vanguard().eclState();
-        const auto& schedule = this->simulator().vanguard().schedule();
+        auto& simulator = this->simulator();
+        auto& eclState = simulator.vanguard().eclState();
+        const auto& schedule = simulator.vanguard().schedule();
         const auto& events = schedule.getEvents();
         const auto& timeMap = schedule.getTimeMap();
+        int episodeIdx = simulator.episodeIndex();
 
-        // The first thing to do in the morning of an episode is update update the
-        // eclState and the deck if they need to be changed.
-        int nextEpisodeIdx = simulator.episodeIndex();
-
-        if (enableExperiments && this->gridView().comm().rank() == 0) {
-            boost::posix_time::ptime curDateTime =
-                boost::posix_time::from_time_t(timeMap.getStartTime(nextEpisodeIdx+1));
-            std::cout << "Report step " << nextEpisodeIdx + 2
-                      << "/" << timeMap.numTimesteps()
-                      << " at day " << timeMap.getTimePassedUntil(nextEpisodeIdx+1)/(24*3600)
-                      << "/" << timeMap.getTotalTime()/(24*3600)
-                      << ", date = " << curDateTime.date()
-                      << "\n ";
-        }
-
-        if (nextEpisodeIdx > 0 &&
-            events.hasEvent(Opm::ScheduleEvents::GEO_MODIFIER, nextEpisodeIdx))
-        {
+        if (episodeIdx >= 0 && events.hasEvent(Opm::ScheduleEvents::GEO_MODIFIER, episodeIdx)) {
             // bring the contents of the keywords to the current state of the SCHEDULE
-            // section
+            // section.
             //
             // TODO (?): make grid topology changes possible (depending on what exactly
             // has changed, the grid may need be re-created which has some serious
             // implications on e.g., the solution of the simulation.)
-            const auto& miniDeck = schedule.getModifierDeck(nextEpisodeIdx);
+            const auto& miniDeck = schedule.getModifierDeck(episodeIdx);
             eclState.applyModifierDeck(miniDeck);
 
             // re-compute all quantities which may possibly be affected.
@@ -776,44 +767,32 @@ public:
             updatePffDofData_();
         }
 
-        // Opm::TimeMap deals with points in time, so the number of time intervals (i.e.,
-        // report steps) is one less!
-        int numReportSteps = timeMap.size() - 1;
-
-        // start the next episode if there are additional report steps, else finish the
-        // simulation
-        while (nextEpisodeIdx < numReportSteps &&
-               simulator.time() >= timeMap.getTimePassedUntil(nextEpisodeIdx + 1)*(1 - 1e-10))
-        {
-            ++ nextEpisodeIdx;
+        if (enableExperiments && this->gridView().comm().rank() == 0 && episodeIdx >= 0) {
+            // print some useful information in experimental mode. (the production
+            // simulator does this externally.)
+            boost::posix_time::ptime curDateTime =
+                boost::posix_time::from_time_t(timeMap.getStartTime(episodeIdx));
+            std::cout << "Report step " << episodeIdx + 1
+                      << "/" << timeMap.numTimesteps()
+                      << " at day " << timeMap.getTimePassedUntil(episodeIdx)/(24*3600)
+                      << "/" << timeMap.getTotalTime()/(24*3600)
+                      << ", date = " << curDateTime.date()
+                      << "\n ";
         }
 
         // react to TUNING changes
         bool tuningEvent = false;
-        if (nextEpisodeIdx > 0
-            && enableTuning_
-            && events.hasEvent(Opm::ScheduleEvents::TUNING_CHANGE, nextEpisodeIdx))
+        if (episodeIdx > 0 && enableTuning_ && events.hasEvent(Opm::ScheduleEvents::TUNING_CHANGE, episodeIdx))
         {
             const auto& tuning = schedule.getTuning();
-            initialTimeStepSize_ = tuning.getTSINIT(nextEpisodeIdx);
-            maxTimeStepAfterWellEvent_ = tuning.getTMAXWC(nextEpisodeIdx);
-            maxTimeStepSize_ = tuning.getTSMAXZ(nextEpisodeIdx);
-            restartShrinkFactor_ = 1./tuning.getTSFCNV(nextEpisodeIdx);
-            minTimeStepSize_ = tuning.getTSMINZ(nextEpisodeIdx);
+            initialTimeStepSize_ = tuning.getTSINIT(episodeIdx);
+            maxTimeStepAfterWellEvent_ = tuning.getTMAXWC(episodeIdx);
+            maxTimeStepSize_ = tuning.getTSMAXZ(episodeIdx);
+            restartShrinkFactor_ = 1./tuning.getTSFCNV(episodeIdx);
+            minTimeStepSize_ = tuning.getTSMINZ(episodeIdx);
             tuningEvent = true;
         }
 
-        Scalar episodeLength = timeMap.getTimeStepLength(nextEpisodeIdx);
-        simulator.startNextEpisode(episodeLength);
-
-        Scalar dt = limitNextTimeStepSize_(episodeLength);
-        if (nextEpisodeIdx == 0 || tuningEvent)
-            // allow the size of the initial time step to be set via an external parameter
-            // if TUNING is enabled, also limit the time step size after a tuning event to TSINIT
-            dt = std::min(dt, initialTimeStepSize_);
-
-        simulator.setTimeStepSize(dt);
-
         const bool invalidateFromHyst = updateHysteresis_();
         const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
         const bool doInvalidate = invalidateFromHyst || invalidateFromMaxOilSat;
@@ -822,14 +801,19 @@ public:
             updateMaxPolymerAdsorption_();
 
         // set up the wells for the next episode.
-        //
-        // TODO: the first two arguments seem to be unnecessary
-        wellModel_.beginEpisode(this->simulator().vanguard().eclState(),
-                                this->simulator().vanguard().schedule(),
-                                isOnRestart);
+        wellModel_.beginEpisode();
 
+        // set up the aquifers for the next episode.
         aquiferModel_.beginEpisode();
 
+        // set the size of the initial time step of the episode
+        Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
+        if (episodeIdx == 0 || tuningEvent)
+            // allow the size of the initial time step to be set via an external parameter
+            // if TUNING is enabled, also limit the time step size after a tuning event to TSINIT
+            dt = std::min(dt, initialTimeStepSize_);
+        simulator.setTimeStepSize(dt);
+
         if (doInvalidate)
             this->model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
     }
@@ -839,17 +823,18 @@ public:
      */
     void beginTimeStep()
     {
-        int epsiodeIdx = this->simulator().episodeIndex();
-        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        const auto& simulator = this->simulator();
+        int epsiodeIdx = simulator.episodeIndex();
+        const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
         if (drsdtActive_())
             // DRSDT is enabled
             for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx)
-                maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*this->simulator().timeStepSize();
+                maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*simulator.timeStepSize();
 
         if (drvdtActive_())
             // DRVDT is enabled
             for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx)
-                maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*this->simulator().timeStepSize();
+                maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*simulator.timeStepSize();
 
         wellModel_.beginTimeStep();
         aquiferModel_.beginTimeStep();
@@ -896,10 +881,13 @@ public:
             // the right thing (i.e., the mass change inside the whole reservoir must be
             // equivalent to the fluxes over the grid's boundaries plus the source rates
             // specified by the problem)
+            std::cout << "checking conservativeness of solution\n";
             this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
+            std::cout << "solution is sufficiently conservative\n";
         }
 #endif // NDEBUG
 
+        const auto& simulator = this->simulator();
         wellModel_.endTimeStep();
         aquiferModel_.endTimeStep();
         tracerModel_.endTimeStep();
@@ -911,7 +899,7 @@ public:
             const auto& residual = this->model().linearizer().residual();
             for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
                 drift_[globalDofIdx] = residual[globalDofIdx];
-                drift_[globalDofIdx] *= this->simulator().timeStepSize();
+                drift_[globalDofIdx] *= simulator.timeStepSize();
                 if (GET_PROP_VALUE(TypeTag, UseVolumetricResidual))
                     drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx);
             }
@@ -925,15 +913,18 @@ public:
     {
         auto& simulator = this->simulator();
         const auto& schedule = simulator.vanguard().schedule();
+        const auto& timeMap = schedule.getTimeMap();
 
         int episodeIdx = simulator.episodeIndex();
 
-        const auto& timeMap = schedule.getTimeMap();
-        int numReportSteps = timeMap.size() - 1;
-        if (episodeIdx + 1 >= numReportSteps) {
+        // check if we're finished ...
+        if (episodeIdx + 1 >= static_cast<int>(timeMap.numTimesteps())) {
             simulator.setFinished(true);
             return;
         }
+
+        // .. if we're not yet done, start the next episode (report step)
+        simulator.startNextEpisode(timeMap.getTimeStepLength(episodeIdx + 1));
     }
 
     /*!
@@ -945,27 +936,25 @@ public:
      */
     bool shouldWriteOutput() const
     {
-        if (this->simulator().timeStepIndex() < 0)
+        const auto& simulator = this->simulator();
+        if (simulator.timeStepIndex() < 0)
             // always write the initial solution
             return true;
 
         if (EWOMS_GET_PARAM(TypeTag, bool, EnableWriteAllSolutions))
             return true;
 
-        return this->simulator().episodeWillBeOver();
+        return simulator.episodeWillBeOver();
     }
 
     /*!
      * \brief Returns true if an eWoms restart file should be written to disk.
+     *
+     * The EclProblem does not write any restart files using the ad-hoc format, only ones
+     * using the ECL format.
      */
     bool shouldWriteRestartFile() const
-    {
-        unsigned n = EWOMS_GET_PARAM(TypeTag, unsigned, RestartWritingInterval);
-        unsigned i = this->simulator().timeStepIndex();
-        if (i > 0 && (i%n) == 0)
-            return true; // we don't write a restart file for the initial condition
-        return false;
-    }
+    { return false; }
 
     /*!
      * \brief Write the requested quantities of the current solution into the output
@@ -977,10 +966,8 @@ public:
         // write the desired VTK files.
         ParentType::writeOutput(isSubStep, verbose);
 
-        if (!eclWriter_)
-            return;
-
-        eclWriter_->writeOutput(isSubStep);
+        if (enableEclOutput_)
+            eclWriter_->writeOutput(isSubStep);
     }
 
     /*!
@@ -1443,21 +1430,20 @@ public:
      */
     void initialSolutionApplied()
     {
-        const auto& eclState = this->simulator().vanguard().eclState();
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
+        const auto& eclState = vanguard.eclState();
 
         // initialize the wells. Note that this needs to be done after initializing the
         // intrinsic permeabilities and the after applying the initial solution because
         // the well model uses these...
-        wellModel_.init(eclState, this->simulator().vanguard().schedule());
+        wellModel_.init();
 
         // let the object for threshold pressures initialize itself. this is done only at
         // this point, because determining the threshold pressures may require to access
         // the initial solution.
         thresholdPressures_.finishInit();
 
-        // release the memory of the EQUIL grid since it's no longer needed after this point
-        this->simulator().vanguard().releaseEquilGrid();
-
         updateCompositionChangeLimits_();
 
         aquiferModel_.initialSolutionApplied();
@@ -1613,8 +1599,9 @@ public:
 
     bool vapparsActive() const
     {
-        int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0);
-        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        const auto& simulator = this->simulator();
+        int epsiodeIdx = std::max(simulator.episodeIndex(), 0);
+        const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
         return (oilVaporizationControl.getType() == Opm::OilVaporizationEnum::VAPPARS);
     }
 
@@ -1633,7 +1620,8 @@ public:
         if (this->nextTimeStepSize_ > 0.0)
             return this->nextTimeStepSize_;
 
-        int episodeIdx = this->simulator().episodeIndex();
+        const auto& simulator = this->simulator();
+        int episodeIdx = simulator.episodeIndex();
 
         // for the initial episode, we use a fixed time step size
         if (episodeIdx < 0)
@@ -1643,7 +1631,6 @@ public:
         // well the previous time step converged. After that, apply the runtime time
         // stepping constraints.
         const auto& newtonMethod = this->model().newtonMethod();
-        const auto& simulator = this->simulator();
         return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
     }
 
@@ -1745,15 +1732,17 @@ private:
 
     bool drsdtActive_() const
     {
-        int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0);
-        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        const auto& simulator = this->simulator();
+        int epsiodeIdx = std::max(simulator.episodeIndex(), 0);
+        const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
         return (oilVaporizationControl.drsdtActive());
     }
 
     bool drvdtActive_() const
     {
-        int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0);
-        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        const auto& simulator = this->simulator();
+        int epsiodeIdx = std::max(simulator.episodeIndex(), 0);
+        const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
         return (oilVaporizationControl.drvdtActive());
 
     }
@@ -1774,7 +1763,8 @@ private:
 
     void updateElementDepths_()
     {
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
         const auto& gridView = vanguard.gridView();
         const auto& elemMapper = this->elementMapper();;
 
@@ -1796,12 +1786,13 @@ private:
     {
         // update the "last Rs" values for all elements, including the ones in the ghost
         // and overlap regions
-        int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0);
-        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        const auto& simulator = this->simulator();
+        int epsiodeIdx = std::max(simulator.episodeIndex(), 0);
+        const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
 
         if (oilVaporizationControl.drsdtActive()) {
-            ElementContext elemCtx(this->simulator());
-            const auto& vanguard = this->simulator().vanguard();
+            ElementContext elemCtx(simulator);
+            const auto& vanguard = simulator.vanguard();
             auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
             const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
             for (; elemIt != elemEndIt; ++elemIt) {
@@ -1830,8 +1821,8 @@ private:
         // update the "last Rv" values for all elements, including the ones in the ghost
         // and overlap regions
         if (drvdtActive_()) {
-            ElementContext elemCtx(this->simulator());
-            const auto& vanguard = this->simulator().vanguard();
+            ElementContext elemCtx(simulator);
+            const auto& vanguard = simulator.vanguard();
             auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
             const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
             for (; elemIt != elemEndIt; ++elemIt) {
@@ -1856,10 +1847,12 @@ private:
 
     bool updateMaxOilSaturation_()
     {
+        const auto& simulator = this->simulator();
+
         // we use VAPPARS
         if (vapparsActive()) {
-            ElementContext elemCtx(this->simulator());
-            const auto& vanguard = this->simulator().vanguard();
+            ElementContext elemCtx(simulator);
+            const auto& vanguard = simulator.vanguard();
             auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
             const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
             for (; elemIt != elemEndIt; ++elemIt) {
@@ -1887,9 +1880,10 @@ private:
 
     void readRockParameters_()
     {
-        const auto& deck = this->simulator().vanguard().deck();
-        const auto& eclState = this->simulator().vanguard().eclState();
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        const auto& deck = simulator.vanguard().deck();
+        const auto& eclState = simulator.vanguard().eclState();
+        const auto& vanguard = simulator.vanguard();
 
         // the ROCK keyword has not been specified, so we don't need
         // to read rock parameters
@@ -1945,7 +1939,8 @@ private:
 
     void readMaterialParameters_()
     {
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
         const auto& deck = vanguard.deck();
         const auto& eclState = vanguard.eclState();
 
@@ -1980,7 +1975,8 @@ private:
         if (!enableEnergy)
             return;
 
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
         const auto& deck = vanguard.deck();
         const auto& eclState = vanguard.eclState();
 
@@ -1996,7 +1992,8 @@ private:
 
     void updatePorosity_()
     {
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
         const auto& eclState = vanguard.eclState();
         const auto& eclGrid = eclState.getInputGrid();
         const auto& props = eclState.get3DProperties();
@@ -2044,7 +2041,7 @@ private:
             // we define the porosity as the accumulated pore volume divided by the
             // geometric volume of the element. Note that -- in pathetic cases -- it can
             // be larger than 1.0!
-            Scalar dofVolume = this->simulator().model().dofTotalVolume(dofIdx);
+            Scalar dofVolume = simulator.model().dofTotalVolume(dofIdx);
             assert(dofVolume > 0.0);
             porosity_[dofIdx] = poreVolume/dofVolume;
         }
@@ -2052,15 +2049,17 @@ private:
 
     void initFluidSystem_()
     {
-        const auto& deck = this->simulator().vanguard().deck();
-        const auto& eclState = this->simulator().vanguard().eclState();
+        const auto& simulator = this->simulator();
+        const auto& deck = simulator.vanguard().deck();
+        const auto& eclState = simulator.vanguard().eclState();
 
         FluidSystem::initFromDeck(deck, eclState);
    }
 
     void readInitialCondition_()
     {
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
 
         const auto& deck = vanguard.deck();
         if (!deck.hasKeyword("EQUIL"))
@@ -2074,9 +2073,11 @@ private:
 
     void readEquilInitialCondition_()
     {
+        const auto& simulator = this->simulator();
+
         // initial condition corresponds to hydrostatic conditions.
         typedef Ewoms::EclEquilInitializer<TypeTag> EquilInitializer;
-        EquilInitializer equilInitializer(this->simulator(), *materialLawManager_);
+        EquilInitializer equilInitializer(simulator, *materialLawManager_);
 
         size_t numElems = this->model().numGridDof();
         initialFluidStates_.resize(numElems);
@@ -2089,17 +2090,22 @@ private:
     void readEclRestartSolution_()
     {
         // Set the start time of the simulation
-        const auto& schedule = this->simulator().vanguard().schedule();
-        const auto& eclState = this->simulator().vanguard().eclState();
+        auto& simulator = this->simulator();
+        const auto& schedule = simulator.vanguard().schedule();
+        const auto& eclState = simulator.vanguard().eclState();
         const auto& timeMap = schedule.getTimeMap();
         const auto& initconfig = eclState.getInitConfig();
         int episodeIdx = initconfig.getRestartStep() - 1;
 
-        this->simulator().setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0));
-        this->simulator().setTime(timeMap.getTimePassedUntil(episodeIdx));
-        this->simulator().setEpisodeIndex(episodeIdx);
-        this->simulator().setEpisodeLength(timeMap.getTimeStepLength(episodeIdx));
-        this->simulator().setTimeStepSize(eclWriter_->restartTimeStepSize());
+        simulator.setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0));
+        simulator.setTime(timeMap.getTimePassedUntil(episodeIdx));
+
+        simulator.startNextEpisode(simulator.startTime() + simulator.time(),
+                                   timeMap.getTimeStepLength(episodeIdx));
+        simulator.setEpisodeIndex(episodeIdx);
+
+        Scalar dt = std::min(eclWriter_->restartTimeStepSize(), simulator.episodeLength());
+        simulator.setTimeStepSize(dt);
 
         eclWriter_->beginRestart();
 
@@ -2125,7 +2131,7 @@ private:
         for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
             auto& elemFluidState = initialFluidStates_[elemIdx];
             elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
-            eclWriter_->eclOutputModule().initHysteresisParams(this->simulator(), elemIdx);
+            eclWriter_->eclOutputModule().initHysteresisParams(simulator, elemIdx);
             eclWriter_->eclOutputModule().assignToFluidState(elemFluidState, elemIdx);
 
             processRestartSaturations_(elemFluidState);
@@ -2139,17 +2145,17 @@ private:
             // if we need to restart for polymer molecular weight simulation, we need to add related here
         }
 
-        const int epsiodeIdx = this->simulator().episodeIndex();
-        const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
+        const int epsiodeIdx = simulator.episodeIndex();
+        const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
         if (drsdtActive_())
             // DRSDT is enabled
             for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx)
-                maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*this->simulator().timeStepSize();
+                maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*simulator.timeStepSize();
 
         if (drvdtActive_())
             // DRVDT is enabled
             for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx)
-                maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*this->simulator().timeStepSize();
+                maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*simulator.timeStepSize();
 
         if (tracerModel().numTracers() > 0)
             std::cout << "Warning: Restart is not implemented for the tracer model, it will initialize with initial tracer concentration" << std::endl;
@@ -2159,7 +2165,7 @@ private:
         // need to be correct for stuff like boundary conditions.
         auto& sol = this->model().solution(/*timeIdx=*/0);
         const auto& gridView = this->gridView();
-        ElementContext elemCtx(this->simulator());
+        ElementContext elemCtx(simulator);
         auto elemIt = gridView.template begin</*codim=*/0>();
         const auto& elemEndIt = gridView.template end</*codim=*/0>();
         for (; elemIt != elemEndIt; ++elemIt) {
@@ -2183,7 +2189,7 @@ private:
         initialFluidStates_ = tmpInitialFs;
         // make sure that the stuff which needs to be done at the beginning of an episode
         // is run.
-        this->beginEpisode(/*isOnRestart=*/true);
+        this->beginEpisode();
 
         eclWriter_->endRestart();
     }
@@ -2215,7 +2221,8 @@ private:
 
     void readExplicitInitialCondition_()
     {
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
         const auto& eclState = vanguard.eclState();
         const auto& eclProps = eclState.get3DProperties();
 
@@ -2241,7 +2248,7 @@ private:
 
         initialFluidStates_.resize(numDof);
 
-        const auto& cartSize = this->simulator().vanguard().cartesianDimensions();
+        const auto& cartSize = simulator.vanguard().cartesianDimensions();
         size_t numCartesianCells = cartSize[0] * cartSize[1] * cartSize[2];
 
         std::vector<double> waterSaturationData;
@@ -2359,7 +2366,8 @@ private:
 
     void readBlackoilExtentionsInitialConditions_()
     {
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
         const auto& eclState = vanguard.eclState();
         size_t numDof = this->model().numGridDof();
 
@@ -2405,8 +2413,9 @@ private:
 
         // we need to update the hysteresis data for _all_ elements (i.e., not just the
         // interior ones) to avoid desynchronization of the processes in the parallel case!
-        ElementContext elemCtx(this->simulator());
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        ElementContext elemCtx(simulator);
+        const auto& vanguard = simulator.vanguard();
         auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
         const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
         for (; elemIt != elemEndIt; ++elemIt) {
@@ -2425,8 +2434,9 @@ private:
     void updateMaxPolymerAdsorption_()
     {
         // we need to update the max polymer adsoption data for all elements
-        ElementContext elemCtx(this->simulator());
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        ElementContext elemCtx(simulator);
+        const auto& vanguard = simulator.vanguard();
         auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
         const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
         for (; elemIt != elemEndIt; ++elemIt) {
@@ -2444,14 +2454,15 @@ private:
 
     void updatePvtnum_()
     {
-        const auto& eclState = this->simulator().vanguard().eclState();
+        const auto& simulator = this->simulator();
+        const auto& eclState = simulator.vanguard().eclState();
         const auto& eclProps = eclState.get3DProperties();
 
         if (!eclProps.hasDeckIntGridProperty("PVTNUM"))
             return;
 
         const auto& pvtnumData = eclProps.getIntGridProperty("PVTNUM").getData();
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& vanguard = simulator.vanguard();
 
         unsigned numElems = vanguard.gridView().size(/*codim=*/0);
         pvtnum_.resize(numElems);
@@ -2463,14 +2474,15 @@ private:
 
     void updateSatnum_()
     {
-        const auto& eclState = this->simulator().vanguard().eclState();
+        const auto& simulator = this->simulator();
+        const auto& eclState = simulator.vanguard().eclState();
         const auto& eclProps = eclState.get3DProperties();
 
         if (!eclProps.hasDeckIntGridProperty("SATNUM"))
             return;
 
         const auto& satnumData = eclProps.getIntGridProperty("SATNUM").getData();
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& vanguard = simulator.vanguard();
 
         unsigned numElems = vanguard.gridView().size(/*codim=*/0);
         satnum_.resize(numElems);
@@ -2482,14 +2494,15 @@ private:
 
     void updateMiscnum_()
     {
-        const auto& eclState = this->simulator().vanguard().eclState();
+        const auto& simulator = this->simulator();
+        const auto& eclState = simulator.vanguard().eclState();
         const auto& eclProps = eclState.get3DProperties();
 
         if (!eclProps.hasDeckIntGridProperty("MISCNUM"))
             return;
 
         const auto& miscnumData = eclProps.getIntGridProperty("MISCNUM").getData();
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& vanguard = simulator.vanguard();
 
         unsigned numElems = vanguard.gridView().size(/*codim=*/0);
         miscnum_.resize(numElems);
@@ -2501,14 +2514,15 @@ private:
 
     void updatePlmixnum_()
     {
-        const auto& eclState = this->simulator().vanguard().eclState();
+        const auto& simulator = this->simulator();
+        const auto& eclState = simulator.vanguard().eclState();
         const auto& eclProps = eclState.get3DProperties();
 
         if (!eclProps.hasDeckIntGridProperty("PLMIXNUM"))
             return;
 
         const auto& plmixnumData = eclProps.getIntGridProperty("PLMIXNUM").getData();
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& vanguard = simulator.vanguard();
 
         unsigned numElems = vanguard.gridView().size(/*codim=*/0);
         plmixnum_.resize(numElems);
@@ -2551,7 +2565,8 @@ private:
     void readBoundaryConditions_()
     {
         nonTrivialBoundaryConditions_ = false;
-        const auto& vanguard = this->simulator().vanguard();
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
 
         if (vanguard.deck().hasKeyword("BC")) {
             nonTrivialBoundaryConditions_ = true;
@@ -2684,31 +2699,27 @@ private:
     }
 
     // this method applies the runtime constraints specified via the deck and/or command
-    // line parameters for the size of the next time step size.
+    // line parameters for the size of the next time step.
     Scalar limitNextTimeStepSize_(Scalar dtNext) const
     {
-        const auto& events = this->simulator().vanguard().schedule().getEvents();
-        int episodeIdx = this->simulator().episodeIndex();
+        if (!enableExperiments)
+            return dtNext;
+
+        const auto& simulator = this->simulator();
+        const auto& events = simulator.vanguard().schedule().getEvents();
+        int episodeIdx = simulator.episodeIndex();
 
         // first thing in the morning, limit the time step size to the maximum size
         dtNext = std::min(dtNext, maxTimeStepSize_);
 
-        // use at least slightly more than half of the maximum time step size by default
-        if (dtNext < maxTimeStepSize_ && maxTimeStepSize_ < dtNext*2)
-            dtNext = 1.01 * maxTimeStepSize_/2.0;
-
         Scalar remainingEpisodeTime =
-            this->simulator().episodeStartTime()
-            + this->simulator().episodeLength()
-            - this->simulator().time()
-            - this->simulator().timeStepSize();
-
-        // limit the time step size to the remaining time of the episode
-        dtNext = std::min(dtNext, remainingEpisodeTime);
+            simulator.episodeStartTime() + simulator.episodeLength()
+            - (simulator.startTime() + simulator.time());
+        assert(remainingEpisodeTime >= 0.0);
 
         // if we would have a small amount of time left over in the current episode, make
         // two equal time steps instead of a big and a small one
-        if (dtNext < remainingEpisodeTime*(1.0 - 1e-5) && dtNext*1.25 > remainingEpisodeTime)
+        if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
             // note: limiting to the maximum time step size here is probably not strictly
             // necessary, but it should not hurt and is more fool-proof
             dtNext = std::min(maxTimeStepSize_, remainingEpisodeTime/2.0);
@@ -2768,6 +2779,8 @@ private:
 
     EclWellModel wellModel_;
     EclAquiferModel aquiferModel_;
+
+    bool enableEclOutput_;
     std::unique_ptr<EclWriterType> eclWriter_;
 
     PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
diff --git a/ebos/eclthresholdpressure.hh b/ebos/eclthresholdpressure.hh
index e64463829c..7b17aee86f 100644
--- a/ebos/eclthresholdpressure.hh
+++ b/ebos/eclthresholdpressure.hh
@@ -91,10 +91,6 @@ public:
         enableThresholdPressure_ = false;
     }
 
-
-    void setFromRestart(const std::vector<Scalar>& values)
-    { thpres_ = values; }
-
     /*!
      * \brief Actually compute the threshold pressures over a face as a pre-compute step.
      */
@@ -201,10 +197,21 @@ public:
         return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
     }
 
-    const std::vector<Scalar>& data() const {
-        return thpres_;
-    }
+    /*!
+     * \brief Return the raw array with the threshold pressures
+     *
+     * This is used for the restart capability.
+     */
+    const std::vector<Scalar>& data() const
+    { return thpres_; }
 
+    /*!
+     * \brief Set the threshold pressures from a raw array
+     *
+     * This is used for the restart capability.
+     */
+    void setFromRestart(const std::vector<Scalar>& values)
+    { thpres_ = values; }
 
 private:
     // compute the defaults of the threshold pressures using the initial condition
diff --git a/ebos/eclwellmanager.hh b/ebos/eclwellmanager.hh
index 19d64ac1f7..e7609c287e 100644
--- a/ebos/eclwellmanager.hh
+++ b/ebos/eclwellmanager.hh
@@ -101,9 +101,10 @@ public:
      *
      * I.e., well positions, names etc...
      */
-    void init(const Opm::EclipseState& eclState OPM_UNUSED,
-              const Opm::Schedule& deckSchedule)
+    void init()
     {
+        const Opm::Schedule& deckSchedule = simulator_.vanguard().schedule();
+
         // create the wells which intersect with the current process' grid
         for (size_t deckWellIdx = 0; deckWellIdx < deckSchedule.numWells(); ++deckWellIdx)
         {
@@ -130,8 +131,10 @@ public:
      * \brief This should be called the problem before each simulation
      *        episode to adapt the well controls.
      */
-    void beginEpisode(const Opm::EclipseState& eclState, const Opm::Schedule& deckSchedule, bool wasRestarted=false)
+    void beginEpisode(bool wasRestarted=false)
     {
+        const Opm::EclipseState& eclState = simulator_.vanguard().eclState();
+        const Opm::Schedule& deckSchedule = simulator_.vanguard().schedule();
         unsigned episodeIdx = simulator_.episodeIndex();
 
         WellConnectionsMap wellCompMap;
@@ -589,7 +592,7 @@ public:
     void deserialize(Restarter& res OPM_UNUSED)
     {
         // initialize the wells for the current episode
-        beginEpisode(simulator_.vanguard().eclState(), simulator_.vanguard().schedule(), /*wasRestarted=*/true);
+        beginEpisode(/*wasRestarted=*/true);
     }
 
     /*!
diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh
index 7159bae016..4e7d7f1092 100644
--- a/ebos/eclwriter.hh
+++ b/ebos/eclwriter.hh
@@ -160,11 +160,11 @@ public:
         // output using eclWriter if enabled
         Opm::data::Wells localWellData = simulator_.problem().wellModel().wellData();
 
-        int episodeIdx = simulator_.episodeIndex() + 1;
+        int reportStepNum = simulator_.episodeIndex() + 1;
         const auto& gridView = simulator_.vanguard().gridView();
         int numElements = gridView.size(/*codim=*/0);
         bool log = collectToIORank_.isIORank();
-        eclOutputModule_.allocBuffers(numElements, episodeIdx, isSubStep, log);
+        eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log);
 
         ElementContext elemCtx(simulator_);
         ElementIterator elemIt = gridView.template begin</*codim=*/0>();
@@ -184,7 +184,7 @@ public:
 
         // add cell data to perforations for Rft output
         if (!isSubStep)
-            eclOutputModule_.addRftDataToWells(localWellData, episodeIdx);
+            eclOutputModule_.addRftDataToWells(localWellData, reportStepNum);
 
         if (collectToIORank_.isParallel())
             collectToIORank_.collect(localCellData, eclOutputModule_.getBlockData(), localWellData);
@@ -221,7 +221,7 @@ public:
 
             // first, create a tasklet to write the data for the current time step to disk
             auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(*eclIO_,
-                                                                     episodeIdx,
+                                                                     reportStepNum,
                                                                      isSubStep,
                                                                      curTime,
                                                                      restartValue,
@@ -437,7 +437,7 @@ private:
         : public TaskletInterface
     {
         Opm::EclipseIO& eclIO_;
-        int episodeIdx_;
+        int reportStepNum_;
         bool isSubStep_;
         double secondsElapsed_;
         Opm::RestartValue restartValue_;
@@ -447,7 +447,7 @@ private:
         bool writeDoublePrecision_;
 
         explicit EclWriteTasklet(Opm::EclipseIO& eclIO,
-                                 int episodeIdx,
+                                 int reportStepNum,
                                  bool isSubStep,
                                  double secondsElapsed,
                                  Opm::RestartValue restartValue,
@@ -456,7 +456,7 @@ private:
                                  const std::map<std::pair<std::string, int>, double>& blockSummaryValues,
                                  bool writeDoublePrecision)
             : eclIO_(eclIO)
-            , episodeIdx_(episodeIdx)
+            , reportStepNum_(reportStepNum)
             , isSubStep_(isSubStep)
             , secondsElapsed_(secondsElapsed)
             , restartValue_(restartValue)
@@ -469,7 +469,7 @@ private:
         // callback to eclIO serial writeTimeStep method
         void run()
         {
-            eclIO_.writeTimeStep(episodeIdx_,
+            eclIO_.writeTimeStep(reportStepNum_,
                                  isSubStep_,
                                  secondsElapsed_,
                                  restartValue_,
diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp
index 69049822c3..0ec671752a 100644
--- a/opm/autodiff/BlackoilModelEbos.hpp
+++ b/opm/autodiff/BlackoilModelEbos.hpp
@@ -182,7 +182,6 @@ namespace Opm {
         /// \param[in] timer                  simulation timer
         void prepareStep(const SimulatorTimerInterface& timer)
         {
-
             // update the solution variables in ebos
             if ( timer.lastStepFailed() ) {
                 ebosSimulator_.model().updateFailed();
@@ -194,13 +193,8 @@ namespace Opm {
             // know the report step/episode index because of timing dependend data
             // despide the fact that flow uses its own time stepper. (The length of the
             // episode does not matter, though.)
-            Scalar t = timer.simulationTimeElapsed();
-            ebosSimulator_.startNextEpisode(/*episodeStartTime=*/t, /*episodeLength=*/1e30);
-            ebosSimulator_.setEpisodeIndex(timer.reportStepNum());
-            ebosSimulator_.setTime(t);
+            ebosSimulator_.setTime(timer.simulationTimeElapsed());
             ebosSimulator_.setTimeStepSize(timer.currentStepLength());
-            ebosSimulator_.setTimeStepIndex(ebosSimulator_.timeStepIndex() + 1);
-
             ebosSimulator_.problem().beginTimeStep();
 
             unsigned numDof = ebosSimulator_.model().numGridDof();
diff --git a/opm/autodiff/BlackoilWellModel.hpp b/opm/autodiff/BlackoilWellModel.hpp
index bb98bd4c09..59be94b604 100644
--- a/opm/autodiff/BlackoilWellModel.hpp
+++ b/opm/autodiff/BlackoilWellModel.hpp
@@ -112,7 +112,7 @@ namespace Opm {
 
             BlackoilWellModel(Simulator& ebosSimulator);
 
-            void init(const Opm::EclipseState& eclState, const Opm::Schedule& schedule);
+            void init();
 
             /////////////
             // <eWoms auxiliary module stuff>
@@ -153,18 +153,9 @@ namespace Opm {
                 // TODO (?)
             }
 
-            void beginEpisode(const Opm::EclipseState& /* eclState */,
-                              const Opm::Schedule& /* schedule */,
-                              bool isRestart)
+            void beginEpisode()
             {
-                size_t episodeIdx = ebosSimulator_.episodeIndex();
-                // beginEpisode in eclProblem advances the episode index
-                // we don't want this when we are at the beginning of an
-                // restart.
-                if (isRestart)
-                    episodeIdx -= 1;
-
-                beginReportStep(episodeIdx);
+                beginReportStep(ebosSimulator_.episodeIndex());
             }
 
             void beginTimeStep();
diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp
index 6c74685ffc..c2381d339c 100644
--- a/opm/autodiff/BlackoilWellModel_impl.hpp
+++ b/opm/autodiff/BlackoilWellModel_impl.hpp
@@ -38,8 +38,10 @@ namespace Opm {
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
-    init(const Opm::EclipseState& eclState, const Opm::Schedule& /* schedule */)
+    init()
     {
+        const Opm::EclipseState& eclState = ebosSimulator_.vanguard().eclState();
+
         gravity_ = ebosSimulator_.problem().gravity()[2];
 
         extractLegacyCellPvtRegionIndex_();
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
index 2b8434f54a..f6e7b44a6b 100644
--- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
@@ -133,6 +133,8 @@ public:
     {
         failureReport_ = SimulatorReport();
 
+        ebosSimulator_.setEpisodeIndex(-1);
+
         // handle restarts
         std::unique_ptr<RestartValue> restartValues;
         if (isRestart()) {
@@ -164,6 +166,21 @@ public:
 
             double suggestedStepSize = -1.0;
             if (isRestart()) {
+                // Set the start time of the simulation
+                const auto& schedule = ebosSimulator_.vanguard().schedule();
+                const auto& eclState = ebosSimulator_.vanguard().eclState();
+                const auto& timeMap = schedule.getTimeMap();
+                const auto& initconfig = eclState.getInitConfig();
+                int episodeIdx = initconfig.getRestartStep() - 1;
+
+                ebosSimulator_.setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0));
+                ebosSimulator_.setTime(timeMap.getTimePassedUntil(episodeIdx));
+
+                ebosSimulator_.startNextEpisode(ebosSimulator_.startTime() + ebosSimulator_.time(),
+                                                timeMap.getTimeStepLength(episodeIdx));
+                ebosSimulator_.setEpisodeIndex(episodeIdx);
+
+
                 // This is a restart, determine the time step size from the restart data
                 if (restartValues->hasExtra("OPMEXTRA")) {
                     std::vector<double> opmextra = restartValues->getExtra("OPMEXTRA");
@@ -206,12 +223,18 @@ public:
                 Dune::Timer perfTimer;
                 perfTimer.start();
 
+                ebosSimulator_.setEpisodeIndex(-1);
+                ebosSimulator_.setEpisodeLength(0.0);
+                ebosSimulator_.setTimeStepSize(0.0);
+
                 wellModel_().beginReportStep(timer.currentStepNum());
                 ebosSimulator_.problem().writeOutput(false);
 
                 report.output_write_time += perfTimer.stop();
             }
 
+            ebosSimulator_.setEpisodeIndex(timer.currentStepNum());
+
             // Run a multiple steps of the solver depending on the time step control.
             solverTimer.start();
 
@@ -265,6 +288,14 @@ public:
                 }
             }
 
+            // write simulation state at the report stage
+            Dune::Timer perfTimer;
+            perfTimer.start();
+            const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
+            ebosSimulator_.problem().setNextTimeStepSize(nextstep);
+            ebosSimulator_.problem().writeOutput(false);
+            report.output_write_time += perfTimer.stop();
+
             solver->model().endReportStep();
 
             // take time that was used to solve system for this reportStep
@@ -284,14 +315,6 @@ public:
                 }
             }
 
-            // write simulation state at the report stage
-            Dune::Timer perfTimer;
-            perfTimer.start();
-            const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
-            ebosSimulator_.problem().setNextTimeStepSize(nextstep);
-            ebosSimulator_.problem().writeOutput(false);
-            report.output_write_time += perfTimer.stop();
-
             if (terminalOutput_) {
                 std::string msg =
                     "Time step took " + std::to_string(solverTimer.secsSinceStart()) + " seconds; "
diff --git a/tests/equil_base.DATA b/tests/equil_base.DATA
index 740ae08502..91b13c2dbc 100644
--- a/tests/equil_base.DATA
+++ b/tests/equil_base.DATA
@@ -113,3 +113,6 @@ PRESSURE
 SUMMARY
 
 SCHEDULE
+
+TSTEP
+1 /
diff --git a/tests/equil_capillary.DATA b/tests/equil_capillary.DATA
index 4289102e5c..e53401b3ea 100644
--- a/tests/equil_capillary.DATA
+++ b/tests/equil_capillary.DATA
@@ -84,4 +84,6 @@ EQUIL
 
 -------------------------------------
 SCHEDULE
--- empty section
+
+TSTEP
+1 /
diff --git a/tests/equil_capillary_overlap.DATA b/tests/equil_capillary_overlap.DATA
index b59e422028..1c5791afd0 100644
--- a/tests/equil_capillary_overlap.DATA
+++ b/tests/equil_capillary_overlap.DATA
@@ -121,6 +121,9 @@ SEPARATE
 
 SCHEDULE  ======
 
+TSTEP
+1 /
+
 RPTSCHED
 'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=3' 'NEWTON=2' /
 
diff --git a/tests/equil_capillary_swatinit.DATA b/tests/equil_capillary_swatinit.DATA
index f4e693a58d..73a2592865 100644
--- a/tests/equil_capillary_swatinit.DATA
+++ b/tests/equil_capillary_swatinit.DATA
@@ -99,4 +99,7 @@ RPTSOL
 /
 -------------------------------------
 SCHEDULE
--- empty section
+
+TSTEP
+1 /
+
diff --git a/tests/equil_deadfluids.DATA b/tests/equil_deadfluids.DATA
index e0d63cb374..f33f1a2f50 100644
--- a/tests/equil_deadfluids.DATA
+++ b/tests/equil_deadfluids.DATA
@@ -85,4 +85,7 @@ EQUIL
 
 -------------------------------------
 SCHEDULE
--- empty section
+
+TSTEP
+1 /
+
diff --git a/tests/equil_livegas.DATA b/tests/equil_livegas.DATA
index a0676f52c7..4717c0481f 100644
--- a/tests/equil_livegas.DATA
+++ b/tests/equil_livegas.DATA
@@ -124,6 +124,9 @@ SEPARATE
 
 SCHEDULE  ======
 
+TSTEP
+1 /
+
 RPTSCHED
 'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
 
diff --git a/tests/equil_liveoil.DATA b/tests/equil_liveoil.DATA
index e1b417c46f..c64bcbe37f 100644
--- a/tests/equil_liveoil.DATA
+++ b/tests/equil_liveoil.DATA
@@ -133,6 +133,9 @@ SEPARATE
 
 SCHEDULE  ======
 
+TSTEP
+1 /
+
 RPTSCHED
 'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
 
diff --git a/tests/equil_pbvd_and_pdvd.DATA b/tests/equil_pbvd_and_pdvd.DATA
index 59da42d692..78093c16ec 100644
--- a/tests/equil_pbvd_and_pdvd.DATA
+++ b/tests/equil_pbvd_and_pdvd.DATA
@@ -145,6 +145,9 @@ SEPARATE
 
 SCHEDULE  ======
 
+TSTEP
+1 /
+
 RPTSCHED
 'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
 
diff --git a/tests/equil_rsvd_and_rvvd.DATA b/tests/equil_rsvd_and_rvvd.DATA
index 3076524e34..bfc965c259 100644
--- a/tests/equil_rsvd_and_rvvd.DATA
+++ b/tests/equil_rsvd_and_rvvd.DATA
@@ -144,6 +144,9 @@ SEPARATE
 
 SCHEDULE  ======
 
+TSTEP
+1 /
+
 RPTSCHED
 'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
 
diff --git a/tests/test_ecl_output.cc b/tests/test_ecl_output.cc
index 064ff74e24..7e8bb6b239 100644
--- a/tests/test_ecl_output.cc
+++ b/tests/test_ecl_output.cc
@@ -131,6 +131,7 @@ void test_summary()
     simulator->model().applyInitialSolution();
     Opm::data::Wells dw;
     bool substep = false;
+    simulator->startNextEpisode(0.0, 1e30);
     simulator->setEpisodeIndex(0);
     eclWriter->writeOutput(substep);
     simulator->setEpisodeIndex(1);