!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