Square Lid driven not working

Hi Friends,
I am a beginner with LBM. I generated a code in D2Q9 my code is working fine for channel flow problems and I am getting good results for drag coefficients. But when I use the same code for square lid driven cavity in 2-D with north wall moving my results diverge.
Re=1000
Uw=0.1lu/ts
grid = (100,100)

My code for moving boundary is as given below:

  INTEGER  Xmax,Ymax
  REAL*8   Uw,density,f(0:8,1001,1001),ftemp(0:8,1001,1001)

c
c local variables
INTEGER x,y,i
REAL*8 Dxy

c Bounceback at Top boundary
y = Ymax

    DO x = 2, Xmax-1

        Dxy = (ftemp(0,x,y)+ftemp(1,x,y)+ftemp(3,x,y)+
 &  2*(ftemp(2,x,y)+ftemp(6,x,y)+ftemp(5,x,y)))

	   Dxy=Dxy*Uw	
      ftemp(4,x,y)=ftemp(2,x,y)


      ftemp(7,x,y)=ftemp(5,x,y) -Dxy/2. +
 &         .5*(ftemp(1,x,y)-ftemp(3,x,y))


      ftemp(8,x,y)=ftemp(6,x,y) +Dxy/2. +
 &              .5*(ftemp(3,x,y)-ftemp(1,x,y))

    END DO
    END

The other bounceback conditions that I am applying on the other walls is given below:

DO x = 1, Xmax
DO y = 1, Ymax
IF (wall(x,y)) THEN

        f(1,x,y) = ftemp(3,x,y)
        f(2,x,y) = ftemp(4,x,y)
        f(3,x,y) = ftemp(1,x,y)
        f(4,x,y) = ftemp(2,x,y)
        f(5,x,y) = ftemp(7,x,y)
        f(6,x,y) = ftemp(8,x,y)
        f(7,x,y) = ftemp(5,x,y)
        f(8,x,y) = ftemp(6,x,y)

c
END IF
END DO
END DO
I am using full bounceback conditions.
After some time in lid driven cavity the code diverges and gives very large negative values of densities. I am puzzled why it is happening. Kindly help me to sort it out. I am going through various validation processes. I would be grateful to you all.

Regards
Prateek