scalar TZONEmax = T.weightedAverage(mesh.V()).value();

for(label j=0 ; j<=etamax ; j++) // initialize conditional temperature
{
        for(label k = 0 ; k<nf ; k++)
        {
        	TCMC[j+k*(etamax+1) ] = -(TZONEmax - parcels.constProps().T0())*(etaValue[j] - 1)+parcels.constProps().T0();
        }
}

for(label j=0 ; j <= etamax ; j++) //get conditional enthalpy 
{
        for(label k = 0 ; k<nf ; k++)
        {
                scalar Tin = TCMC[ j+k*(etamax+1) ];
                scalarField QiCMCin(Ysize, 0);

                for(label i=0 ; i<Ysize ; i++)
                {
                        QiCMCin[i] = QiCMC[(k*(etamax+1)+j)*Ysize + i];
                }

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