const int Ysize = Y.size();

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

//mixture fraction variance
Info<<"Reading field mfVar\n"<<endl;
volScalarField mfVar
(
    IOobject
    (
        "mfVar",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

//Scalar Dissipation Rate
volScalarField SDR
(
    IOobject
    (
        "SDR",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("SDR", dimEnergy/dimTime/dimEnergy, 0.0)
);

//delta function index
volScalarField deltaftn
(
    IOobject
    (
        "deltaftn",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("deltaftn", dimless, -1.0)
);

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

Switch uni_eta( CMCdict.lookup("uni_eta", false) );
Switch init_start_CMC( CMCdict.lookup("init_start_CMC", true) );

//Switch FG_interaction( CMCdict.lookup("FG_interaction", false));

//scalar Tign(readScalar(CMCdict.lookup("Tign")));
//scalar Kebu(readScalar(CMCdict.lookup("Kebu")));

//scalar nf_p2(readScalar(CMCdict.lookup("nf_p2")));
//scalar nf_p1(readScalar(CMCdict.lookup("nf_p1")));
//scalar nf_main(readScalar(CMCdict.lookup("nf_main")));
//scalar nf = nf_p2 + nf_p1 + nf_main;


scalar nf_p2(readScalar(CMCdict.lookup("nf_p2")));
scalar nf_p1(readScalar(CMCdict.lookup("nf_p1")));
scalar nf_main(readScalar(CMCdict.lookup("nf_main")));
scalar nf = nf_p2 + nf_p1 + nf_main;

scalar etamax( readScalar(CMCdict.lookup("etamax")) );
scalarField etaValue(etamax+1, 0.0);
scalarField deta(etamax+1 ,0.0 );

scalar section(readScalar(CMCdict.lookup("section")));
scalarField MFcut(section+1, 0.0);
scalarField Neta(section+1, 0.0);

scalarField pdf(etamax+1,0.0);
scalarField f(etamax+1, 0.0);
Info<<"Set eta space grid"<<tab<<uni_eta<<endl;

//Uniform eta space grid
if(uni_eta == true)
{
        for(label etacount = 0 ; etacount <= etamax ; etacount++)
        {
		deta[etacount] = 1/etamax;
		if(etacount == 0)
		{
                	etaValue[0] = 0;
		}
		else
		{
                	etaValue[etacount] = etaValue[etacount-1] + deta[etacount];
		}
        }
        //for(label etacount=0 ; etacount <= etamax ; etacount++)
        //{
        //      Info<<etacount<<tab<<etaValue[etacount]<<tab<<deta[etacount]<<endl;
        //}
}

else if(uni_eta == false)
{
        dictionary non_uni_eta_dict = CMCdict.subDict("non_uni_eta_dict");

	MFcut = non_uni_eta_dict.lookup("MFcut");
	Neta = non_uni_eta_dict.lookup("Neta");
		
	int i = 0;
        for(label etacount=1 ; etacount <= etamax ; etacount++)
        {
		if(etacount > Neta[i+1])
		{
			i = i+1;
		}
		if(Neta[i] < etacount && etacount <= Neta[i+1])
		{
			etaValue[etacount] = etaValue[etacount-1] + (MFcut[i+1]-MFcut[i])/(Neta[i+1] - Neta[i]);
		}
        }

        for(label etacount=0 ; etacount <= etamax ; etacount++)
        {
                if(etacount == 0 || etacount == etamax)
                {
                        deta[0] = etaValue[1] - etaValue[0];
                        deta[etamax] = etaValue[etamax] - etaValue[etamax-1];
                }
                else
                {
                        deta[etacount] = ( etaValue[etacount+1] - etaValue[etacount-1] ) / 2;
                }
        }
        //for(label etacount=0 ; etacount <= etamax ; etacount++)
        //{
        //      Info<<etaValue[etacount]<<endl;
        //}
}

scalarField ChemDeltaT( (etamax+1)*nf, 0.0000001); //initial chemical time step;
scalar group(0); //flame group number to add evaporation source (default = 0), Flamegroup.H line# 29
Switch Need_Re_init = false; //To check reinitialization of flame group (default = false), Flamegroup.H line# 35
Switch CMC_on = false; //to check CMC calculation, before injection this shold be false

//Peta field
PtrList<volScalarField> Peta(etamax+1); //Parameter List, number of eta space
for(label etacount = 0 ; etacount <= etamax ; etacount++)
{
        word Petai = "Peta_"+Foam::name(etacount);
        Peta.set (etacount,
        new volScalarField
        (
                IOobject
                (
                Petai,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
                ),
        mesh,
        dimensionedScalar("Peta", dimless, 0.0)
        )
        );
}

wordList FGtypes
(
	mf.boundaryField().size(),
	zeroGradientFvPatchScalarField::typeName
);

forAll(mf.boundaryField(), patchi)
{
	FGtypes[patchi] = mf.boundaryField()[patchi].type();
}

PtrList<volScalarField> F(nf); //flame group field
for(label k = 0 ; k < nf ; k++)
{
	word Fk = "F_"+Foam::name(k);
        F.set (k,
                new volScalarField
                (
                        IOobject
                        (
                        Fk,
                        runTime.timeName(),
                        mesh,
                        IOobject::READ_IF_PRESENT,
                        IOobject::AUTO_WRITE
                        ),
                mesh,
                dimensionedScalar("F", dimless, 0.0),
		FGtypes
		)
	);
}

volScalarField F_total //summation of F[i] for each cell
(
    IOobject
    (
        "F_total",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("F_total", dimless, 0.0)
);

PtrList<volScalarField> Ffrac(nf);
for(label k = 0 ; k < nf ; k++)
{
        word Ffrack = "Ffrac_"+Foam::name(k);
        Ffrac.set (k,
        new volScalarField
        (
                IOobject
                (
                Ffrack,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
                ),
        mesh,
        dimensionedScalar("Ffrac", dimless, 0.0)
        )
        );
}

//conditional SDR
PtrList<volScalarField> CSDRCELL( etamax+1 ); //Parameter List, number of eta space
for(label j = 0 ; j <= etamax ; j++)
{
        word CSDRCELLi = "CSDRCELL_"+Foam::name(j);
        CSDRCELL.set (j,
        new volScalarField
        (
                IOobject
                (
                CSDRCELLi,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
                ),
        mesh,
        dimensionedScalar("CSDRCELL", dimless, 0.0)
        )
        );
}

PtrList<volScalarField> SUMRR(Ysize); //Parameter List, number of eta space
for(label i = 0 ; i < Ysize ; i++)
{
        word SUMRRi = "SUMRR_"+Foam::name(i);
        SUMRR.set (i,
        new volScalarField
        (
                IOobject
                (
                SUMRRi,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
                ),
        mesh,
        dimensionedScalar("SUMRR", dimless, 0.0)
        )
        );
}

//chemistry source term
volScalarField SUMSh
(
    IOobject
    (
        "SUMSh",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("SUMSh", dimEnergy/dimTime/dimVolume, 0.0)
);

// Add eta-space for species Yi using vectorized array method (2D or 3D array -> 1D array)
// Qi(i , j) = Qi[i+j*(imax+1)]
// Qi(i , j , k) = Qi[ (k*(jmax+1)+j)*(imax+1) + i ]
//when i = 0~imax, j=0~jmax, k=0~kmax
// i = number of species
// j = number of eta-space
// k = number of flame group

/////////////////Qi and Qh ZONE value define//////////////////////
scalarField QiCMC( Ysize*(etamax+1)*nf , 0 );
scalarField QhCMC( (etamax+1)*nf , 0 );
scalarField QiCMCNew( Ysize*(etamax+1)*nf , 0 );
scalarField QhCMCNew( (etamax+1)*nf , 0 );
scalarField SDRCMC((etamax+1)*nf,0);

/////////////////////Conditional T field////////////////////////////
scalarField TCMC( (etamax+1)*nf, 0);

/////////////////////Conditional density field////////////////////////////
scalarField rhoCMC( (etamax+1)*nf, 0);

/////////////////////CMC equation source term/////////////////////////
scalarField CMCRR( Ysize*(etamax+1)*nf , 0 );
scalarField CMCSh( (etamax+1)*nf , 0 );
scalar HRR(0);
scalar Pin(p.weightedAverage(mesh.V()).value());
scalar Pin_old(Pin);

scalar o2index(0);
scalar fuelindex(0);

for(label i=0;i<Ysize;i++)
{
        if(Y[i].name()=="O2")
        {
              o2index=i;
        }
        if(Y[i].name()=="C7H16")
        {
              fuelindex=i;
        }
}

scalar outi(0);

OFstream CMCresultsFile
(
	("CMCresults."+ runTime.timeName() +".dat")
);

//make header of Tecplot format
CMCresultsFile<<"TITLE = "<<"\""<<"Results of CMC model calculation (Qi, Qh, T, rho, SDR)"<<"\""<<endl;
CMCresultsFile<<"VARIABLES = "<<"\"eta\"";
CMCresultsFile<<", ";
for(label i = 0 ; i<Ysize ; i++)
{
	CMCresultsFile<<"\""<<Y[i].name()<<"\", ";
}
CMCresultsFile<<"\"Qh\", \"T\", \"RHO\", \"SDR\"";
CMCresultsFile<<endl;
//make header of Tecplot format

volScalarField mf_temp
(
    IOobject
    (
        "mf_temp",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("zero",dimless,0.0)
);



if(init_start_CMC == false)
{
        #include "CMCwrite.H"
}

scalar NewFGindex(0);
scalar NewFGcount(1);