!
!  LBMLID3D.f95
!
!  FUNCTIONS:
!	LBMLID3D        - Entry point of console application.
!	LEEquilibrium	- Linearized entropic equilibrium
!	LQEquilibrium	- Linearized quasi equilibrium
!	Moments			- Conserved moments
!	Trace			- Trace of the pressure tensor
!
!****************************************************************************
!
!  PROGRAM: LBMLID3D, by Pietro ASINARI, 2011
!
!  PURPOSE: Simple LBM code for solving 3D lid driven cavity by E-QE-MRT
! 
! Copyright 2011 Pietro Asinari (Politecnico di Torino, IT)
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
! 
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!****************************************************************************

program lbmlid3d

	implicit none

! Internal Functions
!********************

! Local variables
!********************
	real(SELECTED_REAL_KIND(15,307)),dimension(:,:,:,:),allocatable	:: F,Fp
	real(SELECTED_REAL_KIND(15,307)),dimension(:,:,:),allocatable	:: FWp
	real(SELECTED_REAL_KIND(15,307)),dimension(:,:,:),allocatable	:: rho,ux,uy,uz,P,u,v,w
	real(SELECTED_REAL_KIND(15,307)),dimension(:),allocatable	    :: Estore
	
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19,1:19) :: A,RA
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19) :: Feq,FeqF,FeqW,Fqe,Fge,G,Floc,s,Fdiff,Fnew
	integer,dimension(1:19,1:3)	:: VL
	integer,dimension(1:19)		:: BB

	real(SELECTED_REAL_KIND(15,307)) :: rhoW,uxW,uyW,uzW,uxB,uyB,uzB
	real(SELECTED_REAL_KIND(15,307)) :: E,x,y,z,temp,tempx,tempy,tempz
	real(SELECTED_REAL_KIND(15,307)) :: Sxx,Syy,Szz,Tr,Treq,omega1,omega2
	integer :: i,j,k,t,alpha
	integer :: IO_Stat_Number1 = -1

! Parameters
!********************
	real(SELECTED_REAL_KIND(15,307)),parameter :: Tt =  30.0d0
	real(SELECTED_REAL_KIND(15,307)),parameter :: Lx =  1.0d0
	real(SELECTED_REAL_KIND(15,307)),parameter :: UL =  1.0d0
	real(SELECTED_REAL_KIND(15,307)),parameter :: P0 =  0.0d0
! kinematic viscosity
	real(SELECTED_REAL_KIND(15,307)),parameter :: NU =  0.005d0
! bulk viscosity
	real(SELECTED_REAL_KIND(15,307)),parameter :: XI =  0.050d0
! space discretization
	integer,parameter :: nx = 25
	integer,parameter :: ny = nx
	integer,parameter :: nz = nx

	real(SELECTED_REAL_KIND(15,307)),parameter :: Re   = UL*Lx/NU;
	real(SELECTED_REAL_KIND(15,307)),parameter :: gx   = 0.0d0
	real(SELECTED_REAL_KIND(15,307)),parameter :: gy   = 0.0d0
	real(SELECTED_REAL_KIND(15,307)),parameter :: gz   = 0.0d0
	real(SELECTED_REAL_KIND(15,307)),parameter :: rho0 = 1.0d0

	real(SELECTED_REAL_KIND(15,307)),parameter :: dx = Lx/nx
	real(SELECTED_REAL_KIND(15,307)),parameter :: Kn = dx/Lx
!
! (tau-1/2+theta)/3 = Ma / (Re Kn)
! Diffusive scaling Ma = k Kn
! tau = 1/2+3k/Re-theta
! Stability threshold
! 3k/Re-theta > 0
! not too small
!
!	real(SELECTED_REAL_KIND(15,307)),parameter :: Ma  = 0.01d0*Re*Kn
!	real(SELECTED_REAL_KIND(15,307)),parameter :: Ma  = Kn
	real(SELECTED_REAL_KIND(15,307)),parameter :: Ma  = 0.1d0

	real(SELECTED_REAL_KIND(15,307)),parameter :: c   = UL/Ma
	real(SELECTED_REAL_KIND(15,307)),parameter :: uxL = UL/c
	real(SELECTED_REAL_KIND(15,307)),parameter :: uyL = 0.0d0/c
	real(SELECTED_REAL_KIND(15,307)),parameter :: uzL = 0.0d0/c

	real(SELECTED_REAL_KIND(15,307)),parameter :: dt   = dx/c
	real(SELECTED_REAL_KIND(15,307)),parameter :: lbn  = NU/(dt*c**2)
	real(SELECTED_REAL_KIND(15,307)),parameter :: lbx  = XI/(dt*c**2)
	real(SELECTED_REAL_KIND(15,307)),parameter :: tau1 = 3.d0*lbn+1.0d0/2.0d0
	real(SELECTED_REAL_KIND(15,307)),parameter :: tau2 = 3.d0*lbx+1.0d0/2.0d0
	integer :: nt

! Constants
!********************
	real(SELECTED_REAL_KIND(15,307)),parameter :: pi = 3.141592653589793238d0 
	integer,parameter :: Q = 19

    nt   = int(Tt/dt);

! Lattice D2Q9
	VL(:,1) = (/0,1,0,-1,0,1,-1,-1,1,1,0,-1,0,0,1,0,-1,0,0/) ! x-axis [-]
	VL(:,2) = (/0,0,1,0,-1,1,1,-1,-1,0,1,0,-1,0,0,1,0,-1,0/) ! y-axis [-]
	VL(:,3) = (/0,0,0,0,0,0,0,0,0,1,1,1,1,1,-1,-1,-1,-1,-1/) ! z-axis [-]
	BB(:)   = (/1,4,5,2,3,8,9,6,7,17,18,15,16,19,12,13,10,11,14/) ! bounce back

! Memory Allocation
	allocate(   F(1:nx,1:ny,1:nz,1:19) )
	allocate(  Fp(1:nx,1:ny,1:nz,1:19) )
	allocate( FWp(1:nx,1:nz,1:19) )
	allocate( rho(1:nx,1:ny,1:nz) )
	allocate(  ux(1:nx,1:ny,1:nz) )
	allocate(  uy(1:nx,1:ny,1:nz) )
	allocate(  uz(1:nx,1:ny,1:nz) )
	allocate(   P(1:nx,1:ny,1:nz) )
	allocate(   u(1:nx,1:ny,1:nz) )
	allocate(   v(1:nx,1:ny,1:nz) )
	allocate(   w(1:nx,1:ny,1:nz) )
	allocate( Estore(1:nt) )

! Initialization (by equilibrium)
	do i=1,nx,1
		do j=1,ny,1
			do k=1,nz,1
				rho(i,j,k) = rho0;
				ux(i,j,k)  = 0.0d0;
				uy(i,j,k)  = 0.0d0; 
				uz(i,j,k)  = 0.0d0; 
				call LEEquilibrium(Feq,rho(i,j,k),ux(i,j,k),uy(i,j,k),uz(i,j,k));
				F(i,j,k,:) = Feq(:);
			enddo
		enddo
	enddo
	
    do t=1,nt,1
! Print out
		E = 0.0d0;
		do i=1,nx,1
			do j=1,ny,1
				do k=1,nz,1
					E = E+0.5d0*(ux(i,j,k)**2+uy(i,j,k)**2+uz(i,j,k)**2);
				enddo
			enddo
		enddo
		E = E/(nx*ny*nz)/Ma**2;
		Estore(t) = E;
		write(*,*) t,real(t*dt),E

! Collision step
		do i=1,nx,1
			do j=1,ny,1
				do k=1,nz,1
					call Moments(rho(i,j,k),ux(i,j,k),uy(i,j,k),uz(i,j,k),F(i,j,k,:),VL);
					call LEEquilibrium(Feq,rho(i,j,k),ux(i,j,k),uy(i,j,k),uz(i,j,k));
    				!*************************************************************
                    ! Model (2): Trace
    				!*************************************************************
					!call Trace(Tr,F(i,j,k,:),VL);
					!call LQEquilibrium(Fqe,rho(i,j,k),ux(i,j,k),uy(i,j,k),uz(i,j,k),Tr);
					!Fge(:)  = (tau1/tau2)*Feq(:)+(1.d0-tau1/tau2)*Fqe(:);
					!Fp(i,j,k,:) = F(i,j,k,:)+(Fge(:)-F(i,j,k,:))/(tau1);
    				!*************************************************************
    				! Model (3): Pi_xyy,Pi_xzz,Pi_xxy,Pi_yzz,Pi_xxz,Pi_yyz
    				!*************************************************************
					call NQEquilibrium(Fqe,rho(i,j,k),ux(i,j,k),uy(i,j,k),uz(i,j,k), &
					    & F(i,j,k,:),VL);
					Fge(:)  = (tau1/tau2)*Feq(:)+(1.d0-tau1/tau2)*Fqe(:);
					Fp(i,j,k,:) = F(i,j,k,:)+(Fge(:)-F(i,j,k,:))/(tau1);
				enddo
			enddo
		enddo
		
! Streaming step (inner core)
		do i=2,nx-1,1
			do j=2,ny-1,1
				do k=2,nz-1,1
					do alpha=1,Q,1

                        F(i,j,k,alpha) = Fp(i-VL(alpha,1),j-VL(alpha,2),k-VL(alpha,3),alpha)

					enddo
				enddo
			enddo
		enddo

! Boundary conditions (external layer)
! diagonally 3D lid driven cavity
        uxB = uxL*((2.d0**0.5d0)/2.d0);
        uyB = 0.d0;
        uzB = uxL*((2.d0**0.5d0)/2.d0);
! 3D lid driven cavity
!        uxB = 0.d0;
!        uyB = 0.d0;
!        uzB = uxL*((2.d0**0.5d0)/2.d0);
        
! orthogonal to (i)
	    do j=1,ny,1
			do k=1,nz,1
			    i = 1;
                rhoW = rho(i+1,j,k);
                uxW  = 1.d0/3.d0*ux(i+1,j,k);
                uyW  = 1.d0/3.d0*uy(i+1,j,k);
                uzW  = 1.d0/3.d0*uz(i+1,j,k);
                call LEEquilibrium(FeqW,rhoW,uxW,uyW,uzW);
                call LEEquilibrium(FeqF,rho(i+1,j,k),ux(i+1,j,k),uy(i+1,j,k),uz(i+1,j,k));
			    F(i,j,k,:) = FeqW(:);!+(F(i+1,j,k,:)-FeqF(:));
			    
			    i = nx;
                rhoW = rho(i-1,j,k);
                uxW  = 1.d0/3.d0*ux(i-1,j,k);
                uyW  = 1.d0/3.d0*uy(i-1,j,k);
                uzW  = 1.d0/3.d0*uz(i-1,j,k);
                call LEEquilibrium(FeqW,rhoW,uxW,uyW,uzW);
                call LEEquilibrium(FeqF,rho(i-1,j,k),ux(i-1,j,k),uy(i-1,j,k),uz(i-1,j,k));
			    F(i,j,k,:) = FeqW(:);!+(F(i-1,j,k,:)-FeqF(:));
            enddo
        enddo
			    
! orthogonal to (k)
	    do i=1,nx,1
			do j=1,ny,1
			    k = 1;
                rhoW = rho(i,j,k+1);
                uxW  = 1.d0/3.d0*ux(i,j,k+1);
                uyW  = 1.d0/3.d0*uy(i,j,k+1);
                uzW  = 1.d0/3.d0*uz(i,j,k+1);
                call LEEquilibrium(FeqW,rhoW,uxW,uyW,uzW);
                call LEEquilibrium(FeqF,rho(i,j,k+1),ux(i,j,k+1),uy(i,j,k+1),uz(i,j,k+1));
			    F(i,j,k,:) = FeqW(:);!+(F(i,j,k+1,:)-FeqF(:));
			    
			    k = nz;
                rhoW = rho(i,j,k-1);
                uxW  = 1.d0/3.d0*ux(i,j,k-1);
                uyW  = 1.d0/3.d0*uy(i,j,k-1);
                uzW  = 1.d0/3.d0*uz(i,j,k-1);
                call LEEquilibrium(FeqW,rhoW,uxW,uyW,uzW);
                call LEEquilibrium(FeqF,rho(i,j,k-1),ux(i,j,k-1),uy(i,j,k-1),uz(i,j,k-1));
			    F(i,j,k,:) = FeqW(:);!+(F(i,j,k-1,:)-FeqF(:));
            enddo
        enddo
        
! orthogonal to (j)
	    do i=1,nx,1
			do k=1,nz,1
			    j = 1;
                rhoW = rho(i,j+1,k);
                uxW  = 1.d0/3.d0*ux(i,j+1,k);
                uyW  = 1.d0/3.d0*uy(i,j+1,k);
                uzW  = 1.d0/3.d0*uz(i,j+1,k);
                call LEEquilibrium(FeqW,rhoW,uxW,uyW,uzW);
                call LEEquilibrium(FeqF,rho(i,j+1,k),ux(i,j+1,k),uy(i,j+1,k),uz(i,j+1,k));
			    F(i,j,k,:) = FeqW(:);!+(F(i,j+1,k,:)-FeqF(:));
			    
			    j = ny;
                rhoW = rho0;!rho(i,j-1,k);
                uxW  = uxB+1.d0/3.d0*(ux(i,j-1,k)-uxB);
                uyW  = uyB+1.d0/3.d0*(uy(i,j-1,k)-uyB);
                uzW  = uzB+1.d0/3.d0*(uz(i,j-1,k)-uzB);
                call LEEquilibrium(FeqW,rhoW,uxW,uyW,uzW);
                call LEEquilibrium(FeqF,rho(i,j-1,k),ux(i,j-1,k),uy(i,j-1,k),uz(i,j-1,k));
			    F(i,j,k,:) = FeqW(:);!+(F(i,j-1,k,:)-FeqF(:));
            enddo
        enddo

!	start again

	enddo

! Print out
	write(*,*) 'Re, tau1, tau2: ',Re,tau1,tau2

! Save
	open(UNIT = 1,FILE = 'results3D.dat',IOSTAT = IO_Stat_Number1)
	if (IO_Stat_Number1 == 0) then
		print *,'Saving results3D.dat ...'

		do k=1,nz,1
			do j=1,ny,1
				do i=1,nx,1
					x = real(i)/real(nx)-1.0d0/(2.0d0*real(nx));
					y = real(j)/real(ny)-1.0d0/(2.0d0*real(ny));
					z = real(k)/real(nz)-1.0d0/(2.0d0*real(nz));

					u(i,j,k) = ux(i,j,k)/uxL;
					v(i,j,k) = uy(i,j,k)/uxL;
					w(i,j,k) = uz(i,j,k)/uxL;
					P(i,j,k) = (rho(i,j,k)-rho0)/3.d0/uxL**2;

					write(1,*) x,y,z
					write(1,*) u(i,j,k),v(i,j,k),w(i,j,k),P(i,j,k)
				enddo
			enddo
		enddo

		print *, '... done' 
	endif
	close(UNIT=1)

! Save
	open(UNIT = 1,FILE = 'log.dat',IOSTAT = IO_Stat_Number1)
	if (IO_Stat_Number1 == 0) then
		print *,'Saving results3D.dat ...'

		do t=1,nt,1
			write(1,*) Estore(t)
		enddo

		print *, '... done' 
	endif
	close(UNIT=1)
		
	deallocate(   F )
	deallocate(  Fp )
	deallocate( FWp )
	deallocate( rho )
	deallocate(  ux )
	deallocate(  uy )
	deallocate(  uz )
	deallocate(   P )
	deallocate(   u )
	deallocate(   v )
	deallocate(   w )
	deallocate( Estore )

end program lbmlid3d

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine Moments(rho,ux,uy,uz,F,VL)

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),intent(out) :: rho,ux,uy,uz
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(in) :: F
	integer,dimension(1:19,1:3),intent(in) :: VL

! Local variables
	integer :: a

	rho = 0.0d0;
	ux  = 0.0d0;
	uy  = 0.0d0;
	uz  = 0.0d0;
	do a = 1,19,1
		rho = rho+F(a);
		ux  = ux +F(a)*VL(a,1);
		uy  = uy +F(a)*VL(a,2);
		uz  = uz +F(a)*VL(a,3);
	enddo
	ux = ux/rho;
	uy = uy/rho;
	uz = uz/rho;

end subroutine Moments

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine Equilibrium(Feq,rho,ux,uy,uz)
! Entropic Equilibrium

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(out) :: Feq
	real(SELECTED_REAL_KIND(15,307)),intent(in) :: rho,ux,uy,uz
	real(SELECTED_REAL_KIND(15,307)) :: ss

	Feq(1)  =  rho*( 1.d0/3.d0-1.d0/2.d0*ux**2-1.d0/2.d0*uy**2-1.d0/2.d0*uz**2 );
	Feq(2)  =  rho*( 1.d0/18.d0+1.d0/6.d0*ux+1.d0/6.d0*ux**2-1.d0/12.d0*uy**2-1.d0/12.d0*uz**2 );
	Feq(3)  =  rho*( 1.d0/18.d0+1.d0/6.d0*uy+1.d0/6.d0*uy**2-1.d0/12.d0*ux**2-1.d0/12.d0*uz**2 );
	Feq(4)  =  rho*( 1.d0/18.d0-1.d0/6.d0*ux+1.d0/6.d0*ux**2-1.d0/12.d0*uy**2-1.d0/12.d0*uz**2 );
	Feq(5)  =  rho*( 1.d0/18.d0-1.d0/6.d0*uy+1.d0/6.d0*uy**2-1.d0/12.d0*ux**2-1.d0/12.d0*uz**2 );
	Feq(6)  =  rho*( 1.d0/12.d0*ux**2+(1.d0/12.d0+1.d0/4.d0*uy)*ux+1.d0/36.d0+1.d0/12.d0*uy &
	    & +1.d0/12.d0*uy**2-1.d0/24.d0*uz**2 );
	Feq(7)  =  rho*( 1.d0/12.d0*ux**2+(-1.d0/12.d0-1.d0/4.d0*uy)*ux+1.d0/36.d0+1.d0/12.d0*uy &
	    & +1.d0/12.d0*uy**2-1.d0/24.d0*uz**2 );
	Feq(8)  =  rho*( 1.d0/12.d0*ux**2+(-1.d0/12.d0+1.d0/4.d0*uy)*ux+1.d0/36.d0-1.d0/12.d0*uy &
	    & +1.d0/12.d0*uy**2-1.d0/24.d0*uz**2 );
	Feq(9)  =  rho*( 1.d0/12.d0*ux**2+(1.d0/12.d0-1.d0/4.d0*uy)*ux+1.d0/36.d0-1.d0/12.d0*uy &
	    & +1.d0/12.d0*uy**2-1.d0/24.d0*uz**2 );
	Feq(10) =  rho*( 1.d0/12.d0*ux**2+(1.d0/12.d0+1.d0/4.d0*uz)*ux+1.d0/36.d0+1.d0/12.d0*uz &
	    & +1.d0/12.d0*uz**2-1.d0/24.d0*uy**2 );
	Feq(11) =  rho*( 1.d0/36.d0+1.d0/12.d0*uy+1.d0/12.d0*uz+1.d0/12.d0*uy**2+1.d0/4.d0*uy*uz &
	    & +1.d0/12.d0*uz**2-1.d0/24.d0*ux**2 );
	Feq(12) =  rho*( 1.d0/12.d0*ux**2+(-1.d0/12.d0-1.d0/4.d0*uz)*ux+1.d0/36.d0+1.d0/12.d0*uz &
	    & +1.d0/12.d0*uz**2-1.d0/24.d0*uy**2 );
	Feq(13) =  rho*( 1.d0/36.d0-1.d0/12.d0*uy+1.d0/12.d0*uz+1.d0/12.d0*uy**2-1.d0/4.d0*uy*uz &
	    & +1.d0/12.d0*uz**2-1.d0/24.d0*ux**2 );
	Feq(14) =  rho*( 1.d0/18.d0+1.d0/6.d0*uz+1.d0/6.d0*uz**2-1.d0/12.d0*ux**2-1.d0/12.d0*uy**2 );
	Feq(15) =  rho*( 1.d0/12.d0*ux**2+(1.d0/12.d0-1.d0/4.d0*uz)*ux+1.d0/36.d0-1.d0/12.d0*uz &
	    & +1.d0/12.d0*uz**2-1.d0/24.d0*uy**2 );
	Feq(16) =  rho*( 1.d0/36.d0+1.d0/12.d0*uy-1.d0/12.d0*uz+1.d0/12.d0*uy**2-1.d0/4.d0*uy*uz &
	    & +1.d0/12.d0*uz**2-1.d0/24.d0*ux**2 );
	Feq(17) =  rho*( 1.d0/12.d0*ux**2+(-1.d0/12.d0+1.d0/4.d0*uz)*ux+1.d0/36.d0-1.d0/12.d0*uz &
	    & +1.d0/12.d0*uz**2-1.d0/24.d0*uy**2 );
	Feq(18) =  rho*( 1.d0/36.d0-1.d0/12.d0*uy-1.d0/12.d0*uz+1.d0/12.d0*uy**2+1.d0/4.d0*uy*uz &
	    & +1.d0/12.d0*uz**2-1.d0/24.d0*ux**2 );
	Feq(19) =  rho*( 1.d0/18.d0-1.d0/6.d0*uz+1.d0/6.d0*uz**2-1.d0/12.d0*ux**2-1.d0/12.d0*uy**2 );

end subroutine Equilibrium

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine LEEquilibrium(Feq,rho,ux,uy,uz)
! Entropic Equilibrium

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(out) :: Feq
	real(SELECTED_REAL_KIND(15,307)),intent(in) :: rho,ux,uy,uz
	real(SELECTED_REAL_KIND(15,307)) :: Sxx,Syy,Szz

	Sxx = 1.d0/3.d0+ux**2;
	Syy = 1.d0/3.d0+uy**2;
	Szz = 1.d0/3.d0+uz**2;

	Feq(1)=  rho*( 1.d0-Sxx-Syy-Szz+Sxx*Syy+Syy*Szz+Sxx*Szz ); 
	Feq(2)=  rho*( -1.d0/2.d0*(ux+Sxx)*(-1.d0+Syy+Szz) );
	Feq(3)=  rho*( -1.d0/2.d0*(-1.d0+Sxx+Szz)*(uy+Syy) ); 
	Feq(4)=  rho*( -1.d0/2.d0*(-ux+Sxx)*(-1.d0+Syy+Szz) );
	Feq(5)=  rho*( -1.d0/2.d0*(-1.d0+Sxx+Szz)*(-uy+Syy) );
	Feq(6)=  rho*( 1.d0/4.d0*(ux+Sxx)*(uy+Syy) );
	Feq(7)=  rho*( 1.d0/4.d0*(-ux+Sxx)*(uy+Syy) );
	Feq(8)=  rho*( 1.d0/4.d0*(-ux+Sxx)*(-uy+Syy) );
	Feq(9)=  rho*( 1.d0/4.d0*(ux+Sxx)*(-uy+Syy) );
	Feq(10)= rho*( 1.d0/4.d0*(ux+Sxx)*(uz+Szz) );
	Feq(11)= rho*( 1.d0/4.d0*(uz+Szz)*(uy+Syy) );
	Feq(12)= rho*( 1.d0/4.d0*(-ux+Sxx)*(uz+Szz) );
	Feq(13)= rho*( 1.d0/4.d0*(uz+Szz)*(-uy+Syy) );
	Feq(14)= rho*( -1.d0/2.d0*(uz+Szz)*(Sxx+Syy-1.d0) );
	Feq(15)= rho*( -1.d0/4.d0*(ux+Sxx)*(uz-Szz) );
	Feq(16)= rho*( -1.d0/4.d0*(uz-Szz)*(uy+Syy) );
	Feq(17)= rho*( -1.d0/4.d0*(-ux+Sxx)*(uz-Szz) );
	Feq(18)= rho*( -1.d0/4.d0*(uz-Szz)*(-uy+Syy) );
	Feq(19)= rho*( 1.d0/2.d0*(uz-Szz)*(Sxx+Syy-1.d0) );

end subroutine LEEquilibrium

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine LQEquilibrium(Feq,rho,ux,uy,uz,Tr)
! Quasi Equilibrium

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(out) :: Feq
	real(SELECTED_REAL_KIND(15,307)),intent(in) :: rho,ux,uy,uz,Tr
	real(SELECTED_REAL_KIND(15,307)) :: uu,Sxx,Syy,Szz

	uu  = ux**2+uy**2+uz**2;
	Sxx = Tr/3.d0+ux**2-uu/3.d0;
	Syy = Tr/3.d0+uy**2-uu/3.d0;
	Szz = Tr/3.d0+uz**2-uu/3.d0;

	Feq(1)=  rho*( 1.d0-Sxx-Syy-Szz+Sxx*Syy+Syy*Szz+Sxx*Szz ); 
	Feq(2)=  rho*( -1.d0/2.d0*(ux+Sxx)*(-1.d0+Syy+Szz) );
	Feq(3)=  rho*( -1.d0/2.d0*(-1.d0+Sxx+Szz)*(uy+Syy) ); 
	Feq(4)=  rho*( -1.d0/2.d0*(-ux+Sxx)*(-1.d0+Syy+Szz) );
	Feq(5)=  rho*( -1.d0/2.d0*(-1.d0+Sxx+Szz)*(-uy+Syy) );
	Feq(6)=  rho*( 1.d0/4.d0*(ux+Sxx)*(uy+Syy) );
	Feq(7)=  rho*( 1.d0/4.d0*(-ux+Sxx)*(uy+Syy) );
	Feq(8)=  rho*( 1.d0/4.d0*(-ux+Sxx)*(-uy+Syy) );
	Feq(9)=  rho*( 1.d0/4.d0*(ux+Sxx)*(-uy+Syy) );
	Feq(10)= rho*( 1.d0/4.d0*(ux+Sxx)*(uz+Szz) );
	Feq(11)= rho*( 1.d0/4.d0*(uz+Szz)*(uy+Syy) );
	Feq(12)= rho*( 1.d0/4.d0*(-ux+Sxx)*(uz+Szz) );
	Feq(13)= rho*( 1.d0/4.d0*(uz+Szz)*(-uy+Syy) );
	Feq(14)= rho*( -1.d0/2.d0*(uz+Szz)*(Sxx+Syy-1.d0) );
	Feq(15)= rho*( -1.d0/4.d0*(ux+Sxx)*(uz-Szz) );
	Feq(16)= rho*( -1.d0/4.d0*(uz-Szz)*(uy+Syy) );
	Feq(17)= rho*( -1.d0/4.d0*(-ux+Sxx)*(uz-Szz) );
	Feq(18)= rho*( -1.d0/4.d0*(uz-Szz)*(-uy+Syy) );
	Feq(19)= rho*( 1.d0/2.d0*(uz-Szz)*(Sxx+Syy-1.d0) );

end subroutine LQEquilibrium

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine Trace(Tr,F,VL)

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),intent(out) :: Tr
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(in) :: F
	integer,dimension(1:19,1:3),intent(in) :: VL

! Local variables
	real(SELECTED_REAL_KIND(15,307)) :: rho
	integer :: a

	rho = 0.d0;
	Tr  = 0.d0;
	do a = 1,19,1
		rho = rho+F(a);
		Tr  = Tr +(VL(a,1)**2+VL(a,2)**2+VL(a,3)**2)*F(a);
	enddo
	Tr = Tr/rho;

end subroutine Trace

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine NQEquilibrium(NFqe,rho,ux,uy,uz,F,VL)
! New Quasi Equilibrium

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(out) :: NFqe
	real(SELECTED_REAL_KIND(15,307)),intent(in) :: rho,ux,uy,uz
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(in) :: F
	integer,dimension(1:19,1:3),intent(in) :: VL

	real(SELECTED_REAL_KIND(15,307)) :: mxyy,mxzz,mxxy,myzz,mxxz,myyz
    integer :: a

	mxyy = 0.0d0;
    mxzz = 0.0d0;
    mxxy = 0.0d0;
    myzz = 0.0d0;
    mxxz = 0.0d0;
    myyz = 0.0d0;

	do a = 1,19,1
	    mxyy = mxyy+F(a)*VL(a,1)*VL(a,2)*VL(a,2);
        mxzz = mxzz+F(a)*VL(a,1)*VL(a,3)*VL(a,3);
        mxxy = mxxy+F(a)*VL(a,1)*VL(a,1)*VL(a,2);
        myzz = myzz+F(a)*VL(a,2)*VL(a,3)*VL(a,3);
        mxxz = mxxz+F(a)*VL(a,1)*VL(a,1)*VL(a,3);
        myyz = myyz+F(a)*VL(a,2)*VL(a,2)*VL(a,3);
 	enddo
	
    NFqe(1) = rho*(ux**2*3+uy**2*3+uz**2*3-2)*(-1.0D0/6.0D0);
    
    NFqe(2) = mxyy*(-1.0D0/2.0D0)-mxzz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1&
        &)+rho*ux*(1.0D0/2.0D0)+rho*ux**2*(1.0D0/6.0D0)-rho*uy**2*(1.0D0/1.&
        &2D1)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(3) = mxxy*(-1.0D0/2.0D0)-myzz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1&
        &)+rho*uy*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/6.&
        &0D0)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(4) = mxyy*(1.0D0/2.0D0)+mxzz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1)&
        &-rho*ux*(1.0D0/2.0D0)+rho*ux**2*(1.0D0/6.0D0)-rho*uy**2*(1.0D0/1.2&
        &D1)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(5) = mxxy*(1.0D0/2.0D0)+myzz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1)&
        &-rho*uy*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/6.0&
        &D0)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(6) = mxxy*(1.0D0/4.0D0)+mxyy*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1)&
        &+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.0D0/&
        &2.4D1)+rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(7) = mxxy*(1.0D0/4.0D0)-mxyy*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1)&
        &+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.0D0/&
        &2.4D1)-rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(8) = mxxy*(-1.0D0/4.0D0)-mxyy*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.0D0&
        &/2.4D1)+rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(9) = mxxy*(-1.0D0/4.0D0)+mxyy*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.0D0&
        &/2.4D1)-rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(10) = mxxz*(1.0D0/4.0D0)+mxzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1.0D0&
        &/1.2D1)+rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(11) = myyz*(1.0D0/4.0D0)+myzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1.0D0&
        &/1.2D1)+rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(12) = mxxz*(1.0D0/4.0D0)-mxzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1.0D0&
        &/1.2D1)-rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(13) = myyz*(1.0D0/4.0D0)-myzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1.0D0&
        &/1.2D1)-rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(14) = mxxz*(-1.0D0/2.0D0)-myyz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D&
        &1)+rho*uz*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/1&
        &.2D1)+rho*uz**2*(1.0D0/6.0D0);
        
    NFqe(15) = mxxz*(-1.0D0/4.0D0)+mxzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D&
        &1)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1.0D&
        &0/1.2D1)-rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(16) = myyz*(-1.0D0/4.0D0)+myzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D&
        &1)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1.0D&
        &0/1.2D1)-rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(17) = mxxz*(-1.0D0/4.0D0)-mxzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D&
        &1)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1.0D&
        &0/1.2D1)+rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(18) = myyz*(-1.0D0/4.0D0)-myzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D&
        &1)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1.0D&
        &0/1.2D1)+rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(19) = mxxz*(1.0D0/2.0D0)+myyz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1&
        &)-rho*uz*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/1.&
        &2D1)+rho*uz**2*(1.0D0/6.0D0);
          
end subroutine NQEquilibrium

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine NNQEquilibrium(NFqe,rho,ux,uy,uz,F,VL)
! New Quasi Equilibrium

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(out) :: NFqe
	real(SELECTED_REAL_KIND(15,307)),intent(in) :: rho,ux,uy,uz
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(in) :: F
	integer,dimension(1:19,1:3),intent(in) :: VL

	real(SELECTED_REAL_KIND(15,307)) :: pxx,pyy,pzz
    integer :: a

    pxx = 0.0d0;
    pyy = 0.0d0;
    pzz = 0.0d0;
	do a = 1,19,1
        pxx = pxx+F(a)*VL(a,1)*VL(a,1);
        pyy = pyy+F(a)*VL(a,2)*VL(a,2);
        pzz = pzz+F(a)*VL(a,3)*VL(a,3);
 	enddo
    pxx = pxx/rho; 	
    pyy = pyy/rho; 	
    pzz = pzz/rho; 	

    NFqe(1) = rho*(ux**2*3+uy**2*3+uz**2*3-2)*(-1.0D0/6.0D0);
    
    NFqe(2) = rho*(1.0D0/1.8D1)-pyy*ux*(1.0D0/2.0D0)-pzz*ux*(1.0D0/2.&
        &0D0)+rho*ux*(1.0D0/2.0D0)+rho*ux**2*(1.0D0/6.0D0)-rho*uy**2*(1.0D0&
        &/1.2D1)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(3) = rho*(1.0D0/1.8D1)-pxx*uy*(1.0D0/2.0D0)-pzz*uy*(1.0D0/2.&
        &0D0)+rho*uy*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0&
        &/6.0D0)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(4) = rho*(1.0D0/1.8D1)+pyy*ux*(1.0D0/2.0D0)+pzz*ux*(1.0D0/2.&
        &0D0)-rho*ux*(1.0D0/2.0D0)+rho*ux**2*(1.0D0/6.0D0)-rho*uy**2*(1.0D0&
        &/1.2D1)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(5) = rho*(1.0D0/1.8D1)+pxx*uy*(1.0D0/2.0D0)+pzz*uy*(1.0D0/2.&
        &0D0)-rho*uy*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0&
        &/6.0D0)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(6) = rho*(1.0D0/3.6D1)+pxx*uy*(1.0D0/4.0D0)+pyy*ux*(1.0D0/4.&
        &0D0)+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.&
        &0D0/2.4D1)+rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(7) = rho*(1.0D0/3.6D1)+pxx*uy*(1.0D0/4.0D0)-pyy*ux*(1.0D0/4.&
        &0D0)+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.&
        &0D0/2.4D1)-rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(8) = rho*(1.0D0/3.6D1)-pxx*uy*(1.0D0/4.0D0)-pyy*ux*(1.0D0/4.&
        &0D0)+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.&
        &0D0/2.4D1)+rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(9) = rho*(1.0D0/3.6D1)-pxx*uy*(1.0D0/4.0D0)+pyy*ux*(1.0D0/4.&
        &0D0)+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.&
        &0D0/2.4D1)-rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(10) = rho*(1.0D0/3.6D1)+pxx*uz*(1.0D0/4.0D0)+pzz*ux*(1.0D0/4&
        &.0D0)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1&
        &.0D0/1.2D1)+rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(11) = rho*(1.0D0/3.6D1)+pyy*uz*(1.0D0/4.0D0)+pzz*uy*(1.0D0/4&
        &.0D0)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1&
        &.0D0/1.2D1)+rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(12) = rho*(1.0D0/3.6D1)+pxx*uz*(1.0D0/4.0D0)-pzz*ux*(1.0D0/4&
        &.0D0)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1&
        &.0D0/1.2D1)-rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(13) = rho*(1.0D0/3.6D1)+pyy*uz*(1.0D0/4.0D0)-pzz*uy*(1.0D0/4&
        &.0D0)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1&
        &.0D0/1.2D1)-rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(14) = rho*(1.0D0/1.8D1)-pxx*uz*(1.0D0/2.0D0)-pyy*uz*(1.0D0/2&
        &.0D0)+rho*uz*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D&
        &0/1.2D1)+rho*uz**2*(1.0D0/6.0D0);
        
    NFqe(15) = rho*(1.0D0/3.6D1)-pxx*uz*(1.0D0/4.0D0)+pzz*ux*(1.0D0/4&
        &.0D0)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1&
        &.0D0/1.2D1)-rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(16) = rho*(1.0D0/3.6D1)-pyy*uz*(1.0D0/4.0D0)+pzz*uy*(1.0D0/4&
        &.0D0)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1&
        &.0D0/1.2D1)-rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(17) = rho*(1.0D0/3.6D1)-pxx*uz*(1.0D0/4.0D0)-pzz*ux*(1.0D0/4&
        &.0D0)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1&
        &.0D0/1.2D1)+rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(18) = rho*(1.0D0/3.6D1)-pyy*uz*(1.0D0/4.0D0)-pzz*uy*(1.0D0/4&
        &.0D0)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1&
        &.0D0/1.2D1)+rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(19) = rho*(1.0D0/1.8D1)+pxx*uz*(1.0D0/2.0D0)+pyy*uz*(1.0D0/2&
        &.0D0)-rho*uz*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D&
        &0/1.2D1)+rho*uz**2*(1.0D0/6.0D0);
 	
end subroutine NNQEquilibrium 	

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine NNNQEquilibrium(NFqe,rho,ux,uy,uz,F,VL)
! New Quasi Equilibrium

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(out) :: NFqe
	real(SELECTED_REAL_KIND(15,307)),intent(in) :: rho,ux,uy,uz
	real(SELECTED_REAL_KIND(15,307)),dimension(1:19),intent(in) :: F
	integer,dimension(1:19,1:3),intent(in) :: VL

	real(SELECTED_REAL_KIND(15,307)) :: mxyy,mxzz,mxxy,myzz,mxxz,myyz
	real(SELECTED_REAL_KIND(15,307)) :: pxx,pyy,pzz
    integer :: a

	mxyy = 0.0d0;
    mxzz = 0.0d0;
    mxxy = 0.0d0;
    myzz = 0.0d0;
    mxxz = 0.0d0;
    myyz = 0.0d0;

	do a = 1,19,1
	    mxyy = mxyy+F(a)*VL(a,1)*VL(a,2)*VL(a,2);
        mxzz = mxzz+F(a)*VL(a,1)*VL(a,3)*VL(a,3);
        mxxy = mxxy+F(a)*VL(a,1)*VL(a,1)*VL(a,2);
        myzz = myzz+F(a)*VL(a,2)*VL(a,3)*VL(a,3);
        mxxz = mxxz+F(a)*VL(a,1)*VL(a,1)*VL(a,3);
        myyz = myyz+F(a)*VL(a,2)*VL(a,2)*VL(a,3);
 	enddo
	
! ATTENTION
!    pxx = 0.0d0;
!    pyy = 0.0d0;
!    pzz = 0.0d0;
!	do a = 1,19,1
!        pxx = pxx+F(a)*VL(a,1)*VL(a,1);
!        pyy = pyy+F(a)*VL(a,2)*VL(a,2);
!        pzz = pzz+F(a)*VL(a,3)*VL(a,3);
! 	enddo
!    pxx = pxx/rho; 	
!    pyy = pyy/rho; 	
!    pzz = pzz/rho; 	
!    
!    mxyy = rho*ux*pyy;
!    mxzz = rho*ux*pzz;
!    mxxy = rho*pxx*uy;
!    myzz = rho*uy*pzz;
!    mxxz = rho*pxx*uz;
!    myyz = rho*pyy*uz;
! ATTENTION    
	
    NFqe(1) = rho*(ux**2*3+uy**2*3+uz**2*3-2)*(-1.0D0/6.0D0);
    
    NFqe(2) = mxyy*(-1.0D0/2.0D0)-mxzz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1&
        &)+rho*ux*(1.0D0/2.0D0)+rho*ux**2*(1.0D0/6.0D0)-rho*uy**2*(1.0D0/1.&
        &2D1)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(3) = mxxy*(-1.0D0/2.0D0)-myzz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1&
        &)+rho*uy*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/6.&
        &0D0)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(4) = mxyy*(1.0D0/2.0D0)+mxzz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1)&
        &-rho*ux*(1.0D0/2.0D0)+rho*ux**2*(1.0D0/6.0D0)-rho*uy**2*(1.0D0/1.2&
        &D1)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(5) = mxxy*(1.0D0/2.0D0)+myzz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1)&
        &-rho*uy*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/6.0&
        &D0)-rho*uz**2*(1.0D0/1.2D1);
        
    NFqe(6) = mxxy*(1.0D0/4.0D0)+mxyy*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1)&
        &+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.0D0/&
        &2.4D1)+rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(7) = mxxy*(1.0D0/4.0D0)-mxyy*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1)&
        &+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.0D0/&
        &2.4D1)-rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(8) = mxxy*(-1.0D0/4.0D0)-mxyy*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.0D0&
        &/2.4D1)+rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(9) = mxxy*(-1.0D0/4.0D0)+mxyy*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)+rho*ux**2*(1.0D0/1.2D1)+rho*uy**2*(1.0D0/1.2D1)-rho*uz**2*(1.0D0&
        &/2.4D1)-rho*ux*uy*(1.0D0/4.0D0);
        
    NFqe(10) = mxxz*(1.0D0/4.0D0)+mxzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1.0D0&
        &/1.2D1)+rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(11) = myyz*(1.0D0/4.0D0)+myzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1.0D0&
        &/1.2D1)+rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(12) = mxxz*(1.0D0/4.0D0)-mxzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1.0D0&
        &/1.2D1)-rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(13) = myyz*(1.0D0/4.0D0)-myzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D1&
        &)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1.0D0&
        &/1.2D1)-rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(14) = mxxz*(-1.0D0/2.0D0)-myyz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D&
        &1)+rho*uz*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/1&
        &.2D1)+rho*uz**2*(1.0D0/6.0D0);
        
    NFqe(15) = mxxz*(-1.0D0/4.0D0)+mxzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D&
        &1)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1.0D&
        &0/1.2D1)-rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(16) = myyz*(-1.0D0/4.0D0)+myzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D&
        &1)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1.0D&
        &0/1.2D1)-rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(17) = mxxz*(-1.0D0/4.0D0)-mxzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D&
        &1)+rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/2.4D1)+rho*uz**2*(1.0D&
        &0/1.2D1)+rho*ux*uz*(1.0D0/4.0D0);
        
    NFqe(18) = myyz*(-1.0D0/4.0D0)-myzz*(1.0D0/4.0D0)+rho*(1.0D0/3.6D&
        &1)-rho*ux**2*(1.0D0/2.4D1)+rho*uy**2*(1.0D0/1.2D1)+rho*uz**2*(1.0D&
        &0/1.2D1)+rho*uy*uz*(1.0D0/4.0D0);
        
    NFqe(19) = mxxz*(1.0D0/2.0D0)+myyz*(1.0D0/2.0D0)+rho*(1.0D0/1.8D1&
        &)-rho*uz*(1.0D0/2.0D0)-rho*ux**2*(1.0D0/1.2D1)-rho*uy**2*(1.0D0/1.&
        &2D1)+rho*uz**2*(1.0D0/6.0D0);
          
end subroutine NNNQEquilibrium