玉川研究室で、CReSS を使うには。その3 気象庁GPVデータ篇
海表面温度追加 (2007.11) クリックするとジャンプします
註: この作業は、計算を実行する並列計算機ではなく、
たとえば、192.168.1.12 など。他のPCで行なって下さい。
並列計算機は、あくまでも計算用として使いたいので、GMT などのソフトを 整備していません。

気象データ MANAL

気象庁 からは、 客観解析データが気象業務支援センターを通じて、 販売されている。 ここでは、この客観解析データの内、もっとも解像度の高い、 MANAL をもとに、
CReSS で使用する初期値データ、境界値データを作成する手順を 説明する。

まず、MANAL に限らずこの種のデータは GRIB と言う形式で書かれている。
たとえば、 ここ (京大の坂井亮彦さんのページ) や、 NCEP, NOAA, USA に情報がある。

これを読むには、wgrib というツールを利用するのが簡単である。
wgrib は、 WGRIB in NCEP, NOAA, USA から入手できる。 Cソースコード 。
これを入手し、cc -o wgrib wgrib.c とすると、コンパイルでき、 例えば ./wgrib とすると起動するというように、実行可能となる。

このコマンドの使い方は、多くのページにのっているので、省略するが、
たとえば、 筑波大の hayasaki さんのページ などに詳しい。

とりあえず MANAL のデータを眺めたい人には、研究室の学生用 サンプルのページにある
MANAL を読んでテキスト出力する tcsh スクリプトのサンプル が 役に立つかもしれない。

この wgrib を使うと、データの入っている項目の情報を取り出したり、
選んだデータを、テキスト形式やバイナリ形式で取り出したりできる。
これを利用して、CReSS の前処理プログラム gridata が必要とする バイナリ形式(4バイト実数)で、ベタ書きのデータを作成する。

ただし、CReSS では、最初に地表面を示すデータを入れてやることによって、 地中への内挿処理を避ける事ができるようなので、この機能も使うことにする。
ここで、ちょっとした問題が起こる。それは、MANAL データには、 地上気圧はなくて、海面更正した地上気圧が入っていることである。
これは、気温の低減率 0.5K/100m を使った静力学の式から換算ができるので、 その為のソフトも作った。 (psfc.c)

他にも、相対湿度のように、高度によっては存在しない項目や、 気圧のように、一定値なのでない項目もあるので、
そのような時に、 ただ、決まった値を必要個数出力する他愛のないプログラムも作った。

さらに、相対湿度については、CReSS 用に GPV データを変換する gridata 内の
gridrv 内の chkgpv, chkmxn での検査 (0% - 100% 内)に引っかかる (単なる２進変換誤差)
そこで、そのチェックのために check_under100.c も作った。

また、MANAL データは、MANAL_2005081200.grb のよう( 2005年8月12日00Z ) な ファイル名をしているので、
これを加工して、data.gpv200508120000.bin というファイル名にして出力するようにしてある。

これらの処理を一括して行う以下のスクリプトと、そこで使われている C 言語のプログラムは、以下のようなものである。

ダウンロードの後、適当なディレクトリに置いて、 
cc -o wgrib wgrib.c
cc -o pudding pudding.c
cc -o psfc psfc.c -lm
cc -o check_hgt check_hgt.c
cc -o check_under100 check_under100.c
として、コンパイルしておき、 必要に応じて MANAL_CRESS.tcsh の最初の方にある 
set WGRIB = './wgrib'
set PUDDING = './pudding'
set PSFCPROGRAM  = './psfc'
set CHECKHGT = './check_hgt'
set CHK100 = "./check_under100"
これらの行を絶対パスに書き換える(どちらかというとお勧め)。 更に、 
set TOPO = "/home/tama/JMA_DATA/MANAL_200508/topo/MANAL_TOPO.grb"
の行の、右辺を実際に自分が、MANAL の地形データ MANAL_TOPO.grb を 置いた場所に変更する必要がある。

そののち、 
chmod +x MANAL_CRESS.tcsh
として、MANAL_CRESS.tcsh を実行可能にして、 ./MANAL_CRESS.tcsh MANAL_2005081200.grb などと、MANAL のファイル名を後に 付けて実行すると、
TESTDIR の下に、さまざまな気圧での各要素の２次元バイナリデータが出力され、それらをまとめて、
data.gpv200508120000.bin のような ファイル名にして出力される。このデータが、CReSS 付属のgridata で 処理できるはずである。
もちろん、この際に、気象業務支援センターから購入した CDROM のデータを
必要なだけ、今作業している計算機 (Linux のサーバの方) に転送しておく必要がある。 (地形データ も)

追加、地表面なしバージョン

上と同様な処理ではあるが、地表面の値を作らないで、素直に 1000hPa からデータを
並べるバージョンも作成した。もしかしたら解像度の違いなどで地形がつながっていない場合は
こちらの方が良かったりするかも知れない。
層数が、上の例より１減って 20面になることにも注意が必要である。

追加、終り (2006.10)

もっと下層を減らして、900hPa 以上だけのバージョンも置いておく sfc 1000 950 925 と減って 17層だけになることに注意!

TESTDIR 以下にある各２次元データのクイックルックのために、
CheckBinData.tcsh (シェルスクリプト) と、ダンププログラム dump2dbin.c も用意した。

cc -o dump2dbin dump2dbin.c
chmod +x CheckBinData.tcsh
とした後、 
./CheckBinData.tcsh TESTDIR/RH_900
のように、各２次元データのファイル名を付けて実行すると、 PH_900.ps のような名前で、
ポストスクリプトの画像ができ、 カラーのマップとして、データを見る事ができる。

さて、話を戻して、 MANAL_CRESS.tcsh で出力されたデータは、
高度 HGT
風速U成分 UGRD
風速V成分 VGRD
気圧 P
気温 TMP
相対湿度 RH
である。(単位は m/s Pa K %)
また、高度は、地上と 1000 950 900 850 800 700 600 500 400 300 250 200 150 100 70 50 30 20 10 の計21面である。
地表面なしバージョンだと、20面である。

今回使った MANAL についている説明によると、MANAL はランベルト正角円錐座標系でかかれており、投影基準の緯度は、30°N と 60°N、投影基準経度は 140°E、
投影基準点での格子間隔は 10km 、領域の西から 245番め、北から205番目の格子点が (140°E, 30°N) の点であり、格子数は X軸方向 361 Y軸方向289である。
風は、座標系での格子成分であるが、CReSS の場合はこのままで良いはずである。 図に書いたり、自分で解析したりする際には注意が必要である。

海表面温度データ

CReSS の計算領域に、海がある場合、海表面温度を与える必要もでてくる。
本WWWの例では、定数を与える簡単な方法を示してあるが、領域が広い場合、ちゃんと 海表面温度データを与える方が良いはずである。

海表面温度データとしては、例えば アメリカの海洋気象局 NOAA の NCEP が作成して いるものがあり、 Real-time, global, sea surface temperature analyses から得ることができる。
ここには、RTG_SST_HR と RTG_SST の２種のデータがあり、 それぞれの特徴は、この頁によると以下のようなものである ( 上記 RTG の頁より引用 2007 Nov.)
File Name
 RTG_SST_HR
 RTG_SST
Horizontal Resolution
(Lat/Lon Grid)
 0.083 degree
 0.500 degree
In-Situ Data
 Fixed buoys, drifting buoys, and ships
 Fixed buoys, drifting buoys, and ships
Satellite Data
 NOAA 17 & NOAA 18 AVHRR
 NOAA 17 AVHRR
Satellite Processing
 JCSDA Physical Retrievals
 Navy Retrievals
Implemented
 September 27, 2005
 January 30, 2001
Status
 Operational
 Operational

データは、それぞれ からたどっていくと FTP でダウンロードできる。

これらは、GRIB形式で書かれているので、気象庁のデータと同様の方法で データ変換が可能である。例えば、0.5度解像度の rtg_sst_grb_0.5.20060812 を ダウンロードしたとすると、上記 wgrib を使用して、 

./wgrib -V rtg_sst_grb_0.5.20060812
とすれば、 
rec 1:0:date 2006081200 TMP kpds5=11 kpds6=1 kpds7=0 levels=(0,0) grid=235 sfc anl:
  TMP=Temp. [K]
  timerange 1 P1 0 P2 0 TimeU 2  nx 720 ny 360 GDS grid 0 num_in_ave 0 missing 0
  center 7 subcenter 4 process 44 Table 3 scan: WE:NS winds(N/S) 
  latlon: lat  89.750000 to -89.750000 by 0.500000  nxny 259200
          long 0.250000 to -0.250000 by 0.500000, (720 x 360) scan 0 mode 128 bdsgrid 1
  min/max data 271.34 307.14  num bits 12  BDS_Ref 27134  DecScale 2 BinScale 0
のように内容に関する記述が得られ、CReSS で使用するためには、 
./wgrib rtg_sst_grb_0.5.20060812 | grep TMP | ./wgrib rtg_sst_grb_0.5.20060812 -i -bin -nh -o SST_05_20060812.bin
のようにすれば、SST_05_20060812.bin に、CReSS で読める形式のデータが出力される。 格子数は、nx 720 ny 360 で、(1,1)は、多分 (0.25, -89.75) ではないかと思う。
(Nov. 2007 追加)

安直に確かめるためには、 dump_bin.f90 が使えると思う。
今回は、Fortran95 で書いてみた。使い方は以下のとおり。
(ここで、Fortran95 のコンパイラを gfortran としておく、各自の環境に合わせて 使って欲しい) 

gfortran dump_bin.f90
./a.out SST_05_20060812.bin 720 360 SST_05_20060812.text
これで、i j SST とデータの並んだ安直なテキストファイルができるので GMTでコンターを引く 
 makecpt -Crainbow -T270/330/2 > CPTfile
 pscontour SST_05_20060812.text -CCPTfile -JX10 -R1/720/1/360 -W1 -Ba20 > fig.ps
あるいは 
 makecpt -Crainbow -T270/330/2 > CPTfile
 surface SST_05_20060812.text -GSST_05_20060812.grd -I1 -R1/720/1/360
 grdimage -CCPTfile  -R1/720/1/360 -Ba20  -JX10 SST_05_20060812.grd > GSST_05_20060812.ps
ただし、陸地の上にも適当に補間されたデータが並んでいるので、なかなか どっち向きか判別するのは難しいかもしれない。

追加

研究室内のみアクセス可ですが、192.168.1.12:/work/pub/SST/ に、0.5度SSTデータを ダウンロードしておきました。(2007 Dec.)