*----------------------------------------------------------------------
      PROGRAM BUOYANT CONVECTION
*----------------------------------------------------------------------
*     This code is for the computation of natural convection heated 
*     from a vertical plate at a constant temperature. 
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      Dimension T(101,101), TN(101,101), F(101,101), FN(101,101)
      Dimension W(101,101), WN(101,101), U(101,101), UN(101,101)
      Dimension X(101), Z(101)
      INTEGER(1) kai(2)
      character*30 str0G30
c     'Number of steps ?='
      NS=60000
c     'Number of grids ?='
      NGX=101
      NGZ=101
c     'Time step and Pranstl number'
      DELT=0.0005
      PR=0.7
c     'Width of the fluid layer ?='
      XL=16.0
c     'Height of the plate ?='
      ZL=1.0
*
*
      DELX=XL/(NGX-1)
      DELZ=ZL/(NGZ-1)
      X(1)=0.0
      Z(1)=0.0
*---------------Initial values-------------------------
      DO 15 I=1,NGX
      DO 15 K=1,NGZ
      T(I,K)=0.0
      F(I,K)=0.0
      W(I,K)=0.0
      TN(I,K)=0.0
      FN(I,K)=0.0
      UN(I,K)=0.0
   15 CONTINUE
*---------------Grid values----------------------------
      DO 17 I=1,NGX-1
      X(I+1)=X(I)+DELX
   17 CONTINUE
*
      DO 18 K=1,NGZ-1
      Z(K+1)=Z(K)+DELZ
   18 CONTINUE
*---------------Commence iteration------------------------ 
      DO ICOUNT =1 , NS
  
      DO 10 I=2,NGX-1
      DO 10 K=2,NGZ-1
      TN(I,K)=T(I,K)
     & -(DELT)*(F(I  ,K+1)-F(I  ,K-1))/(2.0*DELZ)
     &        *(T(I+1,K  )-T(I-1,K  ))/(2.0*DELX)
     & -(DELT)*W(I,K)
     &        *(T(I  ,K  )-T(I  ,K-1))/(DELZ)
     & +(DELT)*(T(I+1,K)-2.0*T(I,K)+T(I-1,K))/(DELX**2.0)
*
      WN(I,K)=W(I,K)
     & -(DELT)*(F(I  ,K+1)-F(I  ,K-1))/(2.0*DELZ)
     &        *(W(I+1,K  )-W(I-1,K  ))/(2.0*DELX)
     & -(DELT)*W(I,K)
     &        *(W(I  ,K  )-W(I  ,K-1))/(DELZ)
     & +(DELT)*PR*(W(I+1,K)-2.0*W(I,K)+W(I-1,K))/(DELX**2.0)
     & +(DELT)*PR*T(I,K)
   10 CONTINUE
*
      DO 11 I=2,NGX-1
      DO 11 K=1,NGZ
      FN(I,K)=FN(I-1,K)-(3.0*W(I-1,K)+6.0*W(I,K)-W(I+1,K))/8.0*DELX
   11 CONTINUE
*
      DO 20 I=2,NGX-1
      DO 20 K=2,NGZ-1
      T(I,K)=TN(I,K)
      F(I,K)=FN(I,K)
      W(I,K)=WN(I,K)
   20 CONTINUE
*---------------Boundary values-------------------------
*-----Left rigid wall-----------
      DO 40 K=1,NGZ
      T(1,K)=1.0
      F(1,K)=0.0
      W(1,K)=0.0
*-----Right free boundary-------
      T(NGX,K)=0.0
      F(NGX,K)=F(NGX-1,K)
      W(NGX,K)=0.0
   40 CONTINUE
*-----Bottom free boundary------
      DO 50 I=1,NGX
      T(I,1)=-T(I,2)
      F(I,1)=F(I,2)-(U(I,1)+U(I,2))/2.0*DELZ
      W(I,1)=0.0
*-----Top free boundary---------
      T(I,NGZ)=T(I,NGZ-1)
      F(I,NGZ)=F(I,NGZ-1)+(U(I,NGZ-1)+U(I,NGZ))/2.0*DELZ
      W(I,NGZ)=W(I,NGZ-1)
   50 CONTINUE
*-------------------------------------------------------
      DO 52 I=1,NGX
      U(I,1)=U(I,2)
   52 CONTINUE
      DO 55 I=1,NGX
      DO 55 K=2,NGZ-1
      U(I,K)=(F(I,K+1)-F(I,K-1))/(2.0*DELZ)
   55 CONTINUE
*--------------- Visualization -------------------------
      IF(MOD(ICOUNT,100).EQ.1) THEN
      WRITE(6,*) 'ICOUNT, T(30,30)=', ICOUNT, T(30,30)
      ENDIF
*
      END DO 
*-------------end iteration--------------------------
      kai(1) = 13!改行記号の設定
      kai(2) = 10
      
      open(99 ,file="Res.vtk" , form="unformatted" ,convert="BIG_ENDIAN"
     & ,STATUS="replace",access='stream' )
C     ヘッダを書き込む
C     文字列の書き込みができない。
      write(99) "# vtk DataFile Version 3.0" , kai
      write(99) "vtk output", kai
      write(99) "BINARY", kai
      write(99) "DATASET STRUCTURED_GRID", kai
      write(str0G30,"(I6,I6,I6)") NGX, NGZ , 1
      write(99) "DIMENSIONS " , str0G30, kai
      write(str0G30,"(I6)") NGX*NGZ
      write(99) "POINTS ", str0G30 , " float", kai
C     格子点を記録します。
      do K=1,NGZ
      do I=1,NGX !vtk出力ではループをiから回す必要がある。
      write(99) real(X(I)/20),real(Z(K)),real(0)
      end do
      end do
      
C     じゃんじゃん書き込みます
      write(str0G30,"(a , I8)") "POINT_DATA " , NGX*NGZ
      write(99) str0G30, kai
      write(99) "VECTORS Velocity float"  , kai
      
      do K=1,NGZ
      do I=1,NGX !vtk出力ではループをiから回す必要がある。
      write(99) real(U(I,K))
      write(99) real(W(I,K))
      write(99) real(0)
      end do
      end do
      
      write(99) "SCALARS T float"  , kai 
      write(99) "LOOKUP_TABLE default" , kai
      
      do K=1,NGZ
      do I=1,NGX !vtk出力ではループをiから回す必要がある。
      write(99) real(T(I,K))
      end do
      end do
      
      write(99) "SCALARS F float"  , kai 
      write(99) "LOOKUP_TABLE default" , kai
      
      do K=1,NGZ
      do I=1,NGX !vtk出力ではループをiから回す必要がある。
      write(99) real(F(I,K))
      end do
      end do
      
      close(99)
      STOP
      END
