diff --git a/cpu/Greyscale.cpp b/cpu/Greyscale.cpp
index affe6c2b..a2395e61 100644
--- a/cpu/Greyscale.cpp
+++ b/cpu/Greyscale.cpp
@@ -1566,6 +1566,1039 @@ extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dis
 }
 
 
+extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_MRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,double *Poros,double *Perm, double *Velocity,double rho0,double *Pressure){
+
+	int n, nread;
+	int nr1,nr2,nr3,nr4,nr5,nr6;
+	int nr7,nr8,nr9,nr10;
+	int nr11,nr12,nr13,nr14;
+	double vx,vy,vz,v_mag;
+    double ux,uy,uz,u_mag;
+    double pressure;//defined for this incompressible model
+	// conserved momemnts
+	double rho,jx,jy,jz;
+	// non-conserved moments
+	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
+    double fq;
+	//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
+    double GeoFun;//geometric function from Guo's PRE 66, 036304 (2002)
+    double porosity;
+    double perm;//voxel permeability
+    double c0, c1; //Guo's model parameters
+    double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
+    double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
+    double rlx_setA = rlx;
+    double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
+
+	const double mrt_V1=0.05263157894736842;
+	const double mrt_V2=0.012531328320802;
+	const double mrt_V3=0.04761904761904762;
+	const double mrt_V4=0.004594820384294068;
+	const double mrt_V5=0.01587301587301587;
+	const double mrt_V6=0.0555555555555555555555555;
+	const double mrt_V7=0.02777777777777778;
+	const double mrt_V8=0.08333333333333333;
+	const double mrt_V9=0.003341687552213868;
+	const double mrt_V10=0.003968253968253968;
+	const double mrt_V11=0.01388888888888889;
+	const double mrt_V12=0.04166666666666666;
+
+	for (int n=start; n<finish; n++){
+        //........................................................................
+        //					READ THE DISTRIBUTIONS
+        //		(read from opposite array due to previous swap operation)
+        //........................................................................
+        // q=0
+		fq = dist[n];
+		rho = fq;
+		m1  = -30.0*fq;
+		m2  = 12.0*fq;
+
+		// q=1
+		//nread = neighborList[n]; // neighbor 2 
+		//fq = dist[nread]; // reading the f1 data into register fq		
+		nr1 = neighborList[n]; 
+		fq = dist[nr1]; // reading the f1 data into register fq
+		rho += fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jx = fq;
+		m4 = -4.0*fq;
+		m9 = 2.0*fq;
+		m10 = -4.0*fq;
+
+		// f2 = dist[10*Np+n];
+		//nread = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
+		//fq = dist[nread];  // reading the f2 data into register fq
+		nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist)
+		fq = dist[nr2];  // reading the f2 data into register fq
+		rho += fq;
+		m1 -= 11.0*(fq);
+		m2 -= 4.0*(fq);
+		jx -= fq;
+		m4 += 4.0*(fq);
+		m9 += 2.0*(fq);
+		m10 -= 4.0*(fq);
+
+		// q=3
+		//nread = neighborList[n+2*Np]; // neighbor 4
+		//fq = dist[nread];
+		nr3 = neighborList[n+2*Np]; // neighbor 4
+		fq = dist[nr3];
+		rho += fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jy = fq;
+		m6 = -4.0*fq;
+		m9 -= fq;
+		m10 += 2.0*fq;
+		m11 = fq;
+		m12 = -2.0*fq;
+
+		// q = 4
+		//nread = neighborList[n+3*Np]; // neighbor 3
+		//fq = dist[nread];
+		nr4 = neighborList[n+3*Np]; // neighbor 3
+		fq = dist[nr4];
+		rho+= fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jy -= fq;
+		m6 += 4.0*fq;
+		m9 -= fq;
+		m10 += 2.0*fq;
+		m11 += fq;
+		m12 -= 2.0*fq;
+
+		// q=5
+		//nread = neighborList[n+4*Np];
+		//fq = dist[nread];
+		nr5 = neighborList[n+4*Np];
+		fq = dist[nr5];
+		rho += fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jz = fq;
+		m8 = -4.0*fq;
+		m9 -= fq;
+		m10 += 2.0*fq;
+		m11 -= fq;
+		m12 += 2.0*fq;
+
+
+		// q = 6
+		//nread = neighborList[n+5*Np];
+		//fq = dist[nread];
+		nr6 = neighborList[n+5*Np];
+		fq = dist[nr6];
+		rho+= fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jz -= fq;
+		m8 += 4.0*fq;
+		m9 -= fq;
+		m10 += 2.0*fq;
+		m11 -= fq;
+		m12 += 2.0*fq;
+
+		// q=7
+		//nread = neighborList[n+6*Np];
+		//fq = dist[nread];
+		nr7 = neighborList[n+6*Np];
+		fq = dist[nr7];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx += fq;
+		m4 += fq;
+		jy += fq;
+		m6 += fq;
+		m9  += fq;
+		m10 += fq;
+		m11 += fq;
+		m12 += fq;
+		m13 = fq;
+		m16 = fq;
+		m17 = -fq;
+
+		// q = 8
+		//nread = neighborList[n+7*Np];
+		//fq = dist[nread];
+		nr8 = neighborList[n+7*Np];
+		fq = dist[nr8];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx -= fq;
+		m4 -= fq;
+		jy -= fq;
+		m6 -= fq;
+		m9 += fq;
+		m10 += fq;
+		m11 += fq;
+		m12 += fq;
+		m13 += fq;
+		m16 -= fq;
+		m17 += fq;
+
+		// q=9
+		//nread = neighborList[n+8*Np];
+		//fq = dist[nread];
+		nr9 = neighborList[n+8*Np];
+		fq = dist[nr9];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx += fq;
+		m4 += fq;
+		jy -= fq;
+		m6 -= fq;
+		m9 += fq;
+		m10 += fq;
+		m11 += fq;
+		m12 += fq;
+		m13 -= fq;
+		m16 += fq;
+		m17 += fq;
+
+		// q = 10
+		//nread = neighborList[n+9*Np];
+		//fq = dist[nread];
+		nr10 = neighborList[n+9*Np];
+		fq = dist[nr10];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx -= fq;
+		m4 -= fq;
+		jy += fq;
+		m6 += fq;
+		m9 += fq;
+		m10 += fq;
+		m11 += fq;
+		m12 += fq;
+		m13 -= fq;
+		m16 -= fq;
+		m17 -= fq;
+
+		// q=11
+		//nread = neighborList[n+10*Np];
+		//fq = dist[nread];
+		nr11 = neighborList[n+10*Np];
+		fq = dist[nr11];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx += fq;
+		m4 += fq;
+		jz += fq;
+		m8 += fq;
+		m9 += fq;
+		m10 += fq;
+		m11 -= fq;
+		m12 -= fq;
+		m15 = fq;
+		m16 -= fq;
+		m18 = fq;
+
+		// q=12
+		//nread = neighborList[n+11*Np];
+		//fq = dist[nread];
+		nr12 = neighborList[n+11*Np];
+		fq = dist[nr12];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx -= fq;
+		m4 -= fq;
+		jz -= fq;
+		m8 -= fq;
+		m9 += fq;
+		m10 += fq;
+		m11 -= fq;
+		m12 -= fq;
+		m15 += fq;
+		m16 += fq;
+		m18 -= fq;
+
+		// q=13
+		//nread = neighborList[n+12*Np];
+		//fq = dist[nread];
+		nr13 = neighborList[n+12*Np];
+		fq = dist[nr13];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx += fq;
+		m4 += fq;
+		jz -= fq;
+		m8 -= fq;
+		m9 += fq;
+		m10 += fq;
+		m11 -= fq;
+		m12 -= fq;
+		m15 -= fq;
+		m16 -= fq;
+		m18 -= fq;
+
+		// q=14
+		//nread = neighborList[n+13*Np];
+		//fq = dist[nread];
+		nr14 = neighborList[n+13*Np];
+		fq = dist[nr14];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx -= fq;
+		m4 -= fq;
+		jz += fq;
+		m8 += fq;
+		m9 += fq;
+		m10 += fq;
+		m11 -= fq;
+		m12 -= fq;
+		m15 -= fq;
+		m16 += fq;
+		m18 += fq;
+
+		// q=15
+		nread = neighborList[n+14*Np];
+		fq = dist[nread];
+		//fq = dist[17*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jy += fq;
+		m6 += fq;
+		jz += fq;
+		m8 += fq;
+		m9 -= 2.0*fq;
+		m10 -= 2.0*fq;
+		m14 = fq;
+		m17 += fq;
+		m18 -= fq;
+
+		// q=16
+		nread = neighborList[n+15*Np];
+		fq = dist[nread];
+		//fq = dist[8*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jy -= fq;
+		m6 -= fq;
+		jz -= fq;
+		m8 -= fq;
+		m9 -= 2.0*fq;
+		m10 -= 2.0*fq;
+		m14 += fq;
+		m17 -= fq;
+		m18 += fq;
+
+		// q=17
+		//fq = dist[18*Np+n];
+		nread = neighborList[n+16*Np];
+		fq = dist[nread];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jy += fq;
+		m6 += fq;
+		jz -= fq;
+		m8 -= fq;
+		m9 -= 2.0*fq;
+		m10 -= 2.0*fq;
+		m14 -= fq;
+		m17 += fq;
+		m18 += fq;
+
+		// q=18
+		nread = neighborList[n+17*Np];
+		fq = dist[nread];
+		//fq = dist[9*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jy -= fq;
+		m6 -= fq;
+		jz += fq;
+		m8 += fq;
+		m9 -= 2.0*fq;
+		m10 -= 2.0*fq;
+		m14 -= fq;
+		m17 -= fq;
+		m18 -= fq;
+        //---------------------------------------------------------------------//
+
+        porosity = Poros[n];
+        perm = Perm[n];
+
+        c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
+        if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
+        GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
+        c1 = porosity*0.5*GeoFun/sqrt(perm);
+        if (porosity==1.0) c1 = 0.0;//i.e. apparent pore nodes
+
+        vx = jx/rho0+0.5*porosity*Gx;
+        vy = jy/rho0+0.5*porosity*Gy;
+        vz = jz/rho0+0.5*porosity*Gz;
+        v_mag=sqrt(vx*vx+vy*vy+vz*vz);
+        ux = vx/(c0+sqrt(c0*c0+c1*v_mag));
+        uy = vy/(c0+sqrt(c0*c0+c1*v_mag));
+        uz = vz/(c0+sqrt(c0*c0+c1*v_mag));
+        u_mag=sqrt(ux*ux+uy*uy+uz*uz);
+
+        //Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
+        Fx = rho0*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
+        Fy = rho0*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
+        Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
+        if (porosity==1.0){
+            Fx=rho0*Gx;
+            Fy=rho0*Gy;
+            Fz=rho0*Gz;
+        }
+
+        //Calculate pressure for MRT model
+        pressure=rho/3.f;
+
+        //-------------------- MRT collison where body force has NO higher-order terms -------------//
+		m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) - m1);
+		m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity) - m2);
+        jx = jx + Fx;
+		m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
+                + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
+        jy = jy + Fy;
+		m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
+                + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
+        jz = jz + Fz;
+		m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
+                + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
+		m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) - m9);
+		m10 = m10 + rlx_setA*( - m10);
+        //m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10);
+		m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) - m11);
+		m12 = m12 + rlx_setA*( - m12);
+        //m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12);
+		m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) - m13);
+		m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) - m14);
+		m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) - m15);
+		m16 = m16 + rlx_setB*( - m16);
+		m17 = m17 + rlx_setB*( - m17);
+		m18 = m18 + rlx_setB*( - m18);
+        //.......................................................................................................
+
+
+        //.................inverse transformation......................................................
+        // q=0
+		fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2;
+		dist[n] = fq;
+
+		// q = 1
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10);
+		//nread = neighborList[n+Np];
+		dist[nr2] = fq;
+
+		// q=2
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10);
+		//nread = neighborList[n];
+		dist[nr1] = fq;
+
+		// q = 3
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12);
+		//nread = neighborList[n+3*Np];
+		dist[nr4] = fq;
+
+		// q = 4
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12);
+		//nread = neighborList[n+2*Np];
+		dist[nr3] = fq;
+
+		// q = 5
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11);
+		//nread = neighborList[n+5*Np];
+		dist[nr6] = fq;
+
+		// q = 6
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11);
+		//nread = neighborList[n+4*Np];
+		dist[nr5] = fq;
+
+		// q = 7
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+
+				mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17);
+		//nread = neighborList[n+7*Np];
+		dist[nr8] = fq;
+
+		// q = 8
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+				+mrt_V12*m12+0.25*m13+0.125*(m17-m16);
+		//nread = neighborList[n+6*Np];
+		dist[nr7] = fq;
+
+		// q = 9
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+
+				mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17);
+		//nread = neighborList[n+9*Np];
+		dist[nr10] = fq;
+
+		// q = 10
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+
+				mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17);
+		//nread = neighborList[n+8*Np];
+		dist[nr9] = fq;
+
+		// q = 11
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)
+				+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
+				-mrt_V12*m12+0.25*m15+0.125*(m18-m16);
+		//nread = neighborList[n+11*Np];
+		dist[nr12] = fq;
+
+		// q = 12
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+
+				mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18);
+		//nread = neighborList[n+10*Np];
+		dist[nr11]= fq;
+
+		// q = 13
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)
+				+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
+				-mrt_V12*m12-0.25*m15-0.125*(m16+m18);
+		//nread = neighborList[n+13*Np];
+		dist[nr14] = fq;
+
+		// q= 14
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)
+				+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
+				-mrt_V12*m12-0.25*m15+0.125*(m16+m18);
+		//nread = neighborList[n+12*Np];
+		dist[nr13] = fq;
+
+
+		// q = 15
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)
+				-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18);
+		nread = neighborList[n+15*Np];
+		dist[nread] = fq;
+
+		// q = 16
+		fq =  mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)
+				-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17);
+		nread = neighborList[n+14*Np];
+		dist[nread] = fq;
+
+
+		// q = 17
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)
+				-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18);
+		nread = neighborList[n+17*Np];
+		dist[nread] = fq;
+
+		// q = 18
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)
+				-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18);
+		nread = neighborList[n+16*Np];
+		dist[nread] = fq;
+        //........................................................................
+
+        //Update velocity on device
+        Velocity[0*Np+n] = ux;
+        Velocity[1*Np+n] = uy;
+        Velocity[2*Np+n] = uz;
+        //Update pressure on device
+        Pressure[n] = pressure;
+	}
+}
+
+extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_MRT(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Gx, double Gy, double Gz,double *Poros,double *Perm, double *Velocity,double rho0,double *Pressure){
+
+	int n;
+	double vx,vy,vz,v_mag;
+    double ux,uy,uz,u_mag;
+    double pressure;//defined for this incompressible model
+	// conserved momemnts
+	double rho,jx,jy,jz;
+	// non-conserved moments
+	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
+    double fq;
+	//double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18;
+    double GeoFun;//geometric function from Guo's PRE 66, 036304 (2002)
+    double porosity;
+    double perm;//voxel permeability
+    double c0, c1; //Guo's model parameters
+    double mu_eff = (1.0/rlx_eff-0.5)/3.0;//kinematic viscosity
+    double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz)
+    double rlx_setA = rlx;
+    double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
+
+	const double mrt_V1=0.05263157894736842;
+	const double mrt_V2=0.012531328320802;
+	const double mrt_V3=0.04761904761904762;
+	const double mrt_V4=0.004594820384294068;
+	const double mrt_V5=0.01587301587301587;
+	const double mrt_V6=0.0555555555555555555555555;
+	const double mrt_V7=0.02777777777777778;
+	const double mrt_V8=0.08333333333333333;
+	const double mrt_V9=0.003341687552213868;
+	const double mrt_V10=0.003968253968253968;
+	const double mrt_V11=0.01388888888888889;
+	const double mrt_V12=0.04166666666666666;
+
+
+	for (int n=start; n<finish; n++){
+        //........................................................................
+        //					READ THE DISTRIBUTIONS
+        //		(read from opposite array due to previous swap operation)
+        //........................................................................
+        // q=0
+		fq = dist[n];
+		rho = fq;
+		m1  = -30.0*fq;
+		m2  = 12.0*fq;
+
+		// q=1
+		fq = dist[2*Np+n];
+		rho += fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jx = fq;
+		m4 = -4.0*fq;
+		m9 = 2.0*fq;
+		m10 = -4.0*fq;
+
+		// f2 = dist[10*Np+n];
+		fq = dist[1*Np+n];
+		rho += fq;
+		m1 -= 11.0*(fq);
+		m2 -= 4.0*(fq);
+		jx -= fq;
+		m4 += 4.0*(fq);
+		m9 += 2.0*(fq);
+		m10 -= 4.0*(fq);
+
+		// q=3
+		fq = dist[4*Np+n];
+		rho += fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jy = fq;
+		m6 = -4.0*fq;
+		m9 -= fq;
+		m10 += 2.0*fq;
+		m11 = fq;
+		m12 = -2.0*fq;
+
+		// q = 4
+		fq = dist[3*Np+n];
+		rho+= fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jy -= fq;
+		m6 += 4.0*fq;
+		m9 -= fq;
+		m10 += 2.0*fq;
+		m11 += fq;
+		m12 -= 2.0*fq;
+
+		// q=5
+		fq = dist[6*Np+n];
+		rho += fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jz = fq;
+		m8 = -4.0*fq;
+		m9 -= fq;
+		m10 += 2.0*fq;
+		m11 -= fq;
+		m12 += 2.0*fq;
+
+		// q = 6
+		fq = dist[5*Np+n];
+		rho+= fq;
+		m1 -= 11.0*fq;
+		m2 -= 4.0*fq;
+		jz -= fq;
+		m8 += 4.0*fq;
+		m9 -= fq;
+		m10 += 2.0*fq;
+		m11 -= fq;
+		m12 += 2.0*fq;
+
+		// q=7
+		fq = dist[8*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx += fq;
+		m4 += fq;
+		jy += fq;
+		m6 += fq;
+		m9  += fq;
+		m10 += fq;
+		m11 += fq;
+		m12 += fq;
+		m13 = fq;
+		m16 = fq;
+		m17 = -fq;
+
+		// q = 8
+		fq = dist[7*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx -= fq;
+		m4 -= fq;
+		jy -= fq;
+		m6 -= fq;
+		m9 += fq;
+		m10 += fq;
+		m11 += fq;
+		m12 += fq;
+		m13 += fq;
+		m16 -= fq;
+		m17 += fq;
+
+		// q=9
+		fq = dist[10*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx += fq;
+		m4 += fq;
+		jy -= fq;
+		m6 -= fq;
+		m9 += fq;
+		m10 += fq;
+		m11 += fq;
+		m12 += fq;
+		m13 -= fq;
+		m16 += fq;
+		m17 += fq;
+
+		// q = 10
+		fq = dist[9*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx -= fq;
+		m4 -= fq;
+		jy += fq;
+		m6 += fq;
+		m9 += fq;
+		m10 += fq;
+		m11 += fq;
+		m12 += fq;
+		m13 -= fq;
+		m16 -= fq;
+		m17 -= fq;
+
+		// q=11
+		fq = dist[12*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx += fq;
+		m4 += fq;
+		jz += fq;
+		m8 += fq;
+		m9 += fq;
+		m10 += fq;
+		m11 -= fq;
+		m12 -= fq;
+		m15 = fq;
+		m16 -= fq;
+		m18 = fq;
+
+		// q=12
+		fq = dist[11*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx -= fq;
+		m4 -= fq;
+		jz -= fq;
+		m8 -= fq;
+		m9 += fq;
+		m10 += fq;
+		m11 -= fq;
+		m12 -= fq;
+		m15 += fq;
+		m16 += fq;
+		m18 -= fq;
+
+		// q=13
+		fq = dist[14*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx += fq;
+		m4 += fq;
+		jz -= fq;
+		m8 -= fq;
+		m9 += fq;
+		m10 += fq;
+		m11 -= fq;
+		m12 -= fq;
+		m15 -= fq;
+		m16 -= fq;
+		m18 -= fq;
+
+		// q=14
+		fq = dist[13*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jx -= fq;
+		m4 -= fq;
+		jz += fq;
+		m8 += fq;
+		m9 += fq;
+		m10 += fq;
+		m11 -= fq;
+		m12 -= fq;
+		m15 -= fq;
+		m16 += fq;
+		m18 += fq;
+
+		// q=15
+		fq = dist[16*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jy += fq;
+		m6 += fq;
+		jz += fq;
+		m8 += fq;
+		m9 -= 2.0*fq;
+		m10 -= 2.0*fq;
+		m14 = fq;
+		m17 += fq;
+		m18 -= fq;
+
+		// q=16
+		fq = dist[15*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jy -= fq;
+		m6 -= fq;
+		jz -= fq;
+		m8 -= fq;
+		m9 -= 2.0*fq;
+		m10 -= 2.0*fq;
+		m14 += fq;
+		m17 -= fq;
+		m18 += fq;
+
+		// q=17
+		fq = dist[18*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jy += fq;
+		m6 += fq;
+		jz -= fq;
+		m8 -= fq;
+		m9 -= 2.0*fq;
+		m10 -= 2.0*fq;
+		m14 -= fq;
+		m17 += fq;
+		m18 += fq;
+
+		// q=18
+		fq = dist[17*Np+n];
+		rho += fq;
+		m1 += 8.0*fq;
+		m2 += fq;
+		jy -= fq;
+		m6 -= fq;
+		jz += fq;
+		m8 += fq;
+		m9 -= 2.0*fq;
+		m10 -= 2.0*fq;
+		m14 -= fq;
+		m17 -= fq;
+		m18 -= fq;
+        //---------------------------------------------------------------------//
+
+        porosity = Poros[n];
+        perm = Perm[n];
+
+        c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm);
+        if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes
+        GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity);
+        c1 = porosity*0.5*GeoFun/sqrt(perm);
+        if (porosity==1.0) c1 = 0.0;//i.e. apparent pore nodes
+
+        vx = jx/rho0+0.5*porosity*Gx;
+        vy = jy/rho0+0.5*porosity*Gy;
+        vz = jz/rho0+0.5*porosity*Gz;
+        v_mag=sqrt(vx*vx+vy*vy+vz*vz);
+        ux = vx/(c0+sqrt(c0*c0+c1*v_mag));
+        uy = vy/(c0+sqrt(c0*c0+c1*v_mag));
+        uz = vz/(c0+sqrt(c0*c0+c1*v_mag));
+        u_mag=sqrt(ux*ux+uy*uy+uz*uz);
+
+        //Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
+        Fx = rho0*(-porosity*mu_eff/perm*ux - porosity*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx);
+        Fy = rho0*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy);
+        Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz);
+        if (porosity==1.0){
+            Fx=rho0*Gx;
+            Fy=rho0*Gy;
+            Fz=rho0*Gz;
+        }
+
+        //Calculate pressure for Incompressible-MRT model
+        pressure=rho/3.f;
+
+        //-------------------- IMRT collison where body force has NO higher-order terms -------------//
+		m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) - m1);
+		m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity) - m2);
+        jx = jx + Fx;
+		m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
+                + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
+        jy = jy + Fy;
+		m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
+                + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
+        jz = jz + Fz;
+		m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
+                + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
+		m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) - m9);
+		m10 = m10 + rlx_setA*( - m10);
+        //m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10);
+		m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) - m11);
+		m12 = m12 + rlx_setA*( - m12);
+        //m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12);
+		m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) - m13);
+		m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) - m14);
+		m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) - m15);
+		m16 = m16 + rlx_setB*( - m16);
+		m17 = m17 + rlx_setB*( - m17);
+		m18 = m18 + rlx_setB*( - m18);
+        //.......................................................................................................
+
+        //.................inverse transformation......................................................
+        // q=0
+		fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2;
+		dist[n] = fq;
+
+		// q = 1
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10);
+		dist[1*Np+n] = fq;
+
+		// q=2
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10);
+		dist[2*Np+n] = fq;
+
+		// q = 3
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12);
+		dist[3*Np+n] = fq;
+
+		// q = 4
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12);
+		dist[4*Np+n] = fq;
+
+		// q = 5
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11);
+		dist[5*Np+n] = fq;
+
+		// q = 6
+		fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11);
+		dist[6*Np+n] = fq;
+
+		// q = 7
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+
+				mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17);
+		dist[7*Np+n] = fq;
+
+
+		// q = 8
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11
+				+mrt_V12*m12+0.25*m13+0.125*(m17-m16);
+		dist[8*Np+n] = fq;
+
+		// q = 9
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+
+				mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17);
+		dist[9*Np+n] = fq;
+
+		// q = 10
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+
+				mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17);
+		dist[10*Np+n] = fq;
+
+
+		// q = 11
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)
+				+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
+				-mrt_V12*m12+0.25*m15+0.125*(m18-m16);
+		dist[11*Np+n] = fq;
+
+		// q = 12
+		fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+
+				mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18);
+		dist[12*Np+n] = fq;
+
+		// q = 13
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)
+				+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
+				-mrt_V12*m12-0.25*m15-0.125*(m16+m18);
+		dist[13*Np+n] = fq;
+
+		// q= 14
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)
+				+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11
+				-mrt_V12*m12-0.25*m15+0.125*(m16+m18);
+
+		dist[14*Np+n] = fq;
+
+		// q = 15
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)
+				-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18);
+		dist[15*Np+n] = fq;
+
+		// q = 16
+		fq =  mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)
+				-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17);
+		dist[16*Np+n] = fq;
+
+
+		// q = 17
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)
+				-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18);
+		dist[17*Np+n] = fq;
+
+		// q = 18
+		fq = mrt_V1*rho+mrt_V9*m1
+				+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)
+				-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18);
+		dist[18*Np+n] = fq;
+        //........................................................................
+
+        //Update velocity on device
+        Velocity[0*Np+n] = ux;
+        Velocity[1*Np+n] = uy;
+        Velocity[2*Np+n] = uz;
+        //Update pressure on device
+        Pressure[n] = pressure;
+	}
+}
+
+
 extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den)
 {
 	int n;