Lattice BGK Model for Simulating NPDEs

Big Hi for you all,

I write a FORTRAN code for solving a one-dimensional nonlinear partial differential equation:
ut+ ?uux +?unux ? ?uxx + ?uxxx = F(u), like Burgers’ eq. and so on based on papers like this:

http://www.scichina.com:8083/sciGe/EN/article/downloadArticleFile.do?attachType=PDF&id=413111

! Calculate Initial Distribution Functions for 2000 nodes ( U(Initial)=U(Exact),F(Initial)=Feq(Initial) )
! Feq & H (Amending Function) are defined in paper
DO I=1,2000
FEQ(0,I)=U(I)
FEQ(1,I)=(-1/6.)LAMBDAU(I)
FEQ(2,I)=(1/6.)LAMBDAU(I)
FEQ(3,I)=(1/12.)LAMBDAU(I)
FEQ(4,I)=(-1/12.)LAMBDAU(I)
F(0,I)=FEQ(0,I)
F(1,I)=FEQ(1,I)
F(2,I)=FEQ(2,I)
F(3,I)=FEQ(3,I)
F(4,I)=FEQ(4,I)
END DO
! use Periodic Boundary Condition (This problem has periodic boundaries)
DO I=1,2000
IR1=MOD(I,2000)+1
IR2=MOD(I,2000)+2
IL1=2000-MOD(2001+1-I,2000)
IL2=2000-MOD(2002+1-I,2000)
F(0,I)=F(0,I)
F(1,IR1)=F(1,I)
F(2,IL1)=F(2,I)
F(3,IR2)=F(3,I)
F(4,IL2)=F(4,I)
END DO
! Start Iteration J (Marching in Time & Updating Feq, H, …)
DO J=1,500
DO I=1,2000
U(I)=F(0,I)+F(1,I)+F(2,I)+F(3,I)+F(4,I)
! DEFINE FEQ FROM U
FEQ(0,I)=U(I)
FEQ(1,I)=(-1/6.)LAMBDAU(I)
FEQ(2,I)=(1/6.)LAMBDAU(I)
FEQ(3,I)=(1/12.)LAMBDAU(I)
FEQ(4,I)=(-1/12.)LAMBDAU(I)
! Define H FROM U
H(0,I)=(6*(U(I)2))-(ETAU(I))
H(1,I)=((1./2)ETAU(I))+((((1./2)LAMBDA1-4))(U(I)2)+((1./2)LAMBDA2(U(I)(N+1))))
H(1,I)=((1./2)ETAU(I))-((((1./2)LAMBDA1-4))(U(I)2)-((1./2)LAMBDA2(U(I)(N+1))))
H(3,I)=(U(I)**2)
H(4,I)=(U(I)**2)
END DO
! Streaming & Collision (lattice Boltzmann Equation with amending Function)
DO I=1,2000
IR1=MOD(I,2000)+1
IR2=MOD(I,2000)+2
IL1=2000-MOD(2001+1-I,2000)
IL2=2000-MOD(2002+1-I,2000)
F(0,I)=F(0,I)+((-(1/TAU))
(F(0,I)-FEQ(0,I)))+((DELTA_T
2)H(0,I))
F(1,IR1)=F(1,I)+((-(1/TAU))
(F(1,I)-FEQ(1,I)))+((DELTA_T2)H(1,I))
F(2,IL1)=F(2,I)+((-(1/TAU))
(F(2,I)-FEQ(2,I)))+((DELTA_T
2)H(2,I))
F(3,IR2)=F(3,I)+((-(1/TAU))
(F(3,I)-FEQ(3,I)))+((DELTA_T2)H(3,I))
F(4,IL2)=F(4,I)+((-(1/TAU))
(F(4,I)-FEQ(4,I)))+((DELTA_T
2)*H(4,I))
END DO

IF (MOD(j,50)==0)THEN
   WRITE(100,*)"Zone"
 DO I=1,2000
    WRITE(100,*)I,U(I)
 END DO
 END IF

END DO

I appreciate if you can help me to improove the code and get good results.