..................の間がプログラム本体です
program 15
..................
C サブルーチンの例
C  行列とベクトルを与えるとその積を返す
C  (program 7 をサブルーチンを用いて書き換えたもの)
      program main
      real*8 a(2,2),x(2),y(2)

C  行列に値をいれる(値自身は意味ないです)
      a(1,1)=0.5
      a(1,2)=1.0
      a(2,1)=0.8
      a(2,2)=2.0

C  ベクトルに値をいれる
      x(1)=1.0
      x(2)=2.0

C 積を求めるサブルーチンを呼び出す
      call product(a,x,y)

C 出力 
      write(*,*)"matrix a:"
      do i=1,2
        write(*,*)(a(i,j),j=1,2)
      enddo
      write(*,*)"vector x:"
      write(*,*)(x(i),i=1,2)
      write(*,*)"vector y=a*x:"
      write(*,*)(y(i),i=1,2)
      end

C  行列とベクトルの積を求めるサブルーチン
      subroutine product(a,x,y)
      real*8 a(2,2),x(2),y(2)
    
C  ベクトルyを初期化(値を0に)
      do i=1,2
        y(i)=0
      enddo

C a*xによりyを求める
      do i=1,2
       do j=1,2
         y(i)=y(i)+a(i,j)*x(j)
       enddo
      enddo        

      end
..................
program 16a
..................
C サブルーチンの応用例
C 最低次のsimplectic integrator (対称分解)によって
C 一次元調和振動子の運動の初期値問題を解く。
      program main
      real*8 delta_t,q,p

C 時間刻み、ステップ数、初期の座標と運動量を入力
      read(*,*)delta_t,istep,q,p

C 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ
      do i=1,istep
         call simplectic(q,p,delta_t)

C 各ステップごとに、時刻、座標、運動量、エネルギーを出力 
         write(*,*)i*delta_t,q,p,(p*p+q*q)/2
      enddo
      end

C 1ステップ分の計算をするサブルーチン
      subroutine simplectic(q, p, delta_t)
      real*8 q,p,delta_t

      q = q + p * delta_t/2
      p = p - q * delta_t
      q = q + p * delta_t/2
      end
..................
program 16b
..................
C 16aと同じ問題を、配列をサブルーチンに渡す方式で書き直すとこうなる
      program main
      real*8 delta_t,qp(2)

C 時間刻み、ステップ数、初期の座標と運動量を入力
      read(*,*)delta_t,istep,qp(1),qp(2)

C 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ
      do i=1,istep
         call simplectic(qp,delta_t)

C 各ステップごとに、時刻、座標、運動量、エネルギーを出力 
         write(*,*)i*delta_t,qp(1),qp(2),(qp(1)*qp(1)+qp(2)*qp(2))/2
      enddo
      end

C 1ステップ分の計算をするサブルーチン
      subroutine simplectic(qp,delta_t)
      real*8 qp(2),delta_t

      qp(1) = qp(1) + qp(2) * delta_t/2
      qp(2) = qp(2) - qp(1) * delta_t
      qp(1) = qp(1) + qp(2) * delta_t/2
      end
..................
program 17
..................
C 16より少しだけ汎用のプログラム
      program main
      real*8 delta_t,q,p

C 時間刻み、ステップ数、初期の座標と運動量を入力
      read(*,*)delta_t,istep,q,p

C 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ
      do i=1,istep
         call simplectic(q,p,delta_t)

C 各ステップごとに、時刻、座標、運動量、エネルギーを出力 
         write(*,*)i*delta_t,q,p,(p*p+q*q)/2
      enddo
      end

C 1ステップ分の計算をするサブルーチン
      subroutine simplectic(q, p, delta_t)
      real*8 q,p,delta_t

      q = q + hq(q,p) * delta_t/2
      p = p + hp(q,p) * delta_t
      q = q + hq(q,p) * delta_t/2
      end

C  simplectic integratorの座標部分
C (調和振動子の場合はほとんど自明)
      function hq(q,p)
      real*8 q,p
      hq = p
      end

C  simplectic integratorの運動量部分
C (調和振動子の場合はほとんど自明)
      function hp(q,p)
      real*8 q,p
      hp = -q
      end
..................
program 18
..................
C 16と同じ問題を四次のRunge-Kuttaで解くプログラム
      program main
      real*8 delta_t,q,p

C 時間刻み、ステップ数、初期の座標と運動量を入力
      read(*,*)delta_t,istep,q,p

C 1ステップ分を計算する関数をステップ数だけ繰り返して呼ぶ
      do i=1,istep
         call rungekutta(q,p,delta_t)

C 各ステップごとに、時刻、座標、運動量、エネルギーを出力 
         write(*,*)i*delta_t,q,p,(p*p+q*q)/2
      enddo
      end

C 1ステップ分の計算をするサブルーチン
      subroutine rungekutta(x, v, delta_t)
      real*8 x,v,delta_t
      real*8 k1x,k1v,k2x,k2v,k3x,k3v,k4x,k4v
      call onestep(x, v, k1x, k1v,delta_t)
      call onestep(x+k1x/2, v+k1v/2, k2x, k2v, delta_t)
      call onestep(x+k2x/2, v+k2v/2, k3x, k3v, delta_t)
      call onestep(x+k3x, v+k3v, k4x, k4v, delta_t)
      x = x + k1x/6 + k2x/3 + k3x/3 + k4x/6
      v = v + k1v/6 + k2v/3 + k3v/3 + k4v/6
      end

C rungekuttaの一段分(調和振動子の場合はほとんど自明)
      subroutine onestep(x,v,xx,vv,delta_t)
      real*8 x,v,xx,vv,delta_t
      xx = v*delta_t
      vv = -x*delta_t
      end