program wave
      parameter( lmax=400 )
      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)
      real dx, dy, dt
      real xmin, xmax, ymin, ymax
      real pi, gamma

*-- 出力ファイルをダンプ (保存) するための命令 ---
      CALL SWpSET( 'LDUMP', .TRUE. )
      CALL SWpSET( 'LWAIT', .FALSE. )
      CALL SWpSET( 'LWAIT1', .FALSE. )

*-- 定数の定義
      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*sqrt(real(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

*-- 常微分方程式の計算 ---
*-- 初期計算は初速度がゼロであることから, 1 ステップ先が 1 ステップ前
*-- と同じ値となることを用いている ---
      do k=-mmay,mmay
         do j=-mmax,mmax
            b(2,j,k)=b(1,j,k)
         end do
      end do

*-- ここでは, 2 次の リープフロッグ法で計算 ---
      do l=2,lmax-1
         do k=-mmay,mmay
            do j=-mmax,mmax
               b(l+1,j,k)=(2.0-(j**2+k**2)*(gamma**2)*(dt**2))*b(l,j,k)
     &                    -b(l-1,j,k)
            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,20).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