program wave
      parameter( lmax=2000 )
      parameter( nmax=100 )
      parameter( nmay=100 )
      parameter( mmax=50 )
      parameter( mmay=50 )
      real x(0:nmax-1), y(0:nmay-1)
      real z(0:nmax-1,0:nmay-1)
      complex b(lmax,-mmax:mmax,-mmay:mmay)
      complex a(-mmax:mmax,-mmay:mmay)
      complex k1, k2, k3, k4
      real dx, dy, dt
      real xmin, xmax, ymin, ymax
      real pi, gamma

*-- 定数の定義
      pi=3.14159
      gamma=10.0

*-- 空間の定義 ---
      xmin=-pi
      xmax=pi
      ymin=-pi
      ymax=pi
      dx=(xmax-xmin)/(nmax-1)
      dy=(ymax-ymin)/(nmay-1)

*-- 時間刻み幅の設定 ---
      dt=1.0/(2.0*(mmax**2+mmay**2)*gamma)
*-- 空間座標の設定 ---
      do j=0,nmax-1
         x(j)=xmin+dx*j
      end do
      do k=0,nmay-1
         y(k)=ymin+dy*k
      end do

*-- 初期値の設定 ---
      do k=0,nmay-1
         do j=0,nmax-1
            z(j,k)=exp(-0.5*(x(j)**2+y(k)**2))
         end do
      end do

*-- 二重フーリエ変換 ---
      CALL DDFT(nmax-1,nmay-1,2*mmax,2*mmay,z(0:nmax-1,0:nmay-1),
     & a(-mmax:mmax,-mmay:mmay))

*-- スペクトルの時間変数への代入 ---
      do k=-mmay,mmay
         do j=-mmax,mmax
            b(1,j,k)=a(j,k)
         end do
      end do

*-- 常微分方程式の計算 ---
*-- ここでは, 4 次の RK で計算 ---
      do l=1,lmax-1
         do k=-mmay,mmay
            do j=-mmax,mmax
               k1=-dt*(j**2+k**2)*gamma*b(l,j,k)
               k2=-dt*(j**2+k**2)*gamma*(b(l,j,k)+0.5*k1*dt)
               k3=-dt*(j**2+k**2)*gamma*(b(l,j,k)+0.5*k2*dt)
               k4=-dt*(j**2+k**2)*gamma*(b(l,j,k)+k3*dt)
               b(l+1,j,k)=b(l,j,k)+(1.0/6.0)*(k1+2.0*k2+2.0*k3+k4)
            end do
         end do
      end do

*-- グラフ ----
      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
      READ (*,*) IWS

      CALL GROPN( IWS )
      do l=1,lmax
      if(mod(l-1,100).eq.0)then
      do k=-mmay,mmay
         do j=-mmax,mmax
            a(j,k)=b(l,j,k)
         end do
      end do
      CALL DIDFT(nmax-1,nmay-1,2*mmax,2*mmay,a(-mmax:mmax,-mmay:mmay)
     & ,z(0:nmax-1,0:nmay-1))
      CALL GRFRM

      CALL GRSWND( XMIN, XMAX, YMIN, YMAX )
      CALL GRSVPT(  0.2,  0.8,  0.2,  0.8 )
      CALL GRSTRN( 1 )
      CALL GRSTRF

      CALL USDAXS
      CALL UDCNTR( z, nmax, nmax, nmay )
      end if
      end do

      CALL GRCLS

      stop
      end