--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rungekutta.f90	Tue Jul 23 02:30:57 2024 +0900
@@ -0,0 +1,100 @@
+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