PtrList<volScalarField> Fsource(nf); //spray source term fields for an each flame group
for(label k = 0 ; k < nf ; k++)
{
        word Fsourcek = "Fsource_"+Foam::name(k);
        Fsource.set (k,
        new volScalarField
        (
                IOobject
                (
                Fsourcek,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
                ),
        mesh,
        dimensionedScalar("Fsource", dimless, 0.0)
        )
        );
}

        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;

    Fsource[group] = Sevap;

    for (label k=0; k<=group ; k++)
    {
        volScalarField Fk = F[k];
        volScalarField& Fsourcek = Fsource[k];

        solve
        (
                fvm::ddt(rho,Fk)
                //+ fvm::div(phi,Fk) 
		+ mvConvection->fvmDiv(phi, Fk) 
		//+ fvm::div(phi, Fk, "div(phi,Fk)")
                - fvm::laplacian(1.47*turbulence->mut(), Fk) //let 1/Sc = 1.47
                //- fvm::laplacian(1.47*turbulence->muEff(), Fk) //let 1/Sc = 1.47
                ==
                Fsourcek
                ,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];
    }