Info<< "Reading transportProperties\n" << endl;

    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    );

    autoPtr<phaseModel> phase1 = phaseModel::New
    (
        mesh,
        transportProperties,
        "1"
    );

    autoPtr<phaseModel> phase2 = phaseModel::New
    (
        mesh,
        transportProperties,
        "2"
    );

    volVectorField& U1 = phase1->U();
    surfaceScalarField& phi1 = phase1->phi();
    const dimensionedScalar& nu1 = phase1->nu();
    const dimensionedScalar& k1 = phase1->kappa();
    const dimensionedScalar& Cp1 = phase1->Cp();
    dimensionedScalar rho10
    (
        "rho0",
        dimDensity,
        transportProperties.subDict("phase1").lookup("rho0")
    );
    dimensionedScalar R1
    (
        "R",
        dimensionSet(0, 2, -2, -1, 0),
        transportProperties.subDict("phase1").lookup("R")
    );

    volVectorField& U2 = phase2->U();
    surfaceScalarField& phi2 = phase2->phi();
    const dimensionedScalar& nu2 = phase2->nu();
    const dimensionedScalar& k2 = phase2->kappa();
    const dimensionedScalar& Cp2 = phase2->Cp();
    dimensionedScalar rho20
    (
        "rho0",
        dimDensity,
        transportProperties.subDict("phase2").lookup("rho0")
    );
    dimensionedScalar R2
    (
        "R",
        dimensionSet(0, 2, -2, -1, 0),
        transportProperties.subDict("phase2").lookup("R")
    );

    dimensionedScalar pMin
    (
        "pMin",
        dimPressure,
        transportProperties.lookup("pMin")
    );


    Info<< "Reading field T1" << endl;
    volScalarField T1
    (
        IOobject
        (
            "T1",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info<< "Reading field T2" << endl;
    volScalarField T2
    (
        IOobject
        (
            "T2",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    volScalarField psi1
    (
        IOobject
        (
            "psi1",
            runTime.timeName(),
            mesh
        ),
        1.0/(R1*T1)
    );

    volScalarField psi2
    (
        IOobject
        (
            "psi2",
            runTime.timeName(),
            mesh
        ),
        1.0/(R2*T2)
    );
    psi2.oldTime();

    Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    volScalarField rho1("rho1", rho10 + psi1*p);
    volScalarField rho2("rho2", rho20 + psi2*p);


    Info<< "Reading field alpha1\n" << endl;
    volScalarField alpha1
    (
        IOobject
        (
            "alpha1",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    volScalarField alpha2
    (
        IOobject
        (
            "alpha2",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        scalar(1) - alpha1
    );

    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        alpha1*U1 + alpha2*U2
    );

    dimensionedScalar Cvm
    (
        "Cvm",
        dimless,
        transportProperties.lookup("Cvm")
    );

    dimensionedScalar Cl
    (
        "Cl",
        dimless,
        transportProperties.lookup("Cl")
    );

    dimensionedScalar Ct
    (
        "Ct",
        dimless,
        transportProperties.lookup("Ct")
    );

    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::interpolate(alpha1)*phi1 + fvc::interpolate(alpha2)*phi2
    );

    volScalarField rho
    (
        IOobject
        (
            "rho",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        alpha1*rho1 + alpha2*rho2
    );

    #include "createRASTurbulence.H"

    Info<< "Calculating field DDtU1 and DDtU2\n" << endl;

    volVectorField DDtU1
    (
        fvc::ddt(U1)
      + fvc::div(phi1, U1)
      - fvc::div(phi1)*U1
    );

    volVectorField DDtU2
    (
        fvc::ddt(U2)
      + fvc::div(phi2, U2)
      - fvc::div(phi2)*U2
    );


    Info<< "Calculating field g.h\n" << endl;
    volScalarField gh("gh", g & mesh.C());

    IOdictionary interfacialProperties
    (
        IOobject
        (
            "interfacialProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    );

    autoPtr<dragModel> drag1 = dragModel::New
    (
        interfacialProperties,
        alpha1,
        phase1,
        phase2
    );

    autoPtr<dragModel> drag2 = dragModel::New
    (
        interfacialProperties,
        alpha2,
        phase2,
        phase1
    );

    autoPtr<heatTransferModel> heatTransfer1 = heatTransferModel::New
    (
        interfacialProperties,
        alpha1,
        phase1,
        phase2
    );

    autoPtr<heatTransferModel> heatTransfer2 = heatTransferModel::New
    (
        interfacialProperties,
        alpha2,
        phase2,
        phase1
    );

    word dispersedPhase(interfacialProperties.lookup("dispersedPhase"));

    if
    (
        !(
            dispersedPhase == "1"
         || dispersedPhase == "2"
         || dispersedPhase == "both"
        )
    )
    {
        FatalErrorIn(args.executable())
            << "invalid dispersedPhase " << dispersedPhase
            << exit(FatalError);
    }

    Info << "dispersedPhase is " << dispersedPhase << endl;

    scalar minInterfaceAlpha
    (
        readScalar
        (
            interfacialProperties.lookup("minInterfaceAlpha")
        )
    );

    kineticTheoryModel kineticTheory
    (
        phase1,
        U2,
        alpha1,
        drag1
    );

    surfaceScalarField rAU1f
    (
        IOobject
        (
            "rAU1f",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar("zero", dimensionSet(0, 0, 1, 0, 0), 0.0)
    );

    surfaceScalarField ppMagf
    (
        IOobject
        (
            "ppMagf",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
    );


    label pRefCell = 0;
    scalar pRefValue = 0.0;
    setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);


    volScalarField dgdt
    (
        pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
    );


    Info<< "Creating field dpdt\n" << endl;
    volScalarField dpdt("dpdt", fvc::ddt(p));

    Info<< "Creating field kinetic energy K\n" << endl;
    volScalarField K1("K1", 0.5*magSqr(U1));
    volScalarField K2("K2", 0.5*magSqr(U2));