Info<< "\nCreating thermophysical model for Fluid mesh\n" << endl;
#include "readGravitationalAcceleration.H"

autoPtr<combustionModels::rhoCombustionModel> reaction
(
    combustionModels::rhoCombustionModel::New(mesh)
);

rhoReactionThermo& thermo = reaction->thermo();
thermo.validate(args.executable(), "h", "e");

basicSpecieMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();

word inertSpecie(thermo.lookup("inertSpecie"));

volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh
    ),
    thermo.rho()
);

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


volScalarField& p = thermo.p();
volScalarField& T = thermo.T();

#include "compressibleCreatePhi.H"
dimensionedScalar rhoMax
(
    dimensionedScalar::lookupOrDefault
    (
        "rhoMax",
        simple.dict(),
        dimDensity,
        GREAT
    )
);

dimensionedScalar rhoMin
(
    dimensionedScalar::lookupOrDefault
    (
        "rhoMin",
        simple.dict(),
        dimDensity,
        0
    )
);

Info << "Creating turbulence model.\n" << nl;
autoPtr<compressible::turbulenceModel> turbulence
(
    compressible::turbulenceModel::New
    (
        rho,
        U,
        phi,
        thermo
    )
);

// Set the turbulence into the reaction model
reaction->setTurbulence(turbulence());


Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
    IOobject
    (
        "dpdt",
        runTime.timeName(),
        mesh
    ),
    mesh,
    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
);

Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));


multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;

forAll(Y, i)
{
    fields.add(Y[i]);
}
fields.add(thermo.he());

/***************************************************************/
//SLFM-related Fields 
Info<<"Creating field mf, mfVar, SDR\n"<<endl;

volScalarField mf //mixture fraction
(
    IOobject
    (
        "mf",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

volScalarField mfVar //mixture fraction variance
(
    IOobject
    (
        "mfVar",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

volScalarField SDR //scalar dissipation rate
(
    IOobject
    (
        "SDR",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("SDR", dimEnergy/dimTime/dimEnergy, 0.0)
);

volScalarField jlc //eta-index, lower
(
    IOobject
    (
        "jlc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("jlc", dimless, 0.0)
);

volScalarField jhc //eta-index, higher
(
    IOobject
    (
        "jhc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("jhc", dimless, 0.0)
);

volScalarField jfc //eta-factor
(
    IOobject
    (
        "jfc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("jfc", dimless, 0.0)
);

volScalarField vlc //var-index, lower
(
    IOobject
    (
        "vlc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("vlc", dimless, 0.0)
);

volScalarField vhc //var-index, higher
(
    IOobject
    (
        "vhc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("vhc", dimless, 0.0)
);

volScalarField vfc //var-factor
(
    IOobject
    (
        "vfc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("vfc", dimless, 0.0)
);

volScalarField nlc //sdr-index, lower
(
    IOobject
    (
        "nlc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("nlc", dimless, 0.0)
);

volScalarField nhc //sdr-index, higher
(
    IOobject
    (
        "nhc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("nhc", dimless, 0.0)
);

volScalarField nfc //sdr-factor
(
    IOobject
    (
        "nfc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("nfc", dimless, 0.0)
);

//SLFM property read

IOdictionary SLFMdict //SLFM properties dictionary
(
    IOobject
    (
        "SLFMdict",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    )
);

const label Ysize = Y.size();
const wordList postSpecieNames(SLFMdict.lookup("postSpecieNames"));
labelList postSpecieIndices(postSpecieNames.size(),-1);
forAll(postSpecieNames, i)
{
    const word &spName (postSpecieNames[i]);
    const label spI (composition.species()[spName]);
    postSpecieIndices[i] = spI;
}
scalar stoiMF(readScalar(SLFMdict.lookup("stoiMF")));
scalar stretchFac(readScalar(SLFMdict.lookup("stretchFactor")));
scalar lowerN(readScalar(SLFMdict.lookup("lowerN")));
scalar upperN(readScalar(SLFMdict.lookup("upperN")));
scalar Sc(readScalar(SLFMdict.lookup("Sc")));
const label etamax(lowerN+upperN+1);
word outletName(SLFMdict.lookup("outletName"));

labelList nIntervals (SLFMdict.lookup("detailedN"));
scalarField etaGrid(sum(nIntervals)+1, 0.0);

{
    scalarList etas (SLFMdict.lookup("detailedEta"));

    forAll(nIntervals, i)
    {
        const label itv = nIntervals[i];
        const labelList::subList prev(nIntervals, i);
        const label baseI = sum(prev);

        const scalar low = etas[i];
        const scalar upp = etas[i+1];
        const scalar delta = (upp - low) / scalar(itv);

        for (label j = 0; j < itv; j++)
        {
            etaGrid[baseI+j] = low + delta * j;
        }
    }
    etaGrid.last() = etas.last();
}

Info << etaGrid << endl;

scalarField etaValue(etamax,0);
scalarField etaIndex(etamax,0);
scalarField deta_h(etamax,0);
scalarField deta_l(etamax,0);
scalarField deta_m(etamax,0);

#include "setEtaSpace.H"

Info<<"\nNumber of eta-point = "<<etamax<<endl;


scalar NVar(readScalar(SLFMdict.lookup("NVar")));
scalarField varValue(NVar+1,0);
scalarField varIndex(NVar+1,0);

#include "setVarSpace.H"



PtrList<volScalarField> Neta(etamax);

for(label j=0 ; j<etamax ; j++)
{


    word Netai = "Neta_"+Foam::name(j);
    Neta.set
    (
        j,
        new volScalarField
        (
            IOobject
            (
                Netai,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar(Netai, dimless, 0.0)
        )
    );
}

tmp<fv::convectionScheme<scalar> > mvConvection
(
    fv::convectionScheme<scalar>::New
    (
        mesh,
        fields,
        phi,
        mesh.divScheme("div(phi,Yi_h)")
    )
);

#include "createMRF.H"