﻿! 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)
!
! の下解く. 
!
! 時間積分法は 3 位の アダムス-バッシュフォース スキーム.
! 1 ステップ目はホインの 3 次公式で計算.
! 空間微分は中心差分で離散化する.
!
      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_bb(:)    !振幅変数[t = (n-2)Δ t]
       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
       real(8), allocatable :: U_K3(:)   !増分 3
     
       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_bb(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))
       allocate(U_k3(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
!-----時間積分開始-------------------------------------
!------ホインの3次公式---------------------------------------
!1 ステップ目と 2 ステップ目はホインの 3 次公式で計算

       do t_lupe = 1,2
         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)/3.0d0) &
           & -(U_b(n_x)+U_k1(n_x)/3.0d0))/dx*0.5d0)
         U_k2(2:n_x)= &
           & dt*(-iryuu*((U_b(3:n_x+1) + U_k1(3:n_x+1)/3.0d0) &
           & -(U_b(1:n_x-1)+U_k1(1:n_x-1)/3.0d0))/dx*0.5d0)

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

         U_a(1:n_x) = U_b(1:n_x) + (U_k1(1:n_x)+3.0d0*U_k3(1:n_x))*0.25d0
         U_a(n_x+1)=U_a(1)
         select case (t_lupe)
           case (1)
             U_bb = U_k1
             U_b = U_a
           case (2)
             U_c=U_a
         end select

       end do

!-------2段アダムス-バッシュフォース-------------
!2 ステップ目はアダムス-バッシュフォース公式で計算

       do t_lupe= 1,n_t
!既に計算しているU_k1を使いまわす
          U_k3(1) = -dt*iryuu*((U_c(2) - U_c(n_x))/dx*0.5d0) 
          U_k3(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) + (23.0d0*U_k3(1:n_x)&
            & -16.0d0*U_k1(1:n_x)+5.0d0*U_bb(1:n_x))/12.0d0
          U_a(n_x+1) = U_a(1)

          U_b=U_c
          U_c=U_a
          U_bb = U_k1
          U_k1 = U_k3
         

!---------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_bb,U_b,U_c,U_a,U_k1,U_k2,U_k3,x_yoko)
       stop

      end program sabun
