!----------------------------------------------------------------------
!     Copyright (c) 2004 YAMADA YUKIKO. All rights reserved.
!----------------------------------------------------------------------
!
! Solving equatorial beta plane non-linear shallow water equations
!
! history: 2004/09/07 YAMADA, Y. 
!
!

  program equatorial_beta_shallow_water_eigenvalue
    use gt4_history
    use gt4f90io
    implicit none

 !---- 空間解像度設定 ----
    integer, parameter ::  lm= 41          ! 南北格子 (波数) の設定(Y)
!    integer, parameter ::  lm= 81          ! 南北格子 (波数) の設定(Y)
    integer, parameter ::  km= 40          ! 東西波数の設定(X)
    integer, parameter ::  mm= 2          ! 鉛直波数の設定(X)

 !---- 固有値固有ベクトル ----
    integer, parameter :: mat_m=(3*lm-1)*mm   ! 格子点の設定(Y)
    complex            :: za(mat_m,mat_m) ! 入力: 複素行列
    real               :: ai(mat_m,mat_m) ! 実行列 (実部)
    real               :: ar(mat_m,mat_m) ! 実行列 (虚部)
    real               :: br(mat_m,mat_m) ! 実行列
    integer, parameter :: mode=0          ! 入力: 平衡化: 1=省略する
    complex            :: ze(mat_m)       ! 出力: 固有値
    real               :: eig(km,mat_m)   ! 出力: 固有値
    real               :: growth(km,mat_m)   ! 出力: 固有値
    complex            :: zev(mat_m,mat_m)     ! 出力: 固有ベクトル
    real               :: xy_U(lm,mm,mat_m) 
    real               :: xy_H(lm,mm,mat_m) 
    real               :: xy_V(lm,mm,mat_m) 
    real               :: eta(mm)               ! 加熱モード
    real               :: vw(mat_m),ivw(mat_m) ! 作業配列
    real               :: ip(mat_m),is(mat_m)  ! 作業配列
    integer            :: l,m,k,work           ! 作業配列
    integer            :: work0,work1,work2    ! 作業配列
    integer            :: icon                 ! 出力: コンディションコード
    ! icon=0:エラーなし, icon=10000:N=1である
    ! icon=20000:固有値固有ベクトル求められず
    ! icon=30000: N<1, または K<N であった

 !---- 座標変数など ----
    real                :: xy_K              ! 東西波数
    real                :: xy_Y(lm)          ! 領域範囲(Y 方向) 
    real,parameter      :: dk = 0.1          ! 波数幅
    real, parameter :: ymin=-5.0, ymax=5.0  ! 領域範囲(Y 方向) 
    real, parameter :: del_Y=(ymax-ymin)/(lm-1)   ! 格子幅 (Y 方向)
    real, parameter :: del_Z = 1.0
!    real, parameter :: del_Z = 1.73

 !---- 物理パラメター ----
    real(8), parameter :: Beta= 1.0              ! 赤道βパラメタ
!    real(8), parameter :: Nu  = 0.02             ! 粘性拡散係数
    real(8), parameter :: Nu  = 0.0             ! 粘性拡散係数
    real(8), parameter :: eta1(mm)  = (1.0,1.0)        ! 加熱第一モード
    real(8), parameter :: eta2(mm)  = (-1.0,1.0)       ! 加熱第二モード

 !---- 複素数 ----
    complex, parameter :: img = (0.0,1.0)

    do l = 1,lm
       xy_Y(l) = ymin + (l-1)*del_Y 
    end do

!    eta = 2.0 * ( eta1 + eta2 )
    eta(1) =  1.0 * 0.5
    eta(2) = - 2.0 * 0.5
!    eta = 0.0

 !------------------------------------------------------------
    call output_gtool4_init                       ! ヒストリー初期化

    ! xy_K でループを回して出力

    do k = 1, km
       
       xy_K  = k * dk

       !       a*x = b*\lambda*x     という行列式を
       ! (b-1)*a*x = \lambda*x       と変形して, 固有値問題を解く
       ! 行列式 a_real, a_img, b の設定      

       call gyouretsu

       eig(k,:) = - ze * img        ! 分散関係
       growth(k,:) = - ze           ! 成長率
       write(*,*)  ze
       
       if (xy_K == 0.5) then 

          zev = - zev + img*zev

          do m = 1,mm

             work = (m-1)*(3*lm-1)

             xy_U(1,m,:) =   zev(work+1,:)
             xy_H(1,m,:) =   zev(work+2,:)
             xy_V(1,m,:) = - zev(work+3,:) /2.0  * img
             
             xy_U(lm,m,:) =   zev(work+3*lm-2,:)
             xy_H(lm,m,:) =   zev(work+3*lm-1,:)
             xy_V(lm,m,:) = - zev(work+3*lm-3,:) /2.0 * img
             
             do l = 2,lm-1
                xy_V(l,m,:) = - ( zev(work+3*l-3,:) + zev(work+3*l,:) )/2.0*img
                xy_U(l,m,:) =   zev(work+3*l-2,:)
                xy_H(l,m,:) =   zev(work+3*l-1,:)
             end do

          end do

       endif
       
    end do

       
    call output_gtool4
    call output_gtool4_close

    stop
    
  contains

 !----------------------- 行列の生成  ----------------------------
    subroutine gyouretsu

      ar(:,:) = 0.0
      ai(:,:) = 0.0
      br(:,:) = 0.0

      !---------------------
      ! ar, ai の計算
      do m = 1,mm
         
         do l = 2,lm

            work = (m-1)*(3*lm-1) + 3*l
            
            if (l /= lm) then
               ar(work-2,work  ) = - 1.0/2.0 * Beta * xy_Y(l)
               ar(work-1,work  ) =   1.0 / del_Y + eta(mm) / (del_Z**2)

!               ar(work-1,work-1) =   eta(mm) / (del_Z**2)
            end if

            ar(work-3,work-5) =   Beta/4.0*( xy_Y(l-1)+xy_Y(l) )         
            ar(work-3,work-4) = - 1.0/ del_Y 
            ar(work-3,work-2) =   Beta/4.0*( xy_Y(l-1)+xy_Y(l) )
            ar(work-3,work-1) =   1.0/ del_Y 
            ar(work-2,work-3) = - 1.0/2.0 * Beta * xy_Y(l)
            ai(work-2,work-1) =   xy_K 
            ar(work-1,work-3) = - 1.0 / del_Y + eta(mm) / (del_Z**2)
            ai(work-1,work-2) =   xy_K        + eta(mm) / (del_Z**2)

!            ar(work-1,work-1) =   eta(mm) / (del_Z**2)
         end do
         
         work = (m-1)*(3*lm-1)
         
         ai(work+1,work+2) =   xy_K 
         ar(work+1,work+3) = - 1.0/2.0 * Beta * xy_Y(1) 
         ai(work+2,work+1) =   xy_K          + eta(mm) / (del_Z**2)
         ar(work+2,work+3) =   1.0 / del_Y   + eta(mm) / (del_Z**2)

!         ar(work+2,work+2) =   eta(mm) / (del_Z**2)

      end do

      ! 拡散項の計算
      do m = 1,mm

         if (m==1) then
            
            work0 = (m-1)*(3*lm-1)
            work2 =  m   *(3*lm-1)
         
         
            ar(work0+1,work0+4) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
            ar(work0+1,work0+1) =   2.0 * ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
            ar(work0+2,work0+5) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
            ar(work0+2,work0+2) =   2.0 * ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
         
            ar(work0+1,work2+4) =   Nu / (del_Y**2) / (del_Z**2)
            ar(work0+1,work2+1) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
            ar(work0+2,work2+5) =   Nu / (del_Y**2) / (del_Z**2)
            ar(work0+2,work2+2) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)

            do l = 2,lm
               
               work0 = (m-1)*(3*lm-1) + 3*l
               work2 =  m   *(3*lm-1) + 3*l
               
               if (l /= lm) then
                  ar(work0-3,work0  ) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work0+1) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-1,work0+2) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work0-2) =   2.0 * (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)
                  ar(work0-1,work0-1) =   2.0 * (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)

                  ar(work0-3,work2  ) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work2+1) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-1,work2+2) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work2-2) = - (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)
                  ar(work0-1,work2-1) = - (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)

               end if
            
               if (l == lm) then
                  ar(work0-2,work0-2) =   2.0 * ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)
                  ar(work0-1,work0-1) =   2.0 * ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)

                  ar(work0-2,work2-2) = - ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)
                  ar(work0-1,work2-1) = - ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)

               end if

               ar(work0-3,work0-6) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
               ar(work0-3,work0-3) =   2.0 * ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                    & / (del_Z**2)
               ar(work0-2,work0-5) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
               ar(work0-1,work0-4) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)

               ar(work0-3,work2-6) =   Nu / (del_Y**2) / (del_Z**2)
               ar(work0-3,work2-3) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                    & / (del_Z**2)
               ar(work0-2,work2-5) =   Nu / (del_Y**2) / (del_Z**2)
               ar(work0-1,work2-4) =   Nu / (del_Y**2) / (del_Z**2)

            end do

         elseif (m==mm) then
            
            work0 = (m-1)*(3*lm-1)
            work1 = (m-2)*(3*lm-1)

            ar(work0+1,work0+4) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
            ar(work0+1,work0+1) =   2.0 * ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
            ar(work0+2,work0+5) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
            ar(work0+2,work0+2) =   2.0 * ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
         
         
            ar(work0+1,work1+4) =   Nu / (del_Y**2) / (del_Z**2)
            ar(work0+1,work1+1) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
            ar(work0+2,work1+5) =   Nu / (del_Y**2) / (del_Z**2)
            ar(work0+2,work1+2) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)

            do l = 2,lm
               
               work0 = (m-1)*(3*lm-1) + 3*l
               work1 = (m-2)*(3*lm-1) + 3*l
               
               if (l /= lm) then
                  ar(work0-3,work0  ) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work0+1) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-1,work0+2) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work0-2) =   2.0 * (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)
                  ar(work0-1,work0-1) =   2.0 * (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)

                  ar(work0-3,work1  ) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work1+1) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-1,work1+2) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work1-2) = - (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)
                  ar(work0-1,work1-1) = - (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)

               end if
            
               if (l == lm) then
                  ar(work0-2,work0-2) =   2.0 * ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)
                  ar(work0-1,work0-1) =   2.0 * ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)

                  ar(work0-2,work1-2) = - ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)
                  ar(work0-1,work1-1) = - ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)

               end if

               ar(work0-3,work0-6) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
               ar(work0-3,work0-3) =   2.0 * ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                    & / (del_Z**2)
               ar(work0-2,work0-5) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
               ar(work0-1,work0-4) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)

               ar(work0-3,work1-6) =   Nu / (del_Y**2) / (del_Z**2)
               ar(work0-3,work1-3) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                    & / (del_Z**2)
               ar(work0-2,work1-5) =   Nu / (del_Y**2) / (del_Z**2)
               ar(work0-1,work1-4) =   Nu / (del_Y**2) / (del_Z**2)

            end do

         else

            work0 = (m-1)*(3*lm-1)
            work1 = (m-2)*(3*lm-1)
            work2 =  m   *(3*lm-1)
         
            ar(work0+1,work0+4) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
            ar(work0+1,work0+1) =   2.0 * ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
            ar(work0+2,work0+5) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
            ar(work0+2,work0+2) =   2.0 * ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
         
         
            ar(work0+1,work1+4) =   Nu / (del_Y**2) / (del_Z**2)
            ar(work0+1,work1+1) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
            ar(work0+2,work1+5) =   Nu / (del_Y**2) / (del_Z**2)
            ar(work0+2,work1+2) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)

            ar(work0+1,work2+4) =   Nu / (del_Y**2) / (del_Z**2)
            ar(work0+1,work2+1) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)
            ar(work0+2,work2+5) =   Nu / (del_Y**2) / (del_Z**2)
            ar(work0+2,work2+2) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                 & / (del_Z**2)

            do l = 2,lm
               
               work0 = (m-1)*(3*lm-1) + 3*l
               work1 = (m-2)*(3*lm-1) + 3*l
               work2 =  m   *(3*lm-1) + 3*l
               
               if (l /= lm) then
                  ar(work0-3,work0  ) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work0+1) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-1,work0+2) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work0-2) =   2.0 * (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)
                  ar(work0-1,work0-1) =   2.0 * (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)

                  ar(work0-3,work1  ) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work1+1) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-1,work1+2) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work1-2) = - (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)
                  ar(work0-1,work1-1) = - (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)

                  ar(work0-3,work2  ) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work2+1) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-1,work2+2) =   Nu / (del_Y**2) / (del_Z**2)
                  ar(work0-2,work2-2) = - (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)
                  ar(work0-1,work2-1) = - (2.0*Nu/(del_Y**2) + Nu*(xy_K**2))&
                       & / (del_Z**2)

               end if
            
               if (l == lm) then
                  ar(work0-2,work0-2) =   2.0 * ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)
                  ar(work0-1,work0-1) =   2.0 * ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)

                  ar(work0-2,work1-2) = - ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)
                  ar(work0-1,work1-1) = - ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)

                  ar(work0-2,work2-2) = - ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)
                  ar(work0-1,work2-1) = - ( Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                       & / (del_Z**2)

               end if

               ar(work0-3,work0-6) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
               ar(work0-3,work0-3) =   2.0 * ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                    & / (del_Z**2)
               ar(work0-2,work0-5) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)
               ar(work0-1,work0-4) = - 2.0 * Nu / (del_Y**2) / (del_Z**2)

               ar(work0-3,work1-6) =   Nu / (del_Y**2) / (del_Z**2)
               ar(work0-3,work1-3) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                    & / (del_Z**2)
               ar(work0-2,work1-5) =   Nu / (del_Y**2) / (del_Z**2)
               ar(work0-1,work1-4) =   Nu / (del_Y**2) / (del_Z**2)

               ar(work0-3,work2-6) =   Nu / (del_Y**2) / (del_Z**2)
               ar(work0-3,work2-3) = - ( 2.0*Nu/(del_Y**2) + Nu*(xy_K**2) ) &
                    & / (del_Z**2)
               ar(work0-2,work2-5) =   Nu / (del_Y**2) / (del_Z**2)
               ar(work0-1,work2-4) =   Nu / (del_Y**2) / (del_Z**2)

            end do

         end if

      end do

      !---------------------
      ! br の計算
      !  m = 2,3
      do m = 1,mm

         if (m==1) then

            br(1,1         ) =   1.0
!            br(2,2         ) =   1.0 /(del_Z**2)
            br(2,2         ) =   2.0 /(del_Z**2)
            br(2,(3*lm-1)+2) = - 1.0 /(del_Z**2)
            do l = 2,lm
               br(3*l-3,3*l-3)          =   1.0
               br(3*l-2,3*l-2)          =   1.0
!               br(3*l-1,3*l-1)          =   1.0 / (del_Z**2)
               br(3*l-1,3*l-1)          =   2.0 / (del_Z**2)
               br(3*l-1,(3*lm-1)+3*l-1) = - 1.0 / (del_Z**2)
            end do

         elseif (m==mm) then
            
            work0 = (m-1)*(3*lm-1)
            work1 = (m-2)*(3*lm-1)
            
            br(work0+1,work0+1) =   1.0
!            br(work0+2,work0+2) =   1.0 /(del_Z**2)
            br(work0+2,work0+2) =   2.0 /(del_Z**2)
            br(work0+2,work1+2) = - 1.0 /(del_Z**2)

            do l = 2,lm
               
               work0 = (m-1)*(3*lm-1) + 3*l
               work1 = (m-2)*(3*lm-1) + 3*l
               
               br(work0-3,work0-3) =   1.0
               br(work0-2,work0-2) =   1.0
!               br(work0-1,work0-1) =   1.0 /(del_Z**2)
               br(work0-1,work0-1) =   2.0 /(del_Z**2)
               br(work0-1,work1-1) = - 1.0 /(del_Z**2)
               
            end do

         else

            write(*,*) m

            work0 = (m-1)*(3*lm-1)
            work1 = (m-2)*(3*lm-1)
            work2 =  m   *(3*lm-1)
            
            br(work0+1,work0+1) =   1.0
            br(work0+2,work0+2) =   2.0 /(del_Z**2)
            br(work0+2,work1+2) = - 1.0 /(del_Z**2)
            br(work0+2,work2+2) = - 1.0 /(del_Z**2)
            
            do l = 2,lm
               
               work0 = (m-1)*(3*lm-1) + 3*l
               work1 = (m-2)*(3*lm-1) + 3*l
               work2 =  m   *(3*lm-1) + 3*l
               
               br(work0-3,work0-3) =   1.0
               br(work0-2,work0-2) =   1.0
               br(work0-1,work0-1) =   2.0 /(del_Z**2)
               br(work0-1,work1-1) = - 1.0 /(del_Z**2)
               br(work0-1,work2-1) = - 1.0 /(del_Z**2)
               
            end do

         endif

      end do

      ! ---------------
      ! (b-1)*a の計算
      call ALU(br,mat_m,mat_m,0.0,ip,is,vw,icon) ! LU 分解
      call LUIV(br,mat_m,mat_m,ip,icon)          ! 逆行列
      ! 行列のかけ算
      call MGGM(br,mat_m,ar,mat_m,ar,mat_m,mat_m,mat_m,mat_m,vw,icon)
      call MGGM(br,mat_m,ai,mat_m,ai,mat_m,mat_m,mat_m,mat_m,vw,icon)
      
      za = ar + img*ai
      
      !---- SSL2 で複素行列の固有値固有ベクトルを求める ----
      
      CALL CEIG2 (za,mat_m,mat_m,mode,ze,zev,vw,ivw,icon)
       
    end subroutine gyouretsu


 !----------------------- gtool4 設定 ----------------------------
    subroutine output_gtool4_init

      real                :: cood_eig(mat_m)    ! 固有関数の数
      real                :: cood_k(km)         ! 東西波数
      real                :: cood_z(mm)         ! 鉛直層

      do l = 1,mat_m
         cood_eig(l) = l
      end do
      do m = 1,mm
         cood_z(m) = m
      end do
      do k = 1,km
         cood_k(k) = k * dk
      end do

      call HistoryCreate( &                             ! ヒストリー作成
           file='wavecisk.nc',                                 &
           title='wave-cisk model',                        &
           source='wavecisk-eigen.f90',                    &
           institution='GFD_Dennou Club aqua-planet experiment project',  &
           dims=(/'k','n','l','m'/), dimsizes=(/km,mat_m,lm,mm/),   &
           longnames=(/'wave-number ','eigen-number','y           ','z           '/),    &
           units=(/'1','1','1','1'/),                              &
           origin=0.0, interval=real( dk )                 &
                         ) 

      call HistoryPut('k',cood_k)                              ! 変数出力
      call HistoryPut('n',cood_eig)                           ! 変数出力
      call HistoryPut('l',xy_Y)                           ! 変数出力
      call HistoryPut('m',cood_z)                           ! 変数出力
      call HistoryAddVariable( &                                 ! 変数定義
           varname='eigen_val', dims=(/'k','n'/), & 
           longname='eigen_value', units='1', xtype='single')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='growth', dims=(/'k','n'/), & 
           longname='growth_rate', units='1', xtype='single')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='U', dims=(/'l','m','n'/), & 
           longname='U_keq1', units='1', xtype='double')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='V', dims=(/'l','m','n'/), & 
           longname='V_keq1', units='1', xtype='double')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='H', dims=(/'l','m','n'/), & 
           longname='H_keq1', units='1', xtype='double')
    end subroutine output_gtool4_init
    
    subroutine output_gtool4
      call HistoryPut('eigen_val',eig )
      call HistoryPut('growth',growth )
      call HistoryPut('U',xy_U )
      call HistoryPut('V',xy_V )
      call HistoryPut('H',xy_H )
    end subroutine output_gtool4
    
    subroutine output_gtool4_close
      call HistoryClose
    end subroutine output_gtool4_close

 end program equatorial_beta_shallow_water_eigenvalue
  





