      program Thompson_1970_AB
      implicit NONE
! originally produced by Mr. Nanamichi
! 2015.12.18 modified by T. IWAYAMA
! 2015.12.30 時間発展のバグを修正
! 2015.12.31 境界条件のバグを修正 Re=100ではKarman渦発生せず.
! 2016.01.01 omegaの境界条件を再考
! 2016.02.01 収束判定(Mizushima and Takemoto)
! 2016.11.05-11.17 2D turbulence
! 2016.11.17- Boussinesq fluids, convection
! 2016.11.18 Poisson方程式の解法 : SOR法完成 絶対誤差で収束性判定
! 2016.11.19 Dennou Sapmleを再現
! http://www.gfd-dennou.org/library/spmodel/2d-channel-esc/boussinesq/sample-tt/f90/bsncnv-tt-1.f90
! 2016.11.20 熱対流用のプログラムを改変
! 2016.12.27 Thompsonの条件にE, K, u_bar, Energyを合わせる
! 2016.12.28 RKG法に改良
! 2016.12.31 水平粘性を導入
! 2017.01.01- 境界上での温度方程式
!***********● 変数の定義 ●***********
      integer i, j, t, t1
      integer Nx, Nz                              !x, z方向の格子点数
      integer t_max, t1_max                       !
      real*8  Lx, Lz                              !計算領域のx, y方向の長さ
      real*8  dx, dz                              !格子点間隔
      real*8  dt                                  !時間差分間隔
      real*8  cap_G, cap_E_h, cap_E_v, gamma
      real*8  cap_K_h, cap_K_v          !
!***********● 定数の定義 ●***********
      parameter (Nx=10, Nz=5)                  !格子点数
!***********● 配列の定義 ●***********
      real*8 x(0:Nx+1), z(0:Nz+1)     ! x,y座標
      real*8 psi(0:Nx+1,-1:Nz+2)       ! 流れ関数の配列
      real*8 omg(0:Nx+1, 0:Nz+1)       ! 渦度の配列
      real*8 tmp(0:Nx+1,-1:Nz+2)       ! 温度の配列
      real*8 adv_omg(Nx,Nz), adv_tmp(Nx,0:Nz+1)
c
      real*8 ENRG_Z, ENRG_C, ENRG     ! エネルギー, エンストロフィー
      real*8 cap_H(1:Nx, 0:Nz+1)        ! 熱強制の関数
      real*8 work_time, work_1, work_2, work_3 
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* WORK SPACE FOR RUNGE-KUTTA METHOD
c     W_  means work space
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      REAL*8   WK_psi(Nx, Nz), 
     +         WX_psi(Nx, Nz), 
     +         WR_psi(Nx, Nz), 
     +         WQ_psi(Nx, Nz)
      REAL*8   WK_tmp(Nx, 0:Nz+1), 
     +         WX_tmp(Nx, 0:Nz+1), 
     +         WR_tmp(Nx, 0:Nz+1), 
     +         WQ_tmp(Nx, 0:Nz+1)
c constans for RKG method
      real*8   c1, c2, unity, zero, half, pi
      parameter (
     &            C1=1.0D+0-DSQRT(0.5D+0),
     &            C2=1.0D+0+DSQRT(0.5D+0),
     &            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
      t1_max=1000
      dt=1.d-3                        !時間差分間隔

      open(26, file='energy_x10_z05_RKG.dat', status='unknown')

      cap_G=1.D+0
      cap_E_h=0.0d-2
      cap_K_h=0.0d-2
      cap_E_v=1.0d-2
      cap_K_v=1.0d-2
      gamma=1.D+0
      work_time=dble(0)
!**************************************
      do i=0, Nx+1                             !x座標の設定
           x(i)=dx*dble(i)
      end do                                    
C
      do j=0, Nz+1                             !z座標の設定
           z(j)=dz*dble(j)-Lz
      end do                                    
! initialize
      do j=1,Nz
       do i=1, Nx
         WQ_psi(i,j)=dble(0)
       end do
      end do
      do j=0,Nz+1
       do i=1, Nx
         WQ_tmp(i,j)=dble(0)
       end do
      end do
     
!-initial condition******************************************
      t=0
      call initial_condition(psi, omg, tmp, Nx, Nz, Lx, Lz)
      call energy(psi, omg, tmp, Nx, Nz, ENRG_Z, ENRG_C,
     &           ENRG, gamma, dx, dz)
       write(26,*) work_time, ENRG_Z, ENRG_C, ENRG
 !     call output_data(psi, omg, tmp, t, Nx, Nz, dx, dz, Lx, Lz)
! heating function
      do j=0,Nz+1
       do i=1, Nx
         cap_H(i,j)=(dble(1)+z(j)/pi)*dcos(x(i))
       end do
      end do

      do 100 t=1, t_max
       do 200 t1=1, t1_max
         work_time=work_time+dt
         if (work_time .LT. 100.) then 
           cap_E_v=1.0D+0-0.99*work_time/100.
         else
           cap_E_v=0.01
         end if
         cap_K_v=cap_E_v


         if (work_time .GE. 130. .and. work_time .LT. 130.+dt) then 
           do j=0,Nz+1
             do i=0,Nx+1
               psi(i,j)=psi(i,j)-dble(0.2)*dsin(z(j))
             end do
           end do
          end if
         call bc_psi(psi,Nx,Nz)           

! STEP1
      CALL EVLTN_rkg(omg, psi, tmp, WK_psi, WK_tmp, Nx, Nz, 
     +                half*dt, cap_E_h, cap_K_h, cap_E_v, cap_K_v, 
     +               cap_G, cap_H, gamma, dx, dz)          
      CALL NEXT_RKG(omg, WK_psi, WX_psi, WR_psi, WQ_psi, 
     +      Nx, Nz, unity, unity)
      CALL NEXT_RKG_tmp(tmp, WK_tmp, WX_tmp, WR_tmp, WQ_tmp, 
     +      Nx, Nz, unity, unity)
      call bc_omg(omg, Nx, Nz)
      call bc_tmp(tmp, Nx, Nz)
      call poisson(psi, omg, Nx, Nz, dx, dz)
! STEP2
      CALL EVLTN_rkg(omg, psi, tmp, WK_psi, WK_tmp, Nx, Nz, 
     +                dt, cap_E_h, cap_K_h, cap_E_v, cap_K_v, 
     +               cap_G, cap_H, gamma, dx, dz)          
      CALL NEXT_RKG(omg, WK_psi, WX_psi, WR_psi, WQ_psi, 
     +      Nx, Nz, C1, C1)
      CALL NEXT_RKG_tmp(tmp, WK_tmp, WX_tmp, WR_tmp, WQ_tmp, 
     +      Nx, Nz, C1, C1)
      call bc_omg(omg, Nx, Nz)
      call bc_tmp(tmp, Nx, Nz)
      call poisson(psi, omg, Nx, Nz, dx, dz)
! STEP3
      CALL EVLTN_rkg(omg, psi, tmp, WK_psi, WK_tmp, Nx, Nz, 
     +                dt, cap_E_h, cap_K_h, cap_E_v, cap_K_v, 
     +               cap_G, cap_H, gamma, dx, dz)          
      CALL NEXT_RKG(omg, WK_psi, WX_psi, WR_psi, WQ_psi, 
     +      Nx, Nz, C2, C2)
      CALL NEXT_RKG_tmp(tmp, WK_tmp, WX_tmp, WR_tmp, WQ_tmp, 
     +      Nx, Nz, C2, C2)
      call bc_omg(omg, Nx, Nz)
      call bc_tmp(tmp, Nx, Nz)
      call poisson(psi, omg, Nx, Nz, dx, dz)
! STEP4
      CALL EVLTN_rkg(omg, psi, tmp, WK_psi, WK_tmp, Nx, Nz, 
     +                half*dt, cap_E_h, cap_K_h, cap_E_v, cap_K_v, 
     +               cap_G, cap_H, gamma, dx, dz)          
      CALL NEXT_RKG(omg, WK_psi, WX_psi, WR_psi, WQ_psi, 
     +      Nx, Nz, unity/dble(3), unity)
      CALL NEXT_RKG_tmp(tmp, WK_tmp, WX_tmp, WR_tmp, WQ_tmp, 
     +      Nx, Nz, unity/dble(3), unity)
      call bc_omg(omg, Nx, Nz)
      call bc_tmp(tmp, Nx, Nz)
      call poisson(psi, omg, Nx, Nz, dx, dz)

 200    continue

        call energy(psi, omg, tmp, Nx, Nz, ENRG_Z, ENRG_C,
     &           ENRG, gamma, dx, dz)
        write(26,*) work_time, ENRG_Z, ENRG_C, ENRG
!        call output_data(psi, omg, tmp, t, Nx, Nz, dx, dz, Lx, Lz)
!--------------------------------------------------------------------
!      work_1=dble(0)
!      work_2=dble(0)
!      work_3=dble(0)
!      call arakawa_jacobian(psi, omg, Nx, Nz, dx, dz, adv_omg)
!      call arakawa_jacobian_tmp(tmp, omg, Nx, Nz, dx, dz, adv_tmp)
!      do j=1,Nz
!        do i=1, Nx
!          work_1=work_1+psi(i,j)*adv_omg(i,j)
!          work_2=work_2+omg(i,j)*adv_omg(i,j)
!        end do
!      end do
!      do j=1,Nz
!        do i=1, Nx
!          work_3=work_3+psi(i,j)*adv_tmp(i,j)
!        end do
!      end do
!
!       write(6,*) work_time, work_1, work_2, work_3
!---------------------------------------------------------------------
 100  continue
!
      stop
      end
!------------------------------------------------------------------------------
      subroutine bc_psi(Gdata, Nx, Nz)
! periodic boundary condition for x
! field vanishes at the wall
      integer i, j, k
      integer Nx, Nz                
      real*8  dt                           
!***********● 配列の定義 ●***********
      real*8 Gdata(0:Nx+1,-1:Nz+2)         

! at Wall
      do i=1,Nx
!        Gdata(i,0)=Gdata(i,1)
!        Gdata(i,Nz+1)=Gdata(i,Nz)
        Gdata(i,0)=dble(0)
        Gdata(i,Nz+1)=dble(0)
      end do
!  in Wall 
      do i=1,Nx
!        Gdata(i,0)=Gdata(i,1)
!        Gdata(i,Nz+1)=Gdata(i,Nz)
        Gdata(i,-1)=-Gdata(i,1)
        Gdata(i,Nz+2)=-Gdata(i,Nz)
      end do

! cyclic
      do j=-1,Nz+2
        Gdata(0,j)=Gdata(Nx,j)
        Gdata(Nx+1,j)=Gdata(1,j)
      end do
C
!        Gdata(0,0)=Gdata(Nx,1)
!        Gdata(Nx+1,0)=Gdata(1,1)
!        Gdata(0,Nz+1)=Gdata(Nx,Nz)
!        Gdata(Nx+1,Nz+1)=Gdata(1,Nz)
!
      Return
      End
!------------------------------------------------------------------------------
      subroutine bc_omg(Gdata, Nx, Nz)
! periodic boundary condition for x
! field vanishes at the wall
      integer i, j, k
      integer Nx, Nz                
      real*8  dt                           
!***********● 配列の定義 ●***********
      real*8 Gdata(0:Nx+1,0:Nz+1)         

! at Wall
      do i=1,Nx
        Gdata(i,0)=dble(0)
        Gdata(i,Nz+1)=dble(0)
      end do
! cyclic
      do j=0,Nz+1
        Gdata(0,j)=Gdata(Nx,j)
        Gdata(Nx+1,j)=Gdata(1,j)
      end do

      Return
      End
!------------------------------------------------------------------------------
      subroutine bc_tmp(Gdata, Nx, Nz)
! periodic boundary condition for x
! no-gradient at the wall
      integer i, j, k
      integer Nx, Nz                
!***********● 配列の定義 ●***********
      real*8 Gdata(0:Nx+1,-1:Nz+2)         

! in Wall
      do i=1,Nx
        Gdata(i,-1)=Gdata(i,1)
        Gdata(i,Nz+2)=Gdata(i,Nz)
      end do
! cyclic
      do j=-1,Nz+2
        Gdata(0,j)=Gdata(Nx,j)
        Gdata(Nx+1,j)=Gdata(1,j)
      end do
C

      Return
      End
!------------------------------------------------------------------------------
      subroutine initial_condition(psi, omg, tmp, Nx, Nz, Lx, Lz)
!     初期値を指定するsurboutine
!***********● 変数の定義 ●***********
      integer i, j, k
      integer Nx, Nz                
      real*8  Lx, Lz                         
      real*8  dx, dz                       
!***********● 配列の定義 ●***********
      real*8 psi(0:Nx+1,-1:Nz+2)         
      real*8 omg(0:Nx+1,0:Nz+1)       
      real*8 tmp(0:Nx+1,-1:Nz+2)       
      real*8 work(1:Nx,1:Nz)       
      real*8 x(0:Nx+1), z(0:Nz+1)
!      
      real*8 ENRG_Z, ENRG_C, pi
!************● 変数の値 ●************
      dx=Lx/dble(Nx)                  
      dz=Lz/dble(Nz+1)
      pi=dble(acos(-1.))                  
!**************************************
      do i=0, Nx+1
         x(i)=dble(i)*dx
      end do
      do j=0, Nz+1
         z(j)=dble(j)*dz-Lz
      end do

      do j=0,Nz+1
        do i=0,Nx+1
          psi(i,j)=-dble(0.2)*dsin(z(j))
          omg(i,j)=dble(0)
          tmp(i,j)=dble(0)
        end do
      end do

      call bc_psi(psi, Nx, Nz)
      call bc_tmp(tmp, Nx, Nz)

      call laplacian_1dz(psi, Nx, Nz, dx, dz, omg)

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

      Return
      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 
!*******************************************************************************
      subroutine output_data(psi, omg, tmp, 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       psi(0:Nx+1, -1:Nz+2)
      real*8       omg(0:Nx+1, 0:Nz+1)
      real*8       tmp(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(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,Nz+1
        do 15 i=0,Nx
          write(7,*) dble(i)*dx, dble(j)*dz-Lz, psi(i,j)
          write(8,*) dble(i)*dx, dble(j)*dz-Lz, omg(i,j)
          write(9,*) dble(i)*dx, dble(j)*dz-Lz, tmp(i,j)
   15 continue

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


      return
      end
C
!***************************************
      subroutine arakawa_jacobian(psi, omg, Nx, Nz, dx, dz,
     &                            Jacobian)

      integer   i, j, Nx, Nz
      real*8    dx, dz
      real*8    psi(0:Nx+1,-1:Nz+2), omg(0:Nx+1,0:Nz+1)
      real*8    Jacobian_1, Jacobian_2, Jacobian_3, 
     &          Jacobian(1:Nx,1:Nz),
     &          factor

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

      do j=1,Nz                     !内部領域のみ
        do i=1,Nx
! \psi_x \omg_y - \psi_y \omg_x
         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
! (\psi \omg_y)_x - (\psi \omg_x)_y
         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
! (\omg \psi_x)_y - (\omg \psi_y)_x
         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
!***************************************
      subroutine arakawa_jacobian_tmp(psi, tmp, Nx, Nz, dx, dz,
     &                            Jacobian)

      integer   i, j, Nx, Nz
      real*8    dx, dz
      real*8    psi(0:Nx+1,-1:Nz+2), tmp(0:Nx+1,-1:Nz+2)
      real*8    Jacobian_1, Jacobian_2, Jacobian_3, 
     &          Jacobian(1:Nx,0:Nz+1),
     &          factor

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

      do j=0,Nz+1                     !内部領域のみ
        do i=1,Nx
! \psi_x \tmp_y - \psi_y \tmp_x
         Jacobian_1=((psi(i+1,j)-psi(i-1,j))*
     &   (tmp(i,j+1)-tmp(i,j-1))-(psi(i,j+1)-psi(i,j-1))*
     &   (tmp(i+1,j)-tmp(i-1,j)))*factor
! (\psi \tmp_y)_x - (\psi \tmp_x)_y
         Jacobian_2=((psi(i+1,j))*
     &   (tmp(i+1,j+1)-tmp(i+1,j-1))-(psi(i-1,j))*
     &   (tmp(i-1,j+1)-tmp(i-1,j-1))-(psi(i,j+1))*
     &   (tmp(i+1,j+1)-tmp(i-1,j+1))+(psi(i,j-1))*
     &   (tmp(i+1,j-1)-tmp(i-1,j-1)))*factor
! (\tmp \psi_x)_y - (\tmp \psi_y)_x
         Jacobian_3=((tmp(i,j+1))*
     &   (psi(i+1,j+1)-psi(i-1,j+1))-(tmp(i,j-1))*
     &   (psi(i+1,j-1)-psi(i-1,j-1))-(tmp(i+1,j))*
     &   (psi(i+1,j+1)-psi(i+1,j-1))+(tmp(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
!**************************************************
      subroutine laplacian_1dz(A, Nx, Nz, dx, dz, Lap_A)
! 2nd derivative of a quantity A with respect to z   
      integer   i, j, Nx, Nz
      real*8    dx, dz
      real*8    A(0:Nx+1, -1:Nz+2)
      real*8    Lap_A(0:Nx+1, 0:Nz+1)
     
      do j=1, Nz
        do i=1, Nx
        Lap_A(i,j)=(A(i,j+1)+A(i,j-1)-dble(2)*A(i,j))/(dz**2)
        end do
      end do  

      call bc_omg(Lap_A, Nx, Nz)
      return
      end
C
!**************************************************
      subroutine laplacian_1dz_tmp(A, Nx, Nz, dx, dz, Lap_A)
! 2nd derivative of a quantity A with respect to z   
      integer   i, j, Nx, Nz
      real*8    dx, dz
      real*8    A(0:Nx+1, -1:Nz+2)
      real*8    Lap_A(1:Nx, 0:Nz+1)
     
      do j=0, Nz+1
        do i=1, Nx
        Lap_A(i,j)=(A(i,j+1)+A(i,j-1)-dble(2)*A(i,j))/(dz**2)
        end do
      end do  

      return
      end
C
!**************************************************
      subroutine laplacian_1dz_omg(A, Nx, Nz, dx, dz, Lap_A)
! 2nd derivative of a quantity A with respect to z   
      integer   i, j, Nx, Nz
      real*8    dx, dz
      real*8    A(0:Nx+1, 0:Nz+1)
      real*8    Lap_A(1:Nx, 1:Nz)
     
      do j=1, Nz
        do i=1, Nx
        Lap_A(i,j)=(A(i,j+1)+A(i,j-1)-dble(2)*A(i,j))/(dz**2)
        end do
      end do  

      return
      end
C
!**************************************************
      subroutine laplacian_1dx_omg(A, Nx, Nz, dx, dz, Lap_A)
! 2nd derivative of a quantity A with respect to x   
      integer   i, j, Nx, Nz
      real*8    dx, dz
      real*8    A(0:Nx+1, 0:Nz+1)
      real*8    Lap_A(1:Nx, 1:Nz)
     
      do j=1, Nz
        do i=1, Nx
        Lap_A(i,j)=(A(i+1,j)+A(i-1,j)-dble(2)*A(i,j))/(dx**2)
        end do
      end do  

      return
      end
C
!**************************************************
      subroutine laplacian_1dx_tmp(A, Nx, Nz, dx, dz, Lap_A)
! 2nd derivative of a quantity A with respect to x   
      integer   i, j, Nx, Nz
      real*8    dx, dz
      real*8    A(0:Nx+1, -1:Nz+2)
      real*8    Lap_A(1:Nx, 0:Nz+1)
     
      do j=0, Nz+1
        do i=1, Nx
        Lap_A(i,j)=(A(i+1,j)+A(i-1,j)-dble(2)*A(i,j))/(dx**2)
        end do
      end do  

      return
      end
C
!*******************************************************************************
       subroutine poisson(psi, omg, Nx, Nz, dx, dz)
c
! Solver of the Poisson eqn. psi_{yy}=omega
c
       integer   c, i, j, Nx, Nz, c_max, i_poit, j_poit
       real*8    psi(0:Nx+1, -1:Nz+2), omg(0:Nx+1, 0:Nz+1)
       real*8    old_psi(0:Nx+1, -1:Nz+2)
       real*8    dx, dz
       real*8    w 
       real*8    epsilon, error, pres, value
      

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

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

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

           psi(i,j)=(psi(i,j+1)+psi(i,j-1)-omg(i,j)*(dz)**2)/dble(2)

           error=psi(i,j)-old_psi(i,j)
!           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_psi(psi, Nx, Nz)
!        write(95,*) c, epsilon
           if (epsilon .LT. pres) then
           exit
           end if
C
   90 continue

! 91   write(16,*) c, epsilon

 95   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
!-------------------------------------------------------------------
      subroutine EVLTN_rkg(omg, psi, tmp, WK_psi, WK_tmp, Nx, Nz, 
     +           dt, cap_E_h, cap_K_h, cap_E_v, cap_K_v, cap_G, cap_H, 
     +           gamma, dx, dz)

      integer Nx, Nz, i, j
      real*8  omg(0:Nx+1, 0:Nz+1)
      real*8  psi(0:Nx+1, -1:Nz+2), tmp(0:Nx+1, -1:Nz+2)
      real*8  WK_psi(Nx, Nz), WK_tmp(Nx, 0:Nz+1)
      real*8  dt, cap_E_h, cap_K_h, cap_E_v, cap_K_v, cap_G, gamma
      real*8  dx, dz, cap_H(Nx, 0:Nz+1)
      real*8  adv_omg(Nx, Nz), adv_tmp(Nx, 0:Nz+1)
      real*8  omg_xx(Nx, Nz), tmp_xx(Nx, 0:Nz+1)
      real*8  omg_zz(Nx, Nz), tmp_zz(Nx, 0:Nz+1)

      call arakawa_jacobian(psi, omg, Nx, Nz, dx, dz, adv_omg)
      call arakawa_jacobian_tmp(psi, tmp, Nx, Nz, dx, dz, adv_tmp)
      call laplacian_1dx_omg(omg, Nx, Nz, dx, dz, omg_xx)
      call laplacian_1dx_tmp(tmp, Nx, Nz, dx, dz, tmp_xx)
      call laplacian_1dz_omg(omg, Nx, Nz, dx, dz, omg_zz)
      call laplacian_1dz_tmp(tmp, Nx, Nz, dx, dz, tmp_zz)

      cap_E_h=cap_E_v
      cap_K_h=cap_K_v

      do j=1, Nz
        do i=1, Nx
          WK_psi(i,j)=-adv_omg(i,j)
     &                  +cap_E_h*omg_xx(i,j)
     &                  +cap_E_v*omg_zz(i,j)
     &                  +cap_G*(tmp(i+1,j)-tmp(i-1,j))/(dble(2)*dx)
          WK_psi(i,j)=WK_psi(i,j)*dt
        end do
      end do

      do j=0, Nz+1
        do i=1, Nx
          WK_tmp(i,j)=-adv_tmp(i,j)+cap_H(i,j)
     &                  +cap_K_h*tmp_xx(i,j)
     &                  +cap_K_v*tmp_zz(i,j)
     &                  -gamma*(psi(i+1,j)-psi(i-1,j))/(dble(2)*dx)
          WK_tmp(i,j)=WK_tmp(i,j)*dt
        end do
      end do

      return
      end
!------------------------------------------------------------------
      subroutine NEXT_RKG(A, WK, WX, WR, WQ, 
     +      Nx, Nz, Cx, Cq)

      REAL*8   CX, CQ
      INTEGER  Nx, Nz
      REAL*8   A(0:Nx+1,0:Nz+1)
      REAL*8   WK(Nx, Nz), WX(Nx, Nz), WR(Nx, Nz), WQ(Nx, Nz)
*
      DO j=1, Nz
       DO i=1, Nx
         WX(i,j)=A(i,j)+CX*(WK(i,j)-WQ(i,j))
         WR(i,j)=WX(i,j)-A(i,j)
         WQ(i,j)=WQ(i,j)+3.0D+0*WR(i,j)-CQ*WK(i,j)
         A(i,j)=WX(i,j)
       end do
      end do

      RETURN
      END
!---------------------------------------------------------------------
      subroutine NEXT_RKG_tmp(A, WK, WX, WR, WQ, 
     +      Nx, Nz, Cx, Cq)

      REAL*8   CX, CQ
      INTEGER  Nx, Nz
      REAL*8   A(0:Nx+1,-1:Nz+2)
      REAL*8   WK(Nx, 0:Nz+1), WX(Nx, 0:Nz+1), 
     +         WR(Nx, 0:Nz+1), WQ(Nx, 0:Nz+1)
*
      DO j=0, Nz+1
       DO i=1, Nx
         WX(i,j)=A(i,j)+CX*(WK(i,j)-WQ(i,j))
         WR(i,j)=WX(i,j)-A(i,j)
         WQ(i,j)=WQ(i,j)+3.0D+0*WR(i,j)-CQ*WK(i,j)
         A(i,j)=WX(i,j)
       end do
      end do

      RETURN
      END
!------------------------------------------------------------------
      subroutine NEXT_AB(A, WK, WX, Nx, Nz, C1, C2)

      REAL*8   C1, C2
      INTEGER  Nx, Nz
      REAL*8   A(0:Nx+1,0:Nz+1)
      REAL*8   WK(Nx, Nz), WX(Nx, Nz)
*
      DO j=1, Nz
       DO i=1, Nx
         A(i,j)=A(i,j)+C1*WK(i,j)-C2*WX(i,j)
         WX(i,j)=WK(i,j)
       end do
      end do

      RETURN
      END
c--------------------------------------------------------
      subroutine NEXT_AB_tmp(A, WK, WX, Nx, Nz, C1, C2)

      REAL*8   C1, C2
      INTEGER  Nx, Nz
      REAL*8   A(0:Nx+1,-1:Nz+2)
      REAL*8   WK(Nx, 0:Nz+1), WX(Nx, 0:Nz+1)
*
      DO j=0, Nz+1
       DO i=1, Nx
         A(i,j)=A(i,j)+C1*WK(i,j)-C2*WX(i,j)
         WX(i,j)=WK(i,j)
       end do
      end do

      RETURN
      END

