forAll(rho, celli)
        {
                if(F_total[celli] > 0)
                {
			scalarField Ytemp(Ysize, 0.0);
			
                        for(label i=0 ; i < Ysize ; i++) //reaction rate integration over an eta-space for species i
                        {
        			for(label k=0 ; k<=group ; k++)
        			{
					for(label j=0 ; j<= etamax ; j++)
					{
						pdf[j] = Peta[j][celli];
						f[j] = QiCMC[(k*(etamax+1)+j)*Ysize + i];
					}
					Ytemp[i] += Ffrac[k][celli]*integration(deltaftn[celli], MFcut, Neta, pdf, f);
                		}
        		}

                        if( Ytemp[o2index] > Y[o2index][celli] || Ytemp[fuelindex] > Y[fuelindex][celli])
                        {
                                SUMSh[celli] = 0;
                                //keep Yi value
                        }
                        else
                        {
                                scalarField ReactionRate(Ysize, 0.0);
                                scalar tShtemp(0);

                                for(label i = 0 ; i<Ysize ; i++)
                                {
                                        ReactionRate[i] = rho[celli] * (Ytemp[i] - Y[i][celli]) / runTime.deltaT().value();
                                        tShtemp -= chemistry->calculateShCMC( ReactionRate , i);
                                }
                                SUMSh[celli] = tShtemp;

                                for(label i = 0 ; i <Ysize ; i++)
                                {
                                        Y[i][celli] = Ytemp[i];
                                }

                        }
                }
                else
                {
                        SUMSh[celli] = 0;
                }
        }