!----------------------------------------------------------------------
!     Copyright (c) 2004 YAMADA YUKIKO. All rights reserved.
!----------------------------------------------------------------------
!
! Solving equatorial beta plane linear shallow water equations
!
! command: spmfrt shallow-eigen.f90 -lssl2
!
! history: 2004/09/07 YAMADA, Y. 
!          2004/10/25 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 :: mat_m=3*lm-1     ! 行列の大きさの設定(Y)
    complex            :: za(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,mat_m) 
    real               :: xy_H(lm,mat_m) 
    real               :: xy_V(lm,mat_m) 
    real               :: vw(mat_m),ivw(mat_m) ! 作業配列
    integer            :: l,k                  ! 作業配列
    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=-3.0, ymax=3.0  ! 領域範囲(Y 方向) 
!    real, parameter :: ymin=-5.0, ymax=5.0  ! 領域範囲(Y 方向) 
    real, parameter :: del_Y=(ymax-ymin)/(lm-1)   ! 格子幅 (Y 方向)

 !---- 物理パラメター ----
    real(8), parameter :: Grav= 1.0              ! 重力加速度
    real(8), parameter :: Beta= 1.0              ! 赤道βパラメタ
    real(8), parameter :: H0  = 1.0              ! 平均の深さ
!    complex, parameter :: H0  = (2.3, 0.5)
!    real(8), parameter :: Nu  = 0.02             ! 粘性拡散係数
    real(8), parameter :: Nu  = 0.0             ! 粘性拡散係数
    real(8), parameter :: Rd  = 0.0             ! レイリー摩擦

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

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

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


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


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

          zev = - zev + img*zev

          xy_U(1,:) = zev(1,:)
          xy_H(1,:) = zev(2,:)
          xy_V(1,:) = - zev(3,:) /2.0  * img

          xy_U(lm,:) = zev(3*lm-2,:)
          xy_H(lm,:) = zev(3*lm-1,:)
          xy_V(lm,:) = - zev(3*lm-3,:) /2.0 * img

          do l = 2,lm-1
             xy_V(l,:) = - ( zev(3*l-3,:) + zev(3*l,:) ) /2.0 * img
             xy_U(l,:) = zev(3*l-2,:)
             xy_H(l,:) = zev(3*l-1,:)
          end do
       endif
       
    end do
    
    call output_gtool4
    call output_gtool4_close

    stop
    
  contains

 !----------------------- 行列の生成  ----------------------------

    subroutine gyouretsu
 
      za(:,:)=0.0
      
      do l = 2,lm
        
         if (l /= lm) then
            za(3*l-2,3*l  ) = - 1.0/2.0 * Beta * xy_Y(l)
            za(3*l-1,3*l  ) =   Grav *  H0 / del_Y

            za(3*l-3,3*l  ) = - Nu / (del_Y**2)
            za(3*l-2,3*l+1) = - Nu / (del_Y**2)
            za(3*l-1,3*l+2) = - Nu / (del_Y**2)
            za(3*l-2,3*l-2) =   2.0*Nu / (del_Y**2) + Nu * xy_K**2 + Rd
            za(3*l-1,3*l-1) =   2.0*Nu / (del_Y**2) + Nu * xy_K**2 + Rd
         end if

         if (l == lm) then
!            za(3*l-2,3*l-2) =   Nu / (del_Y**2) + Nu * xy_K**2 + Rd
!            za(3*l-1,3*l-1) =   Nu / (del_Y**2) + Nu * xy_K**2 + Rd
            za(3*l-2,3*l-2) =   3.0 * Nu / (del_Y**2) + Nu * xy_K**2 + Rd
            za(3*l-1,3*l-1) =   3.0 * Nu / (del_Y**2) + Nu * xy_K**2 + Rd
         end if

         za(3*l-3,3*l-5) =   Beta/4.0*( xy_Y(l-1)+xy_Y(l) )         
         za(3*l-3,3*l-4) = - 1.0/ del_Y 
         za(3*l-3,3*l-2) =   Beta/4.0*( xy_Y(l-1)+xy_Y(l) )
         za(3*l-3,3*l-1) =   1.0/ del_Y 
         za(3*l-2,3*l-3) = - 1.0/2.0 * Beta * xy_Y(l)
         za(3*l-2,3*l-1) =   img * xy_K 
         za(3*l-1,3*l-3) = - Grav *  H0 / del_Y
         za(3*l-1,3*l-2) =   img * xy_K  * Grav * H0 

         za(3*l-3,3*l-6) = - Nu / (del_Y**2)
         za(3*l-3,3*l-3) =   2.0*Nu / (del_Y**2) + Nu * xy_K**2 + Rd
         za(3*l-2,3*l-5) = - Nu / (del_Y**2)
         za(3*l-1,3*l-4) = - Nu / (del_Y**2)
            
      end do

      za(1,2) =   img * xy_K 
      za(1,3) = - 1.0/2.0 * Beta * xy_Y(1) 
      za(2,1) =   img * xy_K  * Grav * H0
      za(2,3) =  Grav *  H0 / del_Y  

      za(1,4) = - Nu / (del_Y**2)
!      za(1,1) =   Nu / (del_Y**2) + Nu * xy_K**2  + Rd
      za(1,1) =   3.0 * Nu / (del_Y**2) + Nu * xy_K**2  + Rd
      za(2,5) = - Nu / (del_Y**2)
!      za(2,2) =   Nu / (del_Y**2) + Nu * xy_K**2  + Rd
      za(2,2) =   3.0 * Nu / (del_Y**2) + Nu * xy_K**2  + Rd

!          write(*,*) za(5,:)
      
      !---- 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)         ! 東西波数

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

      call HistoryCreate( &                             ! ヒストリー作成
           file='shallow.nc',                              &
           title='equatorial beta plane shallow water model',  &
           source='shallow-eigen.f90',                    &
           institution='GFD_Dennou Club aqua-planet experiment project',  &
           dims=(/'k','n','m'/), dimsizes=(/km,mat_m,lm/),   &
           longnames=(/'wave-number ','eigen-number','y           '/),    &
           units=(/'1','1','1'/),                              &
           origin=0.0, interval=real( dk )                 &
                         ) 
      call HistoryPut('k',cood_k)                              ! 変数出力
      call HistoryPut('n',cood_eig)                           ! 変数出力
      call HistoryPut('m',xy_Y)                           ! 変数出力
      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=(/'m','n'/), & 
           longname='U_keq1', units='1', xtype='double')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='V', dims=(/'m','n'/), & 
           longname='V_keq1', units='1', xtype='double')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='H', dims=(/'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


