      program explicit
      parameter( nmax=500 )
      parameter( mmax=50 )
      parameter( VXMIN=-0.4, VXMAX=0.4, VYMIN=-0.3, VYMAX=0.3 )
      parameter( VZMIN=-0.2, VZMAX=0.2 )
      real u(nmax+1,mmax+1,mmax+1), v(mmax+1,mmax+1)
      real x(mmax+1), y(mmax+1)
      real XP(3), YP(3), ZP(3)
*-- 出力ファイルをダンプ (保存) するための命令 ---
      CALL SWpSET( 'LDUMP', .TRUE. )
      CALL SWpSET( 'LWAIT', .FALSE. )
      CALL SWpSET( 'LWAIT1', .FALSE. )

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

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

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

*-- 時間ステップの計算 ---
      do 100 i=2,nmax
*-- 境界条件の設定 ---
      do 50 j=1,mmax+1
         u(i+1,j,1)=0.0
         u(i+1,j,mmax+1)=0.0
         v(j,1)=u(i+1,j,1)
         v(j,mmax+1)=u(i+1,j,mmax+1)
  50  continue
      do 60 j=2,mmax
         u(i+1,1,j)=0.0
         u(i+1,mmax+1,j)=0.0
         v(1,j)=u(i+1,1,j)
         v(mmax+1,j)=u(i+1,mmax+1,j)
  60  continue
         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)+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))
               v(j,k)=u(i+1,j,k)
  80  continue
  70  continue
*-- t=100 ごとで描画する設定 ---
      if(mod(i,10).eq.0)then
*-- X-Y PLANE ----
      CALL SGFRM
      CALL SGSWND(  xmin,  xmax,  ymin,  ymax )
      CALL SGSVPT( VXMIN, VXMAX, VYMIN, VYMAX )
      CALL SGSTRN( 1 )
      CALL SGSTRF

      CALL SCSPLN( 1, 2, VZMIN)
      CALL SCSEYE( -1.1, -1.1, 1.2 )
      CALL SCSOBJ(  0.0,  0.0, 0.0 )
      CALL SCSPRJ
      CALL UXAXDV( 'B', 0.05, 0.1 )
      CALL UXAXDV( 'T', 0.05, 0.1 )
      CALL UXSTTL( 'B', 'X-AXIS', 0.5 )
      CALL UYAXDV( 'L', 0.05, 0.1 )
      CALL UYAXDV( 'R', 0.05, 0.1 )
      CALL UYSTTL( 'L', 'LATITUDE', 0. )

*-- X-Z PLANE ----

      CALL SGSWND(  xmin,  xmax,  -1.0,  1.0 )
      CALL SGSVPT( VXMIN, VXMAX, VZMIN, VZMAX )
      CALL SGSTRN( 1 )
      CALL SGSTRF

      CALL SCSPLN( 1, 3, VYMAX)
      CALL SCSPRJ

      CALL UZINIT
      CALL UXAXDV( 'T', 0.05, 0.1 )
      CALL UZLSET( 'LABELXB', .FALSE. )
      CALL UXAXDV( 'B', 0.05, 0.1 )
      CALL UYAXDV( 'L', 0.2, 1.0 )
      CALL UYAXDV( 'R', 0.2, 1.0 )
      CALL UYSTTL( 'L', 'Y-AXIS', 0. )

*-- Y-Z PLANE ----

      CALL SGSWND(  ymin,  ymax,  -1.0,  1.0 )
      CALL SGSVPT( VYMIN, VYMAX, VZMIN, VZMAX )
      CALL SGSTRN( 1 )
      CALL SGSTRF

      CALL SCSPLN( 2, 3, VXMAX)
      CALL SCSPRJ

      CALL UZINIT
      CALL UXAXDV( 'T', 0.05, 0.1 )
      CALL UXAXDV( 'B', 0.05, 0.1 )
      CALL UZLSET( 'LABELYL', .FALSE. )
      CALL UYAXDV( 'L', 0.2, 1.0 )
      CALL UYAXDV( 'R', 0.2, 1.0 )

*---------------- 3-D ------------------

      CALL SCSVPT(VXMIN, VXMAX, VYMIN, VYMAX, VZMIN, VZMAX)
      CALL SCSWND( xmin,  xmax,  ymin,  ymax,  -1.0,  1.0)
      CALL SCSLOG(.FALSE., .FALSE., .FALSE.)
      CALL SCSTRN(1)
      CALL SCSTRF

      CALL SZL3OP(1)
C      CALL SZT3OP(2999,4999)

      DO 140 J=mmax, 1, -1
      DO 140 IM=mmax, 1, -1
        XP(1) = x(IM)
        YP(1) = y(J)
        ZP(1) = v(IM,J)
        XP(2) = x(IM)
        YP(2) = y(J+1)
        ZP(2) = v(IM,J+1)
        XP(3) = x(IM+1)
        YP(3) = y(J+1)
        ZP(3) = v(IM+1,J+1)

C        CALL SZT3ZU(XP, YP, ZP)
        CALL SZL3ZU(3, XP, YP, ZP)

        XP(1) = x(IM)
        YP(1) = y(J)
        ZP(1) = v(IM+1,J+1)
        XP(2) = x(IM+1)
        YP(2) = y(J)
        ZP(2) = v(IM+1,J)
        XP(3) = x(IM)
        YP(3) = y(J)
        ZP(3) = v(IM,J)

C        CALL SZT3ZU(XP, YP, ZP)
        CALL SZL3ZU(3, XP, YP, ZP)

  140 CONTINUE

      CALL SZL3CL
C      CALL SZT3CL
      end if
  100  continue
  111  CALL SGCLS
      stop
      end