
{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");

    surfaceScalarField phir
    (
        IOobject
        (
            "phir",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf() //phase fraction * flux/cellsurface * cell face unit interface normal flux 
/* copied from original interMixingFoam but should probably be min(interface.cAlpha()*mag(phi/mesh.magSf()), max(mag(phi/mesh.magSf())) * interface.nHatf())*/
    );
    
    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
    {
        // Create the limiter to be used for all phase-fractions
        scalarField allLambda(mesh.nFaces(), 1.0);

        // Split the limiter into a surfaceScalarField
        slicedSurfaceScalarField lambda
        (
            IOobject
            (
                "lambda",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            mesh,
            dimless,
            allLambda,
            false   // Use slices for the couples
        );
        // Create the complete convection flux for alpha1
        surfaceScalarField phiAlpha1
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, alpha2, alpharScheme),
                alpha1,
                alpharScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, alpha3, alpharScheme),
                alpha1,
                alpharScheme
            )
        );
        // Create the complete flux for alpha2 
        surfaceScalarField phiAlpha2
        (
            fvc::flux
            (
                phi,
                alpha2,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(phir, alpha1, alpharScheme),
                alpha2,
                alpharScheme
            )
        );
        // Create the bounded (upwind) flux for alpha2
        surfaceScalarField phiAlpha2BD
        (
            upwind<scalar>(mesh, phi).flux(alpha2)
        );
        // Calculate the flux correction for alpha2
        phiAlpha2 -= phiAlpha2BD;
        // Calculate the limiter for alpha2
        MULES::limiter
        (
            allLambda,
            1.0/runTime.deltaT().value(),
            geometricOneField(),
            alpha2,
            phiAlpha2BD,
            phiAlpha2,
            zeroField(),
            zeroField(),
            1,
            0,
            3
        );

        // Construct the limited fluxes
        phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;

        // Create the complete flux for alpha3
        surfaceScalarField phiAlpha3
        (
            fvc::flux
            (
                phi,
                alpha3,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(phir, alpha1, alpharScheme),
                alpha3,
                alpharScheme
            )
        );
        // Create the bounded (upwind) flux for alpha3
        surfaceScalarField phiAlpha3BD
        (
            upwind<scalar>(mesh, phi).flux(alpha3)
        );
        // Calculate the flux correction for alpha3
        phiAlpha3 -= phiAlpha3BD;
        // Further limit the limiter for alpha3
        MULES::limiter
        (
            allLambda,
            1.0/runTime.deltaT().value(),
            geometricOneField(),
            alpha3,
            phiAlpha3BD,
            phiAlpha3,
            zeroField(),
            zeroField(),
            1,
            0,
            3
        );
        // Construct the limited fluxes
        phiAlpha3 = phiAlpha3BD + lambda*phiAlpha3;
	volScalarField Dc23(D23*max(alpha3, scalar(0))*pos(alpha2));
        volScalarField Dc32(D23*max(alpha2, scalar(0))*pos(alpha3));
	
		// Add the diffusive flux for alpha3->alpha2
		// Solve for alpha2
		fvScalarMatrix alpha2Eqn
		(
		    fvm::ddt(alpha2)
		  + fvc::div(phiAlpha2) 
		//see multiphaseEulerFoam=>multiphaseSystem: surfaceScalarField& phiAlphaCorr = fvc::flux(phi_, phase1, "div(phi," + alpha1.name() + ')' ) += fvc::flux(-fvc::flux(-phir, phase2, phirScheme),
		  - fvm::laplacian(Dc32 + Dc23 , alpha2)
		);
		alpha2Eqn.solve();
		// Add the diffusive flux for alpha2->alpha3
		
		alpha2 = max(alpha2, scalar(0.0));
		fvScalarMatrix alpha3Eqn
		(
		    fvm::ddt(alpha3)
		  + fvc::div(phiAlpha3)
		  - fvm::laplacian(Dc23 + Dc32  , alpha3)//
		);
		alpha3Eqn.solve();
		alpha3 = min(alpha3, 1.0 - alpha2); //was called after alpha3 = max(alpha3, scalar(0.0)); before, but caused sometimes instability though keeping mass
		alpha3 = max(alpha3, scalar(0.0));
		// Info<< "min alpha3 volume fraction = " << min(alpha3).value();
		// Construct the complete mass flux */
		rhoPhi =
		      phiAlpha3 * rho3
		    + phiAlpha2 * rho2
		    + phiAlpha1 * rho1; 
	
        alpha1 = 1.0 - alpha2 - alpha3;
    }

    Info<< "Mud phase volume fraction* = "
        << alpha2.weightedAverage(mesh.V()).value()
        << "  Min(alpha2) = " << min(alpha2).value()
        << "  Max(alpha2) = " << max(alpha2).value()
        << endl;

    Info<< "Granular phase volume fraction* = "
        << alpha3.weightedAverage(mesh.V()).value()
        << "  Min(alpha3) = " << min(alpha3).value()
        << "  Max(alpha3) = " << max(alpha3).value()
        << endl;
	

	

}
