program analysis_spml

  use ee_module

  implicit none

  integer :: i, j
  integer :: k, l
  integer :: a
  integer, parameter :: b=300
  real(8) :: x, y
  real(8) :: t
  real(8) :: nu
  real(8) :: eta, k_eta
  
  real(8), allocatable :: yx_zeta(:,:)
  real(8), allocatable :: yx_psi(:,:)
  real(8), allocatable :: yx_F(:,:)
  
  real(8), allocatable :: ee_zeta(:,:)
  real(8), allocatable :: ee_psi(:,:)
  real(8), allocatable :: ee_F(:,:)

  real(8), allocatable :: E(:)

  real(8), parameter  :: pi = 3.141592653589793D0 ! 円周率

  character(4) :: data_number
  
  integer :: km=21                               ! X 方向切断全波数
  integer :: lm=21                               ! Y 方向切断全波数
  integer :: im=64                               ! X 方向格子点数 (>3*km)
  integer :: jm=64                               ! Y 方向格子点数 (>3*lm)

  integer :: kf                                  ! 強制波数 

  real(8) :: XLength                             ! 領域の大きさ(X 方向)
  real(8) :: YLength                             ! 領域の大きさ(Y 方向)
 
  allocate(yx_zeta(0:jm-1,0:im-1))
  allocate(yx_psi(0:jm-1,0:im-1))
  allocate(yx_F(0:jm-1,0:im-1))
  
  allocate(ee_zeta(-lm:lm,-km:km))
  allocate(ee_psi(-lm:lm,-km:km))
  allocate(ee_F(-lm:lm,-km:km))

  allocate(E(0:km))

  nu=2.5d-4/dble((im/64)**2)
  XLength=dble(2)*pi
  YLength=dble(2)*pi

  call ee_Initial(im,jm,km,lm,0.0D0,XLength,0.0D0,YLength)

  a=1
  do a=1, b

     call counter(a, data_number)
     
     open(50,file='pre_spectrum_'//data_number//'.dat')
     do j=0, jm-1
        do i=0, im-1
           read(50,*) t, x, y, yx_zeta(j,i), yx_psi(j,i), yx_F(j,i)
        end do
     end do
     close(50)

     ee_zeta=ee_yx(yx_zeta)

!----- エネルギースペクトルの計算 -----     
     do kf=0, km
        E(kf)=0.0d0
     end do
          
     do l=-lm, lm
        do k=-km, km
           kf=nint(sqrt(dble(k**2+l**2)))
           E(kf)=E(kf)+(ee_zeta(l,k)**2+ee_zeta(-l,-k)**2)/dble(2*(kf**2))
        end do
     end do
     

!----- 散逸規格化 -----
     eta=0.0d0
     do kf=1, km      
        eta=eta+(dble(kf)**dble(4))*E(kf)
     end do
     eta=dble(-2)*nu*eta
     k_eta=(nu**dble(-0.5))*(eta**dble(1/6))

     open(54,file='scaled_energy_spectrum_'//data_number//'.dat')
     do kf=1, km
        write(54,*) kf, E(kf)*(eta**dble(-2/3))*(k_eta**dble(3)), k_eta
     end do
     close(54)

  end do
  
end program analysis_spml

!----- datanumber -----
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
