! Sample program for gtool_history/gtool5 and deepconv
! 
! Solving 2-D non-divergent (barotropic) fluid system
!     d\zeta/dt  = \nu\nabla^2\zeta
!
program diffuse

  !---- モジュール読み込み ----  
  use dc_types,      only: STRING, DP
  use gtool_historyauto, only: HistoryAutoPut, HistoryAutoAddVariable
  ! MPI
  use mpi_wrapper, only : MPIWrapperInit, MPIWrapperFinalize
  ! 引数処理
  use argset,     only: argset_init
  ! 格子点
  use gridset,    only: gridset_init, imin, imax, jmin, jmax, kmin, kmax, nx, ny, nz
  ! 軸
  use axesset,    only: axesset_init, x_X, z_Z
  ! ファイル情報
  use fileset,    only: fileset_init
  ! 時刻
  use timeset,    only: timeset_init, DelTime, EndTime, StartTime
  ! 微分演算
  use xyz_deriv_module, only: pyz_dx_xyz, xqz_dy_xyz, xyr_dz_xyz, &
    &                         xyz_dx_pyz, xyz_dy_xqz, xyz_dz_xyr
  ! 境界条件
  use xyz_bc_module, only: BoundaryXCyc_xyz, BoundaryZCyc_xyz          
  ! ファイル出力
  use HistoryFileIO, only: HistoryFileAutoOpen, HistoryFileAutoClose

  !---- 暗黙の型宣言禁止 ----  
  implicit none

  !---- 変数 ----
  real(DP), allocatable :: xyz_ZetaN(:,:,:)
  real(DP), allocatable :: xyz_ZetaA(:,:,:)
  real(DP)              :: TimeA, TimeN
  
 !---- 物理パラメター ----
  real(DP), parameter :: nu=1.0         ! 粘性係数
  real(DP), parameter :: sigma=0.1      ! 初期分布の大きさ
  real(DP), parameter :: x1=5.0d-1      ! 初期分布 X 座標
  real(DP), parameter :: z1=5.0d-1      ! 初期分布 Z 座標

  !---- 入力ファイル -----------
  character(STRING) :: cfgfile

  !---- 作業用 -----------
  integer :: i, k

 !---------------- 初期化 ---------------------
  call MPIWrapperInit
  call argset_init( cfgfile )  !(out)
  call fileset_init( cfgfile ) !(in)
  call timeset_init( cfgfile ) !(in)
  call gridset_init( cfgfile ) !(in)
  call axesset_init( cfgfile ) !(in)
  call HistoryFileAutoOpen (cfgfile) !(in)

  call HistoryAutoAddVariable(  &
    & varname='zeta',           &  
    & dims=(/'x','y','z','t'/), &
    & longname='vorticity',     &
    & units=' ',                &        
    & xtype='double' )

  TimeN = StartTime
  TimeA = StartTime + DelTime

 !---------------- 変数の allocate ---------------------
  allocate(xyz_ZetaN(imin:imax,jmin:jmax,kmin:kmax))
  allocate(xyz_ZetaA(imin:imax,jmin:jmax,kmin:kmax))

 !------------------- 初期値設定 ----------------------
  do k = imin, imax
    do i = kmin, kmax
      xyz_ZetaN(i,:,k) = dexp(-((x_X(i)-x1)**2 + (z_Z(k)-z1)**2)/ (2*sigma**2))
    end do
  end do

 !------------------- 時間積分 ----------------------
  do while (TimeN <= EndTime)    
    
    xyz_ZetaA = xyz_ZetaN                 &
      &        + DelTime * nu * (            &
      &              xyz_dx_pyz(pyz_dx_xyz(xyz_ZetaN)) &
      &            + xyz_dy_xqz(xqz_dy_xyz(xyz_ZetaN)) &
      &            + xyz_dz_xyr(xyr_dz_xyz(xyz_ZetaN)) &
      &          )                                   ! Euler 法による時間積分
    call BoundaryXcyc_xyz( xyz_ZetaA ) 
    call BoundaryZCyc_xyz( xyz_ZetaA ) 

    ! ループ処理
    xyz_ZetaN = xyz_ZetaA
    TimeA = TimeA + DelTime
    TimeN = TimeA

    ! 出力
    call HistoryAutoPut(TimeN, 'zeta', xyz_ZetaN(1:nx,1:ny,1:nz))    
  end do

  !------------------- 終了処理 ----------------------
  call HistoryFileAutoClose
  call MPIWrapperFinalize
  stop

end program diffuse
