Fortranで遊ぶ

投稿日: 2021年5月31日 2022年6月7日
投稿者: 老猫こてつ.f90
カテゴリー: Fortran

Fortranは、1954年、IBMのジョン・ワーナー・バッカスによって考案された科学技術計算向けのプログラミング言語である。今やレガシー言語になってしまったが、気候モデルや、構造力学といった大規模な計算において、現在でもスーパーコンピュータ上で使われており、数値計算分野においてはずば抜けたプログラミング言語なのである。また、PythonやPHPのようなスクリプト言語ではなく、コンパイルして実行する。

Fortranには、その長い歴史においていくつかの種類が存在する。Fortran、FORTRAN77、Fortran90/95、Fortran2003、Fortran2008、Fortran2018などである。Fortran2003からはオブジェクト指向が取り入れられ、2008では並行計算が取り入れられた。

Fortranは、Java、C#、C/C++のように型宣言が必要なプログラミング言語である。使用する変数名はすべて、型を指定しなければならない。PHPやPythonのように、自由に使えるわけではなく、同じコード内でaという変数に数値を入れたり、文字列を入れたりすることができない。どちらがいいかは個人の好みだが、PHPやPythonでは意識しないと、バグを発生させる要因になりかねない。

Fortranでは、数値を扱う場合、単精度、倍精度、整数があり、混在させた式を扱う場合は、思った計算結果にならないという状況が発生する場合がある。数値計算において使用する変数および数値は型を統一しなければならない。また、倍精度型の変数に単精度の数値を代入すると面倒なことになりかねない。倍精度の計算であれば、

implicit none ! 暗黙の型宣言禁止
real(8) :: x  ! 変数xを倍精度型の変数として宣言

x = 1.0d0 ! 整数の数値であっても、小数点を付ける
x = 2.0d-6 ! 指数表記でも明確に分ける。2.0e-6と2.0d-6では異なる。

というようにする。厳格なのである。Pythonも科学技術計算ライブラリがあるなら、Fortranのように、型宣言する、しない、とすべきだろう。科学技術計算の多くは小数点第３桁までで十分なことが多いので、倍精度にすればいいのではないのである。

さて、現在、フリーで使えるFortranとしてgfortranがある。FORTRAN77をサポートするが、基本的にはFortran90/95である。Linux mintには、入っていなかったのでインストールした。

sudo apt install gfortran

これでインストールされたので、あとはVisual Studio Codeのエクテンションとして、

Modern Fortran
FORTRAN Intellisense

をインストールして終わりである。設定を調べると、コンパイラまでのパスを設定しろとあるが、別に設定しなくともコンパイルやインテリセンスは動作する。

Fortranのコードの書き方については、

Fortran演習 (地球惑星物理学演習)

を読んだだけである。PHPやJavascript、Pythonなど、いずれかのプログラミング言語の基本的なことを理解していれば、後述する他言語からのコード書き換えぐらいはできるようになるし、とりあえずコードを書くことができる。わからなければ、該当するところを読めばいいだけである。とりたてて深く学習しなかった。Fortranをはじめてから２日程度である。

いつものようにHello Worldを出力してみる。ファイル名を、hello.f90とする。hello.fとかだと、書き方が変わる。f90は大文字、小文字を問わない自由形式、fは固定形式で、インデント９文字、すべて大文字で書く。ま、小文字で書いてもエラーにならないけどね。

! はじめてのFortranのコード
program hello
    write(*, *) 'Hello World'
    stop
end program hello

コードを書いたら、コンパイルする。

gfortran -o hello hello.f90

コンパイル時に-oオプションを付けて出力先のファイルをhelloとした。オプションを付けなければ、a.outというファイル名になる。実行は、

./hello

とする。コンパイラなので実行時に、python hello.pyなどとしなくてもよい。パーミッションも自動的に実行権限が付与される。

「マイコン宇宙講座」に書かれているコード（N-BASIC）をFortranで書き直してみた。グレゴリオ暦の日付から修正ユリウス日を求めるものである。

メイン１つ、サブルーチン３つで構成している。それぞれコンパイルするのではなく、

gfortran -o m23 m23.f90 mjd.f90 julian.f90 jdate.f90

としてコンパイルする。メインからサブルーチン、サブルーチンから別のサブルーチンを呼び出しているので、一括でコンパイルしないとエラーになる。本来なら、サブルーチンではなく関数にするといいのかもしれないが、コードの中に書かなければいけない。これだと利用において汎用性に欠けるためにサブルーチンにしている。

メイン m23.f90

program m23
    implicit none

    real(8) :: jd, dy, dt
    real(8) :: year, month, day, dtime
    real(8) :: hour, minute, second

    write(*,*)
    write(*, fmt='(a)', advance='no') 'DATE AND TIME(JST) ? '
    read(*,*) dy, dt

    call ModifiedJulian(dy, dt, jd, year, month, day, hour, minute, second)
    dtime = hour / 24.0d0 + minute / 1440.0d0 + second / 86400.0d0
    write(*, *)
    write(*, '(" INPUT DATE = ", i5, " ", i2, " ", f8.5, " UT")') int(year), int(month), day + dtime - 0.375
    write(*, '("        MJD = ", f12.5, " UT")') jd
    write(*, *)

    call Jdate(jd, year, month, day)
    write(*, '("RETURN DATE = ", i5, " ", i2, " ", f8.5, " UT")') int(year), int(month), day

    call Julian(year, month, day, jd)
    write(*, '("        MJD = ", f12.5, " UT")') jd - 2400000.5
    write(*, *)

    stop
end program m23

サブルーチン mjd.f90

! 修正ユリウス日を求める
subroutine ModifiedJulian(dy, dt, jd, year, month, day, hour, minute, second)
    implicit none

    real(8) :: dy, dt
    real(8) :: jd, year, month, day
    real(8) :: hour, minute, second
    real(8) :: d1, d2
    integer :: s1

    s1 = sgn(dy)
    dy = abs(dy)

    ! 年月日および時刻を抽出する
    year = int(dy / 10000.0d0)
    d1 = dy - 10000.0d0 * year
    month = int(d1 / 100.0d0)
    day = d1 - 100.0d0 * month
    hour = int(dt / 10000.0d0)
    d2 = dt - 10000.0d0 * hour
    minute= int(d2 / 100.0d0)
    second = d2 - 100.0d0 * minute

    if (s1 < 0.0d0) then
        year = year * s1
    else
        if (year < 100.0d0) then
            year = year + 1900.0d0
        end if
    end if

    call julian(year, month, day, jd)

    jd = jd + hour / 24.0d0 + minute / 1440.0d0 + second / 86400.0d0
    jd = jd - 0.375d0

    jd = jd - 2400000.5d0

    return
    contains
    integer function sgn(x) result(y)
        implicit none

        real(8) :: x

        if (x < 0) then
            y = -1
        else
            y = 1
        end if

        return
    end function sgn
end subroutine ModifiedJulian

サブルーチン julian.f90

! ユリウス日を求める
subroutine Julian(year, month, day, jd)
    implicit none

    real(8) :: p5, jd
    real(8) :: year, month, day

    p5 = year + (month - 1.0d0) / 12.0d0 + day / 365.25d0
    if (month < 3) then
        month = month + 12.0d0
        year = year - 1.0d0
    end if

    if (year < 0.0d0) then
        jd = int(year * 365.25d0) + int(30.59d0 * (month - 2.0d0)) + day + 1721085.5d0
    end if

    if (p5 >= 0.0d0) then
        jd = int(year * 365.25d0) + int(30.59d0 * (month - 2.0d0)) + day + 1721086.5d0
    end if

    if (p5 >= 1582.78d0) then
        jd = int(year * 365.25d0) + int(year / 400.0d0) - int(year / 100.0d0) + int(30.59d0 * (month - 2.0d0)) + day + 1721088.5d0
    end if
    
    year = year * sgn(year)
    if (month > 12.0d0) then
        month = month - 12.0d0
        year = year + 1.0d0
    end if

    return
    contains
    integer function sgn(x) result(y)
        implicit none

        real(8) :: x

        if (x < 0) then
            y = -1
        else
            y = 1
        end if

        return
    end function sgn
end subroutine Julian

サブルーチン jdate.f90

! ユリウス日から年月日を求める
subroutine Jdate(jd, year, month, day)
    implicit none

    integer :: T(12) = (/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)
    real(8) :: jd, jj
    real(8) :: year, month, day
    real(8) :: r1, r2
    integer :: m

    jj = jd
    year = int(2.7379093d-3 * jd + 1858.877d0)
    month = 1.0d0
    day = 0.0d0

    call Julian(year, month, day, jd)

    r2 = jj - (jd - 2400000.5d0)
    if (mod(year, 4.0d0) == 0 .and. mod(year, 100.0d0) /= 0 .or. mod(year, 400.0d0) == 0) then
        T(2) = 29
    end if

    r1 = 0.0d0
    do m = 1, 12
        if (int(r2) - r1 - T(m) <= 0) then
            exit
        end if
        r1 = r1 + T(m)
    end do

    month = m
    day = r2 - r1
    T(2) = 28
    jd = jj

    if (month == 13) then
        year = year + 1.0d0
        month = month - 12.0d0
    end if

    return
end subroutine Jdate

実行結果は次のとおりになる。

DATE AND TIME(JST) ? 20210601,090000 

  INPUT DATE =  2021  6  1.00000 UT
         MJD =  59366.00000 UT

 RETURN DATE =  2021  6  1.00000 UT
         MJD =  59366.00000 UT

老猫こてつ.f90