!----- 渦度方程式の数値計算 -----
program forced_vor_eq
  
  implicit none

!----- 変数の定義 -----  
  integer :: i, j, k_t, l
  integer, parameter :: n=64, nn=n**2
  integer :: k_ell=8
  integer :: ir
  real(8) :: t, dt, t_max, t_out
  real(8) :: x, dx
  real(8) :: y, dy
  real(8) :: pi, Lx, Ly
  real(8) :: nu
  real(8) :: R=0.5d0
  real(8), parameter :: c1=dble(3)/dble(2)
  real(8), parameter :: c2=dble(1)/dble(2)
  
!----- 配列の定義 -----  
  real(8) :: zeta(0:n+1,0:n+1)
  real(8) :: psi(0:n+1,0:n+1)
  
  real(8) :: d2zeta_x(1:n,1:n)
  real(8) :: d2zeta_y(1:n,1:n)

  real(8) :: jacobian(1:n,1:n)

  real(8) :: dtzeta(1:n,1:n)
  real(8) :: dtzeta_old(1:n,1:n)
  
  real(8) :: F(0:n+1,0:n+1)
  
  real(8) :: e ! エネルギー
  real(8) :: z ! エンストロフィー
  real(8) :: kk
  
  integer :: t_step, k_out
  character(4) :: data_number

!----- 各 patameter の設定 -----
  t_max=600.0d0
  dt=4.0d-2
  t_out=2.0d0
  t_step=int(t_max/dt) !---全時間 step 数---15000
  k_out=int(t_out/dt)  !---50
  
  pi=acos(-1.0d0)
  Lx=dble(2)*pi
  Ly=dble(2)*pi
  dx=Lx/dble(n)
  dy=Ly/dble(n)

  nu=2.5d-4/dble((n/64)**2)  

!----- 初期条件(F) -----
  k_t=0
  l=0
  ir=1
  t=dble(0)
  e=dble(0)
  z=dble(0)
  
  do j=1, n
     do i=1, n
        F(i,j)=0.0d0
     end do
  end do
  call bc_cyclic(n, F)

  call force(n, nn, k_ell, ir, dx, dy, R, F)

  do j=1, n
     do i=1, n
        zeta(i,j)=0.0d0
     end do
  end do    
  call bc_cyclic(n, zeta)

  open(30,file='iteration.dat')

  call poisson(n, k_t, dx, dy, zeta, psi)

  do j=1, n
     do i=1, n
        e=e-zeta(i,j)*psi(i,j)/dble(2)
        z=z+zeta(i,j)**dble(2)/dble(2)
     end do
  end do

  e=e/dble(Lx*Ly)
  z=z/dble(Lx*Ly)
  kk=z/e
  
  call counter(l, data_number)
  
  open(10,file='fort_'//data_number//'.dat')
  do j=0, n
     y=dble(j)*dy
     do i=0, n
        x=dble(i)*dx
        write(10,*) t, x, y, zeta(i,j), psi(i,j), F(i,j)
     end do
     write(10,*) 
  end do
  close(10)

  open(40,file='pre_spectrum_'//data_number//'.dat')
  do j=0, n-1
     y=dble(j)*dy
     do i=0, n-1
        x=dble(i)*dx
        write(40,*) t, x, y, zeta(i,j), psi(i,j), F(i,j)
     end do
  end do
  close(40)

  open(20,file='conv.dat')
  write(20,*) t, e, z, kk
  
!----- 時間発展の計算(Adams-Bashforth 法) -----
  call second_deriv(n, dx, dy, zeta, d2zeta_x, d2zeta_y)
  call arakawa_jacobian(n, dx, dy, zeta, psi, jacobian)
  do j=1, n
     do i=1, n
        dtzeta(i,j)=nu*(d2zeta_x(i,j)+d2zeta_y(i,j))-jacobian(i,j) &
             +F(i,j)
     end do
  end do

  do j=1, n
     do i=1, n
        zeta(i,j)=zeta(i,j)+dtzeta(i,j)*dt
        dtzeta_old(i,j)=dtzeta(i,j)
     end do
  end do
  
  call bc_cyclic(n, zeta)
  call poisson(n, k_t, dx, dy, zeta, psi)
  
  do k_t=1, t_step
     t=t+dt
     e=dble(0)
     z=dble(0)

     call force(n, nn, k_ell, ir, dx, dy, R, F)
     
     call second_deriv(n, dx, dy, zeta, d2zeta_x, d2zeta_y)
     call arakawa_jacobian(n, dx, dy, zeta, psi, jacobian)
     
     do j=1, n
        do i=1, n
           dtzeta(i,j)=nu*(d2zeta_x(i,j)+d2zeta_y(i,j))-jacobian(i,j) &
                +F(i,j)
        end do
     end do

     do j=1, n
        do i=1, n
           zeta(i,j)=zeta(i,j)+(c1*dtzeta(i,j)-c2*dtzeta_old(i,j))*dt
           dtzeta_old(i,j)=dtzeta(i,j)
        end do
     end do

     call bc_cyclic(n, zeta)
     call poisson(n, k_t, dx, dy, zeta, psi)

!----- data 出力 -----
     if ( mod(k_t, k_out) == 0 ) then
        l=l+1

        do j=1, n
           do i=1, n
              e=e-zeta(i,j)*psi(i,j)/dble(2)
              z=z+zeta(i,j)**dble(2)/dble(2)
           end do
        end do

        e=e/dble(Lx*Ly)
        z=z/dble(Lx*Ly)
        kk=z/e
        
        call counter(l, data_number)

        open(10,file='fort_'//data_number//'.dat')
        do j=0, n
           y=dble(j)*dy
           do i=0, n
              x=dble(i)*dx
              write(10,*) t, x, y, zeta(i,j), psi(i,j), F(i,j)
           end do
           write(10,*) 
        end do
        close(10)

        open(40,file='pre_spectrum_'//data_number//'.dat')
        do j=0, n-1
           y=dble(j)*dy
           do i=0, n-1
              x=dble(i)*dx
              write(40,*) t, x, y, zeta(i,j), psi(i,j), F(i,j)
           end do
        end do
        close(40)
 
        write(20,*) t, e, z, kk
        
     end if
  end do

  close(20)
  close(30)
  
end program forced_vor_eq

!----- poisson 方程式の数値計算 -----
subroutine poisson(n, k_t, dx, dy, zeta, psi)

  implicit none

  integer, intent(in) :: n
  integer, intent(in) :: k_t
  real(8), intent(in) :: dx, dy
  real(8), intent(in) :: zeta(0:n+1,0:n+1)
  real(8), intent(inout) :: psi(0:n+1,0:n+1)

  integer :: i, j, k
  real(8) :: x, y

  real(8) :: psi_old(0:n+1,0:n+1)
  real(8) :: e(1:n,1:n)
  real(8) :: e_max

  integer, parameter :: n_max=100000000
  
!----- psi に仮の初期値を与える -----
  if ( k_t == 0 ) then
     
     do j=1, n
         do i=1, n
            psi(i,j)=dble(0)
        end do
     end do

     call bc_cyclic(n, psi)
    
  end if
  
  !----- 反復計算 -----
  do k=1, n_max
     e_max=dble(0)
     
     do j=1, n
        do i=1, n
           psi_old(i,j)=psi(i,j)
           psi(i,j)=((psi(i+1,j)+psi(i-1,j))*(dy**2) &
                +(psi(i,j+1)+psi(i,j-1))*(dx**2) &
                -zeta(i,j)*(dx**2)*(dy**2)) &
                /(dble(2)*((dx**2)+(dy**2)))

           e(i,j)=abs(psi(i,j)-psi_old(i,j))
           e_max=max(e_max, e(i,j))
        end do        
     end do

     call bc_cyclic(n, psi)

     if ( e_max < 1.0d-8 ) then
        write(30,*) k
        exit
     end if

  end do
  
end subroutine poisson

!----- 空間一階微分(ζ,ψ)-----
subroutine central_deriv(n, dx, dy, u, d1u_x, d1u_y)

  implicit none
  
  integer, intent(in) :: n
  real(8), intent(in) :: dx, dy
  real(8), intent(in) :: u(0:n+1,0:n+1)
  real(8), intent(out) :: d1u_x(1:n,1:n)
  real(8), intent(out) :: d1u_y(1:n,1:n)
  integer :: i, j

  do j=1, n
     do i=1, n
        d1u_x(i,j)=(u(i+1,j)-u(i-1,j))/(dble(2)*dx)
        d1u_y(i,j)=(u(i,j+1)-u(i,j-1))/(dble(2)*dy)
     end do
  end do

end subroutine central_deriv

!----- 空間二階微分(ζ) -----
subroutine second_deriv(n, dx, dy, u, d2u_x, d2u_y)

  implicit none

  integer, intent(in) :: n
  real(8), intent(in) :: dx, dy
  real(8), intent(in) :: u(0:n+1,0:n+1)
  real(8), intent(out) :: d2u_x(1:n,1:n)
  real(8), intent(out) :: d2u_y(1:n,1:n)
  integer :: i, j

  do j=1, n
     do i=1, n
        d2u_x(i,j)=(u(i+1,j)+u(i-1,j)-dble(2)*u(i,j))/(dx**2)
        d2u_y(i,j)=(u(i,j+1)+u(i,j-1)-dble(2)*u(i,j))/(dy**2)
     end do
  end do

end subroutine second_deriv

!----- Arakwa-Jacobian Method -----
subroutine arakawa_jacobian(n, dx, dy, omg, psi, jacobian)

  implicit none
  
  integer, intent(in) :: n
  real(8), intent(in) :: dx, dy
  real(8), intent(in) :: omg(0:n+1,0:n+1)
  real(8), intent(in) :: psi(0:n+1,0:n+1)
  real(8), intent(out) :: jacobian(1:n,1:n)

  integer :: i, j
  real(8) :: Jacobian_1
  real(8) :: Jacobian_2
  real(8) :: Jacobian_3
  real(8) :: factor
  
  factor=dble(1)/(dble(4)*dx*dy)

  do j=1, n   
     do i=1, n
        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)
     end do
  end do

end subroutine arakawa_jacobian

!----- 二重周期境界条件 -----
subroutine bc_cyclic(n, A)

  implicit none

  integer, intent(in) :: n
  real(8), intent(out) :: A(0:n+1,0:n+1)
  integer :: i, j

  do i=1, n
     A(i,0)=A(i,n)
     A(i,n+1)=A(i,1)
  end do

  do j=1, n
     A(0,j)=A(n,j)
     A(n+1,j)=A(1,j)
  end do

  A(0,0)=A(n,n)
  A(n+1,0)=A(1,0)
  A(0,n+1)=A(0,1)
  A(n+1,n+1)=A(1,1)
  
end subroutine bc_cyclic

!----- 強制関数の計算 -----
subroutine force(n, nn, k_ell, ir, dx, dy, R, F)
  
  implicit none

  integer, intent(in) :: n, nn
  integer, intent(in) :: k_ell
  integer, intent(inout) :: ir
  real(8), intent(in) :: dx, dy
  real(8) ,intent(in) :: R
  real(8), intent(inout) :: F(0:n+1,0:n+1)

  integer :: i, j
  integer :: k_x, k_y
  integer, parameter :: m=1
  real(8) :: x, y
  real(8) :: pi
  real(8) :: phi(1:m)          !--- 0-1 の乱数
  real(8) :: theta(1:m)        !--- 0-2*pi の乱数
  real(8), parameter :: A=1.0d0
  complex(8), parameter :: z=(0.0d0,1.0d0)

  pi=acos(-1.0d0)
  
  call urand1(m, phi, ir)
  theta(m)=dble(2)*pi*phi(m)
  
  call urand1(m, phi, ir)
  k_x=nint((2.0d0*phi(m)-1.0d0)*dble(k_ell))

  call urand1(m, phi, ir)
  k_y=nint((dble(k_ell)**dble(2)-dble(k_x)**dble(2))**dble(0.5))
  
  if ( phi(m) < 0.5d0 ) then
     k_y=-k_y
  end if
  
  do j=1, n
     y=dble(j)*dy
     do i=1, n
        x=dble(i)*dx
        F(i,j)=R*F(i,j) &
             +A*((dble(1)-R**dble(2))**dble(0.5))* &
             sin(dble(k_x)*x+dble(k_y)*y+theta(m))
     end do
  end do

  call bc_cyclic(n, F)

end subroutine force

!----- 乱数出力(0-1) -----
subroutine urand1(n, x, ir)

  implicit none
  
  integer, intent(in) :: n
  real(8), intent(out) :: x(1:n)
  integer, intent(inout) :: ir

  integer :: i
  integer, parameter :: m = 1664501, lambda = 1229, mu = 351750
  real(8), parameter :: invm = 1.0d0 / m

  if( n .le. 0) then
     
     write(6,*) '(subr.urand1) parameter error. n = ', n
     write(6,*) 'return with no further calculation.'

     return

  endif

  if( ir .lt. 0 .or. ir .ge. m) then

     write(6,*) '(subr.urand1) warning. ir = ', ir

  endif

  do i = 1, n
     ir = mod( lambda * ir + mu, m)
     x(i) = ir * invm
  end do
  
end subroutine urand1

!----- データ出力 -----
subroutine counter(i, data_number)

  implicit none

  integer, intent(in) :: i
  character(4), intent(out) :: data_number
  integer :: o_4, o_3, o_2, o_1
  character(10) :: i_data
  
  i_data='0123456789'

  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)
  
end subroutine counter

    
   
    
