for(label j = 0 ; j<=etamax ;j++)
{
	IFstream ffFile
	(
		"line.dat"
	);

	scalar pz_old(0);
	scalarField zz(4,0), nn(4,0), dd(4,0);
	for(label e = 0 ; e<107 ; e++)
	{
		scalar px, py, pz, ga;
		label ci;

		ffFile>>px>>py>>pz>>ga;
		ci = mesh.findCell(point(px, py, pz));
		scalar drr = pz - pz_old;
        	//Info<<ci<<tab<<px<<tab<<py<<tab<<pz<<endl;
		scalarField ss(4,0);
		for(label k=0 ; k<nf ; k++)
		{
			//sum += (F[k][ci]/mf[ci])*SDRCMC[j+k*(etamax+1)];
			ss[0] += (F[k][ci]/mf[ci])*TCMC[ j+k*(etamax+1) ]; //T
			ss[1] += (F[k][ci]/mf[ci])*QiCMC[(k*(etamax+1)+j)*Ysize + 14]; //CO
                        ss[2] += (F[k][ci]/mf[ci])*QiCMC[(k*(etamax+1)+j)*Ysize + 0]; //H2
                        ss[3] += (F[k][ci]/mf[ci])*QiCMC[(k*(etamax+1)+j)*Ysize + 4]; //OH
                        //ss[0] += (F[k][ci]/mf[ci])*Gij[(k*(etamax+1)+j)*N_Geqn + 5]; //Tvar
                        //ss[1] += (F[k][ci]/mf[ci])*Gij[(k*(etamax+1)+j)*N_Geqn + 2]; //COvar
                        //ss[2] += (F[k][ci]/mf[ci])*Gij[(k*(etamax+1)+j)*N_Geqn + 3]; //OHvar
                        //ss[0] += (F[k][ci]/mf[ci])*Gij[(k*(etamax+1)+j)*N_Geqn + 10]; //T O2 covar
                        //ss[1] += (F[k][ci]/mf[ci])*Gij[(k*(etamax+1)+j)*N_Geqn + 11]; //T CO covar
                        //ss[2] += (F[k][ci]/mf[ci])*Gij[(k*(etamax+1)+j)*N_Geqn + 12]; //T OH covar
		}
		nn[0] += rho[ci] * Peta[j][ci] * ss[0] * pz * drr ;
		dd[0] += rho[ci] * Peta[j][ci] * pz * drr;
                nn[1] += rho[ci] * Peta[j][ci] * ss[1] * pz * drr ;
                dd[1] += rho[ci] * Peta[j][ci] * pz * drr;
                nn[2] += rho[ci] * Peta[j][ci] * ss[2] * pz * drr ;
                dd[2] += rho[ci] * Peta[j][ci] * pz * drr;
                nn[3] += rho[ci] * Peta[j][ci] * ss[3] * pz * drr ;
                dd[3] += rho[ci] * Peta[j][ci] * pz * drr;
		pz_old = pz;	
	}
Info<<etaValue[j]<<tab;
	for(label iii= 0 ; iii<4 ; iii++)
	{
		if(dd[iii] <= 0)
		{
			zz[iii] = 0.0;
		}
		else
		{
			zz[iii] = nn[iii] / dd[iii];
		}
		Info<<zz[iii]<<tab;
	}
Info<<endl;
}
/*
postFile<<runTime.value()<<tab;
scalarField tSDR(10, 0);
scalar ll(0), hh(0);
scalar ff(0), dd(0);
for(label j=0; j <=etamax ; j++)
{
	if(etaValue[j] < 0.351)
	{
		ll = j;
	}
}
hh = ll+1;
if(NewFGindex < 0.5)
{
	for(label k=group ; k>=0 ; k--)
	{
		ff = (0.351 - etaValue[ll])/(etaValue[hh]-etaValue[ll]);
		dd = SDRCMC[hh+k*(etamax+1)] - SDRCMC[ll+k*(etamax+1)];
		tSDR[k] = SDRCMC[ll+k*(etamax+1)] + ff*dd;
		postFile<<tSDR[k]<<tab;
	}
	postFile<<endl;
}
else if(NewFGindex > 0.5 && NewFGindex < 10.5)
{
	for(label k=group ; k>NewFGindex ; k--)
	{
                ff = (0.351 - etaValue[ll])/(etaValue[hh]-etaValue[ll]);
                dd = SDRCMC[hh+(k-NewFGindex)*(etamax+1)] - SDRCMC[ll+(k-NewFGindex)*(etamax+1)];
                tSDR[k] = SDRCMC[ll+(k-NewFGindex)*(etamax+1)] + ff*dd;
	}
	for(label k=NewFGindex ; k >=0 ; k--)
	{
		ff = (0.351 - etaValue[ll])/(etaValue[hh]-etaValue[ll]);
                dd = SDRCMC[hh+0*(etamax+1)] - SDRCMC[ll+0*(etamax+1)];
                tSDR[k] = SDRCMC[ll+0*(etamax+1)] + ff*dd;
	}
        for(label k=group ; k>=0 ; k--)
        {
                postFile<<tSDR[k]<<tab;
        }
        postFile<<endl;
}
else
{
	abort();//nothing
}
*/