Macで台風のシミュレーションを回してみる

# はじめに
この記事は初めてMacで気象データを扱う方向けに書いています。

また、こちらはScaleユーザーズガイドの3.2 現実大気実験を参考に作成しています。
合わせてご覧ください。
![wind_01_anim.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/514848/0378cc36-e00b-1cb5-f6be-b17438d44924.gif)
2015年台風15号のシミュレーション


# シミュレーション環境のセットアップ
セットアップは別記事にまとめましたので、こちらを参考にしてください。

- [Macで気象シミュレーションをするためのセットアップ](https://qiita.com/Tori-chan/items/0ef2b904c442edc75d4d)

# 入力データの取得

## scaleデータベースの準備
[こちら](https://scale.riken.jp/archives/scale_database.tar.gz)のサイトからSCALE-RM用に形式を変換した標高・土地利用区分のデータを取得します。
※リンクが切れている場合は[こちら](http://scale.riken.jp/download/index.html)の下の方にある`Topo and landuse data`のデータベースのリンクをクリックしてください

```
$ tar -zxvf scale_database.tar.gz
```
ダウンロードしたtarを解凍します。(できない場合はクリックで開きます)

データベースへのパスを設定します。

```terminal:~/.bash_profile
$ export SCALE_DB="${path_to_directory_of_scale_database}/scale_database"
# 例えば、ホームディレクトリにデータベースファイルをおいた場合は
$ export SCALE_DB="/Users/[ユーザー名]/scale_database"
```

またディレクトリ間の移動を楽にするために

```terminal:~/.bash_profile
export Tutorial_DIR="/Users/[ユーザー名]/scale-5.3.6/scale-rm/test/tutorial"
```
とすると便利です。`.bash_profile`を書き換えた後は`source ~/.bash_profile`をするかターミナルを立ち上げ直してください。


## 実データの取得
FNLdata（全球再解析データ）をダウンロードできるサイトがあるのでそこからダウンロードします。[こちら](https://rda.ucar.edu/datasets/ds083.2/)にアクセスしてデータを取得します。初回はユーザー登録が必要です。

- 登録が終わった後メインページのメニューバーから`Data Access > Web File Listing("Union of Available Products"の横) > Faceted Browse`を選択します。

- ダウンロードしたい日時(UTCで指定することに注意)を入力して、チェック項目はスルーし、一番下の「continue」をクリックします。
    - ノートPCで行う場合はデータ量が多いと計算に時間がかかるため期間を短くするのを推奨します。
    - 今回は2015年8月23日18時から2015年8月26日18時までの72時間のデータをダウンロードしました。

- ファイル⼀覧（.grib1 または.grib2 で終わるファイル名の書かれた表）が表⽰されるので、`filename`の列をクリックして全てのファイルをダウンロードします。


次に、

```
$ mkdir ${Tutorial_DIR}/real/tools/FNL_input/grib[1or2]/[再解析データの年]

(例) 
$ mkdir ${Tutorial_DIR}/real/tools/FNL_input/grib2/2015
```
として、解析データ年に対応したディレクトリを作っておきます。
今回は2015年8月23日18時から2015年8月26日18時までのデータを扱っているので2015のディレクトリを作成しました。

このディレクトリに先ほどダウンロードしたファイルを入れましょう。

その後、scale が読めるデータ（grd 形式）に修正する作業を⾏います。
grib1 形式のデータに対しても、convert_FNL-grib2grads.sh ファイルを使ってください。

```
$ cd ${Tutorial_DIR}/real/tools/
$ sh convert_FNL-grib2grads.sh [計算したいはじめの時刻] [計算したい終わりの時刻] FNL_input FNL_output
(例) sh convert_FNL-grib2grads.sh 2015082318 2015082618 FNL_input FNL_output
$ ls FNL_output/*/*
```
とコマンドを入力して(.grd)データができれば成功です。

# 実験セットの準備

一個前のディレクトリ(`scale-5.3.6/scale-rm/test/tutorial/real/`)に戻ります

```
$ cd ${Tutorial_DIR}/real/
$ vi USER.sh
```
として、USER.shを編集します。

```terminal:USER.sh
#!/bin/bash
  
#################################################
#
# &PARAM_COMM_CARTESC_NEST(ALL configs)
# &PARAM_TIME             (ALL configs)
# &PARAM_FILE_HISTORY     (run config)
# &PARAM_ATMOS_BOUNDARY   (run config)
# &PARAM_ATMOS_REFSTATE   (run config)
#
#################################################

NUM_DOMAIN=1 # set number of domains

RUN_DATE_YEAR=2015 #ダウンロードしたデータの開始年
RUN_DATE_MON=8 #ダウンロードしたデータの開始月
RUN_DATE_DAY=23 #ダウンロードしたデータの開始日
RUN_DATE_HOUR=18 #ダウンロードしたデータの開始時
RUN_DATE_MIN=0
RUN_DATE_SEC=0
RUN_DATE_MSEC=0

TIME_DURATION="72"[開始時から終了時までの時間]
#例えば、2015/8/23[18時]~2015/8/26[18時]なら72にする

...

# もし2005年5月31日以前のデータや2016年5月11日以降のデータを扱いたい場合は書き換える必要があります
# 表記については${Tutorial_DIR}/real/dataのファイルを参照してください
BASENAME_ORG="namelist.grads_boundary.FNL.2005053112-2016051106"
```

と修正し、そのディレクトリで`make`を行います

```
$ make

bash config/mkconf.sh


#################################################
#     SCALE-RM Configuration File Generator     #
#################################################

START DATE: 2015/08/23 - 18:00:00.000
```
コマンドを打つと上記のように表示されます。`START DATE:`部分がダウンロードしたファイルのスタート時間と一致していればOKです。
その後、real以下のディレクトリに「expreriment」というディレクトリができます。




## 地形データの作成
次に地形データを作成します。

```
$ cd ${Tutorial_DIR}/real/experiment/pp/
または、
$ cd scale-5.3.6/scale-rm/test/tutorial/real/experiment/pp/
```
で移動し、下記のコマンドを実行します。

```
$ mpirun -n 4 --oversubscribe ./scale-rm_pp pp.d01.conf
```
本環境では並列計算する領域の数が、物理コア数を超えているため`--oversubscribe`をつける必要があります。
(逆にクアッドコア以上のスペックがあるマシンの場合は、oversubscribeをつけなくても大丈夫です。)

また、実行には少し時間がかかり、実行時にアクセスの許可を求められる画面が出ることがあるので、その時は許可を押してください。


実行が終わった後、エラーメッセージが出ることもありますが、`ls`コマンドで`topo_d01.pe######.nc`と`landuse_d01.pe######.nc`ファイルができていればmpirunには成功しています。

また`pp_LOG_d01.pe000000`をいうログファイルに対して、以下のコマンドを打ち、

```
$ tail pp_LOG_d01.pe000000
...中略...
+++++ finalize MPI...
+++++ MPI is peacefully finalized
```
と出ていれば正常にデータが作成されています。


逆にこのようなログが出ていない場合はUSER.shなどで失敗している可能性が高いため、`experiment`のディレクトリを消して、USER.shを編集するところから再び始めてください。


## 初期値/境界値データの作成

```
$ cd ${Tutorial_DIR}/real/experiment/init
または
$ cd scale-5.3.6/scale-rm/test/tutorial/real/experiment/init
```
でディレクトリを移動してください。


```
$ cp ../../data/gradsinput-link_FNL.sh ./
$ vi gradsinput-link_FNL.sh
```
で`set start date`、`set end date`の部分をダウンロードしたファイルの時刻に修正します。
次に先ほどダウンロードしたファイルに初期値や境界値にリンクを貼ります。

```
$ sh gradsinput-link_FNL.sh
...
./ATM_00000.grd -> ../../tools/FNL_output/201508/FNL_ATM_2015082318.grd
./LND_00000.grd -> ../../tools/FNL_output/201508/FNL_LND_2015082318.grd
./SFC_00000.grd -> ../../tools/FNL_output/201508/FNL_SFC_2015082318.grd
./ATM_00001.grd -> ../../tools/FNL_output/201508/FNL_ATM_2015082400.grd
./LND_00001.grd -> ../../tools/FNL_output/201508/FNL_LND_2015082400.grd
./SFC_00001.grd -> ../../tools/FNL_output/201508/FNL_SFC_2015082400.grd
... 略(ダウンロードしたgrdファイルの数だけリンクが張られます)
```
このような感じのリンクが作成されていれば成功です。

次にnamelistにもリンクを貼ります。`mpirun`で初期値・境界値データの作成を行います。
リンクについては2016年よりも前のデータを使うと仮定して貼っていますが、もし2016年以降のデータを利用したい場合は`${Tutorial_DIR}/real/data`にあるファイル名を参考に置き換えてください。

```
$ ln -s ../../data/namelist.grads_boundary.FNL.2005053112-2016051106 ./
$ mpirun -n 4 --oversubscribe ./scale-rm_init init.d01.conf
```
mpirunが終了し、以下のファイルができていれば成功です。
boundary_d01.pe000000.nc
boundary_d01.pe000001.nc
boundary_d01.pe000002.nc
boundary_d01.pe000003.nc
init_d01_20150823-180000.000.pe000000.nc
init_d01_20150823-180000.000.pe000001.nc
init_d01_20150823-180000.000.pe000002.nc
init_d01_20150823-180000.000.pe000003.nc
init_LOG_d01.pe000000

また、次のコマンドで実行が正常に終了したかを確認することができます。

```
$ tail init_LOG_d01.pe000000 
...中略...
+++++ finalize MPI...
+++++ MPI is peacefully finalized
```

## シミュレーションの実行
```
$ cd ${Tutorial_DIR}/real/experiment/run
または
$ cd scale-5.3.6/scale-rm/test/tutorial/real/experiment/run
```
で移動して、次のコマンドでシミュレーションを実行します。

```
$ mpirun -n 4 --oversubscribe ./scale-rm run.d01.conf >& log &
```
データ数にもよりますが、すごく時間がかかるため(MacbookAir(Mojave)で4時間ほど)バックグラウンド実行でコマンドを打っています。

シミュレーションが進んでいるかどうかは、

```
$ tail -n 50 LOG_d01.pe000000
```
でログが書き換わっているかどうかで確認できます。
実行が終了した場合は、

```
+++++ finalize MPI...
+++++ MPI is peacefully finalized
```
と出力され、以下のファイルができたら成功です。
history_d01.pe000000.nc
history_d01.pe000001.nc
history_d01.pe000002.nc
history_d01.pe000003.nc


# 結果の描画
最後に作成したデータから画像を生成します。

## GrADSバイナリに変換

```
$ cd ${Tutorial_DIR}/real/experiment/net2g
または、
$ cd scale-5.3.6/scale-rm/test/tutorial/real/experiment/net2g
```
で移動し、2次元変数と3次元変数のnetCDF形式をGrADSで表示できるフォーマットに変更するために`mpirun`を行います。

```
$ mpirun -n 4 --oversubscribe ./net2g net2g.2D.d01.conf
$ mpirun -n 4 --oversubscribe ./net2g net2g.3D.d01.conf
```
`+++ MPI COMM: Corrective Finalize`と標準出力されていれば変換が成功しています。
成功すると、以下のファイルが生成されます。
MSLP_d01z-2d.ctl
MSLP_d01z-2d.grd
MSLP_d01z-2d_lccr.ctl
PREC_d01z-2d.ctl
PREC_d01z-2d.grd
PREC_d01z-2d_lccr.ctl
PRES_d01z-3d.ctl
PRES_d01z-3d.grd
PRES_d01z-3d_lccr.ctl
Umet_d01z-3d.ctl
Umet_d01z-3d.grd
Umet_d01z-3d_lccr.ctl
Vmet_d01z-3d.ctl
Vmet_d01z-3d.grd
Vmet_d01z-3d_lccr.ctl

## 計算結果の確認

```
$ cp ../../data/checkfig_real.gs ./
$ grads -blc checkfig_real.gs
```
を実行し、下記の名前のpngができていれば成功です。
real_mslp.png
real_prec.png
real_wind.png


また、gifファイルを作成したい場合は

```
$ touch wind_anim.gs
$ vi wind_anim.gs
```
でファイルを開き、以下の内容をコピペして`${Tutorial_DIR}/real/experiment/net2g`ディレクトリに保存してください。

```terminal:wind_anim.gs
'!mkdir wind'
'set display color white'
'set mproj nps'
'set mpdset worldmap'
'set xlint 20'
'set ylint 20'
'set parea 1 9.036 1 7.5'
'set gxout grfill'
'c'
**
** WIND
**
'set parea 1 9 1 7.5'
'open U_d01z-3d_lccr.ctl'
'open V_d01z-3d_lccr.ctl'
'set grads off'
***'set lev 500'
'set lev 850'
count=1
while(count<=73)
'set grads off'
'set t 'count
'color 5 45 5'
'd sqrt(u*u+v.2*v.2)'
'cbarn'
***'set clevs 5 10 15 20 25 30 35 40 45 50'
'set arrscl 0.5 60'
'd u;skip(v.2,8)'
'q dims'
temp=sublin(result,5)
tim=subwrd(temp,6)
'draw title Wind (m/s) @ 850hPa on 'tim
if(count<=9)
'printim ./wind/wind_01-00'count'.png'
else
if(count<=99)
'printim ./wind/wind_01-0'count'.png'
else
'printim ./wind/wind_01-'count'.png'
endif
endif
count=count+1
'c'
'!rm -f tmp'
endwhile
'close 2'
'close 1'
'!convert -loop 1 -delay 40 ./wind/wind_01-*.png ./wind_01_anim.gif'
'quit'
```

その後、

```
$ grads
で選択画面が出るため、エンターキーを押す

$ run wind_anim.gs
```
をすると下記のようなpngファイルと最初にあげたようなgifが生成されます。
![wind_01-009.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/514848/78077865-8ba7-5065-0ecd-faa17a796458.png)


# 参考
SCALE ユーザーズガイド:
https://scale.riken.jp/documents/

wind_anim.gs:
東大三浦研・灘本さんのスクリプトを改変させていただきました。