はじめての Fortran プログラミング vscode 版

葛西真寿　弘前大学大学院理工学研究科

Jupyter Fortran Kernel を使った Fortran 90/95 プログラミング入門である

はじめての Fortran プログラミング – 弘前大学 Homw Sweet Home

の Visual Studio Code 版。

HEROIC 2021 推奨実行環境の確認

gfortran および Visual Studio Code を使います。情報基盤センターPC実習室の Ubuntu, Windows 10, macOS にはインストール済みです。

自分のパソコンにインストールする場合，センターでは Windows では Chocolatey，macOS のでは Homebrew というパッケージ管理システムを利用してインストールしていますので，以下を参考にしてください。

拡張機能のインストールと設定

Visual Studio Code をインストールしたままでは，メニューも日本語化されていません。プログラミング用テキストエディタとして統合環境的に活用するため，いくつか拡張機能をインストールし，簡単な設定を行います。

以下のページを参考に，機能拡張のインストールと設定（settings.json の編集）を行ってください。何やら難しいことが書いていそうですが，基本的にターミナル画面にコピペすれば作業完了です。

Visual Studio Code の拡張機能と設定

参考：Jupyter Fortran kernel ではだめなのか

弘大 JupyterHub では Fortran kernel も使えますし，以下のように Jupyter Notebook の Fortran kernel の使用例もまとめてあります。

はじめての Fortran プログラミング

Jupyter Notebook 環境そのものは（特に初心者にとって）大変便利なのですが，残念なことに，本稿執筆時点での jupyter-fortran-kernel は，read 文に対応していないという致命的な弱点があります。（キーボードからの入力を読み取れず，エラーとなってしまう。）

ですので，Fortran プログラミングの教育用の環境としては，Jupyter Notebook ではなく，Visual Studio Code の Modern Fortran とCode Runner 機能拡張を利用した環境を使ってみることにします。

最初のプログラミングセッション

Visual Studio Code を起動します。
アプリのアイコンをダブルクリックして起動します。
また，ターミナル（Windows では「PowerShell」，Ubuntu では「端末」，macOS では「ターミナル」）のプロンプトに code と入力してエンターキーを押しても起動します。
「ファイル」メニューから「新規ファイル」を選択します。

「Untitled-1」などというファイル名で新規ファイルができますので，（まだ何も書いていませんが）「ファイル」メニューから「保存」または「名前を付けて保存… 」を選択し，「ai1.f90」というファイル名で保存します。

拡張子は「.f90」としています。
Visual Studio Code は，この拡張子をみて，あぁこれから書く内容は Fortran 90/95 のプログラム文なんだなぁと認識し，いろいろなお節介をやいてくれます。

プログラムの実行

以下のプログラム文を一字一句間違いなく打ち込みます。コピペしてもいいです。

ai1.f90

program ai1
  boku = 2.5
  kanojo = 0
  love = boku * kanojo
  print *, love
end program

右上の ▷ マークをクリックすると，まだファイルが保存されていない場合は自動的に保存して，ウインドウの下にターミナル画面が開き，自動でコンパイルと実行を行います。

下のターミナルウインドウに，コンパイル済みのコマンド ./ai1 を打ち込んで，もう一度実行することもできます。

implicit none にする動機付け

上のプログラム ai1.f90 は無事コンパイルできてエラーなく実行されたが，結果が 0 というのは，あまりに寂しい。

以下のように，少し変更したプログラムを ai2.f90 として保存し，右上の ▷ マークをクリックして実行すると…

ai2.f90

program ai2
  watashi = 2.5
  kareshi = 0.4
  love = watashi * kareshi
  print *, love
end program

結果はやはり 0 となってしまう。おかしい… 2.5 に 0.4 をかければ 1.0 になるはずだが…

暗黙の型宣言という歴史的慣習

歴史的慣習により，FORTRAN では，黙っていれば i, j, k, l, m, n で始まる変数は整数，それ以外のアルファベットで始まる変数は実数として扱われてきた。

上記の例題 program ai2 では，kareshi = 0.4 は実は kareshi = 0 のこと。

このような思いもかけないトラブル（アイにトラブルはつきもの？）を防ぐ意味も兼ねて，最近の教科書では，この暗黙の型宣言を無効にするように，プログラム文の冒頭に

implicit none

と書くことを推奨している。

コーディングスタイル

ここでは，以下のサイトなどを参考にして，それらしいスタイルで書いてみます。

参考：

ai.f90

program ai
  ! ! から始まる行はコメント行。実行に影響しません。
  ! 暗黙の型宣言は無効にして，
  ! プログラムで使う変数は全て宣言します。
  implicit none
  real :: watashi, kareshi, love
  watashi = 2.5
  kareshi = 0.4
  love = watashi * kareshi
  print *, love
end program ai

上のプログラム ai.f90 を例に，現代的視点から推奨されている書き方を以下に例示します。

Fortran 主プログラムは，program ... ではじまります。
Fortran は大文字と小文字を区別しません。ここでは原則として全て小文字で書くことにします。
! からはじまる行はコメント行です。コメント行はそれを読む人間のための注釈であり，プログラムの実行時には無視されます。
5行目でまず暗黙の型宣言を無効にします。implicit none と宣言したら，以降はプログラム文内で使用する変数は，全て以下の例のように宣言する必要があります。
次にこれから使う変数の型を real（実数）として宣言します。
7行目では，数値 2.5 を変数 watashi に代入しています。人一倍どころか，人より 2.5 倍の思いがあることを示しています。
8行目では，数値 0.4 を変数 kareshi に代入しています。
9行目では，変数 watashi に代入された数値と，変数 kareshi に代入された数値をかけた答えを，変数 love に代入しています。
10行目では，変数 love の値を出力しています。標準出力（画面）への出力は print に統一します。
主プログラムの最後は end program ... で終わります。

実行結果は以下のとおりです。めでたく 1.00000000 と出力されています。

エンゲル係数を求めるプログラム例

エンゲル係数とは，家計の消費支出に占める飲食費（食料費，食費）の割合（パーセント単位）です。

以下の engel.f90 は消費支出 110,000円，飲食費 37,000円の人のエンゲル係数を求めるプログラム例です。コピペして engel.f90 として保存し，実行してください。

engel.f90

program engel
  ! エンゲル係数を求める
  implicit none
  ! 消費支出 spending, 飲食費 food, エンゲル係数 eng
  real :: spending, food, eng
  spending = 110000
  food = 37000
  eng = food * 100 / spending
  print *, eng
end program engel

エンゲル係数の値が 33.6363640 と出力されています。

データの型

実数型

real

上の例では，real（実数）として変数 spending, food, eng を宣言していました。 real は real(4) と同等であり，real で宣言される実数が表示できる最大値・最小値は有効桁数は以下のようにして確認できます。コピペして real.f90 として保存し，実行してみてください。

real.f90

program real
  implicit none
  real :: r
  print *, '1 / 3 = ', 1.0 / 3.0
  print *, 'real の最大値は', huge(r)
  print *, 'real の最小値は', tiny(r)
end program real

実行結果は

 1 / 3 =   0.333333343    
 real の最大値は   3.40282347E+38
 real の最小値は   1.17549435E-38

real(8), real(16)

real では表せない非常に大きい実数や，より有効桁数の多い実数を使う場合は，以下のように real(8) や real(16) と宣言します。表示できる最大値・最小値は有効桁数は以下のようにして確認できます。コピペして real816.f90 として保存し，実行してみてください。

real816.f90

program real816
  implicit none
  real(8)  :: r8
  real(16) :: r16
  print *, '1 / 3 = ', 1.0_8 / 3.0_8
  print *, 'real(8) の最大値は', huge(r8)
  print *, 'real(8) の最小値は', tiny(r8)
  print *, ' '
  print *, '1 / 3 = ', 1.0_16 / 3.0_16
  print *, 'real(16) の最大値は', huge(r16)
  print *, 'real(16) の最小値は', tiny(r16)  
end program real816

実行結果は以下のようになります。

 1 / 3 =   0.33333333333333331     
 real(8) の最大値は   1.7976931348623157E+308
 real(8) の最小値は   2.2250738585072014E-308
  
 1 / 3 =   0.333333333333333333333333333333333317      
 real(16) の最大値は   1.18973149535723176508575932662800702E+4932
 real(16) の最小値は   3.36210314311209350626267781732175260E-4932
$

参考：selected_real_kind

例えば有効桁数 12 桁の実数を表したい場合，real, real(8), real(16) のどれを使えばいいか調べるのが面倒な場合は

integer, parameter :: dp = selected_real_kind(12)
  real(dp)...

のようにします。以下のプログラムをコピペし，printau1.f90 として保存し，実行してみてください。

printau1.f90

program printau1
  implicit none
  ! 天文単位は 12 桁の定義定数（誤差無し）
  integer, parameter :: dp = selected_real_kind(12)
  real(dp), parameter :: au1 = 149597870700.0_dp
  print *, "1 天文単位は", au1, " m です。"
end program printau1

実行結果は以下のようになります。

 1 天文単位は   149597870700.00000       m です。

整数型

integer

整数型は integer です。表示できる最大値は有効桁数は以下のようにして確認できます。

int.f90

program int
  implicit none
  integer :: i
  print *, '1 / 3 = ', 1 / 3
  print *, 'integer の最大値は', huge(i)
end program int

実行結果は以下のようになります。

 1 / 3 =            0
 integer の最大値は  2147483647

integer(8), integer(16)

より桁数の多い整数を扱う場合は以下のように integer(8) や integer(16) を使います。表示できる最大値は有効桁数は以下のようにして確認できます。

int816.f90

program int816
  implicit none
  integer(8)  :: i8
  integer(16) :: i16 
  print *, 'integer(8) の最大値は ', huge(i8)
  print *, 'integer(16) の最大値は', huge(i16)
end program int816

実行結果は，

 integer(8) の最大値は   9223372036854775807
 integer(16) の最大値は  170141183460469231731687303715884105727

参考：selected_int_kind

例えば有効桁数 12 桁（以上）の整数を表したい場合，integer, integer(8), integer(16) のどれを使えばいいか調べるのが面倒な場合は

integer, parameter :: ik = selected_int_kind(12) 
integer(ik)...

のようにします。

printau.f90

program printau
  implicit none
  ! 天文単位は 12 桁の定義定数（誤差無し）
  integer, parameter :: ik = selected_int_kind(12)
  integer(ik), parameter :: au = 149597870700_ik
  print *, "1 天文単位は", au, " m です。"
end program printau

実行結果は，

 1 天文単位は         149597870700  m です。

文字型

文字型は以下のように character で宣言し，文字の最大数を括弧で書きます。

変数 moji に文字列を代入する場合は，以下のように ' または " で囲みます。

printchar.f90

program printchar
  implicit none
  ! 文字型の変数 moji（長さ12）の宣言
  character(12) :: moji
  moji = 'engel keisu '
  print *, moji
end program printchar

式と演算

主な演算子

足し算は +，引き算は -，掛け算は *，割り算は / です。

優先的に計算したい箇所は丸括弧 ( ) でグループ化します。以下の例を参照：

yusen.f90

program yusen
  print *, 2 +  50 - 5 * 6 /2
  print *, 2 + (50 - 5 * 6)/2
end program yusen

// は文字列の連結演算に使います。以下の例を参照。

charadd.f90

program charadd
  implicit none
  character(6) :: x0, y0
  x0 = 'engel '
  y0 = 'keisu '
  print *, x0 // y0
end program charadd

文字列と数値を直接連結することはできません。

文字列と数値を連結などしなくても，print 文は，カンマ (,) で区切ることによって複数の項目（変数）を表示することができます。

engel2.f90

program engel2
  implicit none
  real :: spending, food, eng
  spending = 110000
  food = 37000
  eng = food * 100 / spending
  print *, 'エンゲル係数は ', eng, 'です。'
end program engel2

実行結果：

 エンゲル係数は    33.6363640     です。

除算 / の例を以下に示します。

divide.f90

program divide
  print *, '17 / 3      = ', 17 / 3
  print *, '17.0 / 3    = ', 17.0 / 3
  print *, '17.0_8 / 3  = ', 17.0_8 / 3
  print *, '17.0_16 / 3 = ', 17.0_16 / 3
end program divide

実行結果：

 17 / 3      =            5
 17.0 / 3    =    5.66666651    
 17.0_8 / 3  =    5.6666666666666670     
 17.0_16 / 3 =    5.66666666666666666666666666666666692

2行目では，整数を整数で割るので，小数分を切り捨てた整数値を返します。
3行目では，17.0 と最後に小数点をつけたので，実数型（real）の値を返します。この例では，小数点以下6桁くらいまでの精度です。
4行目は 17.0_8 で real(8) の実数としています。この例では，小数点以下15桁くらいまでの精度です。
5行目では，17.0_16 で real(16) の実数としています。この例では，小数点以下33桁くらいまでの精度です。

べき乗は ** です。

以下の例では $2^4=2\times 2\times 2\times 2$ と $2^{1/2}=\sqrt{2}$ の値を表示します。

power.f90

program power
  print *, '2**4   =', 2**4
  print *, '2**0.5 =', 2**0.5
end program power

 2**4   =          16
 2**0.5 =   1.41421354

べき乗に限らず，あまりにも桁数が多くなる計算すると，基本整数型（integer）の範囲を超えてしまい，おかしな値を出力することがあるので注意。

以下の例では， $2^{31}$ あたりで答えが変になっています。

oor.f90

program oor
  ! out of range
  print *, '2**29 = ', 2**29
  print *, '2**30 = ', 2**30
  print *, '2**31 = ', 2**30 * 2
  print *, '2**32 = ', 2**30 * 2 * 2
end program oor

実行結果：

 2**29 =    536870912
 2**30 =   1073741824
 2**31 =  -2147483648
 2**32 =            0

非常に大きい（または極めて小さい）数値を扱う場合や，real の有効桁数約6桁では足りない場合，例えば有効桁数を15（以上）とする場合は，以下のようにします。

上の例で整数型の場合に正しい答えが出なかった $2^{31}, 2^{32}$ なども，以下の例のように有効桁数を多くとることによって精度よく出力されます。

precision.f90

program precision 
  implicit none
  ! 15桁（以上）の有効桁数を持つ実数・整数の定義用
  integer, parameter :: dp = selected_real_kind(p=15)
  integer, parameter :: ik = selected_int_kind(15)
  real(dp) :: sqrtwo
  sqrtwo = sqrt(2.0_dp)
  print *, 'sqrt(2.0)    = ', sqrt(2.0)
  print *, 'sqrt(2.0_dp) = ', sqrtwo
  print *, '2**31 = ', 2_ik**31
  print *, '2**32 = ', 2_ik**32
end program precision

実行結果：

 sqrt(2.0)    =    1.41421354    
 sqrt(2.0_dp) =    1.4142135623730951     
 2**31 =            2147483648
 2**32 =            4294967296

練習 1

1 天文単位（1 au）の距離を光速 c で進むのに要する時間は（何時間）何分何秒か。ただし，

1 天文単位は 12 桁の定義定数（誤差無し）で 149597870700 m
光速 c は 9 桁の定義定数（誤差無し）で 299792458 m/s である。

参考：

対話型プログラム

上記の例では，プログラムの中で消費支出や飲食費を決めていましたが，今度はプログラムを実行する人が入力できるようにしてみます。

以下のプログラム例では，read *, spending でキーボードからの入力値を変数 spending に入れます。

engelread.f90

program engelread
  implicit none
  real :: spending, food, eng
  print *, "消費支出? "
  read *, spending
  print *, "飲食費? "
  read *, food
  eng = food * 100 / spending
  print *, "エンゲル係数は", eng
end program engelread

実行結果：

 消費支出? 
110000     <----- キーボードから入力
 飲食費? 
37000      <----- キーボードから入力
 エンゲル係数は   33.6363640

関数の定義と呼び出し

（これまでの例では1回しか計算していませんが）よく使う計算を「関数」として定義する例です。

以下の例では，エンゲル係数を計算する関数 calc() を定義し，主プログラムでその関数を呼び出しています。

calc が実数型であることは，主プログラム内でも，関数副プログラム内でも real として宣言する必要があります。

engelcalc.f90

program engelcalc
  implicit none
  real :: spending, food, eng, calc
  spending = 110000
  food = 37000
  eng = calc(spending, food)
  print *, "エンゲル係数は", eng
end program engelcalc

function calc(spending, food)
  ! 消費支出 spending と食費 food から
  ! エンゲル係数を計算
  implicit none
  real :: calc, spending, food
  calc = food * 100 / spending
end function calc

実行結果：

 エンゲル係数は   33.6363640

練習 2

消費支出を入力すると，エンゲル係数が 20 ~ 40 となる飲食費を表示するプログラムをつくりなさい。

組み込み関数

関数は上記の例のように自分で定義することもできますが，あらかじめ用意されている「組み込み関数」を使うこともできます。

以下のプログラムでは，組み込み関数の例として，

小数点以下を切り捨てる floor()
小数点以下を四捨五入する nint()

を使用しています。ちなみに，

小数点以下を切り上げる関数は ceiling()
sqrt(x) は平方根，
三角関数は sin(x), cos(x), tan(x)
指数関数は exp(x), 自然対数は log(x), 常用対数は log10(x)
int(x) は実数を整数に変換し，real(n) は整数を実数に変換します。 dble(x) は倍精度実数に変換するのでした。

func.f90

program func
  implicit none
  real :: spending, food, eng
  spending = 110000
  food = 37000
  eng = food * 100 / spending
  print *, 'エンゲル係数は', eng
  print *, '切り捨てでは  ', floor(eng)
  print *, '四捨五入では  ', nint(eng)
end

実行結果：

 エンゲル係数は   33.6363640    
 切り捨てでは            33
 四捨五入では            34

Fortran の組み込み関数

Fortran で使える組み込み関数については，インターネット検索すると，いろいろ見つかります。たとえば，以下の JAMSTEC のサイトが参考になります。

JAMSTEC | 海洋研究開発機構 | ジャムステック
- Fortran90 コース
  - 4.1 組み込み関数

練習 3

小数点以下を切り上げたエンゲル係数の値を表示させるように上記のプログラムを変更しなさい。

条件分岐

条件分岐とは，例えばエンゲル係数の値によって実行文を変えることです。例を示します。

以下のプログラムでは，エンゲル係数の値が 40 以上だと，「高い値です。」と表示します。

engelif.f90

program engelif
  implicit none
  real :: spending, food, eng
  print *, "消費支出? "
  read *, spending
  print *, "飲食費? "
  read *, food
  eng = food * 100 / spending
  if (eng >= 40) then
    print *, "エンゲル係数は", eng, "です。高い値です。"
  else
    print *, "エンゲル係数は", eng, "です。"
  end if
end program engelif

実行結果：

 消費支出? 
110000     <----- キーボードからの入力
 飲食費? 
45000      <----- キーボードからの入力
 エンゲル係数は   40.9090919     です。高い値です。

条件分岐の if 文の書式は以下のようになっています。

if (条件式1) then
  条件式1 が満たされた場合に実行する文
else if (条件式2) then
  条件式2 が満たされた場合に実行する文
else if (条件式3) then
  条件式3 が満たされた場合に実行する文
...
else
  それ以外の場合に実行する文
end if

関係演算子

上の例では，eng >= 40 つまり eng が 40 以上のとき，という関係演算子を使いました。他の例は以下のとおりです。

a == b    a と b が等しいとき
a /= b    a と b が等しくないとき
a < b    a が b より小さいとき
a <= b    a が b 以下のとき
a > b    a が b より大きいとき
a >= b    a が b 以上のとき

練習 4

上記のプログラム engelif.f90 の条件分岐を，エンゲル係数が

40 以上の場合には「高い値です。」
40 未満 20 以上の場合には「正常値です。」
20 未満の場合には「低い値です。」

と表示させるように変更し，engelif2.f90 として保存し，実行してみなさい。

繰り返し処理

以下の例では，1行目で代入された消費支出に対してエンゲル係数を計算します。

1回だけ計算して表示するのではなく，飲食費の初期値（ここでは 30000 円）から，エンゲル係数を計算して表示し，条件式（ここでは 50000 円以下）が成り立つ場合に，5000 円ずつ増やしながら計算と表示を続けます。

表示の際に，実数を整数に変換する関数 int や小数点以下を四捨五入する関数 nint を使っています。

engelloop1.f90

program engelloop1
  implicit none
  real :: spending, food, eng
  spending = 110000
  food = 30000
  print *, "消費支出", int(spending), &
            "円に対するエンゲル係数の値"
  do while (food <= 50000)
    eng = food * 100 / spending
    print *, "  飲食費", int(food), "円の場合: ", nint(eng)
    food = food + 5000
  end do
end program engelloop1

6行目の最後に & をつけ，7行目が継続行であることを示しています。1行が長くなるときに使います。

実行結果：

 消費支出      110000 円に対するエンゲル係数の値
   飲食費       30000 円の場合:           27
   飲食費       35000 円の場合:           32
   飲食費       40000 円の場合:           36
   飲食費       45000 円の場合:           41
   飲食費       50000 円の場合:           45

do while 文の書式は以下の通りです。

do while (条件式)
    条件式が成り立つ場合に実行される文
end do

この繰り返し処理の do while 文の中に出てくる

food = food + 5000

という書き方は，最初はびっくりするかも知れません。

数学の等式では，両辺が等しくなることはないので意味のない式ですが，Fortran では，= は左辺と右辺が等しいことを示す「等号」ではなく，右辺で計算した数値を左辺に「代入」することを表します。したがってこの式は，

変数 food の現在の値に 5000 を足した値をあらためて変数 food に代入する

という意味になります。

同様の内容を do 文で書いた例が以下のプログラムです。

engelloop2.f90

program engelloop2
  implicit none
  integer :: i
  real :: spending, food, eng
  spending = 110000
  print *, "消費支出", int(spending), &
           "円に対するエンゲル係数の値"
  do i = 30000, 50000, 5000
    food = float(i)
    eng = food * 100 / spending
    print *, "  飲食費", int(food), "円の場合: ", nint(eng)
  end do
end program engelloop2

実行結果：

 消費支出      110000 円に対するエンゲル係数の値
   飲食費       30000 円の場合:           27
   飲食費       35000 円の場合:           32
   飲食費       40000 円の場合:           36
   飲食費       45000 円の場合:           41
   飲食費       50000 円の場合:           45

do 文の書式は以下の通りです。

do ivar = istart, iend, istep
    ivar が istart から iend の間，実行される文
end do

ivar は整数変数で，その初期値は istart。文が実行されるごとに，ivar の値は istep だけ加算され，iend の値になるとループが終了します。istep を省略した場合は 1 とみなします。

上記の例では，念のため，基本整数型（integer）の i を実数型（real）に変換する組み込み関数 float() を使っています。float() は real() と同じです。

練習 5

以下の値を求めるプログラムを作成し，実行してみなさい。

$\sum_{n=1}^{10} n^2$

配列とファイル

配列

成分が5個の配列を宣言し，各成分に数値を入れたり表示したりする例です。

array.f90

program array
  implicit none
  integer :: a(5), b(5), i
  ! a の各成分を設定
  a = [5, 4, 3, 2, 1]
  print *, '全成分を一挙に表示'
  print *, a
  print *, '縦に表示'
  do i = 1, 5
    print *, a(i)
  end do
  ! b の全成分を 0 に
  b = 0
  ! 第3成分の値を設定
  b(3) = 3
  print *, '1つだけ 0 でない成分をもつ配列'
  print *, b
end program array

実行結果：

 全成分を一挙に表示
           5           4           3           2           1
 縦に表示
           5
           4
           3
           2
           1
 1つだけ 0 でない成分をもつ配列
           0           0           3           0           0

ファイルから配列へ読み込み

以下のようなテキストファイル mydata.txt を用意します。

$ cat mydata.txt
1
2
3
4
5
6

以下のプログラムでは，ファイル mydata.txt の数値を成分が6個の配列 a に読み込み，成分の総和，平均，最大値，最小値を出力します。

readdata1.f90

program readdata1
  implicit none
  integer :: row
  real :: a(6)
  open(10, file = 'mydata.txt')
  read(10, *) a
  do row = 1, 6
    print *, a(row)
  end do
  ! 成分の和は sum(a)
  print *, '成分の総和は', sum(a)
  ! 成分の個数は size(a)
  print *, '成分の平均は', sum(a)/size(a)
  print *, '最大値は    ', maxval(a)
  print *, '最小値は    ', minval(a)
end program readdata1

実行結果：

   1.00000000    
   2.00000000    
   3.00000000    
   4.00000000    
   5.00000000    
   6.00000000    
 成分の総和は   21.0000000    
 成分の平均は   3.50000000    
 最大値は       6.00000000    
 最小値は       1.00000000

次は，以下のような「6行2列」のデータが入ったテキストファイル mydata2.txt を読み込む例です。

$ cat mydata2.txt
1  1
2  4
3  9
4 16
5 25
6 36

readdata2ng.f90

program readdata2ng
  implicit none
  integer :: row, col
  integer :: b(6, 2)
  open(10, file = 'mydata2.txt')
  read(10, *) b
  do row = 1, 6
    print *, (b(row, col), col = 1, 2)
  end do
end program readdata2ng

実行結果：

           1           4
           1          16
           2           5
           4          25
           3           6
           9          36

上の結果からわかるように，dimension(6, 2) の配列 b を用意して read(10, *) b とするだけでは，期待したようには読み込まれません。

ちょっとトリッキーですが，dimension(2, 6) の配列に読み込ませて，その転置をとってみました。以下の例を参照。

readdata2.f90

program readdata2
  implicit none
  integer :: row, col
  integer :: b(2, 6), c(6, 2)
  open(10, file = 'mydata2.txt')
  ! あえて dimension(2, 6) の配列 b に読み込む
  read(10, *) b
  ! b の転置行列を c に代入
  c = transpose(b)
  do row = 1, 6
    print *, (c(row, col), col = 1, 2)
  end do
end program readdata2

実行結果：

           1           1
           2           4
           3           9
           4          16
           5          25
           6          36

ファイルへ書き込み

ファイルの書き込み例です。

writedata.f90

program writedata
  implicit none
  integer :: i
  real :: x
  open(12, file = 'myoutput.txt')
  do i = 1, 10
    x = i * 0.1
    write(12, *) x, sin(x)
  end do
end program writedata

書き込まれたファイルの内容を表示：

$ cat myoutput.txt
  0.100000001       9.98334214E-02
  0.200000003      0.198669329    
  0.300000012      0.295520216    
  0.400000006      0.389418334    
  0.500000000      0.479425550    
  0.600000024      0.564642489    
  0.699999988      0.644217670    
  0.800000012      0.717356086    
  0.900000036      0.783326924    
   1.00000000      0.841470957    
$

書式指定

画面への出力の際，print *, ... などと書いていたが，この * の部分について。

* は（特に書式をしないで）おまかせで出力させる。書式を設定したいときには * の場所に以下のように書きます。

浮動小数点数の指定

実数を表示する際，全体の文字数と小数点以下の桁数を指定して表示させる例です。

例えば '(f8.5)' は全体で 8 桁の表示範囲をとり，小数点以下は 5 桁表示します。

formatf.f90

program formatf
  implicit none
  character(10) :: ten = '1234567890'
  ! 15桁（以上）の有効桁数を持つ実数の定義
  integer, parameter :: dp = selected_real_kind(p=15)
  real(dp) :: sqrtwo
  sqrtwo = sqrt(2.0E0_dp)
  print '(2a10, a)', ten, ten, ' 桁数'
  print '(f8.5)', sqrtwo
  print '(f18.15)', sqrtwo
  ! 入れ子にする場合は" と ' で使い分ける
  print '("2の平方根は", f18.15, " です。")', sqrtwo
end program formatf

実行結果：

12345678901234567890 桁数
 1.41421
 1.414213562373095
2の平方根は 1.414213562373095 です。

練習 6

engel2.f90 の出力を，エンゲル係数の値を（四捨五入して）小数点以下1桁まで表示するように変更しなさい。

指数表式の指定

非常に大きい数や小さい数を表示するときには， $1.9891\times 10^{30}$ や $6.67430\times 10^{-11}$ などのような指数表示をします。

例えば '(e9.2)' では，全体で 9 桁の表示分を確保し，小数点以下は 2 桁表示します。

このままだと整数部分が 0 で始まるので，'(1pe9.2)' としてみます。

formate.f90

program formate
  implicit none
  character(10) :: ten
  ! 光速は9桁の定義定数（誤差無し）
  integer, parameter :: c = 299792458
  ! 9桁（以上）の有効桁数を持つ実数の定義
  integer, parameter :: dp = selected_real_kind(p=9)
  real(dp) :: val
  val = c
  ten = '1234567890'
  print '(2a10, a)',  ten, ten, ' 桁数'
  print '(e9.2)', val
  print '(1pe9.2)', val
  print '(1pe15.8)', val
  ! 入れ子にする場合は" と ' で使い分ける
  print '("光速は", 1pe15.8, " m/s です。")', val
end program formate

実行結果：

12345678901234567890 桁数
 0.30E+09
 3.00E+08
 2.99792458E+08
光速は 2.99792458E+08 m/s です。

文字型の書式指定

formata.f90

program formata
  implicit none
  character(10) :: ten
  ten = '1234567890'
  print '(a, a)', ten, ten
  print '(2a)', ten, ten
  print '(2a5)', ten, ten
  print '(a6,a4)', ten, ten
end program formata

実行結果：

12345678901234567890
12345678901234567890
1234512345
1234561234

整数の書式指定

formati.f90

program formati
  implicit none
  integer :: v(4), m(4, 4)
  integer :: row, col
  ! v の成分
  v = [1, 23, 456, 7890]
  do row = 1, 4
    print '(i5)', v(row)
  end do
  ! 4行4列の成分を全て 0 に
  m = 0
  ! 対角成分
  m(1,1) = 1
  m(2,2) = 23
  m(3,3) = 456
  m(4,4) = 7890
  ! 書式指定で小綺麗に表示
  do row = 1, 4
    print '(4i5)', (m(row, col), col = 1, 4)
  end do
end program formati

実行結果：

    1
   23
  456
 7890
    1    0    0    0
    0   23    0    0
    0    0  456    0
    0    0    0 7890