!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  T2PACK/ISPACK-F90 テストプログラム
!
!履歴  2002/01/25  竹広真一
!
program t2diffusion

  use t2pack
  use gtool_history
  implicit none

 !---- 空間解像度設定 ----
  integer, parameter  :: im=32, jm=32       ! 格子点の設定(X,Y)
  integer, parameter  :: km=10, lm=20       ! 切断波数の設定(X,Y)

 !---- 座標変数など ----
  real(8), parameter  :: xmin=0,  xmax=1    ! X 座標範囲
  real(8), parameter  :: ymin=0, ymax=1     ! Y 座標範囲

 !---- 変数 ----
  real(8), dimension(0:jm,0:im-1) :: gg_zeta     ! 格子点データ
  real(8), dimension(-km:km,0:lm) :: et_zeta     ! スペクトルデータ
  real(8), dimension(0:im-1,2)    :: g_bvalues   ! 境界値
  real(8), dimension(-km:km,2)    :: e_bvalues   ! 境界値

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

 !---- 物理パラメター ----
  !character(len=2), parameter :: bcond='DD'     ! 境界条件
  !character(len=2), parameter :: bcond='ND'     ! 境界条件
  character(len=2), parameter :: bcond='DN'     ! 境界条件
  !character(len=2), parameter :: bcond='NN'     ! 境界条件
  real(8), parameter :: nu=1.0                  ! 粘性係数
  real(8), parameter :: x1=(xmin+xmax)/2.0      ! 初期分布 X 座標
  real(8), parameter :: y1=(xmin+xmax)/2.0      ! 初期分布 Y 座標
  real(8), parameter :: sigma=0.1               ! 初期分布の大きさ

 !---- その他 ----
  real(8), parameter  :: pi=3.1415926535897932385D0
  integer :: it


 !---------------- 座標値の設定 ---------------------
  call t2initial(im,jm,km,lm,xmin,xmax,ymin,ymax)

 !------------------- 初期値設定 ----------------------
  gg_zeta=  exp(-((gg_x-x1)**2+(gg_y-y1)**2)/(2*sigma**2))
  et_zeta = et_gg(gg_zeta)

 !------------------- 境界値設定 ----------------------

  g_bvalues(:,1) = 0.1*sin(2*pi*g_x)
  g_bvalues(:,2) = 0

  e_bvalues(:,1) = e_g(g_bvalues(:,1))
  e_bvalues(:,2) = e_g(g_bvalues(:,2))

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

 !------------------- 時間積分 ----------------------
  do it=1,nt                                         ! Euler 法による時間積分
     et_zeta = et_zeta + dt * nu * et_lapla_et(et_zeta) 
     call etboundaries(et_zeta,e_bvalues,cond=bcond)

     if(mod(it,ndisp) .eq. 0)then                    ! 出力
        gg_zeta = gg_et(et_zeta)
        call output_gtool4
     endif
  enddo

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

contains

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

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

    call HistoryAddVariable( &                                ! 変数定義
           varname='zeta', dims=(/'x','y','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',transpose(gg_zeta))
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program t2diffusion
