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] = QiCMCNew[(k*(etamax+1)+j)*Ysize + i] ;
        }

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