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

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