玉川研究室で、CReSS を使うには。6
計算後の処理

solver.exe による計算が終了し、unite.exe による各PEのデータの
結合処理が終ったら、次にその計算結果を読みだして処理しなければならない。
CReSS の出力データは、指定にもよるが、多くの場合バイナリ形式を選ぶことになるし、
各出力時刻のさまざまなデータを並べたものになっているので、同様のフォーマットの
データの処理が得意な GraDS ででもなければいきなり図を描く事は難しい。
また、望む出力を得るためにはデータを読んで加工する必要もでてくる。

そこで、バイナリ形式で、計算している座標系 (dmpcmp=1) での出力ファイルを
極力汎用性を重んじて切り出す処理を行なうソフトを作成した。
readcress200609.c である。
U = Ubase + U , V = Vbase + V と処理していたが、間違っていた ような気がする、
素直に U = U, V = V とするバージョンを作成した
readcress200802.c (2008 Feb.) CReSS の出力結果から、さまざまな指定をした切り出しが可能である。

使い方は、まず計算時の LOG を見て、何が出力されているのかを確認し
その記述ファイルを Descrition.txt のように作成する。
書き方は、# から始まる行がコメントで無視されることを踏まえて、実際に読んでもらえれば分かると思うが、
最初にモデルの設定情報を書き、その後に、出力変数を出力される順に
「２次元データなのか、３次元データなのか、」
「そのデータを何と呼ぶか」
「データに関するコメント」
の順で書き並べるだけである。サンプルは、雲の計算を行なった状態での全ての出力である。
この記述は、geography のファイルと dmp ファイルと両方について行なう。

コンパイルは簡単で、ただ

cc readcress200802.c -o readcress200802
とするだけである (作成は FreeBSD 6.1 上で行なった)。
実行は、./readcress200802 とすれば出て来る以下のメッセージにあるように
**** ERROR MESSAGE ****
Program for reading CReSS data
USAGE:
for 3d valiables
./readcress dmpdata geodata descdata 3 x  xnum item
./readcress dmpdata geodata descdata 3 xy xnum ynum item
./readcress dmpdata geodata descdata 3 xyz xnum ynum znum item
for 2d valiables
./readcress dmpdata geodata descdata 2 item
./readcress dmpdata geodata descdata 2 x xnum item
    x canbe y or z, also y , z
    xnum ... : from 2 to  xdim-2
    item: in itemlist
３次元変数、２次元変数のどちらに対しても、指定した格子番号での、
断面図、あるいはプロファイル、点の値を標準出力に出力する。
パターンが xyz の時には３次元データの３つの格子を指定し、その点での値を出力し、
パターンが、xy の時には、Xの格子番号、Yの格子番号の順で指定して、Z 方向の分布を 出力する。
この時、パターンを yx と指定したら、Yの格子番号、Xの格子番号 の順として 処理する。
xz, zx, yz, zy も同様である。
出力の際に、U V P PT のように、ベースと変動に分かれている物に対しては、
両者の和をとり実際の物理値として出力するようになっている。
(ただし、変数の名前を利用しているので、「記述ファイル」の U V P Pt の名前は 変更してはならない。)
パターンや、変数名は、大文字小文字を区別しないで処理する。
出力は、いつでも
経度 緯度 高度 値
の順に並んでいる。

(例1)

./readcress ../CReSS2.1m/takayama2006a.dmp000450.united.bin ../CReSS2.1m/takayama2006a.geography.united.bin Descrition.txt 3 z 2 U > outfile
とすると、
ヒストリデータ ../CReSS2.1m/takayama2006a.dmp000450.united.bin
地形データ ../CReSS2.1m/takayama2006a.geography.united.bin
記述データ Descrition.txt
から、３次元データである U を 格子番号 z=2 つまり最下層で水平に切り取った 断面のデータを outfile に出力する。

これを、例えば awk などで必要なカラムだけ取り出して、GMT などで図化すれば良い。
(例1)の場合は、xy 方向の分布なので (本当は、地形に沿った座標で z も変化しているので美しくない例なのだが)
awk '{print $1, $2, $4}' < outfile > xyU
とすると、xyU に、経度 緯度 U の形で並んだデータセットができるので、
GMT のsurface makecpt grdimage psscale などを使って、グラフを書く事ができる。

(例2)

./readcress ../CReSS2.1m/takayama2006a.dmp000450.united.bin ../CReSS2.1m/takayama2006a.geography.united.bin Descrition.txt 3 xy 10 10 PT > outfile
とすると、格子番号で x=10 y=10 の場所での温位の鉛直分布が得られる。
先程と同様に 
awk '{print $3, $4}' < outfile > zPT
で zPT に、高度 温位 のデータが得られるので、xgraph や gnuplot あるいは psxy などで図化できる。
gnuplot なら、plot "outfile" using 3,4 with lines でいきなり書くこともできる。

(例3)
時系列を書くには、３次元データなら３つの格子番号を、２次元データなら２つを指定し、
点の値を出力するようにして、それを foreach のループで回すと良い。tcsh を使っている。
もし bash なら do loop があるはずだが、tcsh と打って下を使うこともできる。
rm outfile
foreach t (../CReSS2.1m/takayama2006a.dmp*united.bin)
set s=`echo $t:t | sed 's/takayama2006a.dmp//' | sed 's/.united.bin//'`
printf "%s " $s >> outfile
./readcress $t ../CReSS2.1m/takayama2006a.geography.united.bin Descrition.txt 3 xyz 10 10 2 P >> outfile
end
時刻の情報はヒストリデータ(dmpファイル)のファイル名にあるので、
set と printf の行で作り出しているところが、ちょっとテクニカルである。
その結果、時刻 経度 緯度 高度 気圧 と言ったデータが出て来るので、
gnuplot で plot "outfile" using 1:4 とする、などして図化すれば良い。

実際に gnuplot を使う場合には、タイトルやラベルの設定、 出力する画像データの形式、ファイル名などを指定する必要がある。
(例3)の場合、gnuplot 中で
set title "time sequnce of pressure"
set xlabel "time (sec)"
set ylabel "Pressure (hPa)"
plot "outfile" using 1:($4/100) notitle wi li
などとして、画面上で図を整えた後、 
set term post eps color "Helvetica" 18
set output "output.eps"
replot
save "outfile.gnuplot"
として、出力形式、ファイルを指定し、描画した後、gnuplot のスクリプトを 保存しておくと良い。
保存したスクリプトは、term と output の指定の場所がコメントアウトされているので
そこを変更するなど、必要な編集をした後、gnuplot outfile.gnuplot と言う風に
非対話的に使用することができ便利である。
蛇足であるが、ちょっと線が細めなので、vi output.eps して、 gnulinewidth を 5 から 15 くらいするのが、筆者の好みである。

風速のベクトル図のために

風速のベクトル図を書くのにGMTでは、「経度 緯度 風速 風向」の順に ならんだデータが必要である。
(風向は、ベクトルの方向で書くので、気象学の常識とは 違うので注意!)。
このようなデータの変換には、上の readcress200802 と awk を組み合わせた下のようなスクリプトが役立つ。
このスクリプトを 書いて、変数を状況に合わせて設定すれば、「経度 緯度 風向 風速」 で 並んだデータセットができるので、
これを使ってGMT の psxy の -Sv オプションを 使って表示できるはずである。(2008.1) 
#!/bin/tcsh 
set PROG = readcress200802
set DATADIR = '../CReSS2.2m'
set GEO  = takayama2007a.geography.united.bin
set DMP  = takayama2007a.dmp031800.united.bin
set DESC = Descrition.txt
set OUT  = UDIRout


set DMP  = takayama2007a.dmp031800.united.bin

./$PROG $DATADIR/$DMP $DATADIR/$GEO $DESC 3 z 2 U > Uout
./$PROG $DATADIR/$DMP $DATADIR/$GEO $DESC 3 z 2 V > Vout

cat -n Uout | sort -b >  Uout1
cat -n Vout | sort -b >  Vout1

join Uout1 Vout1 | sort -n  > UVout1
# awk '{print $2, $3, $5, $9}' < UVout1 > UVout2
awk '{print $2, $3,(180.0/3.1415926535) *atan2($9, $5), sqrt($5*$5+$9*$9)}' < UVout1 > $OUT

# head Uout1 Vout1 UVout1 UVout2 UDIRout
rm Uout Vout Uout1 Vout1 UVout1

次に地上気温の処理について

#!/bin/tcsh 
set PROG = readcress200802
set DATADIR = '../CReSS2.2m'
set GEO  = takayama2007a.geography.united.bin
set DMP  = takayama2007a.dmp031800.united.bin
set DESC = Descrition.txt
set OUT  = Tsurface
set POUT = P15out
set PTOUT = Pt15out
set PPTOUT = PPtOUT


set DMP  = takayama2007a.dmp031800.united.bin

./$PROG $DATADIR/$DMP $DATADIR/$GEO $DESC 2 P15 > $POUT
./$PROG $DATADIR/$DMP $DATADIR/$GEO $DESC 2 Pt15 > $PTOUT

cat -n $POUT | sort -b >  ${POUT}.1
cat -n $PTOUT | sort -b >  ${PTOUT}.1

join ${POUT}.1 ${PTOUT}.1 | sort -n  > $PPTOUT
awk '{T=$9 / ((100000.0/$5)**(287.05/1005.0)) ;  print $2, $3, (T-273.15)}' < $PPTOUT > $OUT

# head  $POUT $PTOUT ${POUT}.1 ${PTOUT}.1 $PPTOUT
rm $POUT $PTOUT ${POUT}.1 ${PTOUT}.1 $PPTOUT

水平断面

CReSS では、地形に沿った座標系で計算が行なわれ、今回は地形に沿った座標系での 出力であるので、
よくある、高度 XXXX m での断面図とか、XXXX hPa での断面図を作るのに、１仕事必要である。
そこで、計算結果を線形内挿して、水平方向の断面データを出力するプログラムも作成した。
Horizontal_Cross.c である。
readcress 同様、これも U V の処理がおかしかった。
 修正版を Horizontal_Cross2.c におく(2008.2)

これも、コンパイルは簡単で、ただ cc -o Horizontal_Cross Horizontal_Cross.c とするだけである。
使い方は、やはり同様に、何もオプションを付けずに起動すると下のように出て来る。

./Horizontal_Cross item number unit descritionfile dmpfile geofile
           item:   shuld be 3d variables in dmpfile
           number: height or pressure
           unit:   m or hPa
item の指定は、上の readcress200802 と同じであるが、dmpfile にある３次元データしか指定できない。
number unit は、例えば 500 m とすると高度500m での断面になり、 700 hPa とすると、700hPa での断面になる。
descritionfile dmpfile geofile は、readcress200802 と同じで、それぞれ 記述ファイル、ヒストリファイル、geography ファイルである。
内挿はもちろん鉛直方向に行なわれ、線形で行なわれています。
精密な図を作りたい時などには、検討が必要かも知れません。
また、内挿でなく外挿になってしまう場合、つまり計算領域の上端より上、
あるいは地表面より下の場合には、欠測を示す -9999 を出力するようにしてあります。
(2006.9 玉川)

Nesting のために

CReSS の計算結果を使って、次の CReSS による計算の境界条件や初期値を作成する (ネスティングなどのため)
ためのデータ変換のソフトを として、作成した。(2008.1 修正) UやVを Ubase + U Vbase+V として 処理していたのは、間違いで U, V そのままで良いのではないかという 気がした (CReSS 2.2m) 。そこで、U V の処理を替えたバージョンを 作成した (2008.2) コンパイルの方法は 
cc -o CReSStoCReSS2 CReSStoCReSS2.c
実行方法は、ただ ./CReSStoCReSS2 とするとメッセージにでるように 
./CReSStoCReSS2 ヒストリファイル名  出力ファイル名 descritionfile
である。descritionfile は、上記の読み込みソフトと同じものでよい。
main 関数の冒頭部分に、読み出すデータが列記されているので、適当に変更すると
違う状況でも使用できると思う。ここでは、雲のあるシミュレーションで、計算した座標系で 全ての出力
で、ヒストリファイルが出力されており、そこから使えるもの全てを取り出す
という想定で書かれている。
出力されるデータは、 ZPH U V P PT W QV QC QR QI QS QG である(はず)。

これを、CReSS に与えるには、ヒストリファイル作成時に CReSS の設定ファイルで
指定していた xdim ydim zdim から 3 引いた値を、次の CReSS に与える GPV データの
格子数として指定する必要があります。もちろん、地図の座標系などもヒストリファイル作成時
の設定を、次のGPVデータの設定として与える必要がありますし、温位、混合比でデータが与えられている
という設定 (pm) に変更する必要もあります。(2006.10 玉川)