NCEP Final analysis、Reynolds SST（High Resolution) のデータを以下のサイトから
ダウンロードしてそれぞれ適当なdirectory 
(以下では/home/maro/DATA/FNL/20150910/, /home/maro/DATA/RTGSST/20150910/以下)
に置いておく. 以下にダウンロードの方法を記す. 

注意: 特定の事例の再現実験では、見たい日時の数日前後の分も含めて、
データをダウンロードして保存しておく. 

FNLのダウンロード方法
NCEP Research Data Archive
一番下のData Accessをクリック
→ Web File Listing (期間によってGRIB1かGRIB2のいずれか,このケースはGRIB2)
→ Faceted Browse 
→ Temporal (Valid Time) Selection で期間(対象期間の数日前〜数日後)を指定 (Parameter Selectionはマークしない)
→一番下の Continueをクリック
→必要なファイルのboxをマークして Downloadをクリック
→datafiles.tarが/home/maro/ダウンロード/以下にダウンロードされるので、
作成したdirectoryに移して展開(tar xvf datafiles.tar). 

SSTのダウンロード方法（FNLのSSTを使う場合は不要）
Verification and archival of the RTG_SST_HR analysis
ftp://polar.ncep.noaa.gov/pub/history/sst/rtg_high_res (from two days ago back to Sept.27, 2005) と書いてあるところをクリック
→必要なファイル(rtg_sst_....) をクリックしてダウンロード. 
→rtg_sst... が/home/maro/ダウンロード/以下にダウンロードされるので、
作成したdirectoryに移す. 

他の客観解析データ、Vtableについては
こちら


[maro@shear7 ]$ cd WPS_V3.5.1                        
[maro@shear7 WPS_V3.5.1]$ mkdir WPS_20150910_FNL              # 適当なdirectoryを作る
[maro@shear7 WPS_V3.5.1]$ cp -rf WPS/*  WPS_20150910_FNL/.    # 実行形式などコピーしてしまう
[maro@shear7 WPS_V3.5.1]$ cd WPS_20150910_FNL/ 
[maro@shear7 WPS_20150910_FNL]$ emacs namelist.wps &         # 期間、領域、格子サイズなどを設定するファイルを編集. 
とりあえず以下のものをコピーして使うこと（geog_data_pathは自分の地形データの場所に変更すること）
 namelist.wps 
編集したnamelist.wpsはどこかに保存しておき、以降はそれをコピーして適宜
編集すると良い. 
 comments on namelist.wps>

注意: start_date (end_date)は、あとでWRFで計算することになると思われる期間より
充分前(後)にとっておき、長めにデータを処理しておくと後で二度手間にならない. 
例えば、9月10日の事例に着目して計算する場合、start_dateは 9月7日、end_dateは9月13日
とかにしておけば良い. 
また、start_date, end_date の時刻は 00時 にしておくと、SSTデータの処理の際
変更しなくてすむので面倒がない. 

*FNLデータのunpack(展開）
namelist.wps中の prefix= 'FILE' とする. 

[maro@shear7 WPS_20150910_FNL]$ ln -sf ungrib/Variable_Tables/Vtable.GFS Vtable 
[maro@shear7 WPS_20150910_FNL]$ ./link_grib.csh  /home/maro/DATA/FNL/20150910/fnl
[maro@shear7 WPS_20150910_FNL]$ ./ungrib.exe
[maro@shear7 WPS_20150910_FNL]$ ls -l 
FILE:yyyy-mm-dd_hh ができていることを確認

*SSTデータのunpack (FNLのSSTを使う場合はこの作業は不要）
[maro@shear7 WPS_20150910_FNL]$ emacs namelist.wps & 
prefix= 'FILE'を prefix= 'SST'に書き換える
[maro@shear7 WPS_20150910_FNL]$ ln -sf ungrib/Variable_Tables/Vtable.SST  Vtable 
[maro@shear7 WPS_20150910_FNL]$ ./link_grib.csh  /home/maro/DATA/RTGSST/20150910/rtg
[maro@shear7 WPS_20150910_FNL]$ ./ungrib.exe
[maro@shear7 WPS_20150910_FNL]$ ls -l 
SST:yyyy-mm-dd_hh ができていることを確認

*地形データの処理
[maro@shear7 WPS_20150910_FNL]$ ./geogrid.exe
[maro@shear7 WPS_20150910_FNL]$ ls -l 
geo_em.d01.nc ... ができていることを確認. 

*WRFの入力データの作成 (FNL、SST、地形データから作成）
[maro@shear7 WPS_20150910_FNL]$ ./metgrid.exe
[maro@shear7 WPS_20150910_FNL]$ ls -l 
met_em.d01.yyyy-mm-dd....ができていることを確認. 
なお、FNLデータのSSTを使う場合はnamelist.wps の
fg_name = 'FILE','SST'は
fg_name = 'FILE' としておくこと。

[maro@shear7 WPS_20150910_FNL]$ ncdump -h  met_em.d01.....(fileはどれでもよい) | more 
として、冒頭に出てくる 
num_metgrid_levels の値  (古いFNLデータの場合は27、最近のデータは32、さらに最近のデータは34になっている)
num_st_layers, num_st_layers の値  (古いFNLデータは2, 最近のデータは4になっている)
を確認しておくこと。

以下でWRFを実行する際の制御ファイル namelist.inputの中の
num_metgrid_levels, num_metgrid_soil_levels 
の値は、それぞれこれらの値に合せること。

注）最近のfnlデータ（2017年後半以降）は上記の方法では処理できないが、WPS ver4.0をコンパイルして
(その歳、 WRF_ver4.0のコンパイルは不要）、できたungrib.exe, ungrib/Variable_Tables/Vtable.GFS
を古いversionのWPSにコピーすれば処理できます. 

geo_em_d0*.nc や met_em_d0*.ncなど netcdf形式の入力データは
read_wrf_nc.f で読み込み、改変することができる。
read_wrf_nc.fのページ 

以下、モデルの地形(標高)を変える方法について説明する。
以下の地形改変用にread_wrf_nc.f を修正したものをコピーして、自分の課題に合わせ、
適宜修正して使うこと。
read_wrf_nc_d01.f

修正した箇所は プログラム中の subroutine USER_CODE (最後にあります) で、
! modified by Kawashima と書いてあります。

座標 A (icent, jcent)を中心に、Aからの距離に応じて地形の標高を下げるよう
にしています。
点Aからの距離が十分離れた所では、元の標高と殆ど変わらないようにしています。

地形を修正する場合は、geo_em_d0*.ncを修正します。
その後で metgrd.exeを行えば、修正した地形を反映した入力データ met_em_d0*.ncができます。

コンパイル (geo_em_d0*nc のあるdirectoryで行う)
[maro@shear7]$ ifort read_wrf_nc_d01.f -L/usr/local/netcdf-3.6.3/lib -lnetcdf -lm -I/usr/local/netcdf-3.6.3/include  -free  -o read_wrf_nc_d01

実行
[maro@shear7]$ ./read_wrf_nc_d01 -EditData  HGT_M  geo_em.d01.nc 
INPUT FILE IS: geo_em.d01.nc

 Attempting to open netCDF file with write access

 CAUTION variable HGT_M      is about to change. Continue? (yes/no)

と出るので yes とすれば、地形の標高(HGT_M)が修正される。

Nestingして実験する際は、read_wrf_nc_d01.fをコピーして適宜書き換え
geo_em_d02.nc, geo_em_d03.nc....も同様に修正すること。
(icent, jcent, dx, dy の値を書き換える必要があります）


同様の方法で、標高以外のデータ、例えば地表面の種類も修正できます。
以下は、陸面を海面に変更するプログラムの例で、上と同様のやり方で、
HGT_Mの他、LAMDMASK, LANDUSEF, LU_INDEXを変更します。
変更する領域などは各自問題に合わせ、適宜修正すること。
read_wrf_nc_d01_landtosea.f

注）最近のFNLデータでは, met_em** のSNOW, SNOWH を0に直すべき場合がある
(wrfの出力のtsk, snowcをチェックすること）

[maro@shear7 ]$ cd /home/maro/WRF_V3.7.1/WRFV3/test
[maro@shear7 test]$ mkdir em_20150910                         # 計算用direcotoryの作成
[maro@shear7 test]$ cp -rf em_real/*  em_20150910/.           # 実行形式など、全部コピーしてしまうと面倒がない.  
[maro@shear7 test]$ cd  em_20150910/
[maro@shear7 em_20150910]$ emacs namelist.input &             # 計算の設定ファイルを編集、とりあえず以下のものをコピーして使うこと. 
namelist.input

注: em_real/以下のdefaultのnamelist.inputは書き換えるべきところが非常に多いので、
一度作った namelist.inputは保存しておき、以降はそれを適宜編集して使うと良い. 
    
namelist.inputの説明(よく変更する箇所のみ)

詳しくはWRF User's Guide Chapter 5 を参照のこと


注: いきなり最初から Nestingして広解像度の大規模な計算をすると計算時間がかかるので、始めは格子の粗い一番外側の領域のみで計算して (namelist.input中で max_dom = 1とする)、その結果を見てみると良い. 
この段階で初期時刻等を変えて計算を実行して、最適な設定を見つけ、その後Nestingした計算をするようにしたほうが効率が良い. 

[maro@shear7 em_20150910]$ ulimit -s 819200000       # システムリソースのstack size の上限変更. 
これをしないと大規模な計算ができません. 現在の状況は　ulimit -a　で調べられます. 

[maro@shear7 em_20150910]$ ln -sf /home/maro/WPS_V3.7.1/WPS_20150910_FNL/met_em.d0* .    # WPSで作った入力データ(met_em.d0...)にリンクする. 
[maro@shear7 em_20150910]$ ./real.exe                # WRF計算用の初期値、境界値などの作成
[maro@shear7 em_20150910]$ ls -l 
wrfbdy_d01 (側面境界値データ, d01のみ),  
wrfinput_d01, wrfinput_d02, ... （初期値+底面境界データ, max_domで指定した領域の数だけできる) 
ができていることを確認. 
海面水温をupdateする設定の場合、
wrflowinp_d01, wrflowinp_d02, ...(時間変化する海面水温データ, max_domで指定した領域の数だけできる) 
もできます. 

[maro@shear7 em_20150910]$ mpiexec -n 16 ./wrf.exe < /dev/null &    # 計算開始 -n 16 は 16並列計算. 

# 各PCの CPU coreの数は $ cat /proc/cpuinfo で確認できます. 

（shear6: 24個,  shear7: 32個,  shear8:32個,  shear9:16個, shear10とshear11: 4個) 


連属して複数の計算を行う場合、以下のようなバッチファイルを用意して実行すれば、一つの計算が
終ったら次の計算をしてくれます．(各実験の計算directoryで./real.exeまではやっておき、境界値
ができているか確認しておくこと）

バッチファイルの例:
batch_mpi.sh

バッチジョブの実行は
[maro@shear7] sh batch_mpi.sh &   　

[maro@shear7 em_20150910]$ ls -l 
計算がされていれば、wrfout_d01....という出力ができます(多重Nestingの場合
wrfout_d01.., wrfout_d02.., wrfout_d03...）. 

ちゃんと計算されているか、1ステップの計算にどのくらい時間がかかっているかチェックするには、ログファイル
rsl.out.00**, rsl.error.00** (16並列計算の場合、**は00〜15)を 
[maro@shear7 em_20150910]$ tail -f rsl.out.0000 
[maro@shear7 em_20150910]$ tail -f rsl.error.0000
などとして見れば分かります. 

注: 例えば、namelist.input中の time_step の設定が長すぎたりすると、
ONE OF THE PROCESSES TERMINATED BADLY: CLEANING UP
APPLICATION TERMINATED WITH THE EXIT STRING: Killed (signal 9)
と出て計算が止まります. 

rsl.out.00** もしくはrsl.error.00** (12並列の場合、**は0〜11のいずれか)をチェックすると
d01 2015-09-08_00:16:40           22 points exceeded cfl=2 in domain d01 at time 2015-09-08_00:16:40 hours
d01 2015-09-08_00:16:40  MAX AT i,j,k:           96         120          43 vert_cfl,w,d(eta)=   2.974685     -0.2250476      3.1409189E-02
などと出ており、CFL条件が満たされてない旨のメッセージが確認できますので、
その場合はtime_stepを短かくして計算すること.

*風速、地形、鉛直格子間隔などの条件によるが、time_step はd01の格子サイズが X kmなら、3X〜5X (sec) 程度にすると良い (これより短いぶんには問題ないですが、その分計算時間がかかるのでそれとの兼ね合いで).