CReSS 初期値・境界値変換プログラム
本ファイルは京都大学生存圏研究所の グローバル大気観測データアーカイブプロジェクト にて学術目的で配信されている気象庁数値予報 GPV データから CReSS の初期値・境界値用大気場データに変換するプログラムの説明である.
GRIB2 データ変換用プログラム gpv_convert
必要ライブラリ・プログラム
- 数値解析ライブラリ STPK.
- 上記リンクからソースを取得し, ビルド, インストールして頂きたい.
- インストール時にエンディアンは "big endian" で統一して頂きたい. 以下では CReSS の統一エンディアンである big endian を前提に進める.
- GRIB2 フォーマット変換プログラム wgrib2
- 上記リンクからソースをダウンロード・ビルド後, 適切なパスの通ったディレクトリにインストールして頂きたい.
変換前データの準備
- グローバル大気観測データアーカイブプロジェクトページの 数値予報 GPV から必要なデータをダウンロードしておく. データの種類等は内部の README 参考.
- 地形データも同ページの こちら から適切なデータをダウンロードする. README 参考.
- 地形データの格子点はそれを利用する大気モデルの格子点と一致しているので, README 内のモデル詳細から地形データの格子点情報は取得できる.
実行までの前準備 (wgrib2 による変換)
- 本プログラムは wgrib2 にて GRIB2 からダイレクトアクセス形式になっていることを前提に処理を行う. よって, 先に wgrib2 にて下処理を行っておく.
以下, 変換したいオリジナル GRIB2 ファイルを "idata", wgrib2 で処理したデータを "gdata", 最終的に変換された CReSS 用のデータを "cfile" として話を進める.
$ wgrib2 idata | grep "anl" | wgrib2 -i idata -no_header -ieee gdata # idata の中の初期値のみ (anl) を取り出して, gdata に変換.
この作業を行うと, 以下のような標準出力がある:
1.1:0:d=2016010221:PRMSL:mean sea level:anl: 1.2:0:d=2016010221:PRES:surface:anl: 1.3:0:d=2016010221:UGRD:10 m above ground:anl: 1.4:0:d=2016010221:VGRD:10 m above ground:anl: 1.5:0:d=2016010221:TMP:1.5 m above ground:anl: 1.6:0:d=2016010221:RH:1.5 m above ground:anl: ...
この出力順序が変換後データ gdata での格納順である.
注意事項 1
- オリジナル GPV データは同じファイルに初期値と予報値が格納されている場合がある (MSM など). しかし, CReSS は 1 時刻 1 ファイルの入力を行うので, この時点で分離する.
なお, 予報値を変換する場合は
$ wgrib2 idata | grep ":3 hour fcst" | wgrib2 -i idata --no_header -ieee gdata
などと grep の引数を変更する. 実はこれ,
$ wgrib2 idata
としたときの出力結果の末尾からデータを取得しているだけである.
- もちろん, "anl" と "3 hour fcst" は別のファイルに出力しなければならない.
注意事項 2
- MSM は地表面データと 3 次元データが分離されているので, CReSS の初期値として地表面データも含める場合, 両方とも上記の処理を行わなければならない.
- このようにファイルが 2 種類に分けられている場合はそれぞれ別のファイルとして出力する.
- 以下, 地表面データは gsdata, 3 次元データは g3data として話を進める.
注意事項 3
変換するファイルが複数時刻に及ぶ場合, シェル上で以下のスクリプトを実行すれば自動的に変換が可能である (bash を想定).
#!/bin/bash for i in [data list] do wgrib2 $i | grep "anl" | wgrib2 -i $i -no_header -ieee "$i".out done
- [data list] には正規表現等を用いて変換するファイルを全て指定するようにする.
- 上の例では, 入力ファイル名の拡張子に .out を追加した出力ファイル名にしている.
最終的に 3 時間ごとのデータを 24 時間分, 地表面と 3 次元大気場分変換し終わると, gdata のファイル数は,
(24 / 3 + 1) * 2 = 18 ファイル
出力されているはずである.
変換プログラムの設定
変換プログラム gpv_convert は実行前に設定ファイル gpv_convert.nml を適切に設定する必要がある. 以下, 設定ファイルの内容である.
&basic uflag = 'u' ! 東西風の変数名 (任意) vflag = 'v' ! 南北風の変数名 (任意) wflag = 'w' ! 鉛直風の変数名 (任意) pflag = 'p' ! 圧力の変数名 (任意) tflag = 't' ! 温度の変数名 (任意) qflag = 'q' ! 相対湿度の変数名 (任意) zflag = 'z' ! ジオポテンシャル高度の変数名 (任意) sflag = 's' ! 追加変数 (現在は使用していない) dord = 'zuvtxqzuvtxqzuvtxqzuvtxqzuvtxqzuvtxqzuvtxqzuvtxqzuvtxqzuvtxqzuvtxqzuvtxqzuvtxzuvtxzuvtxzuvtx' (後述) nx1 = 241 ! 下層側の東西格子点数. (後述) ny1 = 253 ! 下層側の南北格子点数. nx2 = 241 ! 上層側の東西格子点数. ny2 = 253 ! 上層側の南北格子点数. nz = 16 ! 鉛直方向の格子点数. lonmin1 = 120.0 ! 下層側の西端経度 [deg] latmin1 = 22.4 ! 下層側の南端経度 [deg] lonmin2 = 120.0 ! 上層側の西端経度 [deg] latmin2 = 22.4 ! 上層側の南端経度 [deg] lonmax1 = 150.0 ! 下層側の東端経度 [deg] latmax1 = 47.6 ! 下層側の北端経度 [deg] lonmax2 = 150.0 ! 上層側の東端経度 [deg] latmax2 = 47.6 ! 上層側の北端経度 [deg] ipflag = 3 ! 水平方向の統一格子点数フラグ (後述) iplev = 1 ! 統一格子点区分 (後述) rhflag = .false. ! 水蒸気データが相対湿度 (RH) なら .true., 混合比なら .false. ofoot = '.out' ! 変換後のデータにつける拡張子 (任意) rlist = 'list.dat'! 変換する 3 次元大気データのリストファイル (後述) / &vercoord z = 100000.0, & & 97500.0, & & 95000.0, & & 92500.0, & & 90000.0, & & 85000.0, & & 80000.0, & & 70000.0, & & 60000.0, & & 50000.0, & & 40000.0, & & 30000.0, & & 25000.0, & & 20000.0, & & 15000.0, & & 10000.0 ! 各層の高度 (気圧の場合 [Pa], 幾何距離の場合 [m]). z の配列数は上で設定した nz に一致させる. / &terrain terfile = 'TOPO.MSM200409' ! 地形データのリストファイル (地表面情報を入れない場合は空白) nxt = 240 ! 地形データの東西格子点数 nyt = 252 ! 地形データの南北格子点数 lonmint = 120.0 ! 地形データの西端 [deg] latmint = 47.6 ! 地形データの北端 [deg] dlont = 1.25e-1 ! 地形データの東西格子点間隔 [deg] dlatt = -1.0e-1 ! 地形データの南北格子点間隔 [deg] / &surface surfile = 'lists.dat' ! 地表面データ (地表面情報を入れない場合は空白) d2ord = 'xpuvtqxxxx' ! 地表面データの格納順 (後述) nxs = 481 ! 地表面データの東西格子点数 nys = 505 ! 地表面データの南北格子点数 lonmins = 120.0 ! 地表面データの西端 [deg] latmins = 22.4 ! 地表面データの南端 [deg] lonmaxs = 150.0 ! 地表面データの東端 [deg] latmaxs = 47.6 ! 地表面データの北端 [deg] /
設定時の注意点
- 各データの格子点や東西南北端の位置座標は配信データ元の README に記載があるのでそれを参考すること.
- 変数 "dord", "d2ord" は wgrib2 で処理された gdata のうち, 何番目にどの変数が入っているかを示すフラグ.
例 : wgrib2 が実行されたとき, 以下のような標準出力が得られたとする:
1.1:0:d=2016010221:HGT:1000 mb:anl: 1.2:0:d=2016010221:UGRD:1000 mb:anl: 1.3:0:d=2016010221:VGRD:1000 mb:anl: 1.4:0:d=2016010221:TMP:1000 mb:anl: 1.5:0:d=2016010221:VVEL:1000 mb:anl: 1.6:0:d=2016010221:RH:1000 mb:anl: 1.7:0:d=2016010221:HGT:975 mb:anl: 1.8:0:d=2016010221:UGRD:975 mb:anl: 1.9:0:d=2016010221:VGRD:975 mb:anl: 1.10:0:d=2016010221:TMP:975 mb:anl: 1.11:0:d=2016010221:VVEL:975 mb:anl: 1.12:0:d=2016010221:RH:975 mb:anl: 1.13:0:d=2016010221:HGT:950 mb:anl: 1.14:0:d=2016010221:UGRD:950 mb:anl: 1.15:0:d=2016010221:VGRD:950 mb:anl: 1.16:0:d=2016010221:TMP:950 mb:anl: 1.17:0:d=2016010221:VVEL:950 mb:anl: ....
この場合, HGT (ジオポテンシャル高度), UGRD (東西風), VGRD (南北風), TMP (温度), VVEL (鉛直風), RH (相対湿度) を表している. その直後は高度 (気圧) である. ここで出力される先頭から 2 番目の整数が dord, d2ord の先頭からの順番に対応する. たとえば, uflag, vflag, wflag, pflag, tflag, qflag, zflag が上記のように設定されていたとすると, 上の標準出力から "鉛直風情報を除いて" CReSS 用の入力データを作成しようとする場合,
dord = 'zuvtxqzuvtxqzuvtx' ! 'x' はその部分だけ読まないことを表す.
と設定する. データが等圧面で格納されているため, pflag の変数は用いない. その代わり, 変数 z にその圧力情報を記載していく. その圧力情報は上の標準出力で出力される気圧順に記載する.
- d2ord も地表面データについて同様に設定を行えばよい.
- 第一セクションで変数名の末尾に 1 あるいは 2 がついている変数は 1 が 3 次元データの下層の情報, 2 が上層の情報についての設定項目である. 気象庁から提供されているデータはファイル容量節約のため, 下層と上層で格子点数が異なっている (間引かれている). これは配信元の README にも明記されており, 地表面に近い方を下層, 大気上端に近い方を上層とここでは定義し, 水平格子点数や東西南北端はそれぞれの層の情報を設定すること.
ipflag, iplev は CReSS 用に変換されるデータの格子点数, 水平格子間隔をどのデータを基準として作成するかを指定するフラグである. 前項目の通り, 3 次元データが上下層で格子点が異なる. それだけではなく, 地表面データや地形データも追加したデータを作成する場合, それぞれの格子点数, 東西南北端は共通ではないので, CReSS 用データに変換する際に, いずれかの格子点情報で統一する必要がある. ipflag はそのときの基準格子点をどのデータから用いるかを示す:
ipflag : 1 = 3 次元データ下層の水平格子点数, 2 = 3 次元データ上層の水平格子点数, 3 = 地表面データの水平格子点数, 4 = 地形データの水平格子点数.- iplev は ipflag = 1 or 2 の場合, 水平格子点の不連続が起こる dord の位置を指定する. たとえば, dord(200:200) で 199 番目までが下層, 200 番目から上層の場合を指定する. これについても, 配信元の README から不連続が起こる高度を調べ, wgrib2 が標準出力する情報から iplev を指定する.
rlist, surflist は wgrib2 で処理したファイルをリストしたテキストファイル名をそれぞれ指定する. これは ls コマンド等で簡単に行える.
$ ls gdata* > [rlist に指定したファイル名]
- もちろん, surflist を用いる場合は, rlist と surflist の各行に指定したファイルの時刻は合わせておかなければならない.
i
変換プログラムの実行
- 上記の準備が全て整ったら, 実際に変換プログラムを実行させる.
- 自分でいちから設定してもよいが, 同梱した設定ファイルは以下のテンプレートが用意されており, rlist, surflist に指定されたファイルでリストファイルを作成すればそのまま変換できるようになっている:
- gpv_convert_MSM_nosfc.nml : MSM で 3 次元データのみを変換.
- gpv_convert_MSM_sfc.nml : MSM で 地表面データ + 3 次元データを変換.
- gpv_convert_GSM_nosfc.nml : GSM で 3 次元データのみを変換.
- gpv_convert_GSM_sfc.nml : GSM で 地表面データ + 3 次元データを変換.
先の gdata を以下のコマンドで 2 つのリストファイルにする.
$ ls gsdata* > lists.dat # 地表面データ用 $ ls g3data* > list.dat # 3 次元大気用
これで, 地表面データと 3 次元データの処理ファイルリストができる.
[注意] : ここで, プログラムは 2 つのリストファイルを上から順に処理していくが, 各ファイルの同じ行に記載されているデータが同じ時刻のデータであると思い込んで処理している. なので, 同じ行に違う時刻のデータが入ってしまった場合, 不整合なデータが出来上がってしまうことになるので, 充分注意されたい. gdata のファイル名に時刻を使用すればそのような誤りは防ぐことができるはずである.
最後にプログラム gpv_convert を実行する.
$ ./gpv_convert < gpv_convert.nml
変換されたファイル名が順次画面に出力されるはずである.
Finished normally.
と出力されれば, 変換は終了である.
あとは, CReSS のプリプロセスにあうように, ファイル名を変更すれば終了である.
なお, 同梱した .ctl ファイルは同梱の namelist ファイルで変換された CReSS 用のデータのコントロールファイルを用意してあるので, 適宜時刻等を修正の後, grads, gpview 等で確認されたい.
NetCDF データ変換用プログラム gpv_convert_nc
[注意] : 本プログラムは NetCDF 化されたファイルのうち, MSM のみを変換するためのプログラムである.
必要ライブラリ・プログラム
- 数値解析ライブラリ STPK.
- 上記リンクからソースを取得し, ビルド, インストールして頂きたい.
- インストール時にエンディアンは "big endian" で統一して頂きたい. 以下では CReSS の統一エンディアンである big endian を前提に進める.
- NetCDF ラッパーライブラリ gtool5
- 上記リンクからソースをダウンロード・ビルド後, 適切なパスの通ったディレクトリにインストールして頂きたい.
- このライブラリにはもちろん, NetCDF ライブラリが必要である.
変換前データの準備
- グローバル大気観測データアーカイブプロジェクトページの NetCDF化した数値予報GPVデータ から必要なデータをダウンロードしておく. データの種類等は内部の README 参考.
- 地形データも同ページの こちら から適切なデータをダウンロードする. README 参考.
- 地形データの格子点はそれを利用する大気モデルの格子点と一致しているので, README 内のモデル詳細から地形データの格子点情報は取得できる.
変換プログラムの設定
変換プログラム gpv_convert_nc は実行前に設定ファイル gpv_convert_nc.nml を適切に設定する必要がある. 以下, 設定ファイルの内容である.
ただし, 同梱した設定ファイルはこちらで既に MSM 用に設定してある.
&upper ! 3 次元データ設定セクション nx = 241 ! 東西格子点数 ny = 253 ! 南北格子点数 nz = 16 ! 鉛直格子点数 nt = 8 ! 1 ファイルに格納されている初期値の個数 (1 日分 1 ファイルに 3 時間ごとの初期値 = 8 時刻) fname = 'upper.dat' ! 変換する NetCDF ファイルを列挙したリストファイル ctype = 'u', 'v', 'w', 'z', 'temp', 'rh' ! 各変数の名前 (ncdump 等で事前に確認) wflag = 2, 3, 0, 1, 5, 6 ! CReSS 用のデータへの格納順 (固定) off_flag = .true., .true., .false., .false., .true., .true. ! オフセット処理するか typ_flag = 'i', 'i', 'd', 'd', i, 'i' ! 変数の型 axis_name = 'lon', 'lat', 'p', 'time' ! 各座標軸の名前 axis_ord = 0, 0, 4, 0 ! このまま固定 axis_typ = 'r', 'r', 'r', 'r' ! このまま固定 axis_off = .false., .false., .true., .false.! オフセット処理するか / &surface ! 地表面データ設定セクション nxs = 481 ! 東西格子点数 nys = 505 ! 南北格子点数 nts = 24 ! 1 ファイルに格納されている初期値の個数 (1 日分 1 ファイルに 3 時間ごとの初期値 = 8 時刻) sfname = 'surface.dat' ! 変換する NetCDF ファイルを列挙したリストファイル cstype = 'u', 'v', 'sp', 'psea', 'temp', 'rh' ! 各変数の名前 (ncdump 等で事前に確認) wsflag = 2, 3, 4, 0, 5, 6 ! CReSS 用のデータへの格納順 (固定) offs_flag = .true., .true., .true., .true., .true., .true. ! オフセット処理するか typs_flag = 'i', 'i', 'i', 'i', i, 'i' ! 変数の型 axiss_name = 'lon', 'lat' ! 各座標軸の名前 axiss_ord = 0, 0 ! このまま固定 axiss_typ = 'r', 'r' ! このまま固定 axiss_off = .false., .false. ! オフセット処理するか htfile = 'TOPO.MSM_5K_20071121' ! 地形データの名前 ht_ord = 1 ! このまま固定 / &inter yrev = .true. ! If y-direction is reversed, yrev = .true. vert_conv = .true. ! If unit of vert axis is hPa, vert_conv = .true. sflag = .true. ! 地表面データを使うか. interpo_flag = 'surf' ! どちらの格子点にあわせるか. 'surf' = 地表面格子点, 'uppr' = 3 次元データ格子点 /
- 注意点
- 本プログラムは地表面データと 3 次元データを結合して 1 つのファイルにすることができる.
これを行うフラグが sflag である.
これが .true. になっていると, 2 つのデータを結合して 1 つのファイルとして変換する.
- sflag = .true. の場合, 格子点間隔の異なる 2 つのファイルを結合するので, 内挿処理が必要となる.
これらは, convert_nc.nml の "interpo_flag" で設定できる.
本プログラムでは, 2 つの方法で線形内挿を行う.
- 地表面データを 3 次元データに合わせる. "interpo_flag" = 'uppr'
- 3 次元データを地表面データに合わせる. "interpo_flag" = 'surf'
実行前準備
ファイルリストの作成.
3 次元大気データ : UP_0101.nc, UP_0102.nc 2 次元地表面データ : SF_0101.nc, SF_0102.nc
があったとする. このとき, 以下のように, 変換する NetCDF ファイルをリスト化したテキストファイルを用意する.
$ ls UP_*.nc > upper.dat $ ls SF_*.nc > surface.dat
リスト化されたテキストファイル名を convert_nc.nml の指定の場所に記載.
実行方法
$ ./convert_nc < convert_nc.nml
- [注意] : 本プログラムは実行するために 500 MB 程度のメモリを必要とします.