
program stream

  use gtool_history                                   ! モジュール指定

  integer                :: i, j, k
  integer, parameter     :: nx=100                   ! グリッド数
  integer, parameter     :: ny=100                   ! グリッド数
  integer, parameter     :: nt=100                 
  real(8), parameter     :: dx=12.0/nx           ! グリッド間隔
  real(8), parameter     :: dy=12.0/ny           ! グリッド間隔
  real(8), parameter     :: dt=1.0               ! 時間間隔
  real(8), dimension(nx) :: x
  real(8), dimension(ny) :: y
  real(8), dimension(nt) :: t
  real(8)                :: psi(nx, ny, nt)
  !real(8)                :: g(nx, ny, nt)
  !real(8)                :: h(nx, ny, nt)
  real(8), parameter     :: U = 1.0
  real(8), parameter     :: wnum = 0.9
  !complex(8), parameter  :: CI = (0.0,0.187)
  complex(8), parameter  :: CI = (0.4349,0.0)
  complex(8), parameter  :: Im = (0.0,1.0)
  complex(8)             :: A
  real(8), parameter     :: B = 1.0
  !real(8), parameter     :: C = -1.0
  real(8), parameter     :: C = 1.0
  
  do i = 1, nx
     x(i) = -6 + dx*(i-1)
  end do

  do j = 1, ny
     y(j) = -6 + dy*(j-1)
  end do

  do k = 1, nt
     t(k) = dt*(k-1)
  end do

    
  call HistoryCreate( &                             ! ヒストリー作成
    & file='stream.nc', title='Diffusion equation', &
    & source='Sample program of gtool_history/gtool5',   &
    & institution='GFD_Dennou Club davis project',       &
    & dims=(/'x','y','t'/), dimsizes=(/nx,ny,nt/),               &
    & longnames=(/'X-coordinate','Y-coordinate','time        '/),       &
    & units=(/'m','m','s'/) )

  call HistoryPut('x',x)                           ! 次元変数出力
  call HistoryPut('y',y)                           ! 次元変数出力
  call HistoryPut('t',t)                           ! 次元変数出力

  call HistoryAddVariable( &                        ! 変数定義
    & varname='psi', dims=(/'x','y','t'/), &
    & longname='streamfunction', units='/s', xtype='double')


  !A=(2*wnum*(1-CI)-1)*exp(2*wnum)
   A=( (wnum*(1-CI)-1)/(-wnum*(1-CI)) )*exp(2*wnum)
  
  ! caluculation
     do j = 1, ny
        
        if ( y(j) < C ) then
           do k = 1, nt
              do i = 1, nx
                 psi(i, j, k) = aimag(exp(wnum*(y(j)+1))*exp(Im*wnum*(x(i)-CI*t(k)))) 
              end do
           end do
           
           !else if ( y(j) > B ) then
        else
           do k = 1, nt
              do i = 1, nx                 
                 psi(i, j, k) = aimag( A*exp(-wnum*(y(j)-1))*exp(Im*wnum*(x(i)-CI*t(k))) )
              end do
           end do
           
        ! else
        !    do k = 1, nt
        !       do i = 1, nx
        !          psi(i, j, k) = aimag((exp(wnum*(y(j)-1))+ A*exp(-wnum*(y(j)+1)))*exp(Im*wnum*(x(i)-CI*t(k))))
        !       end do
        !    end do
           
        end if
       
  end do
  
  call HistoryPut('psi',psi)
  
  call HistoryClose
  stop
end program stream
