!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  CTPACK/ISPACK-F90 テストプログラム
!      1 次元拡散方程式
!
!履歴  2002/01/24  竹広真一
!
program ctdiff

  use gtool_history
  use ctpack
  use lumatrix

  implicit none
  integer, parameter :: im=8, km=5, nm=4
  real(8), parameter :: xmin=-1.0D0, xmax=1.0D0

  real(8), dimension(0:im) :: g_temp
  real(8), dimension(0:im) :: g_temp_right
  real(8), dimension(0:km) :: t_temp
  real(8), parameter       :: kappa = 1.0

  logical, parameter       :: tfix_top    = .false.        ! 上端境界条件
  logical, parameter       :: tfix_bottom = .false.        ! 下端境界条件
  real(8), parameter       :: tvalue_top    = 0.0
  real(8), parameter       :: tvalue_bottom = 0.0
  real(8), parameter       :: dt=1e-2
  integer, parameter       :: nt=100
  integer, parameter       :: ndisp=10

  real(8), parameter       :: pi=3.1415926535897932385D0
  integer  :: it

  call ctinitial(im,km,xmin,xmax)
  g_temp = sin(pi*g_x/2)
  g_temp_right = g_temp

  t_temp = t_g(g_temp)

  call output_gtool4_init
  call output_gtool4

  do it=1,nt
     t_temp = t_temp + kappa * t_dx_t(t_dx_t(t_temp)) * dt

     if ( (tfix_bottom==.false.).AND.(tfix_top==.false.) )then
        call ctboundaries_NN(t_temp,(/tvalue_top,tvalue_bottom/))
     elseif ( (tfix_bottom==.true.).AND.(tfix_top==.false.) )then
        call ctboundaries_ND(t_temp,(/tvalue_top,tvalue_bottom/))
     elseif ( (tfix_bottom==.false.).AND.(tfix_top==.true.) )then
        call ctboundaries_DN(t_temp,(/tvalue_top,tvalue_bottom/))
     else
        call ctboundaries_DD(t_temp,(/tvalue_top,tvalue_bottom/))
     endif

     if ( mod(it,ndisp) == 0 ) then
        g_temp_right = sin(pi*g_x/2)*exp(-kappa*pi**2/4*it*dt)
        g_temp=g_t(t_temp)
        call output_gtool4
     endif
  enddo

  call output_gtool4_close
  stop

contains

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

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

    call HistoryAddVariable( &                                ! 変数定義
           varname='temp', dims=(/'x','t'/), & 
           longname='temperature', units='1', xtype='double')

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

  subroutine output_gtool4
    write(6,*) 'it = ',it
    call HistoryPut('t',real(dt*it))
    call HistoryPut('temp',g_temp)
    call HistoryPut('temp_right',g_temp_right)
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program ctdiff
