if(nf > 1)
{
        scalar groupTemp = group; // To check change of flame group
        for(label k = 0; k<nf-1 ; k++)
        {
                if(runTime.value() > InitializationTime[k])
                {
                        group = k+1;
                }
        }
        if(groupTemp != group)
        {
                Need_Re_init = true;
        }
}

Info<<"Now, evaporation source term is added to flame group = "<<group<<endl;

for (label k=0; k<=group ; k++)
{

    volScalarField Fk = F[k];
    if (k==group)
    {
	fvScalarMatrix FkEqn
	(
	        fvm::ddt(rho,Fk)
		+ mvConvection->fvmDiv(phi, Fk) 
	        - fvm::laplacian(1.47*turbulence->mut(), Fk) //let 1/Sc = 1.47
	        ==
	        parcels.Srho(Fk)
	);
        FkEqn.solve(mesh.solver("Fk"));//define matrix solver for F[i] at fvSolution
    }
    else
    {
	fvScalarMatrix FkEqn
	(
	        fvm::ddt(rho,Fk)
		+ mvConvection->fvmDiv(phi, Fk) 
	        - fvm::laplacian(1.47*turbulence->mut(), Fk) //let 1/Sc = 1.47
	       
	);
       	FkEqn.solve(mesh.solver("Fk"));//define matrix solver for F[i] at fvSolution
    }
    F[k] = Fk;
}

F_total = 0.0;
for(label k = 0 ; k < nf ; k++)
{
    F_total += F[k];
}