Switch Equilibrium = SLFMdict.lookup("Equilibrium");
word NstFolder = SLFMdict.lookup("NstFolder");
scalarField NstList = SLFMdict.lookup("NstList");
scalarField NstIndex(NstList.size(),0);

forAll(NstList, n)
{
    NstIndex[n] = n;
}

typedef List<scalarField> scalarFieldArray1d;
typedef List<scalarFieldArray1d> scalarFieldArray2d;
typedef List<scalarFieldArray2d> scalarFieldArray3d;

CMC::FlameStructure slfmLibrary(SLFMdict, etaGrid.size(), Y.size());

scalarFieldArray2d& Y_SLFM (slfmLibrary.Y());
scalarFieldArray2d& W_SLFM (slfmLibrary.W());
scalarFieldArray1d& T_SLFM (slfmLibrary.T());
scalarFieldArray1d& Q_SLFM (slfmLibrary.Q());
scalarFieldArray1d& rho_SLFM (slfmLibrary.rho());
scalarFieldArray1d& h_SLFM (slfmLibrary.h());
scalarFieldArray1d& Rgas_SLFM (slfmLibrary.Rgas());

for(label n=0 ; n<NstList.size() ; n++)
{
    scalar N = NstList[n];

    fileName fname = "sdr"+Foam::name(N)+".inp";
    IFstream fin(runTime.constant()/NstFolder/fname);

    string gbg;
    label etamax_SLFM, NoSpecies;

    //INPUT FILE...
    fin.getLine(gbg);

    //No.GRID...
    fin.getLine(gbg);

    //Number of eta-point and species
    fin>>etamax_SLFM>>NoSpecies;

    //blank line
    fin.getLine(gbg);

    //PRESSURE...
    fin.getLine(gbg);

    //1.00...
    fin.getLine(gbg);

    //MIXTURE...
    fin.getLine(gbg);

    //read eta space (from SLFM library)
    scalarField etaValue_SLFM(etamax_SLFM, 0.0);
    for(label j=0 ; j<etamax_SLFM ; j++)
    {
        fin>>etaValue_SLFM[j];
    }

    //blank line
    fin.getLine(gbg);

    //INITIAL...
    fin.getLine(gbg);

    //species loop
    scalarField Y_temp(etamax_SLFM, 0.0);
    scalarField W_temp(etamax_SLFM, 0.0);
    for(label i=0 ; i<NoSpecies ; i++)
    {
        //H2... (species name)
        string spLine;
        fin.getLine(spLine);
        IStringStream spStream(spLine);
        word spName;
        spStream >> spName;
        const label spI (composition.species()[spName]);
        // Info<<spI << spName << endl;

        //read species mass fraction (from SLFM library)
        for(label j=0 ; j<etamax_SLFM ; j++)
        {
            scalar Yj = 0.0;
            fin >> Yj;
            Y_temp[j] = max(0.0, Yj);
        }

        Y_SLFM[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, Y_temp);

        //read reaction rate (from SLFM library)
        for(label j=0 ; j<etamax_SLFM ; j++)
        {
            fin>>W_temp[j];
        }

        W_SLFM[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, W_temp);

        //blank line
        fin.getLine(gbg);
    }

    //read temperature
    {
        //H2... (species name)
        fin.getLine(gbg);

        //read temperature
        for(label j=0 ; j<etamax_SLFM ; j++)
        {
            scalar Yj = 0.0;
            fin >> Yj;
            Y_temp[j] = max(0.0, Yj);
        }

        T_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, Y_temp);

        //read heat source
        for(label j=0 ; j<etamax_SLFM ; j++)
        {
            fin>>W_temp[j];
        }

        Q_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, W_temp);

        //blank line
        fin.getLine(gbg);
    }

    //DENSITY...
    fin.getLine(gbg);

    //read density (from SLFM library)
    scalarField rho_temp(etamax_SLFM, 0.0);
    for(label j=0 ; j<etamax_SLFM ; j++)
    {
        fin>>rho_temp[j];
    }
    //Info<<rho_temp<<endl;

    rho_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, rho_temp);

    //blank line
    fin.getLine(gbg);

    //ENTHALPY...
    fin.getLine(gbg);

    //read enthalpy (from SLFM library)
    scalarField h_temp(etamax_SLFM, 0.0);
    for(label j=0 ; j<etamax_SLFM ; j++)
    {
        fin>>h_temp[j];
    }

    h_SLFM[n] = interpolateXY(etaGrid, etaValue_SLFM, h_temp);

    for(label i=0 ; i<Y.size() ; i++)
    {
        Rgas_SLFM[n] += Y_SLFM[i][n]/composition.W(i); // = 1/meanMW
    }
    Rgas_SLFM[n] = 8314.46*Rgas_SLFM[n]; //universal gas constant = 8314.46 [J/K.Kmole]

}