diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
index 62729c9b0b..dc0fd3769a 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -632,7 +632,7 @@ void Foam::multiphaseSystem::solve()
                 fv::localEulerDdt::localRSubDeltaT(mesh_, nAlphaSubCycles);
         }
 
-        PtrList<volScalarField> alpha0s(phases().size());
+        List<volScalarField*> alphaPtrs(phases().size());
         PtrList<surfaceScalarField> alphaPhiSums(phases().size());
 
         forAll(phases(), phasei)
@@ -640,11 +640,7 @@ void Foam::multiphaseSystem::solve()
             phaseModel& phase = phases()[phasei];
             volScalarField& alpha = phase;
 
-            alpha0s.set
-            (
-                phasei,
-                new volScalarField(alpha.oldTime())
-            );
+            alphaPtrs[phasei] = &alpha;
 
             alphaPhiSums.set
             (
@@ -665,9 +661,9 @@ void Foam::multiphaseSystem::solve()
 
         for
         (
-            subCycleTime alphaSubCycle
+            subCycle<volScalarField, subCycleFields> alphaSubCycle
             (
-                const_cast<Time&>(runTime),
+                alphaPtrs,
                 nAlphaSubCycles
             );
             !(++alphaSubCycle).end();
@@ -686,17 +682,7 @@ void Foam::multiphaseSystem::solve()
             phaseModel& phase = phases()[phasei];
             if (phase.stationary()) continue;
 
-            volScalarField& alpha = phase;
-
             phase.alphaPhiRef() = alphaPhiSums[phasei]/nAlphaSubCycles;
-
-            // Correct the time index of the field
-            // to correspond to the global time
-            alpha.timeIndex() = runTime.timeIndex();
-
-            // Reset the old-time field value
-            alpha.oldTime() = alpha0s[phasei];
-            alpha.oldTime().timeIndex() = runTime.timeIndex();
         }
     }
     else
diff --git a/src/OpenFOAM/algorithms/subCycle/subCycle.H b/src/OpenFOAM/algorithms/subCycle/subCycle.H
index 6cbad550a2..a740086131 100644
--- a/src/OpenFOAM/algorithms/subCycle/subCycle.H
+++ b/src/OpenFOAM/algorithms/subCycle/subCycle.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Class
     Foam::subCycle
 
 Description
-    Perform a subCycleTime on a field
+    Perform a subCycleTime on a field or list of fields.
 
 \*---------------------------------------------------------------------------*/
 
@@ -62,9 +62,12 @@ class subCycleField
 
 public:
 
+    //- Type of the field
+    typedef GeometricField FieldsType;
+
     // Constructors
 
-        //- Construct field and number of sub-cycles
+        //- Construct from field and number of sub-cycles
         subCycleField(GeometricField& gf)
         :
             gf_(gf),
@@ -80,10 +83,15 @@ public:
         gf0_ = gf_0_;
 
         // Correct the time index of the field to correspond to the global time
-        gf_.timeIndex() = gf_.time().timeIndex();
-        gf0_.timeIndex() = gf_.time().timeIndex();
+        gf_.timeIndex() = time().timeIndex();
+        gf0_.timeIndex() = time().timeIndex();
     }
 
+    //- Access to time
+    const Time& time() const
+    {
+        return gf_.time();
+    }
 
     //- Correct the time index of the field to correspond to
     //  the sub-cycling time.
@@ -93,8 +101,98 @@ public:
     //  outer iteration
     void updateTimeIndex()
     {
-        gf_.timeIndex() = gf_.time().timeIndex() + 1;
-        gf0_.timeIndex() = gf_.time().timeIndex() + 1;
+        gf_.timeIndex() = time().timeIndex() + 1;
+        gf0_.timeIndex() = time().timeIndex() + 1;
+    }
+};
+
+
+/*---------------------------------------------------------------------------*\
+                         Class subCycleFields Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class GeometricField>
+class subCycleFields
+{
+    // Private data
+
+        //- List of pointers to the fields being sub-cycled
+        List<GeometricField*> gfPtrs_;
+
+        //- List of pointers to the fields old-time field being sub-cycled
+        //  Needed to avoid calls to oldTime() which may cause
+        //  unexpected updates of the old-time field
+        List<GeometricField*> gf0Ptrs_;
+
+        //- Copy of the "real" old-time value of the fields
+        PtrList<GeometricField> gf_0Ptrs_;
+
+
+public:
+
+    //- Type of the list of fields
+    typedef List<GeometricField*> FieldsType;
+
+    // Constructors
+
+        //- Construct from field list and number of sub-cycles
+        subCycleFields(List<GeometricField*>& gfPtrs)
+        :
+            gfPtrs_(gfPtrs),
+            gf0Ptrs_(gfPtrs.size()),
+            gf_0Ptrs_(gfPtrs.size())
+        {
+            forAll(gfPtrs_, i)
+            {
+                gf0Ptrs_[i] = &gfPtrs_[i]->oldTime();
+
+                gf_0Ptrs_.set
+                (
+                    i,
+                    new volScalarField
+                    (
+                        gf0Ptrs_[i]->name() + "_",
+                        *gf0Ptrs_[i]
+                    )
+                );
+            }
+        }
+
+
+    //- Destructor
+    ~subCycleFields()
+    {
+        forAll(gfPtrs_, i)
+        {
+            // Reset the old-time fields
+            *gf0Ptrs_[i] = gf_0Ptrs_[i];
+
+            // Correct the time index of the fields
+            // to correspond to the global time
+            gfPtrs_[i]->timeIndex() = time().timeIndex();
+            gf0Ptrs_[i]->timeIndex() = time().timeIndex();
+        }
+    }
+
+    //- Access to time
+    const Time& time() const
+    {
+        return gfPtrs_[0]->time();
+    }
+
+    //- Correct the time index of the fields to correspond to
+    //  the sub-cycling time.
+    //
+    //  The time index is incremented to protect the old-time value from
+    //  being updated at the beginning of the time-loop in the case of
+    //  outer iteration
+    void updateTimeIndex()
+    {
+        forAll(gfPtrs_, i)
+        {
+            gfPtrs_[i]->timeIndex() = time().timeIndex() + 1;
+            gf0Ptrs_[i]->timeIndex() = time().timeIndex() + 1;
+        }
     }
 };
 
@@ -103,10 +201,14 @@ public:
                           Class subCycle Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class GeometricField>
+template
+<
+    class GeometricField,
+    template<class> class SubCycleField = subCycleField
+>
 class subCycle
 :
-    public subCycleField<GeometricField>,
+    public SubCycleField<GeometricField>,
     public subCycleTime
 {
     // Private Member Functions
@@ -123,10 +225,14 @@ public:
     // Constructors
 
         //- Construct field and number of sub-cycles
-        subCycle(GeometricField& gf, const label nSubCycles)
+        subCycle
+        (
+            typename SubCycleField<GeometricField>::FieldsType& gf,
+            const label nSubCycles
+        )
         :
-            subCycleField<GeometricField>(gf),
-            subCycleTime(const_cast<Time&>(gf.time()), nSubCycles)
+            SubCycleField<GeometricField>(gf),
+            subCycleTime(const_cast<Time&>(this->time()), nSubCycles)
         {
             // Update the field time index to correspond to the sub-cycle time
             this->updateTimeIndex();