津波の伝播計算を自分がやりたいようにやるには,結局一番大変なのは海底地形データを準備することのようです.一例として,遠地地震(チリ地震等)による津波を計算するため,グローバルな地形データを用意する手順をまとめておきます.
まずは,地形データを入手する必要があります.陸域であれば,SRTMやG-DEMが使えますが,海では,GEBCO (General Bathymetric Chart of the Ocean) を使うことが出来ます.ユーザ登録をする必要はありますが,フリーで使えます(詳細なライセンスはサイトを確認してください.以下でダウンロードした zip にも入っています).
ダウンロードすべきは,"GEBCO_2014 Grid (30 arc-second interval) の User-defined area or global grid | 2D netCDF" です.これにチェックをして,"Add data to basket" を押して,"View basket" に進んで待ちます.分かりやすいインタフェースとは言い難いですが,リクエストしたデータが用意できるまで少し待つ必要があります.データが準備出来たら,ダウンロードすると zip ファイルを取得できます.解凍すると GEBCO_2014_2D.nc が入っていて,これをJAGURSで使える形に変換します.
注意点は以下の通り
という感じになります.grdcut で領域の切り出しをしています.東端を300とすることで東経の連続データとして変換します.grdsample でも-Rを指定できますが,120/300 のように指定をしても,意図したようには機能してくれませんでした.(必要に応じて,30秒グリッドのオリジナルデータを grdsample で解像度を変更しています.ここでは, -I2m で,2分角に変更しています.)最後に,正負を逆転(grdmath の NEG)しつつ, =cf でフォーマットをcf に変換しています.
※ 領域を Baba et al. (2017) に合わせました.
まずは,地形データを入手する必要があります.陸域であれば,SRTMやG-DEMが使えますが,海では,GEBCO (General Bathymetric Chart of the Ocean) を使うことが出来ます.ユーザ登録をする必要はありますが,フリーで使えます(詳細なライセンスはサイトを確認してください.以下でダウンロードした zip にも入っています).
ダウンロードすべきは,"GEBCO_2014 Grid (30 arc-second interval) の User-defined area or global grid | 2D netCDF" です.これにチェックをして,"Add data to basket" を押して,"View basket" に進んで待ちます.分かりやすいインタフェースとは言い難いですが,リクエストしたデータが用意できるまで少し待つ必要があります.データが準備出来たら,ダウンロードすると zip ファイルを取得できます.解凍すると GEBCO_2014_2D.nc が入っていて,これをJAGURSで使える形に変換します.
注意点は以下の通り
- 極が入っていると計算できないので,必要な領域(太平洋を中心として日本とチリが含まれる範囲等)の切り出します.
- この時,日付変更線(東経180度と西経180度)をまたぐので,東経で190度等のように東経で0度から360度までの連続的な座標に変換することも同時に必要です.
- JAGURSは水深が正なので,標高データであるGEBCOのデータを正負反転(-1をかける)させます.
- NetCDFの標準(?)である nf ではなく,JAGURSは cf を前提にで作られているので変換します.
- cfとnfは,grdreformat のオプションで確認することが出来ました.cfは,GMT3の古いフォーマットのfloatで,nfがGMT4のデフォルトのフォーマットのようです.以下は grdreformatの出力の抜粋.
| ID | DESCRIPTION - GMT 3 netCDF legacy formats |
| cf | GMT netCDF format (float) (deprecated) |
----------------------------------------------------------------
| ID | DESCRIPTION - GMT 4 netCDF standard |
| nf | GMT netCDF format (float) (COARDS-compliant) [DEFAULT] |
いくつかやり方はあるはずですが,gmt4を使うなら,
% grdcut GEBCO_2014_2D.nc -Gtmp.grd -R120/300/-60/60
% grdsample tmp.grd -Gtmp2.grd -I2m
% grdmath tmp2.grd NEG = gebco2min.grd=cf
という感じになります.grdcut で領域の切り出しをしています.東端を300とすることで東経の連続データとして変換します.grdsample でも-Rを指定できますが,120/300 のように指定をしても,意図したようには機能してくれませんでした.(必要に応じて,30秒グリッドのオリジナルデータを grdsample で解像度を変更しています.ここでは, -I2m で,2分角に変更しています.)最後に,正負を逆転(grdmath の NEG)しつつ, =cf でフォーマットをcf に変換しています.
※ 領域を Baba et al. (2017) に合わせました.
出来上がった gebco.grd は,JAGRUS のgridfile.dat で指定する地形データとして利用することが出来ます.ちゃんと確認していませんが,どうやらディレクトリは正しく処理できないようなので(たぶん,delimiterの違い?),実行ディレクトリにシンボリックリンクを作成するのが妥当かと思います.
gmt4は,port install gmt4 でインストールしています.
% makecpt -I -Cglobe -Z > gebco.cpt
% pscoast -U -Rgebco2min.grd -JQ150/18 -V -Dc -W1 -B30g30/15g15 -K > gebco2min.ps
% grdimage gebco.grd -J -Rgebco.grd -O -P -Cgebco.cpt >> gebco2min.ps
したのが以下の感じです(90度回転や切り出しは別途しています).
参考として,2分角メッシュでdt=5秒で24時間分の計算を,MacBook Airでやってみました.OMP_NUM_THREADS=2 で2並列にしています.参照→JAGURS on Mac.
[Timer Output]
+-------------------------------+----------+------------+
|Timer region |Called |Elapsed |
| | |Time[s] |
+-------------------------------+----------+------------+
|All | 1| 5663.732|
|getpar | 1| 0.001|
|read_grid_info | 1| 0.000|
|read_bathymetry_gmt_grdhdr | 1| 0.017|
|read_bathymetry_gmt_grd | 1| 1.331|
|read_friction_gmt_grd | 1| 0.167|
|initl_wfld | 1| 0.715|
|wet_or_dry | 1| 0.103|
|boundary_rwg | 1| 0.001|
|maxgrd_init_rwg | 1| 0.148|
|maxgrd_v_init_rwg | 1| 0.126|
|tgs_open_ct | 1| 0.000|
|hrise_rwg | 12| 0.394|
|drise_rwg | 12| 1.451|
|dump_gmt_nl_vel | 97| 0.002|
|dump_gmt_nl_hgt | 97| 55.068|
|tstep_grid_vel | 17281| 1265.025|
|- fxy_rwg | 17281| 1264.911|
|- mapgrids | 34562| 0.004|
|- outsea_rwg | 34562| 7.162|
|- recheck_wod | 34562| 0.003|
|maxgrd_v_check_nl | 17281| 2804.151|
|tstep_grid_hgt | 17281| 657.405|
|- hxy_rwg | 17281| 650.099|
|maxgrd_check_nl | 17281| 370.801|
|error_check | 17281| 271.747|
|maxgrd_write_gmt | 1| 0.558|
|maxgrd_v_write_gmt | 1| 0.528|
+-------------------------------+----------+------------+
MacProで4スレッド計算だと以下の感じです.
[Timer Output]
+-------------------------------+----------+------------+
|Timer region |Called |Elapsed |
| | |Time[s] |
+-------------------------------+----------+------------+
|All | 1| 3255.667|
(略)
[Timer Output]
+-------------------------------+----------+------------+
|Timer region |Called |Elapsed |
| | |Time[s] |
+-------------------------------+----------+------------+
|All | 1| 5663.732|
|getpar | 1| 0.001|
|read_grid_info | 1| 0.000|
|read_bathymetry_gmt_grdhdr | 1| 0.017|
|read_bathymetry_gmt_grd | 1| 1.331|
|read_friction_gmt_grd | 1| 0.167|
|initl_wfld | 1| 0.715|
|wet_or_dry | 1| 0.103|
|boundary_rwg | 1| 0.001|
|maxgrd_init_rwg | 1| 0.148|
|maxgrd_v_init_rwg | 1| 0.126|
|tgs_open_ct | 1| 0.000|
|hrise_rwg | 12| 0.394|
|drise_rwg | 12| 1.451|
|dump_gmt_nl_vel | 97| 0.002|
|dump_gmt_nl_hgt | 97| 55.068|
|tstep_grid_vel | 17281| 1265.025|
|- fxy_rwg | 17281| 1264.911|
|- mapgrids | 34562| 0.004|
|- outsea_rwg | 34562| 7.162|
|- recheck_wod | 34562| 0.003|
|maxgrd_v_check_nl | 17281| 2804.151|
|tstep_grid_hgt | 17281| 657.405|
|- hxy_rwg | 17281| 650.099|
|maxgrd_check_nl | 17281| 370.801|
|error_check | 17281| 271.747|
|maxgrd_write_gmt | 1| 0.558|
|maxgrd_v_write_gmt | 1| 0.528|
+-------------------------------+----------+------------+
MacProで4スレッド計算だと以下の感じです.
[Timer Output]
+-------------------------------+----------+------------+
|Timer region |Called |Elapsed |
| | |Time[s] |
+-------------------------------+----------+------------+
|All | 1| 3255.667|
(略)