program explicit
      parameter( nmax=2000 )
      parameter( mmax=100 )
      parameter( kmax=56, vmin=-1, vmax=1 )
      real u(nmax+1,mmax+1,mmax+1), v(mmax+1,mmax+1)
      real x(mmax+1), y(mmax+1)
      real pi(2,kmax+1)
*-- 出力ファイルをダンプ (保存) するための命令 ---
      CALL SWpSET( 'LDUMP', .TRUE. )
      CALL SWpSET( 'LWAIT', .FALSE. )
      CALL SWpSET( 'LWAIT1', .FALSE. )

*-- 空間領域の設定 ---
      xmin=0.0
      xmax=2.0
      ymin=0.0
      ymax=2.0
*-- 空間きざみ幅の設定 ---
*-- ここでは, 空間分解能をx,y方向で同じ間隔とした ---
      dx=(xmax-xmin)/mmax
      dy=(ymax-ymin)/mmax        ! これは y 座標を与えるのみの定義
*-- 時間ステップの設定 ---
      dt=1.0/nmax
      alpha=(dt/dx)**2
*-- 位置座標の設定 ---
      do 210 j=1,mmax+1
         x(j)=xmin+dx*(j-1)
         y(j)=ymin+dy*(j-1)
  210 continue
*-- 初期条件の設定 ---
      do 10 i=1,mmax+1
      do 20 j=1,mmax+1
         u(1,i,j)=exp(-4.0*((x(i)-1.0)**2+(y(j)-1.0)**2))
         v(i,j)=u(1,i,j)
  20  continue
  10  continue

*-- 初速度の設定 ---
      do 11 i=1,mmax+1
      do 21 j=1,mmax+1
         u(2,i,j)=u(1,i,j)
         v(i,j)=u(2,i,j)
  21  continue
  11  continue

*-- グラフの出力装置の指定 (IWS=4 は GTK 指定) ----
      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
*      IWS=4
      READ (*,*) IWS

*-- 出力装置の呼び出し
      CALL GROPN( IWS )

*-- トーン塗り分けの設定 1 ---
*-- UESTLV は 100 以下しか呼び出せないので, 苦肉の策 ---
         dp=(vmax-vmin)/kmax
         do 110 k=1,kmax
            TLEV1=vmin+(k-1)*dp
            TLEV2=TLEV1+dp
            IPAT=(29+k)*1000+999
            CALL UESTLV( TLEV1, TLEV2, IPAT )
  110  continue

*-- 時間ステップの計算 ---
      do 100 i=2,nmax
         do 70 j=2,mmax
            do 80 k=2,mmax
               u(i+1,j,k)=2.0*u(i,j,k)-u(i-1,j,k)+200.0*alpha
     &                  *(u(i,j+1,k)+u(i,j-1,k)+u(i,j,k+1)+u(i,j,k-1)
     &                  -4.0*u(i,j,k))
  80  continue
  70  continue
*-- 境界条件の設定 ---
      do 50 j=2,mmax
         u(i+1,j,1)=u(i+1,j,2)
         u(i+1,j,mmax+1)=u(i+1,j,mmax)
  50  continue
      do 60 j=2,mmax
         u(i+1,1,j)=u(i+1,2,j)
         u(i+1,mmax+1,j)=u(i+1,mmax,j)
  60  continue
*-- 4隅の境界条件の設定 ---
      u(i+1,1,1)=0.5*(u(i+1,2,1)+u(i+1,1,2))
      u(i+1,mmax+1,1)=0.5*(u(i+1,mmax,1)+u(i+1,mmax+1,2))
      u(i+1,1,mmax+1)=0.5*(u(i+1,1,mmax)+u(i+1,2,mmax+1))
      u(i+1,mmax+1,mmax+1)=0.5*(u(i+1,mmax,mmax+1)+u(i+1,mmax+1,mmax))
      do 55 k=1,mmax+1
         do 65 j=1,mmax+1
            v(j,k)=u(i+1,j,k)
  65  continue
  55  continue
*-- t=100 ごとで描画する設定 ---
      if(mod(i-1,10).eq.0)then
         if(i.lt.200)then
         CALL GRFRM
         CALL GRSWND( xmin, xmax, ymin, ymax )
         CALL GRSVPT(  0.2,  0.8,  0.2,  0.8 )
         CALL USPFIT
         CALL GRSTRF

         CALL USSTTL( 'X-AXIS', '', 'Y-AXIS', '' )
         CALL UETONE( v, mmax+1, mmax+1, mmax+1 )  ! 格子点の個数を与える
         CALL USDAXS
         CALL UDCNTR( v, mmax+1, mmax+1, mmax+1 )
*-- トーンバー（凡例）の設定 2 ---
         CALL GRSWND( 0.0, 1.0, vmin, vmax )
         CALL GRSVPT( 0.85, 0.9, 0.2, 0.8 )
         CALL GRSTRN( 1 )
         CALL GRSTRF

         do 120 k=1,kmax+1
            pi(1,k)=vmin+(k-1)*dp
            pi(2,k)=vmin+(k-1)*dp
  120  continue

         CALL UETONE( pi, 2, 2, kmax+1 )

         CALL SLPVPR( 3 )            ! これは GRSTRN より後ろで宣言される
         CALL UZLSET( 'LABELYR', .TRUE. )
*         CALL UZFACT( 0.8 )
         CALL UYSFMT( '(F4.1)' )
         CALL UYAXDV( 'R', 0.1, 0.2 )
         end if
      end if
  100  continue
  200  CALL GRCLS
      stop
      end