!
!  LBMLID.f95 
!
!  FUNCTIONS:
!	LBMLID			- Entry point of console application.
!	EEquilibrium	- Entropic equilibrium
!	QEquilibrium	- Quasi equilibrium
!	Force			- External force in velocity space
!	Moments			- Conserved moments
!	Trace			- Trace of the pressure tensor
!
!****************************************************************************
!
!  PROGRAM: LBMLID, by Pietro ASINARI, 2011
!
!  PURPOSE: Simple LBM code for solving 2D 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 lbmlid

	implicit none

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

! Local variables
!********************
	real(SELECTED_REAL_KIND(15,307)),dimension(:,:,:),allocatable	:: F,Fp
	real(SELECTED_REAL_KIND(15,307)),dimension(:,:),allocatable		:: rho,ux,uy,P,u,v

	real(SELECTED_REAL_KIND(15,307)),dimension(1:9)		:: Feq,Fqe,Fge,G,Floc
	integer,dimension(1:9,1:2)	:: VL
	integer,dimension(1:9)		:: BB

	real(SELECTED_REAL_KIND(15,307)) :: E,x,y,temp
	real(SELECTED_REAL_KIND(15,307)) :: Sxx,Syy,Tr,Treq
	integer :: i,j,t,alpha,Nin,RR,LL
	integer :: IO_Stat_Number1 = -1

! Parameters
!********************
	real(SELECTED_REAL_KIND(15,307)),parameter :: Tt =  100.d0
	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.001d0
! bulk viscosity
	real(SELECTED_REAL_KIND(15,307)),parameter :: XI =  0.010d0
! space discretization
	integer,parameter :: nx = 30

	real(SELECTED_REAL_KIND(15,307)),parameter :: Re = UL*Lx/NU;

	integer,parameter :: ny   = nx
	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 :: 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.05d0;

	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 :: 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 = 9

! Lattice D2Q9
	VL(:,1) = (/0, 1, 0, -1, 0, 1, -1, -1, 1/) ! x-axis
	VL(:,2) = (/0, 0, 1, 0, -1, 1, 1, -1, -1/) ! y-axis
	BB(:)   = (/1, 4, 5, 2, 3, 8, 9, 6, 7/) ! bounce back

    nt = int(Tt/dt);

! Memory Allocation
	allocate(   F(1:nx,1:ny,1:9) )
	allocate(  Fp(1:nx,1:ny,1:9) )
	allocate( rho(1:nx,1:ny) )
	allocate(  ux(1:nx,1:ny) )
	allocate(  uy(1:nx,1:ny) )
	allocate(   P(1:nx,1:ny) )
	allocate(   u(1:nx,1:ny) )
	allocate(   v(1:nx,1:ny) )

	temp = tau1;

! Initialization (by equilibrium)
	do i=1,nx,1
		do j=1,ny,1
			rho(i,j) = rho0;
			ux(i,j)  = 0.0d0;
			uy(i,j)  = 0.0d0; 
			call LEEquilibrium(Feq,rho(i,j),ux(i,j),uy(i,j));
			F(i,j,:)=Feq(:);
		enddo
	enddo

	do t=1,nt,1

! Print out
		E = 0.0d0;
		do i=1,nx,1
			do j=1,ny,1
				E = E+0.5d0*(ux(i,j)**2+uy(i,j)**2);
			enddo
		enddo
		E = E/(nx*ny)/Ma**2;
		write(*,*) t,real(t*dt),E
	    
! Collision step
		do i=1,nx,1
			do j=1,ny,1
				call Moments(rho(i,j),ux(i,j),uy(i,j),F(i,j,:),VL);
				call LEEquilibrium(Feq,rho(i,j),ux(i,j),uy(i,j));
				!*************************************************************
                ! Model (2)
				!*************************************************************
                ! quasi-equilibrium parametrized by the second-order trace
				!*************************************************************
				call Trace(Tr,F(i,j,:),VL);
				call LQEquilibrium(Fqe,rho(i,j),ux(i,j),uy(i,j),Tr);
				Fge(:) = (tau1/tau2)*Feq(:)+(1.d0-tau1/tau2)*Fqe(:);
				Fp(i,j,:) = F(i,j,:)+(Fge(:)-F(i,j,:))/(tau1);
				!*************************************************************
				! Model (3)
				!*************************************************************
				! quasi-equilibrium parametrized by Pi_xxy and Pi_xyy
				!*************************************************************
				!call NQEquilibrium(Fqe,rho(i,j),ux(i,j),uy(i,j),F(i,j,:),VL);
				!Fge(:) = (tau1/tau2)*Feq(:)+(1.d0-tau1/tau2)*Fqe(:);
				!Fp(i,j,:) = F(i,j,:)+(Fge(:)-F(i,j,:))/(tau1);				
			enddo
		enddo
    
! Streaming step
		do i=1,nx,1
			do j=1,ny,1
				do alpha=1,Q,1

					if ( (j+VL(alpha,2))>ny ) then
						! Upper wall / LID
						! PUSH
						call LEEquilibrium(Feq,rho(i,j),0.d0,0.d0);
						F(i,j,BB(alpha)) = Fp(i,j,alpha)-6.d0*Feq(alpha)*(VL(alpha,1)*uxL+VL(alpha,2)*uyL);
						!						
					elseif ( (j+VL(alpha,2))<1 ) then
						! Lower wall
						F(i,j,BB(alpha)) = Fp(i,j,alpha);
						!
					elseif ( (i+VL(alpha,1))>nx ) then
						! Right wall
						F(i,j,BB(alpha)) = Fp(i,j,alpha);
						!
					elseif ( (i+VL(alpha,1))<1 ) then
						! Left wall
						F(i,j,BB(alpha)) = Fp(i,j,alpha);
						!
					else
						! Bulk
						F(i+VL(alpha,1),j+VL(alpha,2),alpha) = Fp(i,j,alpha);
						!
					endif

				enddo
			enddo
		enddo

!	start again

	enddo

! Print out
	write(*,*) 'Re, tau1, tau2: ',Re,tau1,tau2
	
! Save
	open(UNIT = 1,FILE = 'results.dat',IOSTAT = IO_Stat_Number1)

	if (IO_Stat_Number1 == 0) then
		print *,'Saving results.dat ...'

		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));

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

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

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

	read *, temp

end program lbmlid

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

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

	implicit none

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

	Feq(1)=-2.0d0/9.0d0*rho*(-2.0d0+3.0d0*ux**2.0d0+3.0d0*uy**2);

	Feq(2)=1.0d0/18.0d0*rho*(6.0d0*ux+6.0d0*ux**2+2.0d0-3.0d0*uy**2);
	Feq(3)=-1.0d0/18.0d0*rho*(-6.0d0*uy-6.0d0*uy**2-2.0d0+3.0d0*ux**2);
	Feq(4)=1.0d0/18.0d0*rho*(-6.0d0*ux+6.0d0*ux**2+2.0d0-3.0d0*uy**2);
	Feq(5)=-1.0d0/18.0d0*rho*(6.0d0*uy-6.0d0*uy**2-2.0d0+3.0d0*ux**2);

	Feq(6)=1.0d0/36.0d0*rho*(9.0d0*ux*uy+3.0d0*ux+3.0d0*uy+1.0d0+3.0d0*ux**2+3.0d0*uy**2);
	Feq(7)=1.0d0/36.0d0*rho*(-9.0d0*ux*uy-3.0d0*ux+3.0d0*uy+1.0d0+3.0d0*ux**2+3.0d0*uy**2);
	Feq(8)=1.0d0/36.0d0*rho*(9.0d0*ux*uy-3.0d0*ux-3.0d0*uy+1.0d0+3.0d0*ux**2+3.0d0*uy**2);
	Feq(9)=1.0d0/36.0d0*rho*(-9.0d0*ux*uy+3.0d0*ux-3.0d0*uy+1.0d0+3.0d0*ux**2+3.0d0*uy**2);

end subroutine Equilibrium

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine Force(G,rho,gx,gy)

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),dimension(1:9),intent(out) :: G
	real(SELECTED_REAL_KIND(15,307)),intent(in) :: rho,gx,gy

	G(1)=0.0d0;

	G(2)=1.0d0/2.0d0*rho*gx;
	G(3)=1.0d0/2.0d0*rho*gy;
	G(4)=-1.0d0/2.0d0*rho*gx;
	G(5)=-1.0d0/2.0d0*rho*gy;

	G(6)=0.0d0;
	G(7)=0.0d0;
	G(8)=0.0d0;
	G(9)=0.0d0;

end subroutine Force

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

	implicit none

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

! Local variables
	integer :: a

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

end subroutine Moments

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

	implicit none

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

	Sxx = -1.d0/3.d0+2.d0/3.d0*(3.d0*ux**2+1.d0)**(0.5d0);
	Syy = -1.d0/3.d0+2.d0/3.d0*(3.d0*uy**2+1.d0)**(0.5d0);

	Feq(1) = rho*(1.d0-Sxx)*(1.d0-Syy);

	Feq(2) = 0.5d0*rho*(ux+Sxx)*(1.d0-Syy);
	Feq(3) = 0.5d0*rho*(1.d0-Sxx)*(uy+Syy);
	Feq(4) = 0.5d0*rho*(Sxx-ux)*(1.d0-Syy);
	Feq(5) = 0.5d0*rho*(1.d0-Sxx)*(Syy-uy);

	Feq(6) = 0.25d0*rho*(uy+Syy)*(ux+Sxx);
	Feq(7) = 0.25d0*rho*(uy+Syy)*(Sxx-ux);
	Feq(8) = 0.25d0*rho*(Syy-uy)*(Sxx-ux);
	Feq(9) = 0.25d0*rho*(Syy-uy)*(ux+Sxx);

end subroutine EEquilibrium

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

	implicit none

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

!	Sxx = -1.d0/3.d0+2.d0/3.d0*(3.d0*ux**2+1.d0)**(0.5d0);
!	Syy = -1.d0/3.d0+2.d0/3.d0*(3.d0*uy**2+1.d0)**(0.5d0);

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

	Feq(1) = rho*(1.d0-Sxx)*(1.d0-Syy);

	Feq(2) = 0.5d0*rho*(ux+Sxx)*(1.d0-Syy);
	Feq(3) = 0.5d0*rho*(1.d0-Sxx)*(uy+Syy);
	Feq(4) = 0.5d0*rho*(Sxx-ux)*(1.d0-Syy);
	Feq(5) = 0.5d0*rho*(1.d0-Sxx)*(Syy-uy);

	Feq(6) = 0.25d0*rho*(uy+Syy)*(ux+Sxx);
	Feq(7) = 0.25d0*rho*(uy+Syy)*(Sxx-ux);
	Feq(8) = 0.25d0*rho*(Syy-uy)*(Sxx-ux);
	Feq(9) = 0.25d0*rho*(Syy-uy)*(ux+Sxx);

end subroutine LEEquilibrium

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

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),dimension(1:9),intent(out) :: Feq
	real(SELECTED_REAL_KIND(15,307)),intent(in) :: rho,ux,uy,Tr
	real(SELECTED_REAL_KIND(15,307)) :: a,b,c,p,q,r,D,Sxx,Syy

! Cardano formula
	a = -0.5d0*(ux**2-uy**2);
	b = (2.d0-Tr)*(Tr-(ux**2+uy**2));
	c = -0.5d0*((2.d0-Tr)**2)*(ux**2-uy**2);
	p = -a**2/3.d0+b;
	q = 2.d0*(a**3)/27.d0-a*b/3.d0+c;
	D = (q**2)/4.d0+(p**3)/27.d0;
	r = (-q/2.d0+dsqrt(D))**(1.d0/3.d0);

	Sxx = Tr/2.d0+0.5d0*(r-p/(3.d0*r)-a/3.d0);
	Syy = Tr-Sxx;

	Feq(1) = rho*(1.d0-Sxx)*(1.d0-Syy);

	Feq(2) = 0.5d0*rho*(ux+Sxx)*(1.d0-Syy);
	Feq(3) = 0.5d0*rho*(1.d0-Sxx)*(uy+Syy);
	Feq(4) = 0.5d0*rho*(Sxx-ux)*(1.d0-Syy);
	Feq(5) = 0.5d0*rho*(1.d0-Sxx)*(Syy-uy);

	Feq(6) = 0.25d0*rho*(uy+Syy)*(ux+Sxx);
	Feq(7) = 0.25d0*rho*(uy+Syy)*(Sxx-ux);
	Feq(8) = 0.25d0*rho*(Syy-uy)*(Sxx-ux);
	Feq(9) = 0.25d0*rho*(Syy-uy)*(ux+Sxx);

end subroutine QEquilibrium

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

	implicit none

! Inlet/Outlet
	real(SELECTED_REAL_KIND(15,307)),dimension(1:9),intent(out) :: Feq
	real(SELECTED_REAL_KIND(15,307)),intent(in) :: rho,ux,uy,Tr
	real(SELECTED_REAL_KIND(15,307)) :: a,b,c,p,q,r,D,Sxx,Syy

! Cardano formula
!	a = -0.5d0*(ux**2-uy**2);
!	b = (2.d0-Tr)*(Tr-(ux**2+uy**2));
!	c = -0.5d0*((2.d0-Tr)**2)*(ux**2-uy**2);
!	p = -a**2/3.d0+b;
!	q = 2.d0*(a**3)/27.d0-a*b/3.d0+c;
!	D = (q**2)/4.d0+(p**3)/27.d0;
!	r = (-q/2.d0+dsqrt(D))**(1.d0/3.d0);
!	Sxx = Tr/2.d0+0.5d0*(r-p/(3.d0*r)-a/3.d0);

	Sxx = (Tr+(ux**2-uy**2))/2.d0;
	Syy = Tr-Sxx;

	Feq(1) = rho*(1.d0-Sxx)*(1.d0-Syy);

	Feq(2) = 0.5d0*rho*(ux+Sxx)*(1.d0-Syy);
	Feq(3) = 0.5d0*rho*(1.d0-Sxx)*(uy+Syy);
	Feq(4) = 0.5d0*rho*(Sxx-ux)*(1.d0-Syy);
	Feq(5) = 0.5d0*rho*(1.d0-Sxx)*(Syy-uy);

	Feq(6) = 0.25d0*rho*(uy+Syy)*(ux+Sxx);
	Feq(7) = 0.25d0*rho*(uy+Syy)*(Sxx-ux);
	Feq(8) = 0.25d0*rho*(Syy-uy)*(Sxx-ux);
	Feq(9) = 0.25d0*rho*(Syy-uy)*(ux+Sxx);

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:9),intent(in) :: F
	integer,dimension(1:9,1:2),intent(in) :: VL

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

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

end subroutine Trace

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

	implicit none

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

	real(SELECTED_REAL_KIND(15,307)) :: mxyy,mxxy
    integer :: a

	mxyy = 0.0d0;
    mxxy = 0.0d0;
	do a = 1,9,1
	    mxyy = mxyy+F(a)*VL(a,1)*VL(a,2)*VL(a,2);
        mxxy = mxxy+F(a)*VL(a,1)*VL(a,1)*VL(a,2);
 	enddo

    NFqe(1) = rho*(ux**2*3+uy**2*3-2)*(-2.0D0/9.0D0);
    
    NFqe(2) = mxyy*(-1.0D0/2.0D0)+rho*(1.0D0/9.0D0)+rho*ux*(1.0D0/2.0&
        &D0)+rho*ux**2*(1.0D0/3.0D0)-rho*uy**2*(1.0D0/6.0D0);
        
    NFqe(3) = mxxy*(-1.0D0/2.0D0)+rho*(1.0D0/9.0D0)+rho*uy*(1.0D0/2.0&
        &D0)-rho*ux**2*(1.0D0/6.0D0)+rho*uy**2*(1.0D0/3.0D0);
        
    NFqe(4) = mxyy*(1.0D0/2.0D0)+rho*(1.0D0/9.0D0)-rho*ux*(1.0D0/2.0D&
        &0)+rho*ux**2*(1.0D0/3.0D0)-rho*uy**2*(1.0D0/6.0D0);
        
    NFqe(5) = mxxy*(1.0D0/2.0D0)+rho*(1.0D0/9.0D0)-rho*uy*(1.0D0/2.0D&
        &0)-rho*ux**2*(1.0D0/6.0D0)+rho*uy**2*(1.0D0/3.0D0);
        
    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*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*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*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*ux*uy*(1.0D0&
        &/4.0D0);
      	
end subroutine NQEquilibrium