for(label j = 0; j<=etamax ; j++)
{
        scalarField QiCMCin(Ysize, 0);
        scalar Tin = TCMC[ j+k*(etamax+1) ];
        scalar rhoin = rhoCMC[j+k*(etamax+1)];

        for(label i=0 ; i<Ysize ; i++)
        {
                QiCMCin[i] = QiCMC[(k*(etamax+1)+j)*Ysize + i] ;
	//	if(QiCMCin[i] < 1e-30)
	//	{
	//		QiCMCin[i] = 0.0;
	//	}
        }
	//QiCMCin = 1.0/sum(QiCMCin)*QiCMCin;

        TCMC[ j+k*(etamax+1) ] = chemistry->calculateTCMC(QhCMC[ j+k*(etamax+1)], QiCMCin, Tin, rhoin,Pin);
        rhoCMC[j+k*(etamax+1)] = chemistry->calculateRHOCMC(QiCMCin , TCMC[ j+k*(etamax+1) ], Pin);

//Info<<etaValue[j]<<tab<<TCMC[ j+k*(etamax+1) ]<<endl;
}
