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

rho.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];

    scalar T00 = Ttable[jl][vl][nl]*(1-jfac)+Ttable[jh][vl][nl]*jfac;
    scalar T01 = Ttable[jl][vl][nh]*(1-jfac)+Ttable[jh][vl][nh]*jfac;
    scalar T10 = Ttable[jl][vh][nl]*(1-jfac)+Ttable[jh][vh][nl]*jfac;
    scalar T11 = Ttable[jl][vh][nh]*(1-jfac)+Ttable[jh][vh][nh]*jfac;

    scalar T0 = T00*(1-vfac)+T10*vfac;
    scalar T1 = T01*(1-vfac)+T11*vfac;

    T[cellI] = T0*(1-nfac)+T1*nfac;

    scalar R00 = Rtable[jl][vl][nl]*(1-jfac)+Rtable[jh][vl][nl]*jfac;
    scalar R01 = Rtable[jl][vl][nh]*(1-jfac)+Rtable[jh][vl][nh]*jfac;
    scalar R10 = Rtable[jl][vh][nl]*(1-jfac)+Rtable[jh][vh][nl]*jfac;
    scalar R11 = Rtable[jl][vh][nh]*(1-jfac)+Rtable[jh][vh][nh]*jfac;

    scalar R0 = R00*(1-vfac)+R10*vfac;
    scalar R1 = R01*(1-vfac)+R11*vfac;

    scalar Rgas = R0*(1-nfac)+R1*nfac;
    rho[cellI] = p[cellI]/(Rgas*T[cellI]);
}

//correct processor boundary value of T and rho
T.correctBoundaryConditions();
rho.correctBoundaryConditions();