ERA-interim再解析バイナリをFortranで読み込む方法
日常生活において再解析データのお絵描きをしたい時はGrADSを使えば良いのだが，たまにちょっと複雑な解析をしようとする場合はFortran等で計算した方が楽なことがある。言語はFortranでなくても良いのだが，伝統的にみんな使っていて情報量が多いという事情からFortranでやることにする。またERA-interimでなくても，JRA-55だろうとNCEPだろうとほぼ同じようにしてできる(はず)。
(A) バイナリファイル(.grd,.bin)がなく，GRIBファイル(.grb,.grib)しかない場合
自分で公式からダウンロードする場合，ほぼGRIB形式で得られる。GRIBはそのままFortranでは(簡単には)読めないので，バイナリ形式に変換(バイナリダンプ)する必要がある。GRIBにはGRIB-1形式とGRIB-2形式があり，それぞれwgribとwgrib2という異なるコマンドを使ってダンプする。やり方は早崎先生のサイト(ここやここ)に詳しい。気が向いたらこのブログでも記事を書く。
出力する際にはy方向の格納順(北→南 or 南→北)，z方向の格納順(上層→下層 or 下層→上層)に注意。GrADSでも読めるようにする場合はヘッダ無し(wgrib -nh / wgrib2 -no_header)にする。
(B) バイナリファイル(.grd,.bin)とGrADS用ctlファイルがある場合
とりあえず，ここを見れば良い。ただし「u(lat, lon, lev)」という部分は「u(lon, lat, lev)」の誤りだと考えられる。
気象データのバイナリファイルは時間，変数，高度，緯度，経度の辞書順に格納されている。経度は西→東の順。緯度と高度の格納順は場合によるが，これはそのバイナリファイルに対応するctlファイルを見れば良い。optionsにyrevと書かれていれば北→南，そうでなければ南→北。zrevと書かれていれば上層→下層。そうでなければ下層→上層である。
2010年1月のジオポテンシャルのバイナリ"z.201001.grd"から135E 30Nのジオポテンシャル高度を6時間ごと出力するサンプルコードは以下の通り。
program read_data
implicit none
integer, parameter :: nlon = 240
integer, parameter :: nlat = 121
integer, parameter :: nlev = 37
real(4), parameter :: cnst_g = 9.81
real(4), dimension(nlon,nlat) :: dat
open(11, file="z.201001.grd", access="direct", form="unformatted", recl=4*nlon*nlat)
z = 22
do t = 1, 31*4
read(11, rec=(t-1)*nlev+z) dat
write(*,*) dat(91,41)/cnst_g
end do
close(11)
z=22は500hPa，91が135E，41が30Nに対応する。reclの4*は4ビットバイナリだから。それに合わせてdatもreal(4)としている。reclの長さは必ずしもこうである必要はなく，recl=4*nlon*nlat*nlev*ntimとしてread(11, rec=1)として一括で読むこともできる(メモリを食うが)。あるいはrecl=4，rec=(t-1)*nlev*nlat*nlon+(z-1)*nlat*nlon+(y-1)*nlon+xとして1つずつを読み込むこともできる。
コンパイルする時には，コンパイルオプションをつける必要があり
gfortran -fconvert=big-endian -frecord-marker=4 hoge.f90
ifort -convert big_endian -assume byterecl hoge.f90
バイナリファイルにはビッグエンディアンとリトルエンディアンと呼ばれる2つの種類のものがある。気象データのバイナリは通常ビッグエンディアンが用いられるので，前者のオプションをつける必要がある。
またreclの単位にはwordとbyteがあり，ここではbyteで読み出したいので後者のオプションをつける。gfortranの場合は -frecord-marker=4を指定しなくてもうまくいくが，ifortの場合は -assume byterecl は必須。
