!--------------------------------------------------------------------------
!  Copyright (C) 2017 SPMODEL Development Group. All rights reserved.
!--------------------------------------------------------------------------
!
!Title(表題): Two-dimensional barotropic model : vortex Rossby waves
!             (渦ロスビー波 : 2次元順圧モデル)
!
!       Re-calculations of Montgomery and Kallenbach (1997)
!
!       The time integration is performed with the 2th order 
!       Runge-Kutta scheme. 
!
!History(履歴):
!   2017/08/27  竹広真一, 新規作成
!
! The governing equation(支配方程式): 
!
!      (d/dt + V/s d/d\phi)\zeta - (1/s)(d\eta/ds)(d\psi/d\phi)=0
!      \eta = f+(1/r)d(rV)/d\phi
!      \zeta = \nabla^2 \psi
!
program baro_VorRossy

!== モジュール引用宣言 ================================================
  use w_plane_module 
  use gtool_history, only : GT_HISTORY, HistoryCreate, HistoryPut, HistoryGet, &
                            HistoryAddVariable, HistoryClose, HistoryAddAttr
  use dc_message,  only : MessageNotify
  use dc_string, only : StoA
  use dc_types
  use dc_args
  implicit none

!== 宣言部 ============================================================

  !---- 変数(格子点データ) ----
  real(8), allocatable :: ps_VPhi(:,:)        ! 基本場速度方位角成分
  real(8), allocatable :: ps_Eta(:,:)         ! 基本場渦度
  real(8), allocatable :: ps_SInvDEtaDS(:,:)  ! 基本場渦度経度
  
  real(8), allocatable :: ps_Vor(:,:)         ! 渦度擾乱
  real(8), allocatable :: ps_Vor0(:,:)        ! 渦度擾乱
  real(8), allocatable :: ps_StrFunc(:,:)     ! 流線関数擾乱

  !---- 変数(スペクトルデータ) ----
  real(8), allocatable :: w_Vor(:)            ! 渦度擾乱
  real(8), allocatable :: w_StrFunc(:)        ! 流線関数擾乱


  !---- NAMELIST 変数 ----
  ! NAMELIST 入力用テンポラリファイル名(名前変更禁止)
!!$  type(ARGS) :: arg                          ! コマンド引数解析用
!!$  character(STRING), pointer :: argv(:) => null()
!!$  character(STRING)  :: nmlfile              ! NAMELIST 設定ファイル名

  !-- 格子点・スペクトル --
  integer :: nm=85                              ! X 方向切断全波数
  integer :: im=256                             ! 方位角方向格子点数
  integer :: jm=128                             ! 動径方向格子点数
  real(8) :: Radius=5.0D-1                      ! 変換パラメター(球半径)
!!$  integer :: im=128                             ! 方位角方向格子点数
!!$  integer :: jm=64                              ! 動径方向格子点数
!!$  real(8) :: Radius=5.0D-1                      ! 変換パラメター(球半径)
  namelist /gridset/ nm, im, jm, Radius
  !-- 物理パラメター --
  real(8) :: Alpha
  real(8) :: fcori = 0.0D0
  ! -- 初期値 --
  character(len=100) :: initial_file=''          ! 初期値データファイル名
  ! (空なら内部で初期値を計算)
  real(8)            :: initial_time=0.0d0       ! 初期時刻
  namelist /initial/ initial_file, initial_time  ! 
  ! -- 時間積分 --
  real(8) :: delta_t=1.0d-3                      ! 時間積分刻み
  integer :: nstep=60000                         ! 時間積分ステップ数
!!$  real(8) :: delta_t=6.0d-3                      ! 時間積分刻み
!!$  integer :: nstep=1500                          ! 時間積分ステップ数
  namelist /tint/    delta_t, nstep              ! 
  ! -- ヒストリー出力 --
  character(len=100) :: hst_file='inf2dbaro_VorRossby.nc' ! ヒストリーファイル名
  character(len=100) :: title = &                ! タイトル
       'Linear Vortex Rossby waves'
  integer :: hst_intstep=2000                    ! ヒストリー出力間隔ステップ数
!!$  integer :: hst_intstep=50                     ! ヒストリー出力間隔ステップ数
  namelist /history/    hst_file, title, hst_intstep

  !---- 作業変数 ----
  integer :: it=0                            ! 時間ステップ
  real(8) :: time                            ! モデル内時間
  real(8) :: PI                              ! 円周率

  !---------------- NAMELIST 読み込み ---------------------
!!$  call Debug(arg) ; call Help(arg) ; call Strict(arg)
!!$  call Get(arg, argv)
!!$  if ( size(argv) .le. 0 ) then
!!$     call MessageNotify('M','main',&
!!$          'Usage: plbaro-beta_rn4expcs_betadirft.out [namelist file]')
!!$     call MessageNotify('E','main','There is no argument. &
!!$          & Please set the namelist file name after the execution.')
!!$  else
!!$     nmlfile=argv(1)
!!$     call MessageNotify('M','main','Namelist file is '//trim(nmlfile))
!!$  endif
!!$  deallocate(argv)
!!$  call Close(arg)
!!$ 
!!$  write(6,nml=message) 
!!$  open(10,file=nmlfile,status='OLD')
!!$  read(10,nml=message) ; write(6,nml=message) ; close(10)
!!$
!!$  if (verbose) write(6,nml=gridset) 
!!$  open(10,file=nmlfile,status='OLD')
!!$  read(10,nml=gridset) ; write(6,nml=gridset) ; close(10)

  PI = Atan(1.0D0)*4.0D0
  
  !------------------ 変数の割り付け ---------------------
  allocate(ps_VPhi(0:im-1,jm),ps_Eta(0:im-1,jm),ps_SInvDEtaDS(0:im-1,jm))
  allocate(ps_Vor(0:im-1,jm),ps_Vor0(0:im-1,jm),ps_StrFunc(0:im-1,jm))

  allocate(w_Vor((nm+1)**2),w_StrFunc((nm+1)**2))

  !------------------ 座標値の設定 -----------------------
  call w_plane_initial(nm,im,jm,Radius)

  !------------------- 基本場設定 ----------------------
  ps_VPhi = 2*ps_S/(1+ps_S**2)
  ps_Eta = fcori + ps_DS_w(w_ps(ps_S*ps_VPhi))/ps_S

  ps_SInvDEtaDS = ps_DS_w(w_ps(ps_Eta))/ps_S
  
  !------------------- 初期値設定 ----------------------
  time = initial_time

  Alpha = 1.0D0/(6*3**6)

  ! 波数 1
  ps_Vor = (6*Alpha)**(1.0D0/6)*ps_S*exp(1.0D0/6-Alpha*ps_S**6)*cos(ps_Phi)

  ! 波数 2
  !ps_Vor = ps_S**2*exp(1.0D0/2*(1-ps_S**4))*cos(2*ps_Phi)

  w_Vor      = w_ps(ps_Vor)
  w_Strfunc = w_LaplaInv_ps(ps_Vor)
  ps_Strfunc  = ps_w(w_Strfunc)
 
  !---------------- 時間積分(2nd order Runge-Kutta scheme) -----------------
  call output_history_init

  if ( initial_file == '' ) call output_history    ! 内部で与えた初期値は出力

  do it=1,nstep
     time = initial_time + it * delta_t
     
     ps_Vor0 = ps_Vor
     
     ! 1 段目
     ps_Vor = ps_Vor0 &
            - (  (ps_VPhi/ps_S)*ps_DPhi_w(w_Vor) &
                - ps_SInvDEtaDS*ps_DPhi_w(w_Strfunc) ) * delta_t/2
     
     ! 渦度から流線関数を計算
     w_Strfunc = w_LaplaInv_ps(ps_Vor)
     w_Vor = w_ps(ps_Vor)

     ps_Vor = ps_Vor0 &
            - (  (ps_VPhi/ps_S)*ps_DPhi_w(w_Vor) &
                - ps_SInvDEtaDS*ps_DPhi_w(w_Strfunc) ) * delta_t
     
     ! 渦度から流線関数を計算
     w_Strfunc = w_LaplaInv_ps(ps_Vor)
     w_Vor = w_ps(ps_Vor)

     if(mod(it,hst_intstep) .eq. 0)then                    ! ヒストリー出力
        call output_history
     endif

  end do
  
  call output_history_close

contains
! 以上 メインプログラム 
!-----------------------------------------------------------------------------
! 以下 サブルーチン

!=========================== ヒストリー出力 ============================
  !
  ! ヒストリー出力初期化
  !
  subroutine output_history_init

    call HistoryCreate( &
         file=trim(hst_file), &
         title=trim(title), &
         source='baro_VorRossby.f90 (2017/08/27)', &
         institution='GFD_Dennou Club SPMODEL project',     &
         dims=StoA('phi','rad','t'), dimsizes=(/im,jm,0/),&
         longnames=StoA('Azimuth','Radius','time'),&
         units=(/'1','1','1'/),   &
         origin=real(time), interval=real(hst_intstep*delta_t), &
         xtypes=(/'real'/))

    !---- 座標変数定義, 出力 ----
    call HistoryPut('phi',p_Phi*180.0/PI )                  ! 変数出力
    call HistoryAddattr('phi','topology','circular')        ! 周期属性
    call HistoryAddattr('phi','modulo', 360.0)              ! 周期属性
    call HistoryPut('rad',s_S )                             ! 変数出力

    !---- 物理変数定義 ----
    call HistoryAddVariable( &                                ! 変数定義
           varname='VPhi', dims=StoA('rad'), & 
           longname='Basic azimuthal velocity', units='1', xtype='double')
    call HistoryPut('VPhi',ps_VPhi(0,:))
    call HistoryAddVariable( &                                ! 変数定義
           varname='Eta', dims=StoA('rad'), & 
           longname='Basic absolute vorticity', units='1', xtype='double')
    call HistoryPut('Eta',ps_Eta(0,:))    
    call HistoryAddVariable( &                                ! 変数定義
           varname='Beta', dims=StoA('rad'), & 
           longname='Basic absolute vorticity gradient', units='1', xtype='double')
    call HistoryPut('Beta',ps_SInvDEtaDS(0,:))

    call HistoryAddVariable( &                                ! 変数定義
           varname='vor', dims=StoA('phi','rad','t'), & 
           longname='Vorticity', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='strfunc', dims=StoA('phi','rad','t'), & 
           longname='Stream function', units='1', xtype='double')


    !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ----
    call HistoryAddAttr('phi','+delta_t', delta_t )

  end subroutine output_history_init

  !
  ! ヒストリー出力
  !
  subroutine output_history
    write(6,*) ' History file output at it = ',it, '  time = ', time
    call HistoryPut('t',real(time))

    !---- 物理変数出力 ----
    call HistoryPut('vor',ps_Vor)
    call HistoryPut('strfunc',ps_w(w_Strfunc))

  end subroutine output_history

  !
  ! ヒストリー出力終了
  !
  subroutine output_history_close
    call HistoryClose
  end subroutine output_history_close

end program baro_VorRossy
