!1次元拡散方程式
! \DP{U}{t} = \kappa \DP[2]{U}{x}
! 0 \le x \le 1
! 初期条件 U(x,0) = \sin \pi x
! 境界条件 U(0,t)=U(1,t) = 0
! \kappa = 0.01
! 時間方向はオイラー法, 空間方向は中心差分で解く
! U^{n+1}_{i} = U^{n}_{i} + \kappa \frac{dt}{dx}(U^{n}_{i+1} 
!               - 2U^{n}_{i} + U^{n}_{i-1})

      program div

      implicit none

!-----宣言文---------------------
      real(8), allocatable :: U_a(:)     ! U^{n+1}
      real(8), allocatable :: U_b(:)     ! U^{n}
      real(8), allocatable :: X_cal(:)   ! x

      real, allocatable :: U_sing(:)     ! 数値解
      real, allocatable :: X_sing(:)     ! 横軸
      real(8) :: dx                      ! \Delta x
      real(8) :: dt                      ! \Delta t 
      integer :: x_step                  ! x方向分割数
      integer :: t_step                  ! 時間分割数
      integer :: x_lupe                  ! x作業変数
      integer :: t_lupe                  ! t作業変数
      integer :: IWS                     ! dcl のための作業変数
      character(10) chara                ! 空白enterをするための文字列

      real(8), parameter :: pi = 4.0d0*atan(1.0d0) ! \pi
      real(8), parameter :: kappa = 0.01           ! \kappa
      real(8), parameter :: xmin = 0.0d0           ! xの始端
      real(8), parameter :: xmax = 1.0d0           ! xの終端
      real(8), parameter :: tmin = 0.0d0           ! 開始時間
      real(8), parameter :: tmax = 10.0d0          ! 終了時間
      
!----------- 分割数入力-------------------------
      write(*,*) "Enter in time steps. default is 10."

      read(*,'(a)') chara                      !キーボードからの入力
      if(len_trim(chara) == 0)then

        chara ="10"  !キーボードで空白えんたーの場合デフォルト値10を入れる
        read(chara,*) t_step         !仮文字列からの代入

      else

        read(chara,*) t_step !仮文字列からの代入

      endif
      write(*,*) "Enter in X steps. default is 10 " 
      read(*,'(a)') chara                      !キーボードからの入力
      if(len_trim(chara) == 0)then

        chara ="10"  !キーボードで空白えんたーの場合デフォルト値10を入れる
        read(chara,*) x_step         !仮文字列からの代入

      else

        read(chara,*) x_step  !仮文字列からの代入

      endif

!-----配列の動的振り付け-------------------------
      allocate(U_a(0:x_step))
      allocate(U_b(0:x_step))
      allocate(X_cal(0:x_step))
      allocate(U_sing(0:x_step))
      allocate(X_sing(0:x_step))

!------dt,dxの計算----------------------------
      dx = (abs(xmax)-abs(xmin))/real(x_step)
      dt = (abs(tmax)-abs(tmin))/real(t_step)

      ! 安定条件. この値が1以下だと安定
      write(*,*) "K =", kappa*(dt/(dx**2))*((pi)**2)

!-----dclの自動書き出し-----------------------
      call SWpSET('LDUMP',.TRUE.)
      call swpset('LWAIT',.false.)
      call swpset('lwait1',.false.)

!-----dclの開始-----------------------------
      call sgpwsn
      read(*,*) IWS
      call gropn(IWS)

!-----初期値の計算--------------------------
      do x_lupe = 0,x_step
        ! x の計算
        X_cal(x_lupe) = xmin + dx*x_lupe
        ! U(0,x) = sin(\pi x)
        U_b(x_lupe) = sin(pi*X_cal(x_lupe))

      enddo
      !単精度化
      X_sing=real(X_cal)
      U_sing=real(U_b)

!-------初期値書き出し----------------------------------
      call grfrm
      call grswnd(0.0,1.0,0.0,1.0)
!        call usspnt(x_step+1, X_sing(0:x_step),U_sing(0:x_step))

      call uspfit
      call grstrf
      call usdaxs
      call uulin(x_step+1, X_sing(0:x_step),U_sing(0:x_step))        
!-------境界条件-----------
      U_b(0) = 0.0d0
      U_b(x_step) = 0.0d0
      U_a(0) = 0.0d0
      U_a(x_step) = 0.0d0

!---------------オイラー法---------------------------------
!     時間一ステップごとにxを計算してdclで表示している.
!     時間微分はオイラースキーム
!     空間差分は中心差分を行っている  
      do t_lupe = 1,t_step

        do x_lupe = 0, x_step-2

          ! t=t_lupe*dt でのオイラースキーム 
          U_a(x_lupe+1) = U_b(x_lupe+1) & 
            & + kappa * (dt/(dx**2))*(U_b(x_lupe+2) &
            & - 2.0d0*U_b(x_lupe+1) &
            & + U_b(x_lupe))

        enddo

        !単精度化
        U_sing = real(U_a)

        !-------dcl------------------
        call grfrm
        call grswnd(0.0,1.0,0.0,1.0)
!        call usspnt(x_step+1, X_sing(0:x_step),U_sing(0:x_step))

        call uspfit
        call grstrf
        call usdaxs
        call uulin(x_step+1, X_sing(0:x_step),U_sing(0:x_step))

        ! U^{n} = U^{n-1}
        U_b(:) = U_a(:)        

      enddo

      call grcls

      ! 配列の開放
      deallocate(U_a, U_b, X_cal, X_sing, U_sing)
      stop

      end program div

