module rungekutta
  use,intrinsic :: iso_fortran_env
  private
  !---
  integer,parameter,public :: dpkind=real64 
  real(kind=dpkind),public,parameter :: PI=3.1415926535897932384626433_dpkind
  !---
  type,public :: rk
    real(kind=dpkind) :: time=0e0 !<計算開始からの総経過時間
    real(kind=dpkind) :: dt=1e-4   !<時間刻み
    procedure(rhside),pass,pointer :: fun => null() 
                !<微分方程式の右辺を計算する関数 
    procedure(rksol),pass,pointer :: solver => null()      
                !<微分方程式の解法
    contains
      final :: final_rk !<構造体のファイナライズ処理
      procedure :: setup !<数値積分の初期値、解法をセットする手続き
  end type
  !-------------
  interface
    !---
    ! Right hand side function interface
    function rhside(self,time,wk)
      import :: dpkind,rk
      class(rk),intent(inout) :: self !parameter
      real(kind=dpkind) :: time
      real(kind=dpkind),dimension(:) :: wk
      real(kind=dpkind),dimension(size(wk)) :: rhside
    end function
    !---
    ! Solver interface
    subroutine rksol(self,tnext,wk)
      import :: dpkind,rk
      class(rk),intent(inout) :: self
      real(kind=dpkind),intent(in) :: tnext !output time
      real(kind=dpkind),dimension(:) :: wk
    end subroutine
  end interface
  !============================--
  contains
    !---
    pure elemental subroutine final_rk(self)
      type(rk),intent(inout) :: self
      if(associated(self%fun))nullify(self%fun)
      if(associated(self%solver))nullify(self%solver)
    end subroutine
    !---
    subroutine setup(self,fun,solver,dt)
      class(rk),intent(inout) :: self
      procedure(rhside) :: fun  !微分方程式の右辺を計算する関数
      character(len=*),intent(in) :: solver !ソルバーの選択
      real(kind=dpkind),intent(in),optional :: dt !時間刻み幅 
      !時間刻みのセット
      if(present(dt))then
        self%dt=dt
      endif
      self%fun=>fun
      !ソルバーの選択
      select case(solver)
      case("rk4")
        !ルンゲクッタ法
        self%solver=>rk4
      case("rkf45")
        !ルンゲクッタ・フェルベルグ法(刻み幅自動調整型)
        !self%solver=>rkf45
      case default
        write(error_unit,*) "solver must one of the 'rk4, rkf45'"
      endselect
    end subroutine
    !---
    subroutine rk4(self,tn,wk)
      class(rk),intent(inout) :: self
      real(kind=dpkind),intent(in) :: tn !次に出力したい時間
      real(kind=dpkind),dimension(:) :: wk
      real(kind=dpkind),dimension(size(wk)) :: k1,k2,k3,k4
      real(kind=dpkind) :: dt_tmp
      real(kind=dpkind) :: dthalf
      associate(time=>self%time,dt=>self%dt)
        dt_tmp=dt
        do
          if(time+dt>tn)dt_tmp=tn-time
          dthalf=dt_tmp/2e0
          !!$OMP workshare
          k1=self%fun(time,wk)
          k2=self%fun(time+dthalf,wk+k1*dthalf)
          k3=self%fun(time+dthalf,wk+k2*dthalf)
          k4=self%fun(time+dt_tmp,wk+k3*dt_tmp)
          !状態変数を更新
          wk=wk+dt_tmp*(k1+2e0*k2+2e0*k3+k4)/6e0
          !!$OMP end workshare
          !1ステップ次の解が計算できたので、時間を進める
          time=time+dt_tmp 
          !時間の計算範囲がtnを超えたら計算終了
          if(time>=tn)then
            exit
          endif
        enddo
      end associate
    end subroutine
end module