对重力长波方程组，用a网格，使用蛙跃格式进行编程。我运行后发现计算结果很不稳定，增长很厉害

请问各位有什么想法？

program test

 implicit none
   real, parameter  :: g=9.8
   real             :: delt_x = 1000

   real             :: delt_t = 50
   integer        :: n = 100

   integer        :: j = 100 

   real             :: depth = 36

   integer         :: i, k
   real, allocatable, dimension(:,:) :: u , zeta

!
   allocate( u(n,j), zeta(n,j) )

!初始值

   u(1,:)=0
   u(2,:)=0
 
   zeta(1,:)=1
   zeta(2,:)=1

！加上强迫源
   do i=1,n
       zeta(i,1)=sin(0.5*delt_t*(i-1))
   end do

!  开始计算

   do i=2,n-1
     do k=2, j-1
         u (i+1,k)=u(i-1,k)-g*delt_t/delt_x*(zeta(i,k+1)-zeta(i,k-1))
         zeta(i+1,k)=zeta(i-1,k)-depth*delt_t/delt_x*(u(i,k+1)-u(i,k-1))
     end do
   end do
！输出结果
   open (2,file='result.dat', status='unknown',action='write')
   write(2,*) 'u '
    do i=1,n
     write(2,'(100f15.3)') u(i,:)
    end do
   write(2,*)'zeta'
    do i=1,n
     write(2,'(100f15.3)')zeta(i,:)
    end do
   close (2)
   !
   deallocate (u , zeta)

end program test
   do i=2,n-1
     do k=2, j-1
  u (i+1,k)=u(i-1,k)-g*delt_t/delt_x*(zeta(i,k+1)-zeta(i,k-1))
  zeta(i+1,k)=zeta(i-1,k)-depth*delt_t/delt_x*(u(i,k+1)-u(i,k-1))
  end do
   end do
   !!
   !
   !----output  data
   open (2,file='result.dat', status='unknown',action='write')
   write(2,*) 'u '
    do i=1,n
     write(2,'(100f15.3)') u(i,:)
    end do
   write(2,*)'zeta'
    do i=1,n
     write(2,'(100f15.3)')zeta(i,:)
    end do
   close (2)
   !
   deallocate (u , zeta)   do i=2,n-1
     do k=2, j-1
  u (i+1,k)=u(i-1,k)-g*delt_t/delt_x*(zeta(i,k+1)-zeta(i,k-1))
  zeta(i+1,k)=zeta(i-1,k)-depth*delt_t/delt_x*(u(i,k+1)-u(i,k-1))
  end do
   end do
   !!
   !
   !----output  data
   open (2,file='result.dat', status='unknown',action='write')
   write(2,*) 'u '
    do i=1,n
     write(2,'(100f15.3)') u(i,:)
    end do
   write(2,*)'zeta'
    do i=1,n
     write(2,'(100f15.3)')zeta(i,:)
    end do
   close (2)
   !
   deallocate (u , zeta)   do i=2,n-1
     do k=2, j-1
  u (i+1,k)=u(i-1,k)-g*delt_t/delt_x*(zeta(i,k+1)-zeta(i,k-1))
  zeta(i+1,k)=zeta(i-1,k)-depth*delt_t/delt_x*(u(i,k+1)-u(i,k-1))
  end do
   end do
   !!
   !
   !----output  data
   open (2,file='result.dat', status='unknown',action='write')
   write(2,*) 'u '
    do i=1,n
     write(2,'(100f15.3)') u(i,:)
    end do
   write(2,*)'zeta'
    do i=1,n
     write(2,'(100f15.3)')zeta(i,:)
    end do
   close (2)
   !
   deallocate (u , zeta)