program explicit parameter( nmax=200 ) parameter( mmax=100 ) 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) real alpha, dx, dy, dt *-- 出力ファイルをダンプ (保存) するための命令 --- 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 alpha=(dt/dx) *-- 初期条件の設定 --- 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(-2.0*((x(i)-0.5)**2+(y(j)-0.5)**2)) v(i,j)=u(1,i,j) 30 continue 20 continue *-- グラフの出力装置の指定 (IWS=4 は GTK 指定) ---- WRITE(*,*) ' WORKSTATION ID (I) ? ;' CALL SGPWSN * IWS=4 READ (*,*) IWS *-- 出力装置の呼び出し CALL SGOPN( IWS ) *-- 時間ステップの計算 --- do 100 i=1,nmax do 70 j=2,mmax do 80 k=2,mmax u(i+1,j,k)=u(i,j,k)-u(i,j,k)*alpha & *(2.0*u(i,j,k)-u(i,j-1,k)-u(i,j,k-1)) 80 continue 70 continue *-- 境界条件の設定 --- do 50 j=2,mmax u(i+1,j,1)=u(i+1,j,mmax) u(i+1,j,mmax+1)=u(i+1,j,2) 50 continue do 60 j=2,mmax u(i+1,1,j)=u(i+1,mmax,j) u(i+1,mmax+1,j)=u(i+1,2,j) 60 continue *-- 4隅の境界条件の設定 --- u(i+1,1,1)=0.5*(u(i+1,mmax,1)+u(i+1,1,mmax)) u(i+1,1,mmax+1)=0.5*(u(i+1,mmax,mmax+1)+u(i+1,1,2)) u(i+1,mmax+1,1)=0.5*(u(i+1,2,1)+u(i+1,mmax+1,mmax)) u(i+1,mmax+1,mmax+1)=0.5*(u(i+1,mmax+1,2)+u(i+1,2,mmax+1)) do 65 k=1,mmax+1 do 55 j=1,mmax+1 v(j,k)=u(i+1,j,k) 55 continue 65 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, 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 CALL SGCLS stop end