tmp<fv::convectionScheme<scalar>> mvConvection
(
    fv::convectionScheme<scalar>::New
    (
        mesh,
        fields,
        phi,
        mesh.divScheme("div(phi,ft_wv_n2_b_ha_hau)")
		//mesh.divScheme("div(phi,ft_b_ha_hau)")
    )
);

if (composition.contains("ft"))
{
    volScalarField& ft = composition.Y("ft");
    
    
    fvScalarMatrix ftEqn
    (
        fvm::ddt(rho, ft)
      + mvConvection->fvmDiv(phi, ft)
      - fvm::laplacian(turbulence->alphaEff(), ft)
              
	==
        fvOptions(rho, ft) 
    );
    
    //ftEqn.relax();

    fvOptions.constrain(ftEqn);

    ftEqn.solve();

    fvOptions.correct(ft);
	
	
	//wv(water vapor) 
	//volScalarField& pr = composition.Y("pr");
	
    
	volScalarField& wv = composition.Y("wv");
  	
	fvScalarMatrix wvEqn
    (
        fvm::ddt(rho, wv)
      + mvConvection->fvmDiv(phi, wv)
      - fvm::laplacian(turbulence->alphaEff(), wv)
     ==
        fvOptions(rho, wv) 
    );

    //wvEqn.relax();

    
    fvOptions.constrain(wvEqn);

    wvEqn.solve();
	
	//wv += pr;

    fvOptions.correct(wv);
	
	
    //N2 

	
	volScalarField& n2 = composition.Y("n2");

	fvScalarMatrix n2Eqn
    (
        fvm::ddt(rho, n2)
      + mvConvection->fvmDiv(phi, n2)
      - fvm::laplacian(turbulence->alphaEff(), n2)
     ==
        fvOptions(rho, n2) 
    );

    //n2Eqn.relax();

    fvOptions.constrain(n2Eqn);

    n2Eqn.solve();

    fvOptions.correct(n2); 


	
}

