 1    program NS_eq
      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***********● 変数の定義 ●***********
      integer i, j, k, a, b, c, k1, ic
      integer I_grd, J_grd                        !格子点
      integer I1, I2, J1, J2                      !平板の座標
      real*8  W, L                                !計算領域の横，縦
      real*8  dx, dy                              !格子点間隔
      real*8  dt                                  !時間差分間隔
      real*8  U_enter                             !流入する流体の速さ
      real*8  Jacobian_1, Jacobian_2, Jacobian_3  !各種 Arakawa Jacobian
      real*8  nu                                  !動粘性係数
      real*8  dpsi, Max_dpsi                      !(収束判定用)書き換え前と後のψの変化量, それの最大値
      real*8  error
C***********● 定数の定義 ●***********
      parameter (I_grd=400, J_grd=200)            !格子点数
      parameter (L=20.0d0, W=40.0d0)              !計算領域の縦，横の長さ
      parameter (I1=95, I2=105, J1=95, J2=105)      !平板の座標
C***********● 配列の定義 ●***********
      real*8 psi(0:I_grd,0:J_grd)                 !流れ関数の配列
      real*8 omg(0:I_grd,0:J_grd)                 !渦度の配列
      real*8 non_linear(1:I_grd-1, 1:J_grd-1)     !移流項
      real*8 x(0:I_grd), y(0:J_grd)               !x,y座標
      real*8 old_psi(1:I_grd-1,1:J_grd-1)         !Poisson方程式の収束判定用配列
      real*8 now_domgdt(1:I_grd-1,1:J_grd-1)      !Adams-Bashforth法に使われるk1ステップ目のomgの時間変化量
      real*8 old_domgdt(1:I_grd-1,1:J_grd-1)      !Adams-Bashforth法に使われる(k1-1)ステップ目のomgの時間変化量
      real*8 laplacian_omg(1:I_grd-1,1:J_grd-1)   !渦度のLaplacian

      open(10,file='accuracy.dat', status='unknown')
C************● 変数の値 ●************
      dx=W/dble(I_grd)                            !格子点間隔(dx)
      dy=L/dble(J_grd)                            !格子点間隔(dy)
      dt=0.001d0                                  !時間差分間隔(dt)・安定性を調べなさい!
      U_enter=1.0d0                               !流入する流体の速さ
      nu=1.0d-2                                   !粘性係数
C**************************************
      do 10 i=0,I_grd                             !x座標の設定
           x(i)=dx*dble(i)
   10 continue                                    
C
      do 20 j=0,J_grd                             !x座標の設定
           y(j)=dy*dble(j)
   20 continue                                    !y座標の設定
C-initial condition******************************************
      k=0
      call initial_condition(psi, omg,
     &   I_grd, J_grd, L, W, I1, I2, J1, J2, U_enter)     !initial condition
      call boundary_condition(psi, omg,
     &  I_grd, J_grd, L, W, I1, I2, J1, J2, U_enter)      !Boundary condition for vorticity
      call output_data(psi, omg, k, I_grd, J_grd)
C*******● Temporal evolution (Euler method)●*******
      call arakawa_jacobian(psi, omg, I_grd, J_grd, dx, dy, non_linear) 
                                                  !calculation of advection term
      call laplacian(omg, I_grd, J_grd, dx, dy, Laplacian_omg)
                                                  !calculation of Laplacian of vorticity
      do 130 j=1, J_grd-1
        do 131 i=1, I_grd-1
         if ((j .ge. J1 .and. j .le. J2)  .and. 
     &       (i .ge. I1 .and. i .le. I2)) then 
         goto 131
         else
          now_domgdt(i,j)=non_linear(i,j)+nu*laplacian_omg(i,j)
          omg(i,j)=omg(i,j)+dt*now_domgdt(i,j)
          old_domgdt(i,j)=now_domgdt(i,j)
         endif
 131    continue
 130  continue
C***● Solve the Poisson equation ●***
      call poisson(psi, omg, I_grd, J_grd, dx, dy,
     &     L,W,I1,I2,J1,J2,U_enter, error, ic)
      call boundary_condition(psi, omg, 
     &   I_grd,J_grd,L,W,I1,I2,J1,J2,U_enter)     !BC for vorticity
C*******● Temporal evolution (Adams-Bashforth method)●*******
      do 200 k=1,100
       do 100 k1=1,1000 
         call arakawa_jacobian(psi, omg, I_grd, J_grd, dx, dy, 
     &                         non_linear)
         call laplacian(omg, I_grd, J_grd, dx, dy, Laplacian_omg)

       do 230 j=1, J_grd-1
         do 230 i=1, I_grd-1
           if ((j .ge. J1 .and. j .le. J2) .and. 
     &        (i .ge. I1 .and. i .le. I2)) then 
         goto 231
         else
          now_domgdt(i,j)=non_linear(i,j)+nu*laplacian_omg(i,j)
          omg(i,j)=omg(i,j)+
     &    dt*(dble(1.5)*now_domgdt(i,j)-dble(0.5)*old_domgdt(i,j))
          old_domgdt(i,j)=now_domgdt(i,j)
         endif
 231    continue
 230  continue
C***● Solve the Poisson equation ●***
      call poisson(psi, omg, I_grd, J_grd, dx, dy,
     &     L, W, I1, I2, J1, J2, U_enter, error, ic)
      call boundary_condition(psi, omg,
     &   I_grd, J_grd, L, W, I1, I2, J1, J2, U_enter) 
C**************************************
c           Write(*,*) 'k=, Max_dpsi=' ,k , Max_dpsi
c           Max_dpsi=dble(0)
C
       write(10,*) ic, error
 100   continue
       call output_data(psi, omg, k, I_grd, J_grd)

C
  200 continue
C
      close(10)
      Stop
      End
C------------------------------------------------------------------------------
      subroutine boundary_condition(psi, omg, 
     & I_grd, J_grd, L, W, I1, I2, J1, J2, U_enter)

      integer i, j, k
      integer I_grd, J_grd                
      integer I1, I2, J1, J2              
      real*8 W, L                         
      real*8 dx, dy                       
      real*8 dt                           
      real*8 U_enter                      
C***********● 配列の定義 ●***********
      real*8 psi(0:I_grd,0:J_grd)         
      real*8 omg(0:I_grd,0:J_grd)         
      real*8 x(0:I_grd), y(0:J_grd)       
C************● 変数の値 ●************
      dx=W/dble(I_grd)                    
      dy=L/dble(J_grd)                    
C**************************************
C
      do 5 j=0, J_grd
           y(j)=(dble(j)-dble(J_grd/2))*dy
 5    continue

      do 10 i=0,I_grd
        psi(i,J_grd)=dble(J_grd)*U_enter*dy/dble(2)          !left boundary
                                                  
        psi(i,0)     =dble(J_grd)*U_enter*dy/dble(2)*dble(-1)!right boundary
   10 continue
C
      do 20 j=0,J_grd
        psi(0,j)=U_enter*dy*(dble(j)-dble(J_grd)/dble(2))    !upstream
C       psi(I_grd,j)=2*psi(I_grd-1,j)-psi(I_grd-2,j)         !down stream
        psi(I_grd,j)=psi(I_grd-1,j)
   20 continue
C
      do 30 j=J1,J2
        do 30 i=I1,I2
          psi(i,j)=dble(0)                                    !平板の中
   30 continue
C
      do 40 i=1,I_grd-1
c           omg(i,J_grd)=dble(0)                        
c           omg(i,0)     =dble(0)                       
        omg(i,0)=dble(-2)/(dy**2)*
     &       (psi(i,1)-psi(i,0)-U_enter*(dy))            !left boundary
        omg(i,J_grd)=dble(-2)/(dy**2)*
     &       (psi(i,J_grd-1)-psi(i,J_grd)-U_enter*(-dy)) !right boundary
   40 continue
C
      do 50 j=0,J_grd
           omg(0,j)     =dble(0)                         !upstream
           omg(I_grd,j)=omg(I_grd-1,j)                   !down stream
   50 continue
C
      do 60 j=J1+1,J2-1
           do 60 i=I1+1,I2-1
           omg(i,j)=dble(0)                              !平板の中
   60 continue
C
      do 70 j=J1+1,J2-1
          omg(I1,j)=dble(-2)*psi(I1-1,j)/(dx**2) !平板の左壁
          omg(I2,j)=dble(-2)*psi(I2+1,j)/(dx**2) !平板の右壁
   70 continue
C
      do 80 i=I1+1,I2-1
          omg(i,J1)=dble(-2)*psi(i,J1-1)/(dy**2) !平板の下壁
          omg(i,J2)=dble(-2)*psi(i,J2+1)/(dy**2) !平板の上壁
   80 continue
C
          omg(I1,J1)=dble(-2)*psi(I1-1,J1-1)/(dx**2+dy**2) !
          omg(I1,J2)=dble(-2)*psi(I1-1,J2+1)/(dx**2+dy**2) !
          omg(I2,J1)=dble(-2)*psi(I2+1,J1-1)/(dx**2+dy**2) !
          omg(I2,J2)=dble(-2)*psi(I2+1,J2+1)/(dx**2+dy**2) !


      Return
      End
C------------------------------------------------------------------------------
      subroutine initial_condition(psi, omg,
     &                   I_grd, J_grd, L, W, I1, I2, J1, J2, U_enter)
C***********● 変数の定義 ●***********
      integer i, j, k
      integer I_grd, J_grd                
      integer I1, I2, J1, J2              
      real*8 W, L                         
      real*8 dx, dy                       
      real*8 U_enter                        
C***********● 配列の定義 ●***********
      real*8 psi(0:I_grd,0:J_grd)         
      real*8 omg(0:I_grd,0:J_grd)       
      real*8 x(0:I_grd), y(0:J_grd)     
C************● 変数の値 ●************
      dx=W/dble(I_grd)                  
      dy=L/dble(J_grd)                  
C**************************************
      do 5 j=0, J_grd
           y(j)=(dble(j)-dble(J_grd/2))*dy
 5    continue

      do 10 i=0,I_grd
        do 10 j=0,J_grd
          psi(i,j)=U_enter*y(j)
   10 continue
C
      do 20 i=0,I_grd
        do 20 j=0,J_grd
          omg(i,j)=dble(0)                           
   20 continue
C
      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, k, I_grd, J_grd)
C*********●Graphic(by DCL)●**********
      integer      i, j, k
      integer      I_grd, J_grd
      real*8       psi(0:I_grd, 0:J_grd)
      real*8       omg(0:I_grd, 0:J_grd)
      real         psi_dcl(1:I_grd+1, 1:J_grd+1)
      real         omg_dcl(1:I_grd+1, 1:J_grd+1)
      character*4  data_number
C
      do 10 j=0,J_grd
        do 10 i=0,I_grd
          psi_dcl(i+1,j+1)=real(psi(i,j))
          omg_dcl(i+1,j+1)=real(omg(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')


      do 15 j=1,J_grd+1
        do 15 i=1,I_grd+1
          write(7,*) psi_dcl(i,j)
          write(8,*) omg_dcl(i,j)
   15 continue

      close(7)
      close(8)

      return
      end
C
C***************************************
      subroutine arakawa_jacobian(psi, omg, I_grd, J_grd, dx, dy,
     &                            Jacobian)

      integer   i, j, I_grd, J_grd
      real*8    dx, dy
      real*8    psi(0:I_grd,0:J_grd), omg(0:I_grd,0:J_grd)
      real*8    Jacobian_1, Jacobian_2, Jacobian_3, 
     &          Jacobian(1:I_grd-1,1:J_grd-1),
     &          factor

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

      do j=1,J_grd-1                     !内部領域のみ
        do i=1,I_grd-1
         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, I_grd, J_grd, dx, dy, Laplacian_A)
   
      integer   i, j, I_grd, J_grd
      real*8    dx, dy
      real*8    A(0:I_grd, 0:J_grd)
      real*8    Laplacian_A(1:I_grd-1, 1:J_grd-1)
     
      do j=1, J_grd-1
        do i=1, I_grd-1
        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, I_grd, J_grd, dx, dy,
     &            L, W, I1, I2, J1, J2, U_enter, error, ic)

       integer   c, i, j, I_grd, J_grd,I1,I2,J1,J2, ic
       real*8    psi(0:I_grd, 0:J_grd), omg(0:I_grd, 0:J_grd)
       real*8    old_psi(0:I_grd, 0:J_grd)
       real*8    dx, dy,L,W,U_enter, epsilon, error

       error=dble(0)
C***● ψについてのポアソン方程式 ●***
      do 90 c=1,100
        do 50 i=1,I_grd-1
          do 51 j=1,J_grd-1
         if ((j .ge.j1 .and. j .le. j2) .and. 
     &       (i .ge. i1 .and. i .le. i2)) then 
         goto 51
         endif
  
         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)))
C
C                          If ((i.EQ.I_grd-1).AND.(j.EQ.J_grd-1)) Then
C                               do 40 a=1,I_grd-1
C                                    do 40 b=1,J_grd-1
C                                        dpsi=DABS(psi(a,b)-old_psi(a,b))
C
C                                         If (Max_dpsi .LT. dpsi) Then
C                                              Max_dpsi=dpsi             !計算領域の中で最も大きい変化量
C                                         End if
C  40                          continue                                 !Max_dpsiの値が唯一に決まった
C                               If (Max_dpsi .LT. 0.0000001) Then
C                                    Write(*,*) c ,k ,Max_dpsi
C                                    Max_dpsi=dble(0)
C                                    Go to 95
C                               End if
C                          End if

         epsilon=dabs(old_psi(i,j)-psi(i,j))
         error=max(epsilon,error)
   51     continue
   50   continue
C
       call boundary_condition(psi, omg,
     & I_grd, J_grd, L, W, I1, I2, J1, J2, U_enter)
C
c      if (error .lt. 1.D-5) then 
c      goto 95
c      end if     
   90 continue

c   95 continue
 95   ic=c

      return
      end   
C
C
C************○ 製作日記 ○************
C 1130 製作開始．初期条件，境界条件の入力完了．ωの時間発展も入力した．しかし間違って入力している可能性もあるので，明日は確認から．
C 1201 一応最後まで，組んだ．しかし，ポアソン方程式がうまく解けているのかわからない．収束判定がうまくいっていない．どうしよう．	
C 1204 判定がうまくいってなかったのではなくて，ポアソン方程式の差分式の括弧の位置がおかしかった．次は判定条件をつけてみる．
C 1207 オイラー法→アダムスバッシュホース法．Jacobianのバグを修正．Jacobianの符号も修正．
C 1211 DCLでコンター図が描けるようになった．dclないでは倍精度ではなく単精度みたい．
C 1212 U_enterの値を0.01→1.00
C 1213 nuの値を0.1→0.01(Re:10→100)
C 1216 nuの値を0.01→0.1(Re:100→10)
C 1222 出力をデータファイルにするように改善．計算領域を横に倍増させる．
C 1223 nu=0.05 (Re:20)
C 1226 nu=0.025(Re:40)
C 1228 nu=1/60=0.0166666666666(Re:60)
C 1229 nu=0.0125(Re:80)
C 0110 nu=1/120=0.008333333333333(Re:120)
C 0118 nu=1/140=(Re:140)
C 0122 nu=1/160=(Re:160)
C 0205 nu=1/200=(Re=200)
C 0208 nu=1/500
C C  0213 流出境界を2階微分に変更 nu=1/1000
C 0220 nu=
C 0221 poisson方程式の繰り返し 1000→500 異常なし？
C 0222 poisson方程式の繰り返し  500→100 
C 0226 領域2倍
C 0402 領域 10×40 Re:100 
C      注)パラメータを変更するときは, ファイルの先頭から, 全てを置換しなければならない. 
c 2015.12.18 modified by T. IWAYAMA
c 2015.12.30
C**************************************
C
C
