!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
! Sample program for gtool_history/gtool4 and ISPACK   2002/01/25 S.Takehiro
! 
! Solving KdV equation
!     d\zeta/dt  = -\zeta d\zeta/dx + d^3\zeta/dx^3
!
program kdv

  use ftrpack
  use gtool_history
  implicit none

 !---- 空間解像度設定 ----
  integer, parameter :: im=128            ! 格子点の設定(X)
  integer, parameter :: km=42             ! 切断波数の設定(X)

 !---- 変数 ----
  real(8)            :: g_zeta(0:im-1)     ! 格子データ
  real(8)            :: e_zeta(-km:km)     ! スペクトルデータ
  real(8)            :: e_zeta0(-km:km)    ! スペクトルデータ
  real(8)            :: e_zeta1(-km:km)    ! スペクトルデータ

 !---- 座標変数など ----
  real(8), parameter :: xlength=3.0

 !---- 時間積分パラメター ----
  real(8), parameter :: dt=1e-6               ! 時間ステップ間隔
  integer, parameter :: nt=5000, ndisp=100    ! 時間積分数, 表示ステップ

 !---- 物理パラメター ----
  real(8), parameter :: x1=xlength/2.0          ! 初期分布 X 座標
  real(8), parameter :: u1=720.0                ! 初期分布の振幅

  real(8), parameter :: x2=xlength/4.0          ! 初期分布 X 座標
  real(8), parameter :: u2=1440.0               ! 初期分布の振幅

  real(8) :: pi                                 ! 円周率
  integer :: i, j, it                           ! DO 変数

  pi = 4.0D0*atan(1.0D0)

  call ftrinitial(im,km,xlength)                ! スペクトル初期化

 !------------------- 初期値設定 ----------------------
  !g_zeta=  cos(pi*g_x)
  g_zeta= u1*sech((g_x-x1)/sqrt(12/u1))**2 &
         + u2*sech((g_x-x2)/sqrt(12/u2))**2 

  e_zeta = e_g(g_zeta)
  e_zeta1 = e_zeta ; e_zeta0 = e_zeta

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

 !------------------- 時間積分 ----------------------
  do it=1,nt
     e_zeta = e_zeta0 + 2*dt * &                   ! Leap frog 時間積分
          ( -e_g(g_e(e_zeta1)*g_e(e_dx_e(e_zeta1))) &
            - e_dx_e(e_dx_e(e_dx_e(e_zeta1)))  )
     e_zeta0=e_zeta1 ; e_zeta1=e_zeta

     if(mod(it,ndisp) .eq. 0)then                    ! 出力
        g_zeta = g_e(e_zeta)
        call output_gtool4
     endif
  enddo

  call output_gtool4_close                           ! ヒストリー後処理
  stop

contains
  function sech(x)
    double precision :: x(0:im-1)
    double precision :: sech(0:im-1)

    sech = 2/(exp(x)+ exp(-x))

  end function sech

  subroutine output_gtool4_init
    call HistoryCreate( &                                  ! ヒストリー作成
           file='kdv.nc', title='KDV equation model',       &
           source='Sample program of gtool_history/gtool4', &
           institution='GFD_Dennou Club davis project',     &
           dims=(/'x','t'/), dimsizes=(/im,0/),             &
           longnames=(/'X-coordinate','time        '/),     &
           units=(/'1','1'/),                               &
           origin=0.0, interval=real(ndisp*dt) )

    call HistoryPut('x',g_x)                                 ! 変数出力

    call HistoryAddVariable( &                                ! 変数定義
           varname='zeta', dims=(/'x','t'/), & 
           longname='vorticity', units='1', xtype='double')
  end subroutine output_gtool4_init

  subroutine output_gtool4
    write(6,*) 'it = ',it
    call HistoryPut('t',real(it*dt))
    call HistoryPut('zeta',g_zeta)
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program kdv
