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
      complex i
      real dx, dy, dt
      real xmin, xmax, ymin, ymax
      real pi, gammax, gammay
      real sp

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

*-- 定数の定義
      i=(0.0,1.0)
      pi=3.14159
      gammax=10.0
      gammay=5.0
*-- 移流速度の最大値の検索 ---
      if(gammax.gt.gammay)then
         sp=gammax
      else
         sp=gammay
      end if

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

*-- 時間刻み幅の設定 ---
      dt=1.0/(2.0*(mmax+mmay)*sp)
*-- 空間座標の設定 ---
      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

*-- 常微分方程式の計算 ---
*-- 純移流方程式であるため, 初期計算は0.5ステップでオイラー法を使用 ---
      do k=-mmay,mmay
         do j=-mmax,mmax
            b(2,j,k)=(1.0-i*dt*(gammax*j+gammay*k)
     &               *(1.0-0.5*i*(gammax*j+gammay*k)*dt))*b(1,j,k)
         end do
      end do
*-- ここでは, リープフロッグ法で計算 ---
      do l=2,lmax-1
         do k=-mmay,mmay
            do j=-mmax,mmax
               b(l+1,j,k)=b(l-1,j,k)-i*(gammax*j+gammay*k)
     &                    *2.0*dt*b(l,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,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