//force the saving of the old-time values
//important!!

forAll(Y, yi)
{
	Y[yi].oldTime();
}

forAll(rho, cellI)
{
    scalar& Nst = Neta[lowerN][cellI];

    //find Nst index
    if(Nst < NstList.first())
    {
        nlc[cellI] = 0;
        nfc[cellI] = 0;
    }
    else if(Nst > NstList[NstList.size()-2])
    {
        nhc[cellI] = NstList.size()-2; //just before extinction
        nfc[cellI] = 1;
    }
    else
    {
        nlc[cellI] = label( interpolateXY(Nst, NstList, NstIndex) );
        nhc[cellI] = nlc[cellI]+1;
        nfc[cellI] = (Nst-NstList[nlc[cellI]])
                     /(NstList[nhc[cellI]]-NstList[nlc[cellI]]);
    }

    //update Favre mean temp. and rho(=1/RT)
    scalar jl = jlc[cellI], jh = jhc[cellI];
    scalar vl = vlc[cellI], vh = vhc[cellI];
    scalar nl = nlc[cellI], nh = nhc[cellI];
    scalar jfac = jfc[cellI], vfac = vfc[cellI], nfac = nfc[cellI];


    forAll(Y, yi)
    {
        scalar y00 = Ytable[yi][jl][vl][nl]*(1-jfac)+Ytable[yi][jh][vl][nl]*jfac;
        scalar y01 = Ytable[yi][jl][vl][nh]*(1-jfac)+Ytable[yi][jh][vl][nh]*jfac;
        scalar y10 = Ytable[yi][jl][vh][nl]*(1-jfac)+Ytable[yi][jh][vh][nl]*jfac;
        scalar y11 = Ytable[yi][jl][vh][nh]*(1-jfac)+Ytable[yi][jh][vh][nh]*jfac;

        scalar y0 = y00*(1-vfac)+y10*vfac;
        scalar y1 = y01*(1-vfac)+y11*vfac;

        scalar ygas = y0*(1-nfac)+y1*nfac;

        Y[yi][cellI] = ygas;
    }
}

//correct processor boundary value of T and rho
forAll(Y, yi)
{
    Y[yi].correctBoundaryConditions();
}