      program NaturalConvection
      implicit none !自動変数宣言をoff にする。バグが減ります。
!     fortran 自然対流超簡易版
!     移流項: 対流型補間 2次精度中心差分
!     連続の式: HSMAC
!     密度変化モデル: ブシネスク近似
!     出力はvtk Binary出力
C     define constant value
      integer, parameter : : nx = 32 , ny = 64
C     define array variables
      double precision gU(0:(nx+2)*(ny+2)) , gV(0:(nx+2)*(ny+2))
     &,gF(0:(nx+2)*(ny+2)),gFN(0:(nx+2)*(ny+2)) 
     &, gUN(0:(nx+2)*(ny+2)),gVN(0:(nx+2)*(ny+2)) , gP(0:(nx+2)*(ny+2))
C     define parameter variables of double precision
      double precision cTime, dt, dx, dy, Lx, Ly, MaxTime ,Pr ,Ra
     &,vGx, vGy, eps, satSOR, outInterval , NextOutTime
C     define parameter variables of integer
      integer nOut 
      
      INTEGER(1) kai(2)
      character*30 str0G30
C     汎用変数
      double precision conv,diff, press ,thermal,divU,divUN
     &,dP
C     インデックス、カウンター変数
      integer iC,iE, iW, iU, iL, iEL, iUW, i, j,cOut,iEU,
     &pIter,i2,isRed
      character*260 fname
      
      kai(1) = 13
      kai(2) = 10
C     initialization
      MaxTime = 0.1   !計算時間
      eps = 0.001     !圧力修正判定条件
      vGx = 0         !重力ベクトルx成分
      vGy = -1        !重力ベクトルy成分
      Pr = 1          !粘性拡散/熱拡散
      Ra = 50000     !浮力*慣性力/(粘性*熱拡散)Ra数
      dt = 0.00001     !時間刻み
      Lx = 1          !x方向長さ
      Ly = 2          !y方向長さ
      dx = Lx/nx      !x方向刻み
      dy = Ly/ny      !y方向刻み
      nOut = 400       !総出力ファイル数
      cOut = 0
      satSOR = 1.7/((-2/dx/dx-2/dy/dy)*dt) !SOR法にかける加速係数 1.7推奨
      outInterval = MaxTime/DBLE(nOut)   !出力間隔
      NextOutTime = 0   !最後に出力した時間
      
 ! -------変数値を初期化する。好きにするといい-----
      do i=0,nx+1   !これでiが0からnx+1 まで変化するようにループが回る。
      do j=0,ny+1
      iC=i*(ny+2) + j
      gU(iC) = 0
      gV(iC) = 0
      gP(iC) = 0
      gF(iC) = 0
      end do
      end do
!--------------------------------------
C     main program loop
C     ジャンプ先 100がmainloop 990がプログラムエンド
  100 continue
      
C     時間積分
      do i=1,nx
      do j=1,ny
C     indexの初期化 1次元配列に慣れてくれ。アクセスの仕方で見え方は変わるんだよーー。
      iC  = (i+0)*(ny+2) + j+0
      iW  = (i-1)*(ny+2) + j+0
      iE  = (i+1)*(ny+2) + j+0
      iU  = (i+0)*(ny+2) + j+1
      iL  = (i+0)*(ny+2) + j-1
      iEL = (i+1)*(ny+2) + j-1
      iUW = (i-1)*(ny+2) + j+1
C     x方向速度
      conv = ((gU(iC) + gU(iE))/2*(gU(iC)+gU(iE))/2 
     &      - (gU(iC) + gU(iW))/2*(gU(iC)+gU(iW))/2)/dx + 
     &       ((gU(iU) + gU(iC))/2*(gV(iE)+gV(iC))/2 
     &      - (gU(iC) + gU(iL))/2*(gV(iEL)+gV(iL))/2)/dy
      diff = (gU(iE)-2*gU(iC) + gU(iW))/dx/dx 
     &     + (gU(iU)-2*gU(iC) + gU(iL))/dy/dy
      thermal = -(gF(iC) + gF(iE))/2 * vGx
      press  = -(gP(iE)-gP(iC))/dx
      
      gUN(iC) = gU(iC) + dt*(-conv + diff + press + Ra/Pr*thermal)
C     y方向速度
      conv = ((gV(iC) + gV(iE))/2*(gU(iC)+gU(iU))/2 
     &      - (gV(iC) + gV(iW))/2*(gU(iW)+gU(iUW))/2)/dx + 
     &       ((gV(iU) + gV(iC))/2*(gV(iU)+gV(iC))/2 
     &      - (gV(iC) + gV(iL))/2*(gV(iC)+gV(iL))/2)/dy
      diff = (gV(iE)-2*gV(iC) + gV(iW))/dx/dx 
     &     + (gV(iU)-2*gV(iC) + gV(iL))/dy/dy
      thermal = -(gF(iC) + gF(iU))/2 * vGy
      press  = -(gP(iU)-gP(iC))/dx
      
      gVN(iC) = gV(iC) + dt*(-conv + diff + press + Ra/Pr*thermal)
C     エネルギー
      conv = ((gF(iE) + gF(iC))/2*gU(iC) 
     &      - (gF(iC) + gF(iW))/2*gU(iW))/dx + 
     &       ((gF(iU) + gF(iC))/2*gV(iC) 
     &      - (gF(iC) + gF(iL))/2*gV(iL))/dy
      diff = (gF(iE)-2*gF(iC) + gF(iW))/dx/dx 
     &     + (gF(iU)-2*gF(iC) + gF(iL))/dy/dy
      gFN(iC) = gF(iC) + dt*(-conv + 1/Pr*diff)
      end do
      end do 
      
C     圧力修正HSMAC
      divU = eps*2
      pIter = 0
      do while(divU > eps)
!-----------------------------------------------
C     流入出0境界条件 自然対流やっている限りいじるな。
      do j=0,ny+1
      iC = 0*(ny+2) + j
      gUN(iC) = 0
      iC = nx*(ny+2) + j
      gUN(iC) = 0
      end do 
      
      do i=0,nx+1
      iC = i*(ny+2) + 0
      gVN(iC) = 0
      iC = i*(ny+2) + ny
      gVN(iC) = 0
      end do
!-----------------------------------------------
      divU = 0
      do isRed = 0 ,1
      do i2=0,nx/2-1
      do j=1,ny
      i = i2*2 + mod(j+isRed,2) + 1 
C     indexの初期化
      iC = (i+0)*(ny+2) + j+0
      iW = (i-1)*(ny+2) + j+0
      iL = (i+0)*(ny+2) + j-1
      divUN = (gUN(iC)-gUN(iW))/dx + (gVN(iC) - gVN(iL))/dy!速度発散値
      dP = divUN*satSOR           !SOR法
      gP(iC) = gP(iC) + dp
      gUN(iC) = gUN(iC) + dp*dt/dx!速度修正
      gUN(iW) = gUN(iW) - dp*dt/dx
      gVN(iC) = gVN(iC) + dp*dt/dy
      gVN(iL) = gVN(iL) - dp*dt/dy
      
      divUN = abs(divUN)
      if(divUN > divU) then
      divU = divUN
      end if 
      end do
      end do 
      end do
      pIter = pIter + 1
      if(mod(pIter,1000) == 0) then !pIterを1000で割った余りが0ならifを実行
      write(*,"(a ,I6,a, e10.3)") "HSMAC cycle = " , pIter,"div=" , divU
      if(pIter > 100000) goto 990
!     exit !ループを抜ける命令
      end if
      end do
      
!---- 全境界条件を適応する。こちらをいじることで壁の挙動が変わる。----
      do j=0,ny+1
      iC = 0*(ny+2) + j
      iE = 1*(ny+2) + j
      gUN(iC) = 0
      gVN(iC) = 0*2 - gVN(iE)  !すべりなし境界条件
      gFN(iC) = 0.0*2 + gFN(iE)!断熱境界条件
      
      iC = (nx+0)*(ny+2)+ j
      iE = (nx+1)*(ny+2) + j
      gUN(iC) = 0
      gVN(iE) = 0*2 - gVN(iC)
      gFN(iE) = 0.0*2 + gFN(iC)
      end do 
      
      do i=0,nx+1
      iC = i*(ny+2) + 0
      iU = i*(ny+2) + 1
      gUN(iC) = 0*2 - gUN(iU)
      gVN(iC) = 0
      gFN(iC) = 1 + gFN(iU) !等温加熱
     
      iC = i*(ny+2) + ny
      iU = i*(ny+2) + ny+1
      gUN(iU) = 0*2 - gUN(iC)
      gVN(iC) = 0
      gFN(iU) = -1 + gFN(iC)
      end do
!-----------------------------------------------------------
C     値の更新。本当はアドレス入れ替えを行いたいが・・・・・
      do iC=0,(nx+2)*(ny+2)-1
      gU(iC) = gUN(iC)
      gV(iC) = gVN(iC)
      gF(iC) = gFN(iC)
      end do
      
      cTime = cTime+dt   !時間更新
C     気が向いたら出力する VTK binary paraview で読み込めるー　です？
C--------出力部--------------------------------------
      if(cTime > NextOutTime) then
C     ファイル名を生成
      write(fname , "(a,i4.4,a)") "res\result" , cOut , ".vtk" 
      
!     unformatted でバイナリ出力意味し、
!     convert は複数Byte数値の並び順を指定します。通常x86系PCではLittle endianですがvtk binary formatはbig endianにしないと行けません。
!     append は追加書き込みをいみします。ヘッダー部の後方に連続でデータを書き込むためにしようします。
!     stream は指定しないと謎のマーカーをかきこまれます。
!     status replace で新規作成

      open(99 , file=fname , 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)") nx+1, ny+1 , 1
      write(99) "DIMENSIONS " , str0G30, kai
      write(str0G30,"(I6)") (nx+1)*(ny+1)
      write(99) "POINTS ", str0G30 , " float", kai
C     格子点を記録します。
      do j=0,ny
      do i=0,nx !vtk出力ではループをiから回す必要がある。
      write(99) real(i*dx),real(j*dy),real(0)
      end do
      end do
      
C     じゃんじゃん書き込みます
      write(str0G30,"(a , I8)") "POINT_DATA " , (nx+1)*(ny+1)
      write(99) str0G30, kai
      write(99) "VECTORS Velocity float"  , kai
      
      do j=0,ny
      do i=0,nx!vtk出力ではループをiから回す必要がある。
      iC  = (i+0)*(ny+2) + j+0
      iE  = (i+1)*(ny+2) + j+0
      iU  = (i+0)*(ny+2) + j+1
      write(99) real((gU(iC) + gU(iU))/2)
      write(99) real((gV(iC) + gV(iE))/2)
      write(99) real(0)
      end do
      end do
      
      write(99) "SCALARS P float"  , kai 
      write(99) "LOOKUP_TABLE default" , kai
      
      do j=0,ny
      do i=0,nx!vtk出力ではループをiから回す必要がある。
      iC  = (i+0)*(ny+2) + j+0
      iE  = (i+1)*(ny+2) + j+0
      iU  = (i+0)*(ny+2) + j+1
      iEU = (i+1)*(ny+2) + j+1
      write(99) real((gP(iC) + gP(iU) + gP(iE) + gP(iEU))/4)
      end do
      end do
      
      write(99) "SCALARS F float"  , kai
      write(99) "LOOKUP_TABLE default", kai
      
      do j=0,ny
      do i=0,nx!vtk出力ではループをiから回す必要がある。
      iC  = (i+0)*(ny+2) + j+0
      iE  = (i+1)*(ny+2) + j+0
      iU  = (i+0)*(ny+2) + j+1
      iEU = (i+1)*(ny+2) + j+1
      write(99) real((gF(iC) + gF(iU) + gF(iE) + gF(iEU))/4)
      end do
      end do
      
      
      write(*,"(I4,2x,e13.5,f9.5,a,I5)") cOut , cTime , divU , 
     & "HSMAC cycle=", pIter !fortran 書式でググるといろいろわかります。
      NextOutTime = NextOutTime + outInterval !次回の出力予定時間を更新します。
      cOut = cOut + 1
      end if
!----------出力完了------------------------------------------

C     計算を続けるかどうかはんていする。
      if(cTime <= MaxTime) goto 100  !続行。メインループの頭に戻る。
  
  990 end