サンプルデータの中には気象庁提供の MGDSST データが3日分含まれています。 あくまでもサンプルとして載せていますので、実際の解析等に用いる場合には 正しいルートにより入手をして下さい(後日別データに差し替える予定)。
MG では非常に面倒ですが最低でも 5つの定義ファイルを用意する必要が あります。いくつかはすでにデフォルトで用意されています。 以下にそれぞれのファイルの概要を紹介します。 test.namelist 以外の書式については次の節にて説明をします。 test.namelist についてはここでは最低限の箇所の変更のみにとどめ、詳しくは 後程紹介します。
このファイルは Fortran の namelist 文の文法で書かれた実行時に 読み込まれる設定ファイルです。非常にいろいろ出来るので煩雑ですが 通常の使用の際には変更する箇所はそれほどありません。また、いずれ 紹介する util ディレクトリにある work_mg.pl スクリプトを用いると 別の設定ファイルが必要にはなりますが、このファイルを自動生成 してくれます。 近い将来ファイル名を mg.namelist に変えようと思います。
数学定数、物理定数などが定義されたファイルで、デフォルトでは 数学定数ファイル(def/constant_math.def)と物理定数ファイル (def/constant_physics.def) が用意されています。
RPN 式で記述された式が定義されたファイル。それぞれ読み込むデータに 合わせてユーザが準備する必要があります。def ディレクトリの中に サンプルがいくつかありますのでそれを参考にするとよいでしょう。
読み込むデータに関する情報が定義されたファイル。GrADS でいう コントロールファイルに対応します。現在のバージョンでは手作業で 作成するのはかなり面倒なので将来的にはもう少し楽に記述できるように したいと思います。
このファイルはデータの鉛直座標系を変換する際に使用します。
このファイルではよく利用する定数を定義します。1行につき1つの定数を 定義します。書式は
<定数名> <値> = [<コメント>]
とします。= 以降は現在のところ何も処理をしていません。もし単位が存在する 定数の場合、単位を記述しておくとよいでしょう。以下に同梱されている def/constant_physics.def を紹介します。
############################################################################### # PHYSICAL CONSTANT ############################################################################### KAPPA 0.2859 = P0 100000. = Pa EPS 0.622 = GR 9.81 = m s-2 RD 287.04 = J K-1 kg-1 RV 461.50 = J K-1 kg-1 GAMMA 6.5e-3 = OMEGA 7.272e-5 = CPD 1004 = XL 2.5e+6 = GAMMAD 9.7709e-3 =
定義する際注意する点は以下になります。
このファイルではデータを読み込む方法やよく利用する計算式を定義します。 定数ファイルと同様に1行につき1つの式を定義します。書式は
<式名> <RPN式> = [<コメント>]
とします。= 以降は現在のところ何も処理をしていません。もし単位が存在する 定数の場合、単位を記述しておくとよいでしょう。以下に同梱されている def/equation_CReSS.def の一部を紹介します。
############################################################################### # DATA READ SECTION # ############################################################################### PBAR READ(BINARY,"PBAR ",T=T) = Pa PP READ(BINARY,"PP ",T=T) = Pa PTBAR READ(BINARY,"PTBAR ",T=T) = K PTP READ(BINARY,"PTP ",T=T) = K ############################################################################### # DEFINITIONS OF EQUATIONS # ############################################################################### # total pressure P [PBAR] [PP] ADD = Pa # total potential temperature PT [PTBAR] [PTP] ADD = K # temperature T [PT] [P] [P0] DIV [KAPPA] POW MUL = K
定数定義ファイルと同様 # で始まる行はコメント行として扱います。 4〜7行目は外部データから読み込む処理を行います。単純バイナリデータを 読み込む場合
<式名> READ(BINARY,"<変数名> ", T=T) = [<コメント>]
の書式で指定します。<変数名> にはデータ定義ファイルの TotalVariables の項目で書かれた変数名を指定します。通常、データ定義ファイルの変数名と 上記の式名を同じにしておくことをお勧めします。13 行目で圧力 P を 16行目で 温位 PT を定義しています。ここではすでに定義された式 PBAR, PP, PTBAR, PTP を利用し、和(ADD) の演算を行っています。例のようにすでに定義されたものは [ ] で囲むことにより利用することが出来ます。なお、定義式はそれを利用する 式の下にあっても構いません。19行目では温度 T を上で定義された P, PT さらには 物理定数ファイルで定義されている P0、KAPPA を利用し、商(DIV)、べき(POW) 、積(MUL) の演算を行います。
以上のようにして、データファイルよりデータを読み込む方法、物理式を 定義して計算を行う方法を RPN 式で定義をすることが出来ます。
このファイルには読み込むデータに関する情報が定義されています。 非常にたくさん設定する項目がありますので、def 以下に用意されている サンプル(model で始まるファイル)と以下の説明を参考に作成して下さい。 将来的にはもう少し簡単にするもしくは作成支援スクリプトを用意したい と思います 現在 CReSS Ver.2.3 用の支援スクリプトはすでに用意されています ( util/make_CReSS_def2.pl ) このスクリプトの使い方は後ほど説明します。
以下にそれぞれの設定項目を紹介した後、サンプルを紹介します。
読み込むデータのタイプを指定します。現在のバージョンでは単純バイナリのみ 対応していますので BINARY で固定 です。
DataType BINARY
データに含まれている時刻を記述します。現在のバージョンでは1行1時刻で すべて記述する必要があります。簡単なサンプルは
TotalDates 3 2092-08-21_00:00:00.0000 2092-08-21_01:00:00.0000 2092-08-21_02:00:00.0000 /TotalDates
となり、書式は
TotalDates <総時刻数> YYYY-MM-DD_hh:mm:ss.ssss YYYY-MM-DD_hh:mm:ss.ssss : /TotalDates
のように TotalDates の後に総時刻数を記述し、次の行から /TotalDates までの 間に1行1時刻を YYYY-MM-DD_hh:mm:ss.ssss の書式で時刻を記述します。 TotalDates と /TotalDates までの行数が総時刻数と異なる場合、または /TotalDates が存在しない場合には異常終了します。
地形データのように時刻情報のないデータの場合、
TotalDates 1 stationary /TotalDates
のように TotalDates を 1 にして、時刻を指定する部分に stationary と記述して下さい。
unite をする前の CReSS の並列計算のようにファイルが領域分割されているような データを読み込む場合、x 方向の分割数を xsep で y 方向の分割数を ysep で指定 します。通常は両方 1 です。
xsep 2 ysep 2
上の例では領域が4分割(x,yそれぞれ2分割)されているデータであることを 表しています。
読み込むファイルを定義します。1行1ファイルで定義を行います。分割ファイルの 時、非常に長い行になります。書式は
<時刻数> <ファイル名> <時刻番号>
<時刻数>とはファイルに含まれている時刻の数を表します。1つのファイルに 複数の時刻情報が記録されている場合にはその時刻数を、1時刻のみの場合には 1 を記述します。<ファイル名> には実行するディレクトリからの相対パス もしくは絶対パスでファイルを指定します。最後の <時刻番号> には読み込んだデータの時刻を特定します。この数字を先に定義した TotalDates と照らし合わせて例えば 3 ならば TotalDates の3番目の時刻として定義します。 時刻が順番に並んでいて、xsep=1, ysep=1 の場合時刻番号を -1 とすると プログラム内で自動に時刻を定義します。以下例です(これは分割ファイルの例です)。
TotalFiles 12 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00000000.pe0000.bin 1 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00000000.pe0001.bin 1 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00000000.pe0002.bin 1 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00000000.pe0003.bin 1 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00003600.pe0000.bin 2 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00003600.pe0001.bin 2 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00003600.pe0002.bin 2 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00003600.pe0003.bin 2 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00007200.pe0000.bin 3 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00007200.pe0001.bin 3 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00007200.pe0002.bin 3 1 /mnt/e3fs1/CReSS/miroc2092ID01191/miroc2092ID1191.dmp00007200.pe0003.bin 3 /TotalFiles
TotalDates の場合と同様に TotalFiles の後に総ファイル数を指定し、 /TotalFiles の間に上の書式で記述をします。TotalFiles の後のファイル数 と /TotalFiles までの行数が異なる場合、または /TotalFiles が存在しない 場合は異常終了となります。
現在処理していませんが
EndianType BIG
として下さい。little エンディアンで記述されている場合は LITTLE としても 構いません。
現在のバージョンでは 4バイト実数、2バイト整数、1バイト整数のデータを 読み込むことが出来ます。ただし内部ではすべて実数として取り扱っています。 読み込むデータのバイト数に応じて 1, 2, 4 いずれかを指定して下さい。
ByteSize 4 (or 2 or 1 )
FALSE と TRUE が指定できますが現在のところ常に FALSE として下さい。 TRUE にしても動作はしますが実行速度は著しく低下します。
VSLICE FALSE
未定義値を指定します。
Undefined <未定義値>
X(またはLON) 方向の総格子数を指定します。
GridSizeX <x 方向の格子数>
x 方向の領域が分割されている場合 (xsep > 1)、総格子数の後に空白区切りで それぞれの領域における最大格子番号を指定します。例えば総格子数 750 の データを 1〜350, 351〜700 に分割をした場合
GridSizeX 750 375 750
のように指定します。
Y(またはLAT) 方向の総格子数を指定します。
GridSizeY <y 方向の格子数>
y 方向の領域が分割されている場合 (ysep > 1)、総格子数の後に空白区切りで それぞれの領域における最大格子番号を指定します。例えば総格子数 800 の データを 1〜400, 401〜800 に分割をした場合
GridSizeY 800 400 800
のように指定します。
Z(鉛直)方向の総格子数を指定します。
GridSizeZ <z 方向の格子数>
X(または LON)方向の格子間隔を指定します。単位は任意です。
DeltaX <x 方向の格子間隔>
Y(または LAT)方向の格子間隔を指定します。単位は任意です。
DeltaY <y 方向の格子間隔>
時間間隔を指定します。単位は任意です。なお、現在のところこの項目の 処理はおこなっておりません。
DeltaT <時間間隔>
地球の半径(m)を指定します。現在のバージョンでは地球を球として 計算を行っています。
EarthRadius <地球の半径(m)>
水平方向の格子の単位を指定します。ただし、現在のところこの部分は 処理をしていません。
HUNIT <水平格子の単位>
鉛直方向の格子の単位を指定します。ただし、現在のところこの部分は 処理をしていません。
VUNIT <鉛直格子の単位>
水平格子の地図投影法を指定します。
| 0 | 緯度/経度系 |
| 1 | ポーラーステレオ図法 |
| 2 | ランベルト図法 |
| 3 | メルカトル図法 |
MPROJ 0 ( or 1 or 2 or 3 )
MPROJ=1,2,3 の時、地図投影の基準緯度を指定します。
TrueLat1 <地図投影基準緯度>
MPROJ=2 の時、2個目の地図投影の基準緯度を指定します。
TrueLat2 <地図投影基準緯度 #2>
MPROJ=1,2,3 の時、地図投影の基準経度を指定します。
TrueLon <地図投影基準経度>
経度/緯度と格子値の対応をとります。( RefLon, RefLat ) の点が 格子 ( RefX, RefY ) に対応します。RefX, RefY は必ずしも格子点上 (整数値)でなくても構いません
RefLat <基準緯度> RefLon <基準経度> RefX <基準 X 格子番号> RefY <基準 Y 格子番号>
鉛直座標(高度)情報を指定します。2通りの指定方法があります。 高度データが等間隔の場合
Level <bottom> <zinc>
のように z=1 における高度と増分を指定します。これに対して高度データが 等間隔ではない場合、
Level 0 0 1000 925 850 700 500 300 /Level
のように Level の後の <zinc> の項目を 0 として、鉛直総格子数分だけ 順に定義し、最後に /Level を追加します。Level と /Level の間に GridSizeZ で定義した数だけデータがない場合、もしくは /Level が存在しない 場合異常終了します。
変数の記録順序と格子数を定義します。簡単な例を以下に示します。
TotalVariables 2 US 750 800 1 1 U 750 800 70 1 /TotalVariables
まず TotalVariables の後に記録されている変数の数を指定します。そして /TotalVariables まで1行1変数で記述します。変数定義の書式は
<変数名> <x方向格子数> <y方向格子数> <z方向格子数> <1回に読む時刻数>
ここで注意するところは最後の時刻に関する項目で、総時刻数を指定するのではなく、 データを読み込む際に読み込む時刻数を指定します。 通常 1 を指定します。
ByteSize = 1 または 2 の場合切片と傾き情報を与えることにより 圧縮データを復元出来ます。書式は
<変数名> <x方向格子数> <y方向格子数> <z方向格子数> <1回に読む時刻数> <傾き> <切片>
となります。
以上の書式に従って書かれた データ定義ファイルのサンプル を見て、書き方が曖昧なところの理解になればと思います。
ここではサンプルデータを用いて MG の実行の方法を簡単に紹介します。まず MG ディレクトリにてサンプルデータを 展開します。
$ tar zxvf mg_sample_01.tar.gz
これにより以下のファイルが展開されます。
+ def/ | +-- equations_mgdsst.def <-- mgdsst 用式定義ファイル | +-- model_info_mgdsst.def <-- mgdsst 用データ定義ファイル | + data/ | +-- mgdsst_20091026.dat <-- SST データ(2009/10/26) | +-- mgdsst_20091027.dat <-- SST データ(2009/10/27) | +-- mgdsst_20091028.dat <-- SST データ(2009/10/28) | +-- test.namelist <-- namelist ファイル(上書き)
test.namelist を編集して mg を実行することになるのですが、test.namelist の書式に関する詳細は次回以降に紹介します。今回変更する場所は
の5つです。以下いくつか実行サンプルを紹介します。
test.namelist の以下の部分を変更して下さい。
OutFile 'out/sst_%d.grd' OutForm 'FUNC(OUTPUT,GRD,H2D,1,1) =' str_test '[SST] =' stime 1 etime 1
簡単にそれぞれを説明するとOutFileでは出力する ファイル名を指定します。%d が存在する場合、この部分が対応する 日付に置き換わります。 OutFormでは出力するデータの書式を定義します。 詳細は次回以降に説明しますが、この例では GMT の GRD 形式の データで出力をする指定をしています。最後に必ず = をつけて下さい。 str_testでは計算式を RPN にて記述します。 上の例ではすでに式定義ファイルで定義されている SST を利用します。 式の最後には必ず = をつけて下さい。 stimeではデータを読み込む時刻を指定します。 これはデータ定義ファイルの TotalDates の順番を表します。すなわち stime = 2 とした場合、TotalDates の2番目の時刻のデータを読み込みます。 etimeも stime と同様データを読み込む時刻を指定します。 stime と etime が異なる場合、stime の時刻から etime の時刻まで str_test 及び OutForm の処理を繰り返します。
以上の編集が終了後
$ ./mg
とすることにより正常終了した場合には out ディレクトリに sst_2009-10-26_12:00:00.0000.grd というファイルが作成されます。 Windows なファイルシステムでは : をファイル名に用いることが出来ませんので %d を指定した時に出力する時刻を変更することを現在検討しています。
結果を GMT を用いて図にしてみましょう(実際は1行です)。
$ grdcontour out/sst_2009-10-26_12:00:00.0000.grd -JX7i/5i -R1/1440/1/720 -C2 -A10 -BWSNE500g100/200g50 > sample1.ps
test.namelist の以下の部分を変更して下さい。サンプル 1 との変更点は OutForm の H2D の次の 1 を 2 にするだけです。
OutFile 'out/sst_%d.grd' OutForm 'FUNC(OUTPUT,GRD,H2D,2,1) =' str_test '[SST] =' stime 1 etime 1
以上の編集が終了後
$ ./mg
とすることにより正常終了した場合には out ディレクトリに sst_2009-10-26_12:00:00.0000.grd というファイルが作成されます (サンプル 1 の結果を上書きします)。
サンプル 1 では X, Y 情報が格子点値でしたが、OutForm の H2D の次の 数字を 2 にすることにより格子系が等緯度/格子系(MPROJ=0)の時 緯度、経度情報にしてくれます。 結果を GMT を用いて図にしてみましょう(最初の2行は1行です)。
$ grdcontour out/sst_2009-10-26_12:00:00.0000.grd -JN180/8i -R0/360/-90/90 -C2 -A10 -K > sample2.ps $ pscoast -JN -R -W1 -G0/200/0 -A500 -BWSNE60g30/30g15 -O >> sample2.ps
test.namelist の以下の部分を変更して下さい。
OutFile 'out/sst_%d.grd' OutForm 'FUNC(OUTPUT,GRD,H2D,2,1) =' str_test '[SST] 273.15 SUB =' stime 1 etime 3
以上の編集が終了後
$ ./mg
とすることにより正常終了した場合には out ディレクトリに sst_2009-10-26_12:00:00.0000.grd, sst_2009-10-27_12:00:00.0000.grd, sst_2009-10-28_12:00:00.0000.grd, というファイルが作成されます (最初のファイルは上書きします)。
この例では str_test にて読み込んだ SST データから 273.15 を引く演算を 行い、SST の単位を K から摂氏に変換しています。また、etime を 3 にすることにより 3 時刻のデータを一気に変換します。 結果を GMT を用いて図にしてみましょう(最初の2行は1行です)。
$ grdcontour out/sst_2009-10-26_12:00:00.0000.grd -JN180/8i -R0/360/-90/90 -C2 -A10 -K > sample3.ps $ pscoast -JN -R -W1 -G0/200/0 -A500 -BWSNE60g30/30g15 -O >> sample3.ps
test.namelist の以下の部分を変更して下さい。
OutFile 'out/sst_diff.grd' OutForm 'FUNC(OUTPUT,GRD,H2D,2,1) =' str_test '[SSTC] [SSTC:T-2] SUB FUNC(AGRID,2,10,1) =' stime 3 etime 3
以上の編集が終了後
$ ./mg
とすることにより正常終了した場合には out ディレクトリに sst_diff.grd というファイルが作成されます。
サンプル 3 で行った摂氏への単位変換はすでに eqution_mgdsst.def で SSTC という式で定義されています。 stime, etime 共に 3 なので3番目の時刻(2009/10/28) のデータを 読み込みますが、2番目の [SSTC:T-2] のように : の後に T-2 とすることにより2時刻前のデータを読むということが出来ます。 従って上の例では、まず3時刻目の SSTC を読み込み、次に1時刻目の SSTC を読み込み、それらを引くことにより2日間の SST 偏差を 計算します。 最後の FUNC(AGRID,2,10,1) にて平滑化操作を行っていますが 詳しくは次回以降に解説します。 結果を GMT を用いて図にしてみましょう。
$ grdcontour out/sst_diff.grd -JN180/8i -R0/360/-90/90 -C0.2 -A1 -K > sample4.ps $ pscoast -JN -R -W1 -G0/200/0 -A500 -BWSNE60g30/30g15 -O >> sample4.ps