if (nNeSubCycles > 1)
{
    dimensionedScalar totalDeltaT = runTime.deltaT();

    for
    (
        subCycle<volScalarField> neSubCycle(ne, nNeSubCycles);
        !(++neSubCycle).end();
    )
    {
        #include "neEqn.H"
        #include "PhiEqn.H"
        // rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
    }

    // rhoPhi = rhoPhiSum;
}
else
{
    #include "neEqn.H"
    #include "PhiEqn.H"
}