for(label j=0 ; j < etamax ; j++)
{
	QhCMC[ j+k*(etamax+1)] = QhCMCNew[ j+k*(etamax+1)]; //  + (1/rhoCMC[j+k*(etamax+1)]) * (Pin - Pin_old);
}