能書き

基本さえ押さえておけば、あとは知恵と工夫で何でもできる！
・・・と言っても、全てを自分で一から考えるのは不毛です。
定番のやり方はいくつか覚えておきましょう。

能書き２

C/C++やPythonのサンプルコードは検索すればいくらでも出てくるのですが、Fortranで書かれたものはなかなか見つからないですね。

swap

２つの変数『a』と『b』の値を入れ替える

 temp = a
 a = b
 b = temp

1. 『a』の値を一時変数『temp』に待避させる
2. 『b』の値を『a』に代入 3. 『temp』の値(元の『a』の値)を『b』に代入 という手順。

swapサブルーチン

このコードは既出ですが、大事な話なので２回目です。
値を入れ替えるだけの単純作業ですが、コードにすると３行になる上、一時変数を定義する必要があります。
多用するのであれば、サブルーチンにしてしまうと良いでしょう。

 subroutine swap(a,b)
  real(kind(0d0)), intent(inout) :: a, b
  real(kind(0d0)) :: temp
     temp = a
     a = b
     b = temp
 end subroutine

入れ替える変数の型が違ったり、ベクトルだったりする場合は適宜定義を変更して下さい。

swapモジュール

複数のプログラムで使うなら、モジュールにしてしまうと楽です。
swaps.f90
ここでは、総称名とelemental属性を使ってどんな型の変数でも『swap』で入れ替えられるようにしています。

バブルソート

一次元配列(要素数nのベクトル)の中身を小さい順に並べ替える

 do i=1,n-1
    do j=2,n-i+1
       if( vec(j-1) > vec(j) ) then
          temp = vec(j-1)
          vec(j-1) = vec(j)
          vec(j) = temp
       end if
    end do
 end do

最も簡単で基本的なソートアルゴリズム。隣同士を比較して入れ替える、を繰り返す。
Wikipediaのアニメーションを見るとイメージが掴みやすい。
２重ループなので、計算量はO(n²)。
nが大きい場合は他のアルゴリズムも検討すると良い。

バブルソート：サブルーチン

こういう定型コードは、変数名の干渉を避けるためにもサブルーチンにしてしまいましょう。
この例では上のswapモジュールを使っています。

subroutine bubbleSort(vec)
 real(kind(0d0)), intent(inout) :: vec(:)
 integer :: i,j,n
   n = size(vec)
   do i=1,n-1
      do j=2,n-i+1
         if( vec(j-1) > vec(j) ) then
            call swap(vec(j-1),vec(j))
         end if
      end do
   end do
end subroutine

クイックソート：サブルーチン

計算量が平均O(n log n)となる、一般に速いとされているソートアルゴリズム。
このレベルになると、自分で考えて自作する意味はあまり無い。

recursive subroutine quickSort_d(a,left,right)
 real(kind(0d0)), intent(inout) :: a(:)
 integer, intent(in), optional :: left, right
 real(kind(0d0)) :: pivot
 integer :: i, j, imin, jmax
   if(present(left)) then;  imin = left
   else;                    imin = lbound(a,1); end if
   if(present(right)) then; jmax = right
   else;                    jmax = ubound(a,1); end if
   if (imin < jmax) then
      i = imin; j = jmax
      pivot = this__median3_d(a(i), a(i+(j-i)/2), a(j))
      do
         do while( a(i) < pivot ) 
            i = i + 1
         end do
         do while( a(j) > pivot )
            j = j - 1
         end do
         if ( i >= j )  EXIT
         call swap(a(i),a(j))
         i = i + 1
         j = j - 1
      end do
      call quickSort_d(a(:),imin,i-1) !-- sort smaller part
      call quickSort_d(a(:),j+1,jmax) !-- sort larger part
   end if
end subroutine
function this__median3_d(a,b,c) result(med)
 real(kind(0d0)), intent(in) :: a,b,c
 real(kind(0d0)) :: med
   if(a < b) then
      if(c < a) then
         med = a; RETURN
      else if(c < b) then
         med = b; RETURN
      end if
   else if(a < c) then
      med = a; RETURN
   else if(c < b) then
      med = b; RETURN
   end if
   med = c
end function

recursiveは副プログラムを再帰的に呼び出す時に使う属性です。が、使う場面は非常に限られるので「そんなこともできるのか」くらいで大丈夫です。

sortモジュール

これもモジュールにしておきます。
sort.f90
ファイルを開くと"USAGE"のところに使い方のサンプルが書いてあるので参考にして下さい。

ランダムな向きの３次元ベクトルを作る

単位球面上に一様に分布する乱数(ベクトル)を生成する。
単純に極座標表示でθとφを一様乱数にしても、球面上に一様分布しないので注意。

 !----------------------------------
 ! Generate a random vector on the unit sphere
 ! Marsaglia algorithm (1972)
 !----------------------------------
 subroutine rand_sphere(vec)
  real(kind(0d0)), intent(out) :: vec(3)
  real(kind(0d0)) :: r2, r, xx, yy
    r2 = 2.0d0
    do while(r2 > 1.0d0)
       xx = 2.0d0*ran__drnd() - 1.0d0
       yy = 2.0d0*ran__drnd() - 1.0d0
       r2 = xx**2 + yy**2
    end do
    r = sqrt(1.0d0 - r2)
    vec(1) = 2.0d0 * xx * r
    vec(2) = 2.0d0 * yy * r
    vec(3) = 1.0d0 - 2.0d0 * r2
 end subroutine

（『ran__drnd()』は一様乱数。これはranpackを使用）
単位球面上に一様に分布させる方法は色々提案されているが、実用上はこのアルゴリズムが速い。
他のアルゴリズムや注意点を知りたい人は、 Wolframのサイトを見ると良い。