! 1次元線形移流方程式
!
!  dU(x,t)/dt + c dU(x,t)/dx = 0 (c =　1.0)
!
! を, 境界条件は周期条件, 初期条件
!
!  U(x,0) = 1.0 (0 < x < 1.0)
!         = 0.0 (1.0 < x < 5.0)
!
! の下解く. 
!
! 時間積分法は 2 位の アダムス-バッシュフォース スキーム.
! 1 ステップ目はホイン法で計算.
! 空間微分は中心差分で離散化する.
!
      program sabun

!------宣言部------------------------

       implicit none

       real(8), allocatable :: x_yoko(:) !座標値
       real(8) :: xmax = 1.0d0           !描画領域の上限
       real(8) :: xmin = 0.0d0           !描画領域の下限
       real(8) dx                        !空間格子間隔
       real(8) dt                        !時間格子間隔
       real(8) :: time = 1.0d0           !描画の時間
       real(8) pi                        !円周率
       real(8) :: iryuu = 1.0d0             !移流速度
       real(8), allocatable :: U_b(:)    !振幅変数[t = (n-1)Δ t]
       real(8), allocatable :: U_c(:)    !振幅変数[t = (n)Δ t]
       real(8), allocatable :: U_a(:)    !振幅変数[t = (n+1)Δ t]
       real(8), allocatable :: U_k1(:)   !増分1
       real(8), allocatable :: U_k2(:)   !増分2
     
       integer t_lupe                    ! DO ループ変数
       integer j_lupe                    ! DO ループ変数
       integer n_t                       ! 時間格子点数
       integer n_x                       ! 空間格子点数
       integer IWS                       ! dcl 用内部変数
       integer pf                        ! 描画時間の指定用変数
    
!------パラメータ部--------------------------------

       pi = 4.0d0*atan(1.0d0)

       write(*,*) "Enter in cut number n_x"
       read(*,*) n_x

       write(*,*) "Enter in time step n_t"
       read(*,*) n_t

       write(*,*) "Enter in stop time"
       read(*,*) pf

       dx =((xmax-xmin))/dble(n_x+1)
       dt = time/dble(n_t)

       write(*,*) "curant number is ", iryuu*dt/dx

       allocate(x_yoko(1:n_x+1))
       allocate(U_b(1:n_x+1))
       allocate(U_c(1:n_x+1))
       allocate(U_a(1:n_x+1))
       allocate(U_k1(1:n_x+1))
       allocate(U_k2(1:n_x+1))

!      call SWpSET('LDUMP',.TRUE.)
!      call SWpSET('LWAIT',.FALSE.)
!      call SWpSET('LWAIT1',.FALSE.)

!------座標値の設定-----------------------------

       x_yoko = 0.0d0
  
       do j_lupe = 2, n_x+1
          x_yoko(j_lupe) = x_yoko(j_lupe-1) + dx
       end do
 
!------初期値の設定----------------------------
       U_b = 0.0d0       

       U_b = max(0.0d0, cos(2.0d0*pi*x_yoko))

       U_b(n_x + 1) = U_b(1)

!------DCL 初期化-------------------------------
         call sgpwsn
         read(*,*) IWS
         call gropn(IWS) 
       if(pf.eq.0)then
         call grfrm
         call grswnd(real(xmin),real(xmax),-1.2,1.2)
         call uspfit
         call grstrf
         call usdaxs
         call uulin(&
           & n_x+1,real(x_yoko(1:n_x+1)), real(U_b(1:n_x+1)))
       end if
!-----時間積分開始-------------------------------------
!------ホイン法---------------------------------------
!1 ステップ目はホイン法で計算

       U_k1(1)= &
         & dt*(-iryuu*((U_b(2)-U_b(n_x))/dx*0.5d0))
       U_k1(2:n_x)=&
         & dt*(-iryuu*((U_b(3:n_x+1)-U_b(1:n_x-1))/dx*0.5d0))

       U_k2(1)= &
         & dt*(-iryuu*((U_b(2) + U_k1(2)) &
         & -(U_b(n_x)+U_k1(n_x)))/dx*0.5d0)
       U_k2(2:n_x)= &
         & dt*(-iryuu*((U_b(3:n_x+1) + U_k1(3:n_x+1)) &
         & -(U_b(1:n_x-1)+U_k1(1:n_x-1)))/dx*0.5d0)

       U_a(1:n_x) = U_b(1:n_x) + (U_k1(1:n_x)+U_k2(1:n_x))*0.5d0
       U_a(n_x+1)=U_a(1)

       U_c=U_a
!-------2位アダムス-バッシュフォース-------------
!2 ステップ目はアダムス-バッシュフォース公式で計算
       do t_lupe= 1,n_t
!既に計算しているU_k1を使いまわす
          U_k2(1) = -dt*iryuu*((U_c(2) - U_c(n_x))/dx*0.5d0) 
          U_k2(2:n_x) = &
            & -dt*iryuu*((U_c(3:n_x+1)-U_c(1:n_x-1))/dx*0.5d0)

          U_a(1:n_x) = &
            & U_c(1:n_x) + (1.5d0*U_k2(1:n_x) -0.5d0*U_k1(1:n_x))
          U_a(n_x+1) = U_a(1)
          U_b=U_c
          U_c=U_a
          U_k1 = U_k2
!---------DCL による描画-------------------------------------
          if(pf.eq.t_lupe)then
            call grfrm
            call grswnd(real(xmin),real(xmax),-1.2,1.2)
            call uspfit
            call grstrf
            call usdaxs
            call uulin(n_x+1,real(x_yoko(1:n_x+1)),&
              & real(U_a(1:n_x+1)))
          end if
       end do

       call grcls
       deallocate(U_b,U_c,U_a,U_k1,U_k2,x_yoko)
       stop

      end program sabun