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;

scalarFieldArray2d Y_SLFM //species mass frac. [-]
(
    NstList.size(),
    scalarFieldArray1d
    (
        Ysize+1,
        scalarField(etamax+1,0)
    )
);
scalarFieldArray2d W_SLFM //reaction rate [1/sec]
(
    NstList.size(),
    scalarFieldArray1d
    (
        Ysize+1,
        scalarField(etamax+1,0)
    )
);
scalarFieldArray1d rho_SLFM //density [g/cm3]
(
    NstList.size(),
    scalarField(etamax+1, 0)
);
scalarFieldArray1d h_SLFM //enthalpy [erg/gm]
(
    NstList.size(),
    scalarField(etamax+1, 0)
);
scalarFieldArray1d Rgas_SLFM //specific gas constant for meanMW
(
    NstList.size(),
    scalarField(etamax+1, 0)
);

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;
    scalar etamax_SLFM, NoSpecies;

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

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

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

    //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+1, 0.0);
    for(label j=0 ; j<=etamax_SLFM ; j++)
    {
        fin>>etaValue_SLFM[j];
    }

    //blank line
    fin.getLine(gbg);

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

    //species and temperature loop
    scalarField Y_temp(etamax_SLFM+1, 0.0);
    scalarField W_temp(etamax_SLFM+1, 0.0);
    for(label i=0 ; i<=NoSpecies ; i++)
    {
        //H2... (species name)
        fin.getLine(gbg);

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

        Y_SLFM[n][i] = interpolateXY(etaValue, etaValue_SLFM, Y_temp);

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

        W_SLFM[n][i] = interpolateXY(etaValue, etaValue_SLFM, W_temp);

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

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

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

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

    //blank line
    fin.getLine(gbg);

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

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

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

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

}