      program zonal_mean_anal
      implicit NONE
! 2017.01.12 Thompson(1970)のデータ解析
! zonal mean temperatureの計算
! zonal mean zonal velocityの計算
! zonal mean momentum fluxの計算
!***********● 変数の定義 ●***********
      integer i, j, t, t1
      integer Nx, Nz                              !x, z方向の格子点数
      integer t_max
      real*8  Lx, Lz                              !計算領域のx, y方向の長さ
      real*8  dx, dz                              !格子点間隔
      real*8  dt                                  !時間差分間隔
!***********● 定数の定義 ●***********
      parameter (Nx=20, Nz=10)                  !格子点数
!***********● 配列の定義 ●***********
      real*8 x(0:Nx+1), z(0:Nz+1)     ! x,y座標
      real*8 tmp(0:Nx+1,-1:Nz+2)      ! 温度の配列
      real*8 psi(0:Nx+1,-1:Nz+2)      ! 流れ関数の配列
      real*8 tmp_bar(-1:Nz+2)         ! zonal mean temperature 
      real*8 u_bar(-1:Nz+2)           ! zonal mean temperature
      real*8 uw(-1:Nz+2)              ! momentum flux 
c constans 
      real*8   unity, zero, half, pi
      parameter (
     &            unity=1.D+0,
     &            zero=0.D+0, 
     &            half=0.5d+0
     &           )
!************● 変数の値 ●************
      pi=dble(acos(-1.D+0))           ! piの値
      Lx=dble(2)*pi                   ! 計算領域の横の長さ
      Lz=pi                           ! 計算領域の縦の長さ

      dx=Lx/dble(Nx)                  ! 格子点間隔(dx)
      dz=Lz/dble(Nz+1)                ! 格子点間隔(dz)

      t_max=1000
      dt=1.d+0                        !時間差分間隔

      do 100 t=0, t_max
      call input_data(psi, tmp, x, z, t, Nx, Nz, dx, dz, Lx, Lz)
      call zonal_mean(psi, tmp, u_bar, tmp_bar, uw, z, 
     &     Nx, Nz, dx, dz, Lx, Lz)          
      call output_data(u_bar,tmp_bar, uw, z, t, Nz, dz, Lz)
 100  continue
!
      stop
      end
!----------------------------------------------------------------
      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(u_bar, tmp_bar, uw, z, k, Nz, dz, Lz)
      integer      j, k
      integer      Nz
      real*8       dz, Lz
      real*8       tmp_bar(-1:Nz+2), z(0:Nz+1)
      real*8       u_bar(-1:Nz+2), uw(-1:Nz+2)
      character*4  data_number

      call counter(k, data_number)

      open(7,file='./data2/u_bar_'//data_number//'.dat', 
     +          status='unknown')
      open(8,file='./data2/tmp_bar_'//data_number//'.dat', 
     +          status='unknown')
      open(9,file='./data2/mom_flux_'//data_number//'.dat', 
     +          status='unknown')

      do 15 j=0,Nz+1
          write(7,*) z(j), u_bar(j)
          write(8,*) z(j), tmp_bar(j)
          write(9,*) z(j), uw(j)
   15 continue

      close(8)
      close(9)

      return
      end
!*******************************************************************************
      subroutine input_data(psi, tmp, x, z, k, Nx, Nz, dx, dz, Lx, Lz)
!*********●Graphic(by DCL)●**********
      integer      i, j, k
      integer      Nx, Nz
      real*8       dx, dz, Lx, Lz
      real*8       tmp(0:Nx+1, -1:Nz+2), x(0:Nx+1), z(0:Nz+1)
      real*8       psi(0:Nx+1, -1:Nz+2)
      character*4  data_number
C
c      do 10 j=1,Nz
c        do 10 i=1,Nx
c          psi_dcl(i,j)=real(psi(i,j))
c          omg_dcl(i,j)=real(omg(i,j))
c          tmp_dcl(i,j)=real(tmp(i,j))
c   10 continue
C
      call counter(k, data_number)

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


      do 15 j=0,Nz+1
        do 15 i=0,Nx
          read(8,*) x(i), z(j), psi(i,j)
          read(9,*) x(i), z(j), tmp(i,j)
   15 continue

      do j=0, Nz+1
       psi(Nx+1,j)=psi(0,j)
       tmp(Nx+1,j)=tmp(0,j)
      end do
      close(9)

      return
      end
!--------------------------------------------------------------------------
      subroutine energy(psi, omg, tmp, Nx, Nz, ENRG_Z, ENRG_C, 
     &           ENRG, gamma, dx, dz)

! finish 2016.11.07
! 2016.12.27

      integer Nx, Nz, i, j
      real*8  psi(0:Nx+1, -1:Nz+2), omg(0:Nx+1, 0:Nz+1)
      real*8  tmp(0:Nx+1, -1:Nz+2)
      real*8  ENRG_Z, ENRG_C, half, dx, dz, u_bar(0:Nz+1)
      real*8  ENRG, gamma
      real*8  u_p, pi
  
      half=dble(0.5)
      ENRG_Z=dble(0)
      ENRG_C=dble(0)
      ENRG=dble(0)
      pi=dble(acos(-1.))
! zonal component ---------------------------------------------
      do j=0, Nz+1
        u_bar(j)=dble(0)
      end do

      j=0
        do i=1, Nx
           u_bar(j)=u_bar(j)-(psi(i,j+1)-psi(i,j))/(dz)       
        end do
           u_bar(j)=u_bar(j)/dble(Nx)
           ENRG_Z=ENRG_Z+u_bar(j)**2*half*dx*dz*half*dble(Nx)

      do j=1, Nz
        do i=1, Nx
           u_bar(j)=u_bar(j)-(psi(i,j+1)-psi(i,j-1))/(dble(2)*dz)       
        end do
           u_bar(j)=u_bar(j)/dble(Nx)
           ENRG_Z=ENRG_Z+u_bar(j)**2*half*dx*dz*dble(Nx)
      end do

      j=Nz+1
        do i=1, Nx
           u_bar(j)=u_bar(j)-(psi(i,j)-psi(i,j-1))/(dz)       
        end do
           u_bar(j)=u_bar(j)/dble(Nx)
           ENRG_Z=ENRG_Z+u_bar(j)**2*half*dx*dz*half*dble(Nx)

c      do j=0, Nz+1
c      write(10,*) dble(j)*dz-pi, u_bar(j)
c      end do
!--------------------------------------------------------------------
! perturbation component
      j=0
       do i=1, Nx
          u_p=-(psi(i,j+1)-psi(i,j))/(dz)-u_bar(j)
          ENRG_C=ENRG_C+u_p**2*half*dx*dz*half
       end do

      do j=1, Nz
       do i=1, Nx
          u_p=-(psi(i,j+1)-psi(i,j-1))/(dble(2)*dz)-u_bar(j)
          ENRG_C=ENRG_C+u_p**2*half*dx*dz
       end do
      end do  

      j=Nz+1
      do i=1, Nx
          u_p=-(psi(i,j)-psi(i,j-1))/(dz)-u_bar(j)
          ENRG_C=ENRG_C+u_p**2/half*dx*dz*half
       end do
!--------------------------------------------------------------------
! total
      j=0
       do i=1, Nx
          u_p=-(psi(i,j+1)-psi(i,j))/(dz)
          ENRG=ENRG+(u_p**2+tmp(i,j)**2/gamma)*half*dx*dz*half
       end do

      do j=1, Nz
       do i=1, Nx
          u_p=-(psi(i,j+1)-psi(i,j-1))/(dble(2)*dz)
          ENRG=ENRG+(u_p**2+tmp(i,j)**2/gamma)*half*dx*dz
       end do
      end do  

      j=Nz+1
      do i=1, Nx
          u_p=-(psi(i,j)-psi(i,j-1))/(dz)
          ENRG=ENRG+(u_p**2+tmp(i,j)**2/gamma)*half*dx*dz*half
       end do

      return
      end
c-----------------------------------------------------------------------
      subroutine zonal_mean(psi, tmp, u_bar, tmp_bar, uw, z, 
     &           Nx, Nz, dx, dz, Lx, Lz)          
c
      integer Nx, Nz, i, j
      real*8  psi(0:Nx+1, -1:Nz+2), tmp(0:Nx+1, -1:Nz+2)
      real*8  tmp_bar(-1:Nz+2), z(0:Nz+1)
      real*8  u_bar(-1:Nz+2), uw(-1:Nz+2)
      real*8  u_p, w_p
      real*8  dx, dz, Lx, Lz


! zonal mean temperature
      do j=0, Nz+1
        tmp_bar(j)=dble(0)
      end do

      do j=0, Nz+1
        do i=1, Nx
          tmp_bar(j)=tmp_bar(j)+tmp(i,j)
        end do
          tmp_bar(j)=tmp_bar(j)/dble(Nx)+z(j)
      end do

! zonal mean zonal velocity---------------------------------------------
      do j=0, Nz+1
        u_bar(j)=dble(0)
      end do

      j=0
        do i=1, Nx
           u_bar(j)=u_bar(j)-(psi(i,j+1)-psi(i,j))/(dz)       
        end do
           u_bar(j)=u_bar(j)/dble(Nx)

      do j=1, Nz
        do i=1, Nx
           u_bar(j)=u_bar(j)-(psi(i,j+1)-psi(i,j-1))/(dble(2)*dz)       
        end do
           u_bar(j)=u_bar(j)/dble(Nx)
      end do

      j=Nz+1
        do i=1, Nx
           u_bar(j)=u_bar(j)-(psi(i,j)-psi(i,j-1))/(dz)       
        end do
           u_bar(j)=u_bar(j)/dble(Nx)

! zonal mean momentum flux---------------------------------------------
      do j=0, Nz+1
        uw(j)=dble(0)
      end do

      j=0
        do i=1, Nx
           u_p=-(psi(i,j+1)-psi(i,j))/(dz)-u_bar(j)
           w_p=(psi(i+1,j)-psi(i-1,j))/(dble(2)*dx)
           uw(j)=uw(j)+u_p*w_p       
        end do
           uw(j)=uw(j)/dble(Nx)

      do j=1, Nz
        do i=1, Nx
           u_p=-(psi(i,j+1)-psi(i,j-1))/(dble(2)*dz)-u_bar(j)
           w_p=(psi(i+1,j)-psi(i-1,j))/(dble(2)*dx)
           uw(j)=uw(j)+u_p*w_p              
        end do
           uw(j)=uw(j)/dble(Nx)
      end do

      j=Nz+1
        do i=1, Nx
           u_p=-(psi(i,j)-psi(i,j-1))/(dz)-u_bar(j)
           w_p=(psi(i+1,j)-psi(i-1,j))/(dble(2)*dx)       
           uw(j)=uw(j)+u_p*w_p              
        end do
           uw(j)=uw(j)/dble(Nx)

      return
      end
