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

combustionModels::rhoCombustionModel* pRxn = (combustionModels::rhoCombustionModel::New(mesh)).ptr();
autoPtr<combustionModels::rhoChemistryCombustion> combustion
(
    dynamic_cast <combustionModels::rhoChemistryCombustion*> (pRxn)
);
autoPtr<rhoChemistryModel> chemistry = combustion->chem();
rhoReactionThermo& thermo = chemistry->thermo();
thermo.validate(args.executable(), "h", "e");

basicSpecieMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
forAll (Y, i)
{
    Y[i].writeOpt() = IOobject::NO_WRITE;
}

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
combustion->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
    )
);

label Ysize = Y.size();
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;
    Y[spI].writeOpt() = IOobject::AUTO_WRITE;
}
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")));
scalar etamax(lowerN+upperN);
word outletName(SLFMdict.lookup("outletName"));

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

#include "setEtaSpace.H"

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


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

#include "setVarSpace.H"



PtrList<volScalarField> Neta(etamax+1);

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"