grads形式のデータからデータの一部分を取り出す

gradsで、あるシミュレーションデータとあるシミュレーションデータを引き算したい。
でも、格子数が違うからできない！
みたいな経験はないでしょうか。

私の場合は、何かの手違い（初期値作成時だと思われますが）で
格子点が一つだけずれてることがありました。

そんなときに、outputデータからその一部を取り出すことができます。
以下はgrads2dから、viqciのデータのみを取り出す例。
（gradsスクリプトに書いてね）

'open grads2d.ctl'　　　　　　　ctlファイルを開く。
'set fwrite new.grd'　　　　　　new.grdはoutputファイル名
'set gxout fwrite'　　　　　　　gxoutで設定
'set x 1 395'　　　　　　　　　　x方向の格子を1～395まで設定
'set y 1 402'　　　　　　　　　　y方向の格子
'set z 1'　　　　　　　　　　　　　鉛直層数
'set t 1 90'　　　　　　　　　　　t　size
'd viqci'　　　　　　　　　　　　 出力する要素
'disable fwrite'　　　　　　　　出力を閉じる


また、鉛直層が一層だけじゃない場合は、ループを使いましょう。
t sizeについてもループが使えます。
dの前後に↓を追加しましょう。

t=1
tmax=90
while(t<=tmax)
'set t 't''

z=1
zmax=36
while(z<=zmax)
'set z 'z''
 
'd ver'

z=z+1
endwhile

t=t+1
endwhile

また、開くファイルはnetcdfでも可なので、netcdfから取り出したい時にも便利。↓
 
'sdfopen netcdf.nc'

MANALのデータを見よう！

MANALのデータを見る手順。
講座の方は私のディレクトリを参考に。

-------2008.2.3 revised--------

①cloud2/{my name}/MANAL/work/にあるgrib2grads.sh
　　見たい日にちのファイル名に変更する。MANALのデータも同じところにおいておく。
　　シェルの中身はwgribがたくさん書いてあるだけ。
　　shで実行。「.V」「.U」などのファイルが途中に生成される。
　　これらファイルをくっつけて、gradsで読める形に。
　　.drってファイルが出力される。

②あとはctlファイルを作って、gradsで読むだけ。
　　ctlファイルは「grib2ctl.pl」ってコマンドで作る。→できたctlファイルは、いらない部分を消しましょう。
　　あとは、MANALはランベルト、gradsのデフォルトはメルカトルなので、変更が必要な点に注意。

＊現時点ではt size（時間）=1のみ実行できることを確認。

地形を変形する

計算に用いる地形を変形する方法

まずは計算する分の初期値・境界値を作成しておく。
（mkinitとかn2nで作成）

シェルは
＠ Tools/preprocess/topo_edit/script/topo_transform.sh
・output directoryの名前を設定
・この中に、初期値・境界値作成でつくったorg～の地形ファイルにリンクを貼る。
　（original topo fileってとこ）
・nammapのところで変形させる部分の中心緯度を設定
・出力されるファイルは
　org_new（このまま計算に利用するファイル）
　check_zs.grd（地形確認用ファイル。grads形式。zsが元の地形、zsnが変形させた地形）
　ctlファイルはcheck_zs.ctlを参考に。

プログラムは
＠ Tools/preprocess/topo_edit/src/topo_transform.f90
・プログラムを開いて"change"と検索して引っかかったところを変更。
・if ( rad <= 20 ) then zsn(i,j) = 0.2*zs(i,j)
　( rad <= 20 )は中心からの格子点数を表している。
　この範囲の地形を0.2倍している。radの範囲とzsにかける値を適宜変更。
・一部分だけ0.2倍したりすると、断崖絶壁ができたりする恐れがあるため、
　周囲に0.4,0,6…とかけていき、ゆるやかに地形を削る。
＊中心からの格子点数を指定するので、地形は円形に削れます。

雲水などの生成量を出力

NHMで雲水や雲氷の生成量を計算・GrADSで見たい場合

まずは、計算する際のパラメータカード。
NAMKDDで生成量を計算するように設定することが必要。
雲水生成量：pqcw、雲氷生成量：pqci　などです。
namelist一覧で確認のこと。

計算後、pqcwやpqciなどを出力できるように、
nhm2gradsのプログラムを一部追加しました。
＠ Tools/postprocess/nhm2grads/src
プログラム→nhm2grads_new.f90
　　　　　　　　makval_new.f90
＠ Tools/postprocess/nhm2grads
シェル→n2gr_252dx4_new.sh

他のプログラムについては変更しなくて大丈夫です。

緯度・経度によらない鉛直断面を描く

gradsは緯度・経度に沿ってしか断面を見れないので、
斜めに切って見たい時などはこれを使う。


for MANAL
① MANAL/grads/vertical/make_vcs_manal.sh
　　このシェルの必要な部分を書き換える。
　　・ 絵を描きたい日時のMANALファイル名

②MANAL/grads/vertical/src/readmanal.f90
　　readmanal.f90中の必要部分の書き換え
　　・ 断面の中心緯度・経度
　　・ 中心からの距離（単位は"度"）
　　・ 切断する角度
　　＊はじめに定義されているnxiは断面中の格子数である。

③sh make_vcs_manal.sh

④vertical/work_manal中に"vcs.grd"が作成される。
　　これをgradsで開いて見ればよい。


for NHM output
①vertical/fig_make_vcs.sh
　このシェルの必要部分を書き換える。
　・ソースコードやデータの場所など
　・使うプログラム：
　　　252dx4仕様→read_grads3d_252dx4_new.f90
　　　400dx1仕様→read_grads3d_400dx1.f90
　　　↑これらのプログラムには、スタート地点と終了地点の緯度・経度を計算して
　　　　 出力する部分を追加してある。出力する要素もオリジナルより増やしてある。
　　　＊元々のプログラムはread_grads3d.f90。

②/vertical/src/read_grads3d.f90など
　read_grads3d.f90などの必要部分書き換え
　・断面の中心緯度
　・中心からの距離
　・切断する角度
＊nx,ny：水平格子数、nz：鉛直格子数、nt：NHMでの出力ステップ数
　 nv：3-dimensional dataで出力する変数の個数、nv2：2-dimensional dataの同
　 horg：鉛直層の数によって変更しなければならない。fcstのlogに出ているので、
　　　　 それを見るのが手っ取り早い。
　 nxi：断面中の格子数（直径÷格子間隔より大きくする）
　 ＊半径：dis[km]、格子間隔：dr[km]で設定
　 変数の追加・削除→"! read 3-dimensional data"以下にあるループの中身を
　 追加したり消去したりすればよい。

③sh fig_make_vcs.sh
④できたvcs.grdをgradsで見る。

＊gradsで見るときに、始点と終点の位置（緯度・経度）が分からないと不便だったので、
　 簡単な計算式を追加しました。
　 read_grads3d_400dx1.f90などを参照。「add by takuya-s」って書いてあるところです。
　 lon1,lon2,lat1,lat2が始点と終点の緯度・経度です。
　 計算式中の数字（400dx1では115）は、モデルの計算結果における、一度分の格子数
　 （115÷格子間隔）の値を入れればオッケーです。
　 ちなみにlogに出てくるslonとかは、始点の格子点座標なので緯度・経度とは違います。
　 注意ですね。