program explicit parameter( nmax=2000 ) parameter( mmax=10 ) 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/nmax *-- 初期条件の設定 --- do 10 i=2,mmax do 20 j=2,mmax u(1,i,j)=1.0 v(i,j)=u(1,i,j) 20 continue 10 continue do 30 i=1,mmax+1 u(1,i,1)=0.0 u(1,i,mmax+1)=0.0 v(i,1)=u(1,i,1) v(i,mmax+1)=u(1,i,mmax+1) 30 continue do 40 j=2,mmax u(1,1,j)=0.0 u(1,mmax+1,j)=0.0 v(1,j)=u(1,1,j) v(mmax+1,j)=u(1,mmax+1,j) 40 continue *-- グラフの出力装置の指定 (IWS=4 は GTK 指定) ---- WRITE(*,*) ' WORKSTATION ID (I) ? ;' CALL SGPWSN * IWS=4 READ (*,*) IWS *-- 出力装置の呼び出し CALL SGOPN( IWS ) *-- 時間ステップの計算 --- do 100 i=1,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 x(1)=xmin x(mmax+1)=xmax y(1)=ymin y(mmax+1)=ymax do 70 j=2,mmax x(j)=xmin+dx*(j-1) y(j)=ymin+dy*(j-1) do 80 k=2,mmax u(i+1,j,k)=u(i,j,k)+(dt/(dx**2)) & *(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(i.lt.20)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, 2.0 ) 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, 0.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, 0.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, 0.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