scalar TSUM_temp(0);
scalar temp(0);
scalar minimum_mf = 0.005;
RECAL:;
//need to change (use domain integrate function)
forAll(mf , celli)
{
        if(mf[celli] < minimum_mf)
        {
                TSUM_temp += rho[celli] * mesh.V()[celli] * T[celli];
                temp += rho[celli] * mesh.V()[celli];
        }
}

if(TSUM_temp != 0)
{
        TSUM_temp = TSUM_temp / temp;
}
else if(TSUM_temp == 0 || temp == 0)
{
        minimum_mf = minimum_mf * 1.2;
        goto RECAL;
}

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


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

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

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