      program Boussinesq
      implicit NONE
c originally produced by Mr. Nanamichi
c 2015.12.18 modified by T. IWAYAMA
c 2015.12.30 時間発展のバグを修正
c 2015.12.31 境界条件のバグを修正 Re=100ではKarman渦発生せず.
c 2016.01.01 omegaの境界条件を再考
c 2016.02.01 収束判定(Mizushima and Takemoto)
c 2016.11.05-11.17 2D turbulence
c 2016.11.17- Boussinesq fluids, convection
c 2016.11.18 Poisson方程式の解法 : SOR法完成 絶対誤差で収束性判定
c 2016.11.19 Dennou Sapmleを再現
c http://www.gfd-dennou.org/library/spmodel/2d-channel-esc/boussinesq/sample-tt/f90/bsncnv-tt-1.f90
c 201701.20 Nusselt number
C***********● 変数の定義 ●***********
      integer i, j, t, t1
      integer Nx, Ny                            !x, y方向の格子点数
      integer t_max, t1_max                     !
      real*8  Lx, Ly                            !計算領域のx, y方向の長さ
      real*8  dx, dy                            !格子点間隔
      real*8  dt                                !時間差分間隔
      real*8  Nu, Pr, Ra                        !Nusselt数, Prandtle数, Rayleigh数, 温度減率
      real*8  pi                                !円周率
      real*8  C1, C2, UNITY, ZERO               !AB法の定数
C***********● 定数の定義 ●***********
      parameter (Nx=128, Ny=64)                  !格子点数
      parameter (C1=1.5D+0, C2=0.5D+0)
      parameter (UNITY=1.D+0, ZERO=0.D+0)
C***********● 配列の定義 ●***********
      real*8 x(0:Nx+1), y(0:Ny+1)                 !x,y座標
      real*8 psi(0:Nx+1,0:Ny+1)                   !流れ関数の配列
      real*8 omg(0:Nx+1,0:Ny+1)                   !渦度の配列
      real*8 tmp(0:Nx+1,0:Ny+1)                   !温度の配列
      real*8 new_domg(1:Nx,1:Ny)      !AB法に使われるk1ステップ目のomgの時間変化量
      real*8 old_domg(1:Nx,1:Ny)      !AB法に使われる(k1-1)ステップ目のomgの時間変化量
      real*8 new_dtmp(1:Nx,1:Ny)      !AB法に使われるk1ステップ目のtmpの時間変化量
      real*8 old_dtmp(1:Nx,1:Ny)      !AB法に使われる(k1-1)ステップ目のtmpの時間変化量
c
      real*8 work(0:Nx+1,0:Ny+1)      !作業領域
      real*8 ENRG, ENST, time         !エネルギー, エンストロフィー, 

C************● 変数の値 ●************
      pi=dble(acos(-1.D+0))           ! piの値

      Lx=dble(2)                      !計算領域の縦，横の長さ
      Ly=dble(1)                      !計算領域の縦，横の長さ

      dx=Lx/dble(Nx)                  !格子点間隔(dx)
      dy=Ly/dble(Ny+1)                !格子点間隔(dy)
c
      t_max=40                         
      t1_max=2500
      dt=1.D-5                        !時間差分間隔

      Pr=1.0D+0
      Ra=4.D+3
   
      open(26, file='results.dat', status='unknown')
C**************************************
      do 10 i=0, Nx+1                             !x座標の設定
           x(i)=dx*dble(i)
   10 continue                                    
C
      do 20 j=0, Ny+1                             !y座標の設定
           y(j)=dy*dble(j)
   20 continue                                    
C-initial condition ******************************************
      t=0
      time=dble(0)
      call initial_condition(psi, omg, tmp, Nx, Ny, Lx, Ly)
      call energy(psi, omg, Nx, Ny, ENRG, ENST)
      call Nusselt(tmp, Nx, Ny, dy, Nu)

      write(26,*) time, ENRG, ENST, Nu
      call output_data(psi, omg, tmp, t, Nx, Ny, dx, dy)
C*******● Temporal evolution (Euler method)●*******
      do j=1, Ny
        do i=1, Nx
           old_domg(i,j)=dble(0)
           old_dtmp(i,j)=dble(0)
        end do
      end do

      call Governing_Eqns(psi,omg,tmp,new_domg,new_dtmp,
     &     Nx,Ny,dx,dy,Pr,Ra)
      call Next(omg, new_domg, old_domg, Nx, Ny, UNITY, ZERO, dt)
      call Next(tmp, new_dtmp, old_dtmp, Nx, Ny, UNITY, ZERO, dt)
      call bc_channel(omg, Nx, Ny)
      call bc_channel(tmp, Nx, Ny)
C***● Solve the Poisson equation ●***
      call poisson(psi, omg, Nx, Ny, dx, dy)
C*******● Temporal evolution (Adams-Bashforth method)●*******
      do t=1, t_max
       do t1=1, t1_max
      time=time+dt 
      call Governing_Eqns(psi,omg,tmp,new_domg,new_dtmp,
     &     Nx,Ny,dx,dy,Pr,Ra)
      call Next(omg, new_domg, old_domg, Nx, Ny, C1, C2, dt)
      call Next(tmp, new_dtmp, old_dtmp, Nx, Ny, C1, C2, dt)
      call bc_channel(omg, Nx, Ny)
      call bc_channel(tmp, Nx, Ny)
C***● Solve the Poisson equation ●***
      call poisson(psi, omg, Nx, Ny, dx, dy)
C**************************************
       end do

       call output_data(psi, omg, tmp, t, Nx, Ny, dx, dy)
       call energy(psi, omg, Nx, Ny, ENRG, ENST)
       call Nusselt(tmp, Nx, Ny, dy, Nu)

       write(26,*) real(t*t1)*dt, ENRG, ENST, Nu

C
      end do
C
      close(26)
      Stop
      End
C------------------------------------------------------------------------------
      subroutine bc_channel(Gdata, Nx, Ny)
c doubly periodic boundary condition
      integer i, j, k
      integer Nx, Ny                
c      real*8  Lx, Ly                         
      real*8  dx, dy                       
      real*8  dt                           
C***********● 配列の定義 ●***********
      real*8 Gdata(0:Nx+1,0:Ny+1)         
      real*8 x(0:Nx+1), y(0:Ny+1)       
C************● 変数の値 ●************
c      dx=Lx/dble(Nx-1)                    
c      dy=Ly/dble(Ny-1)                    
C**************************************
C
c      do i=0, Nx+1
c           x(j)=dble(i)*dx
c      end do
c
c      do j=0, Ny+1
c           y(j)=dble(j)*dy
c      end do

c at Wall
      do i=1,Nx
        Gdata(i,0)=dble(0)
        Gdata(i,Ny+1)=dble(0)
      end do
C cyclic
      do j=1,Ny
        Gdata(0,j)=Gdata(Nx,j)
        Gdata(Nx+1,j)=Gdata(1,j)
      end do
C
        Gdata(0,0)=dble(0)
        Gdata(Nx+1,0)=dble(0)
        Gdata(0,Ny+1)=dble(0)
        Gdata(Nx+1,Ny+1)=dble(0)

      Return
      End
C------------------------------------------------------------------------------
      subroutine initial_condition(psi, omg, tmp, Nx, Ny, Lx, Ly)
c     初期値を指定するsurboutine
C***********● 変数の定義 ●***********
      integer i, j, k
      integer Nx, Ny                
      real*8  Lx, Ly                         
      real*8  dx, dy                       
C***********● 配列の定義 ●***********
      real*8 psi(0:Nx+1,0:Ny+1)         
      real*8 omg(0:Nx+1,0:Ny+1)       
      real*8 tmp(0:Nx+1,0:Ny+1)       
      real*8 work(1:Nx,1:Ny)       
      real*8 x(0:Nx+1), y(0:Ny+1)
c      
      real*8 ENRG, ENST, pi
C************● 変数の値 ●************
      dx=Lx/dble(Nx)                  
      dy=Ly/dble(Ny+1)
      pi=dble(acos(-1.))                  
C**************************************
      do i=0, Nx+1
         x(i)=dble(i)*dx
      end do
      do j=0, Ny+1
         y(j)=dble(j)*dy
      end do

      do j=0,Ny+1
        do i=0,Nx+1
          psi(i,j)=dble(0)
          omg(i,j)=dble(0)
          tmp(i,j)=dble(0)
        end do
      end do

c      open(20, file='q_0000.dat', status='old')
c      read(20,300) ((x(i), y(j), omg(i,j), i=1,Nx-1), j=1,Ny-1)
c      close(20)

c      open(20, file='p_0000.dat', status='old')
c      read(20,300) ((x(i), y(j), psi(i,j), i=1,Nx), j=1,Ny)
c      close(20)

c      do j=0,Ny+1
c        do i=0,Nx+1
c          psi(i,j)=dble(1)*dcos(x(i)*dble(10))
c          omg(i,j)=-dble(10)**2*
c     &             dcos(x(i)*dble(10))
c        end do
c      end do

c      call bc_channel(psi, Nx, Ny)

c      call laplacian(psi, Nx, Ny, dx, dy, work)

c      do j=1,Ny
c       do i=1, Nx
c        omg(i,j)=work(i,j)
c       end do
c      end do
c      call bc_channel(omg, Nx, Ny)

c normalization of initial condition
c      call energy(psi, omg, Nx, Ny, ENRG, ENST)
c      do j=1,Ny
c        do i=1,Nx
c          psi(i,j)=psi(i,j)/dsqrt(ENRG)*dsqrt(.5D+0)
c        end do
c      end do

c      call bc_channel(psi, Nx, Ny)
c      call laplacian(psi, Nx, Ny, dx, dy, work)

c      do j=1,Ny
c       do i=1, Nx
c        omg(i,j)=work(i,j)
        tmp(Nx/2,Ny/2)=0.01
c       end do
c      end do
c      call bc_channel(omg, Nx, Ny)
c      call bc_channel(tmp, Nx, Ny)
c

c      call poisson(psi, omg, Nx, Ny, dx, dy)

C
c      do 20 i=0,Nx
c        do 20 j=0,Ny
c          omg(i,j)=dble(0)                           
c   20 continue
C

  300 format (3(1X, E12.5))

      Return
      End
C----------------------------------------------------------------
      subroutine counter(i, data_number)
      integer   i, o_4, o_3, o_2, o_1
      character*10 i_data
      data      i_data/ '0123456789' /

      character*4 data_number

       o_4=i/1000
       o_3=(i-o_4*1000)/100
       o_2=(i-o_4*1000-o_3*100)/10
       o_1= i-o_4*1000-o_3*100-o_2*10
       data_number=i_data(o_4+1:o_4+1)//
     +                 i_data(o_3+1:o_3+1)//
     +                 i_data(o_2+1:o_2+1)//
     +                 i_data(o_1+1:o_1+1) 
      
      return
      end 
C*******************************************************************************
      subroutine output_data(psi, omg, tmp, k, Nx, Ny, dx, dy)
C*********●Graphic(by DCL)●**********
      integer      i, j, k
      integer      Nx, Ny
      real*8       dx, dy
      real*8       psi(0:Nx+1, 0:Ny+1)
      real*8       omg(0:Nx+1, 0:Ny+1)
      real*8       tmp(0:Nx+1, 0:Ny+1)
      real         psi_dcl(1:Nx, 0:Ny+1)
      real         omg_dcl(1:Nx, 0:Ny+1)
      real         tmp_dcl(1:Nx, 0:Ny+1)
      character*4  data_number
C
      do 10 j=0,Ny+1
        do 10 i=1,Nx
          psi_dcl(i,j)=real(psi(i,j))
          omg_dcl(i,j)=real(omg(i,j))
          tmp_dcl(i,j)=real(tmp(i,j))
   10 continue
C
      call counter(k, data_number)

      open(7,file='psi_'//data_number//'.dat', 
     +          status='unknown')
      open(8,file='omg_'//data_number//'.dat', 
     +          status='unknown')
      open(9,file='tmp_'//data_number//'.dat', 
     +          status='unknown')


      do 15 j=0,Ny+1
        do 15 i=1,Nx
          write(7,*) dble(i)*dx, dble(j)*dy, psi_dcl(i,j)
          write(8,*) dble(i)*dx, dble(j)*dy, omg_dcl(i,j)
          write(9,*) dble(i)*dx, dble(j)*dy, tmp_dcl(i,j)
   15 continue

      close(7)
      close(8)
      close(9)


      return
      end
C
C***************************************
      subroutine arakawa_jacobian(psi, omg, Nx, Ny, dx, dy,
     &                            Jacobian)

      integer   i, j, Nx, Ny
      real*8    dx, dy
      real*8    psi(0:Nx+1,0:Ny+1), omg(0:Nx+1,0:Ny+1)
      real*8    Jacobian_1, Jacobian_2, Jacobian_3, 
     &          Jacobian(1:Nx,1:Ny),
     &          factor

      factor=dble(1)/(dble(4)*dx*dy)

      do j=1,Ny                     !内部領域のみ
        do i=1,Nx
         Jacobian_1=((psi(i+1,j)-psi(i-1,j))*
     &   (omg(i,j+1)-omg(i,j-1))-(psi(i,j+1)-psi(i,j-1))*
     &   (omg(i+1,j)-omg(i-1,j)))*factor

         Jacobian_2=((psi(i+1,j))*
     &   (omg(i+1,j+1)-omg(i+1,j-1))-(psi(i-1,j))*
     &   (omg(i-1,j+1)-omg(i-1,j-1))-(psi(i,j+1))*
     &   (omg(i+1,j+1)-omg(i-1,j+1))+(psi(i,j-1))*
     &   (omg(i+1,j-1)-omg(i-1,j-1)))*factor

         Jacobian_3=((omg(i,j+1))*
     &   (psi(i+1,j+1)-psi(i-1,j+1))-(omg(i,j-1))*
     &   (psi(i+1,j-1)-psi(i-1,j-1))-(omg(i+1,j))*
     &   (psi(i+1,j+1)-psi(i+1,j-1))+(omg(i-1,j))*
     &   (psi(i-1,j+1)-psi(i-1,j-1)))*factor
          
        Jacobian(i,j)=(Jacobian_1+Jacobian_2+Jacobian_3)/dble(3)
        enddo
      enddo  


      return
      end
C**************************************************
      subroutine laplacian(A, Nx, Ny, dx, dy, Laplacian_A)
   
      integer   i, j, Nx, Ny
      real*8    dx, dy
      real*8    A(0:Nx+1, 0:Ny+1)
      real*8    Laplacian_A(1:Nx, 1:Ny)
     
      do j=1, Ny
        do i=1, Nx
        Laplacian_A(i,j)=((A(i+1,j)+A(i-1,j)-dble(2)*A(i,j))/(dx**2)+
     &                   (A(i,j+1)+A(i,j-1)-dble(2)*A(i,j))/(dy**2))
        end do
      end do  

      return
      end
C
C*******************************************************************************
       subroutine poisson(psi, omg, Nx, Ny, dx, dy)

       integer   c, i, j, Nx, Ny, c_max, i_poit, j_poit
       real*8    psi(0:Nx+1, 0:Ny+1), omg(0:Nx+1, 0:Ny+1)
       real*8    old_psi(0:Nx+1, 0:Ny+1)
       real*8    dx, dy
       real*8    w 
       real*8    epsilon, error, pres, value
      

       c_max=Nx
       w=dble(1.5)
C***● ψについてのポアソン方程式 ●***
      do 90 c=1,c_max*10

        do 50 i=1,Nx
          do 51 j=1,Ny
          epsilon=1.E-30
          pres=1.E-8

           old_psi(i,j)=psi(i,j)

           psi(i,j)=((dy**2)*(psi(i+1,j)+psi(i-1,j))+
     &     (dx**2)*(psi(i,j+1)+psi(i,j-1))-(dx**2)*(dy**2)*(omg(i,j)))/
     &     (dble(2)*((dx**2)+(dy**2)))

           error=psi(i,j)-old_psi(i,j)
c           psi(i,j)=w*psi(i,j)+(dble(1)-w)*old_psi(i,j)
           psi(i,j)=old_psi(i,j)+w*error

           if (dabs(error) .GT. epsilon) then
             epsilon=dabs(error)
             i_point=i
             j_point=j
             value=psi(i,j)
           end if

   51     continue
   50   continue
C
      call bc_channel(psi, Nx, Ny)
c        write(95,*) c, epsilon
           if (epsilon .LT. pres) then
           exit
           end if
C
   90 continue

 91   write(16,*) c, epsilon

 95   return
      end   
c--------------------------------------------------------------------------
      subroutine energy(psi, omg, Nx, Ny, ENRG, ENST)

c finish 2016.11.07

      integer Nx, Ny, i, j
      real*8  psi(0:Nx+1, 0:Ny+1), omg(0:Nx+1, 0:Ny+1)
      real*8  ENRG, ENST, half
  
      half=dble(0.5)
      ENRG=dble(0)
      ENST=dble(0)

      do j=1, Ny
       do i=1, Nx
          ENRG=ENRG-psi(i,j)*omg(i,j)*half
          ENST=ENST+omg(i,j)*omg(i,j)*half
       end do
      end do  

         ENRG=ENRG/dble((Nx)*(Ny))
         ENST=ENST/dble((Nx)*(Ny))

      return
      end

c--------------------------------------------------------------------------
      subroutine Nusselt(tmp, Nx, Ny, dy, Nu)

c finish 2016.11.07

      integer Nx, Ny, i, j
      real*8  tmp(0:Nx+1, 0:Ny+1), t_bar(0:Ny+1)
      real*8  Nu, half, dy
  
      half=dble(0.5)
      Nu=dble(0)

      do j=0,Ny+1
        t_bar(j)=dble(0)
      end do

      do j=0, Ny+1
       do i=1, Nx
          t_bar(j)=t_bar(j)+tmp(i,j)
       end do
          t_bar(j)=t_bar(j)/dble(Nx)
      end do  

      Nu=dble(1)-(t_bar(1)-t_bar(0))/dy

      return
      end
c-------------------------------------------------------------------
      subroutine Governing_Eqns(psi,omg,tmp, domg, dtmp,
     &     Nx,Ny,dx,Dy,Pr,Ra)

      integer Nx, Ny, i, j
      real*8  psi(0:Nx+1, 0:Ny+1), omg(0:Nx+1, 0:Ny+1)
      real*8  tmp(0:Nx+1, 0:Ny+1)
      real*8  domg(Nx, Ny), dtmp(Nx, Ny)
      real*8  adv_omg(Nx, Ny), adv_tmp(Nx, Ny)
      real*8  Laplacian_omg(Nx, Ny), Laplacian_tmp(Nx, Ny)


      real*8  dx, dy, Pr, Ra

C calculation of advection term
      call arakawa_jacobian(psi, omg, Nx, Ny, dx, dy, adv_omg) 
      call arakawa_jacobian(psi, tmp, Nx, Ny, dx, dy, adv_tmp) 
C calculation of Laplacian of vorticity
      call laplacian(omg, Nx, Ny, dx, dy, Laplacian_omg)
      call laplacian(tmp, Nx, Ny, dx, dy, Laplacian_tmp)


      do j=1, Ny
        do i=1, Nx
          domg(i,j)=-adv_omg(i,j)+Pr*laplacian_omg(i,j)
     &                  +Pr*Ra*(tmp(i+1,j)-tmp(i-1,j))/(dble(2)*dx)
c
          dtmp(i,j)=-adv_tmp(i,j)+laplacian_tmp(i,j)
     &                  +(psi(i+1,j)-psi(i-1,j))/(dble(2)*dx)
        end do
      end do

      return
      end
c--------------------------------------------------------------------
      subroutine Next(a, new_da, old_da, Nx, Ny, C1, C2, dt)

      integer  i, j, Nx, Ny
      real*8   a(0:Nx+1, 0:Ny+1), C1, C2, dt
      real*8   new_da(Nx, Ny), old_da(Nx, Ny)

      do j=1, Ny
        do i=1, Nx
          a(i,j)=a(i,j)+dt*(c1*new_da(i,j)-c2*old_da(i,j))
          old_da(i,j)=new_da(i,j)
        end do
      end do
 
      return
      end
