!1次元移流拡散方程式
! \DP{U}{t} + c\DP{U}{x} = \kappa \DP[2]{U}{x}
! c = 1.0, \kappa = 0.01
! 0 \le x \le 1
! 初期条件 U(x,0) = \sin 2\pi x or exp[((x-0.5)/0.25)^2]
! 境界条件 周期境界 U(0,t)=U(1,t)
! 時間:リープフロッグ(アセランフィルタあり,なし), 空間:中心差分
! 時間:オイラー, 空間:中心差分
      
      program advdiv

!-----宣言文---------------------------------------
      implicit none

      real(8), parameter :: pi = 4.0d0*atan(1.0d0) !\pi
      real(8), parameter :: kappa = 0.01d0         !\kappa 拡散係数
      real(8), parameter :: C_ad = 1.0d0           !c 移流速度
      real(8), parameter :: mu = 0.25d0            !フィルタの係数
      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          !終り時間

      real(8), allocatable :: U_a(:)        !U^{n+1}の値
      real(8), allocatable :: U_c(:)        !U^{n}の値
      real(8), allocatable :: U_b(:)        !U^{n-1}の値
      real(8), allocatable :: X_cal(:)      !x の値

      real, allocatable :: U_sing(:)        !dclのための単精度のUの値
      real, allocatable :: X_sing(:)        !dclのための単精度のxの値

      real(8) :: dx                         !\Delta x
      real(8) :: dt                         !\Delta t
      real(8) :: f_i                        !初期条件のための関数

      integer :: x_step                     !空間分割数
      integer :: t_step                     !時間分割数
      integer :: x_lupe                     !空間方向の作業変数
      integer :: t_lupe                     !時間方向の作業変数
      integer :: IWS                        !dclのための作業変数
      integer :: initial                    !初期条件選択をするための変数
      integer :: scheme                     !時間スキーム選択のための変数

      character(10) chara                   !仮入力のための文字列

!-----入力--------------------------
!    キーボード入力が空白の場合デフォルト値を入れるサブルーチンを呼ぶ
!
      write(*,*) "choose initial function. 1=sin(2pi x), &
        & 2=exp[-(x-0.5/0.25)^2] default is 1"

      read(*,'(a1)') chara
      !入力を文字列で取得(空白も取得のため).(a1)は１字複数なら(a)へ
      !初期条件の選択

      call empty(chara,1,initial)
      !仮文字列が空白の場合デフォルト値1を入れる

      select case(initial)
        case(1)
          write(*,*) "You chose sin(2pi x)"
        case(2)
          write(*,*) "You chose exp[-(x-0.5/0.25)^2]"
      end select
      !どちらの初期条件を選んだのかを表示

      write(*,*) "choose time scheme. &
        & 1=euler,2=euler+leapfrog,3=euler+leapfrog+asselin &
        & default is 3"

      read(*,'(a1)') chara
      !時間スキームの選択

      call empty(chara, 3, scheme)
      !仮文字列が空白の場合デフォルト値3を入れる

      select case(scheme)
        case(1)
          write(*,*) "You chose euler"
        case(2)
          write(*,*) "You chose euler + leapfrog"
        case(3)
          write(*,*) "You chose euler + leapfrog + asselinfilter"
      end select
      !どのスキームを選んだか表示

      write(*,*) "Enter in time steps. & 
        & if you want to enter more 10000000000,&
        &  enter in 'a'. and rekey. default is 800."! 500. "

      read(*,'(a)') chara
      !時間分割数の仮入力

      call empty(chara, 800, t_step)
      !仮文字列が空白の場合デフォルト値200を入れる

      write(*,*) "You chose ", t_step

      write(*,*) "Enter in X steps. & 
        & if you want to enter more 10000000000,&
        &  enter in 'a'. and rekey. default is 10."!20 "

      read(*,'(a)') chara
!      !空間分割数の仮文字列入力

      call empty(chara, 10, x_step)
      !仮文字列が空白の場合デフォルト値10を入れる

      write(*,*) "You chose ", x_step
      open(10, file='test.dat', status='replace')
!-----配列の動的振り付け--------------------
!
      allocate(U_a(0:x_step))
      allocate(U_c(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))/dble(x_step)
      dt = (abs(tmax) - abs(tmin))/dble(t_step)

      select case(scheme)
        case(1)
          write(*,*) "Curan=", C_ad*(dt/(2.0d0*dx))
          write(*,*) "K=", kappa*(dt/(dx**2))*((2.0d0*pi)**2)
        case(2:)
          write(*,*) "Curan=", C_ad*(dt/dx)
          write(*,*) "K=", kappa*((2.0d0*dt)/(dx**2))*((2.0d0*pi)**2)
      end select
!-----dclの自動書き出し-----------------------
!
      call SWpSET('LDUMP',.TRUE.)
      call swpset('LWAIT',.false.)
      call swpset('lwait1',.false.)

!----dclの開始--------------------------
!
      call sgpwsn
      read(*,'(a)') chara
      call empty(chara,1,IWS)
      call gropn(IWS)

!-----初期値の計算---------------------
!  初期条件を関数呼び出しにしている. 
!  initial=1のとき\sin(2\pi x)
!  initial=2のとき\exp[(\frac{x -0.5}{0.25})^2]
!
      do x_lupe = 0, x_step

         X_cal(x_lupe) = xmin + dx*dble(x_lupe)

         U_b(x_lupe) = f_i(X_cal(x_lupe), initial)
         !初期条件を関数で代入

      end do

      U_c=U_b

      X_sing=real(X_cal)
      U_sing=real(U_b)

!----初期値書き出し---------------------
!
      call grfrm
      call grswnd(0.0,1.0,-1.0,1.0)

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

!----境界条件-----------------------------
!
      U_b(x_step) = U_b(0)

!----時間積分---------------------------
!    時間1ステップごとにxを計算してdclで表示
!    schemeで選択 1ならオイラー,
!    2ならオイラー+リープフロッグ,3ならオイラー+リープフロッグ+アセラン
!    空間差分は中心差分
!    リープフロッグの場合は一ステップ目はオイラースキームで計算している
!
      do t_lupe = 1, t_step

        if(t_lupe == 1 .or. scheme == 1) then

          U_b = U_c
          U_a(0) = &
            & U_b(0) - C_ad * (dt/(2.0d0*dx)) * &
            & ( U_b(1) - U_b(x_step-1)) + &
            & kappa * (dt/(dx**2)) * &
            & ( U_b(1) - 2.0d0* U_b(0) + U_b(x_step-1))

          !境界条件 -1=I-1 の格子点を考慮した格子点i=0の計算
          U_a(x_step) = &
            & U_b(x_step) - C_ad * (dt/(2.0d0*dx)) * &
            & ( U_b(1) - U_b(x_step-1)) + &
            & kappa * (dt/(dx**2)) * &
            & ( U_b(1) - 2.0d0* U_b(x_step) + U_b(x_step-1))
                   
          !境界条件 I+1=1 の格子点を考慮した格子点i=Iの計算
          do x_lupe = 1, x_step-1

             U_a(x_lupe) = &
               & U_b(x_lupe) - C_ad * (dt/(2.0d0*dx)) * &
               & ( U_b(x_lupe+1) - U_b(x_lupe-1)) + &
               & kappa * (dt/(dx**2)) * &
               & ( U_b(x_lupe+1) - 2.0d0* U_b(x_lupe) + U_b(x_lupe-1))

             !i=1,...I-1 までのオイラースキーム
          end do

          U_a(x_step)=U_a(0) 

          !境界条件 0=I
          U_c = U_a
          U_sing=real(U_a)

          do x_lupe = 0, x_step
            write(10,*) U_sing(x_lupe)
          end do
          !---euler計算dcl------------------------------------
          !
          call grfrm
          call grswnd(0.0,1.0,-1.0,1.0)

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

        else
          ! スキームをリープフロッグを使うことを選んだ場合2ステップ目から

          U_a(0) = &
            & U_b(0) - C_ad * (dt/dx) * &
            & ( U_c(1) - U_c(x_step-1)) + &
            & kappa * ((2.0d0*dt)/(dx**2)) * & 
            & ( U_b(1) - 2.0d0*U_b(0) + U_b(x_step-1)) 
               
          !境界条件 -1=I-1 の格子点を考慮した格子点i=0の計算
          U_a(x_step) = &
            & U_b(x_step) - C_ad * (dt/dx) * &
            & ( U_c(1) - U_c(x_step-1)) + &
            & kappa * ((2.0d0*dt)/(dx**2)) * &
            & ( U_b(1) - 2.0d0*U_b(x_step) + U_b(x_step-1))

          !境界条件 I+1=1 の格子点を考慮した格子点i=Iの計算
          do x_lupe = 0, x_step-2

            U_a(x_lupe + 1) = &
              & U_b(x_lupe + 1) - C_ad * (dt/dx) * &
              & ( U_c(x_lupe + 2) - U_c(x_lupe)) + &
              & kappa * ((2.0d0*dt)/(dx**2)) * &
              & ( U_b(x_lupe + 2) - 2.0d0*U_b(x_lupe + 1) + U_b(x_lupe))
              !i=1,..I-1のリープフロッグスキームの計算
          end do
          U_a(x_step)=U_a(0) 

          !境界条件 0=I 
          if(scheme == 3)then
          !アセランフィルタを使うことを選んだ場合はアセランフィルタを入れる

            do x_lupe = 0, x_step 

              U_c(x_lupe) = &
                & U_c(x_lupe) + 0.5d0 * mu * &
                & ( U_a(x_lupe) - 2.0d0*U_c(x_lupe) + U_b(x_lupe))

              !アセランフィルタ
            end do
          end if
          U_b = U_c
          U_c = U_a
          
          U_sing=real(U_a)

          !---euler+leapfrog(scheme=3のときはasselinも)計算dcl------
          !
          call grfrm
          call grswnd(0.0,1.0,-1.0,1.0)

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

        end if    

      end do

      call grcls

      close(10)
      deallocate(U_a,U_c,U_b,X_cal,U_sing,X_sing)
      !配列の解放
      stop
      end program advdiv

!-----初期条件関数----------------------
!     初期条件を関数として与える.
!     initial=1のときは sin(2\pi x)
!     initial=2のときは e^{(\frac{x-0.5}{0.25})^2}
!
      function f_i(x, n)

      !-----------初期条件の宣言文-------------
      !
      implicit none

      real(8), parameter :: pi = 4.0d0*atan(1.0d0)

      real(8) f_i
      real(8) x

      integer n
      !-----------初期条件------------------
      !
      select case(n)
        case (1)

          f_i = sin(2*pi*x)

          !n=1のときはsin
        case(2)

          f_i = exp(-((x - 0.5d0)/0.25d0)**2)

          !n=2のときはexp
      end select
      return
      end function f_i
!-----入力判断サブルーチン-------------------------
!入力診断のサブルーチン.
!キーボード入力が空白のときはデフォルト値を入れる
!
      subroutine empty(chara,default_step,step)

      !------入力診断の宣言文----------------------
      implicit none

      character(10) chara

      integer default_step
      integer step

      !------入力診断------------------------------
      if(len_trim(chara) == 0)then

        step = default_step
        !空白の場合デフォルト値を代入
      else

        if(chara == "a")then

          read(*,*) step
          !aの場合さらにもう一度代入
        else

          read(chara,*) step
          !仮文字列から実際の変数への代入

        end if     
             
      end if
      return
      end subroutine empty
