海洋上の波浪の数値シミュレーションを行う
波浪モデルとは
・波浪
海洋の波と一言で表しても、その種類はとても多く存在します。おそらく一番最初に頭に浮かぶのは以下に示すような風によって形成される風波（風浪）だと思います。（海の深度や地形由来（地震など）の波もありますが、ここでは一般的なものに絞ります）風波は比較的穏やかな波なので、海を見る機会があれば誰しもが見たことあるでしょう。この風波についてですが、波動論的には風によって発達過程にある波のことを風波と言います。したがって、風が強いほど風波は大きくなります。サーフィンやマリンスポーツの経験者が海上風を必ずチェックしているのは波のコンディションを知るためですね。海上風の風波と似た特徴を持つ波にうねりと呼ばれる波があります。うねりは風波の原因となる風が止んだり、風のない領域まで風波が進んだ場合の波をうねりと言います。波動論的には、風波が発達過程にあるのに対し、うねりは減衰過程であり、ゆっくりと進行します。うねりの特徴としては、波の波長が長いほど位相速度が速いというものがあります。この特徴によって、低緯度で発生した台風由来の風波が強風域を越えてうねりへと変化し、日本沿岸まで到達する現象が説明できます。
この風波とうねりを合わせて波浪といい、波浪モデルとは風（海上風）のデータなどから、風波とうねりを数値シミュレーションするものです。
・波浪モデル
波浪をシミュレーションする際に、大きな問題となるのが現実の空間とモデルの格子間隔の差です。多くの波が相互的に作用し合う波浪の形状を表現するには数10m程度の格子間隔が必要となります。しかし、海洋のシミュレーションをする際に、ある1点のみのシミュレーションを行うということは稀です。したがって、広大な海洋をシミュレーションする際には個々の波動をそれぞれ解析するという手法は現実的ではありません。ここで波動モデルの重要なポイントですが、個々の波動について直接計算をするのではなく、波の性質を応用した波浪スペクトルという値を使用し、その統計量の変化から波の変化を推察するという手法を用います。波浪スペクトルを式で表すと E = (θ, σ) となり、ここで Eは波のエネルギーを示し、θは波の伝播方向 (radian)、σは波の振動数 (1/sec)です。つまり、波をθ とσ を用いることで空間的にエネルギー分布として表現しているのです。したがって、波浪モデルは波浪スペクトルのエネルギーの時間変化を記述する以下のエネルギーバランス方程式が基礎方程式となっています。以下の式において、左辺第1項は波浪エネルギーの時間変化量、左辺第2項の∇（ナブラ）は今議論しているスカラー場にベクトル場のエネルギー値を対応させる演算子を示します。ここでは詳しい説明を省きますが、ベクトル場で記述されている波のエネルギー値は波の群速度 (Cg) によって表すことができますので、このような表現で波の移流を記述しています。（いわゆる移流項とよばれる項）右辺の項は、波浪の発達・減衰・消散に関係する外力について記述している項で、一般的にはソース関数と呼ばれます。Snet は波浪の形状に影響する外力の正味量という意味です。
次にソース関数についてですが、波浪をシミュレーションする上で必要となる要素として大きくは3つの項が考えられています。1つ目の要素としては、海上風によるエネルギーの入力項 Sin、2つ目は砕破などの波の消散・減衰にともなうエネルギー消散項 Sds、3つ目は非線形相互作用によるエネルギー輸送項 Snl です。（他にも海底との摩擦による外力を示す項 Sbot があります）海上風によるエネルギー入力項と砕破などによるエネルギー消散項については先述したように波浪を考慮する上では核となる要素です。3つ目の非線形相互作用によるエネルギー輸送項は、波の発達に関するエネルギー輸送を示しており、波同士の相互作用を記述するので非線形となり、計算量も多く難しいため初期の波浪モデルでは考慮されていませんでした。近年の技術の発達によってこの非線形のエネルギー輸送項の取り扱いが可能となりより精度の高いモデルが開発されています。
最近の研究で広く使われている波浪モデルは第3世代モデルと呼ばれるモデルです。具体的にはNOAA（アメリカ海洋大気庁）が開発・運用するWAVEWATCHや、オランダのデルフト工科大学が開発・運用するSWAN、BODC (イギリス海洋データセンター) が開発・運用するWAMがあります。
The University of Miami Wave Model (UMWM)
今回紹介する波浪モデルはマイアミ大学で開発された波浪モデル (UMWM) です。このモデルは2.5世代の波浪モデルです。上記の通り現在研究の多くで使用されているのは非線形のエネルギー相互作用を考慮した第3世代の波浪モデルですが、UMWMは非線形エネルギー相互作用の項を単純なバルク式で置き換えたモデルになります。バルクとは「おおまかな」という意味であり過去の実験のデータにより経験的にえられた値を定数として導入しているのです。この操作によって煩雑で計算量の多い非線形項の計算量を減らすことができ、モデル全体を軽くすることができています。実はモデルが軽いというのがUMWMの1番の特徴であり、全体のソースコードが短いのでモデル計算で何が行われているのかを把握しやすく、改良も容易です。また、モデル結果の精度もそれなりに高いレベルで計算されるので非常に扱いやすいモデルとなっています。
・UMWMのインストール
それではUMWMを Linux 環境でインストールし実験を行なっていきます。まずインストール (コンパイル) についてですが、こちらのマニュアルに従ってインストールを行います。コンパイルで使用するライブラリはNetCDFです。(詳しくは LOVECLIM 記事を参照)
NetCDFのインストールが終了したら、UMWMのコンパイルを行います。ダウンロードしたUMWMのgzファイルを展開しUMWM-1.0.0ディレクトリを作成します。
UMWM-1.0.0以下にあるソースコードディレクトリ src に移動し作業を行なっていきます。なお、ライブラリやモデルのコンパイルについては上記のLOVECLIM 記事にて詳しく説明していますのでそちらを参照してください。
では、コンパイルにあたっての最初のステップですが /src の Makefile を編集します。
# University of Miami Wave Model Makefile
# Compiler, flags, and libraty paths:
FC = mpiifort -save-temps -DMPI
#FCFFLAGS = -O3 -fastsse
FCFFLAGS = -g
CPPFLAGS =
LD = $(FC)
LDFLAGS = $(FCFFLAGS)
NETCDF = /home/lib/netcdf-4.1.3
NETCDFPATH = -I$(NETCDF)/include -L$(NETCDF)/lib -lnetcdff
NETCDFLIB = $(NETCDF)/lib/libnetcdff90.a
ここで重要なのはコンパイラの指定とNetCDFのパスの設定です。
UMWMのマニュアルに記載がありますが、UMWMはMPI環境での計算をすることができます。MPI環境とはコンピュータ (CPU) を並列化することで計算を短縮するオプションです。なので今回はコンパイラを ifort とするため FC = mpiifort としています。(-DMPIはコンパイラを使用する際にMPIを使いなさいという定義 (define)なのですが、私の環境ではこの定義が必要でしたので...)
次にNetCDFのパスはインストールしたディレクトリを指定してあげます。
以上でMakefileの編集は終了です。
最後に /src で make コマンドを入力しコンパイルします。これでエラーが出なければ /umwm-1.0.0
以下にumwm という実行ファイルが作成されます。エラーが出てしまった場合は、Makefile を修正しmake clean, make と続けてコンパイルをしなおしてください。
また、/tools というディレクトリがありますが、こちらは領域気象モデル (WRF) のアウトプットデータ (グリッドや水深データ) を UMWM で使用する際に使用をします。今回は別の方法を用いて実験を行ったので説明は省略します。
・UMWMの実験の準備
無事にコンパイルが終了したら、実験を始めていきます。
まず、UMWMを動かす上でのインプットデータとして必須なのは海上風のデータです。他にも深さのデータや海面流速のデータをオプションで入力することはできますが、今回は海上風のデータのみをインプットデータとして実験を行います。
海上風のデータとして使用するのは気象庁の予測モデルであるMSMのデータを使用します。MSMの空間解像度と時間解像度は 5 km ごと1 時間 です。MSMはオリジナルデータの間を予報値で埋める構成となっていますので、過去の大気データを再現するにはベストなのです。データはこちらから
ファイル形式は NetCDF 形式（.nc）であるためコントロールファイルを作成すれば（.ctl）簡単に図を作成することができます。こちらも詳しくは上記の LOVECLIM 記事を参照。
それではまず、MSMの u 成分, v 成分の海上風データ（10m）をダウンロードし、図を作って見ます。先ほどのリンクのページで地表面データが欲しいので、解析値を中心に再構成した NetCDF データ（MSM, RSM）でMSMを選び、次のページでMSM-Sのディレクトリを選択します。年のディレクトリの下に、1日ごとの .nc ファイルがまとめられているので、対象の期間の .nc ファイルをダウンロードします。今回の説明では2009年に日本南東の海上を進行した14号（TC0914）を対象としますので、期間を適当な5日間ほどとして、ダウンロードしました。
以下のアニメーションがダウンロードした海上風のデータをGrADSで描いた図になります。
この時に使用したコントロールファイルやスクリプトは省略します。
（スクリプトに関しては下記にて説明）
それでは、データの準備ができましたので UMWM を動かす準備に入ります。
通常、UMWM を動かすには地形データ（格子点間隔）が必要となりますが、今回は MSM のデータを使用するので、格子点が決まってしまっています。なので、ダウンロードしたデータを UMWM のインプット形式のNetCDFファイルに変換する作業を行います。
まず、先ほど紹介したマニュアルの6.4節 Input files の項目を見てもらえば分かりますが、作成したいNetCDFファイルの名前は umwmin_YYYY-MM-DD_hh:mm:ss.nc である必要があります。他にも umwm.gridtopo というファイルが必要となっていますが、前述した通りインプットファイルにこの情報を追加して書き込みます。以下のプログラムが NetCDFファイルを書き換えるプログラムです。
このプラグラムについてですが、形式は .gsとなっておりGrADS script の略です。.gs形式は.gex形式の実行ファイルに比べてできることが多くなります。
ではプログラムを見ていきます。最初に、プログラム全体が 「''」 で囲まれていることについてですが、これはループを指定しています。コマンド全体を囲むことで指定した範囲で再起的な処理を行います。
'reinit'
'open dmsm-s.ctl' MSMファイルを開く.ctlファイル
'q dims'
tt=14905 開始時刻の指定
while (tt<15242) 終了時刻の指定
'set x 340 441' x方向（経度方向）の範囲指定
'set y 140 241' y方向（緯度方向）の範囲指定
'set t 'tt 時間の指定
say tt
'q time'
mytime = subwrd(result,3) ここで組み込み関数を使い、q time で表示した時間の情報から mytime, myhour, myday,
myhour = substr(mytime,1,2) mymonth, myyear にそれぞれを指定しています。
myday = substr(mytime,4,2)
mymonth = substr(mytime,6,3)
myyear = substr(mytime,9,4)
if (mymonth = 'JAN') 月の表示を文字から数字へと変換します。
mymonth = '01'
else
if (mymonth = 'FEB')
mymonth = '02'
else
if (mymonth = 'MAR')
mymonth = '03'
else
if (mymonth = 'APR')
mymonth = '04'
else
if (mymonth = 'MAY')
mymonth = '05'
else
if (mymonth = 'JUN')
mymonth = '06'
else
if (mymonth = 'JUL')
mymonth = '07'
else
if (mymonth = 'AUG')
mymonth = '08'
else
if (mymonth = 'SEP')
mymonth = '09'
else
if (mymonth = 'OCT')
mymonth = '10'
else
if (mymonth = 'NOV')
mymonth = '11'
else (mymonth = 'DEC')
mymonth = '12'
endif
endif
endif
endif
endif
endif
endif
endif
endif
endif
endif
say myhour
say myday
say mymonth
say myyear
'uc=ABS(u*0)' ここからsdfwriteコマンドを使用してNetCDFファイルを作成していきます。
'set sdfwrite -flt tmp1.nc' ucは「x-component of ocean surface currents」のことで、今回の
'sdfwrite uc' の実験では使用しないので 0 というデータを記入します。
**************************
'vc=ABS(v*0)' vcも同様です。ちなみに海表面の流速成分を意味し単位は[m/s]です。
'set sdfwrite -flt tmp2.nc'
'sdfwrite vc'
'!ncks -h -A tmp2.nc tmp1.nc' !はターミナル上でのコマンドです。tmp1.nc にtmp2.nc を加えています。
**************************
'aa = const(maskout(1,100-psea+sp),0,-u)' u10（海上10mの風速）は今回必要な値となるので書き出すのですが、
'u10=u*aa' 陸域のデータを取り除くように工夫をここでしています。
'sdfwrite u10 tmp2.nc' psea（海表面気圧）とsp（地表面気圧）の和を100からひいてマイナスとなった
'!ncks -h -A tmp2.nc tmp1.nc' 地点を陸域として欠損値 0 を、プラスとなった海域に1を与えています。
************************** 最後に MSMでの風速値uをかけて完成です。
'aa = const(maskout(1,100-psea+sp),0,-u)' 同様に v10についても実行します。
'v10=v*aa'
'sdfwrite v10 tmp2.nc'
'!ncks -h -A tmp2.nc tmp1.nc'
**************************
*'aa = const(maskout(1,100-psea+sp),0,-u)'
*'z=100*aa'
*'sdfwrite z tmp2.nc'
*'!ncks -h -A tmp2.nc tmp1.nc'
************************
'!mv tmp1.nc umwmin_'myyear'-'mymonth'-'myday'_'myhour':00:00.nc' ここで、先ほど使用した組み込み関数を用いてUMWMの
*'!rm tmp*.nc' インプット形式の名前に指定します。
*say result
tt=tt+1 最後に時間を 1進めてループの指定をします。
endwhile
'disable fwrite'
say 'finished'
この .gsファイルをGrADS上で runコマンドを使用して実行すれば、前述したUMWMのインプット形式のファイルが1時間毎のデータとして出力されます。
・UMWMの実験
インプットファイルの準備が整ったのでUMWMの実験を開始していきます。
まず、作成したインプットファイルを全て /umwm-1.0.0/input 以下に移動させます。
次に、UMWMの初期設定を行います。inputディレクトリから1つ戻り、umwm-1.0.0以下の /namelists にあるmain.nml を編集します。適当なエディターで開きます。
基本的にはコメントアウトで説明されていますが、重要な箇所を記述します。
&DOMAIN
isGlobal = .FALSE. ! Global (.TRUE.) or regional (.FALSE.)
mm = 282 ! Domain size in x x・y方向（経度・緯度方向）の範囲指定です。
nm = 282 ! Domain size in y inputファイルの大きさと合わせます。
om = 37 ! Number of frequency bins
pm = 24 ! Number of directions
f1 = 0.0313 ! Lowest frequency bin [Hz]
f2 = 2.0 ! Highest frequency bin [Hz]
startTimeStr = '2009-09-18_00:00:00' ! Simulation start time 実験開始時時刻と終了時刻の指定です。
stopTimeStr = '2009-09-23_00:00:00' ! Simulation end time
dtg = 3600 ! Global (I/O) time step [s]
restart = .TRUE. ! Restart from file 実験を途中から始めるための中間ファイルの作成の有無
/ の設定。
&PHYSICS
g = 9.80665 ! Gravitational acceleration [m/s^2]
nu = 0.9E-6 ! Kinematic viscosity of water [m^2/s]
sfct = 0.07 ! Surface tension
kappa = 0.4 ! Von Karman constant
z = 10. ! Height of the input wind speed [m]
gustiness = 0.0 ! Random wind gustiness factor (should be between 0 and 0.2)
dmin = 10. ! Depth limiter [m]
explim = 0.7 ! Exponent limiter (0.69 ~ 100% growth)
sin_fac = 0.11 ! Input factor from following winds
sin_diss1 = 0.10 ! Damping factor from opposing winds
sin_diss2 = 0.01 ! Damping factor from swell overrunning wind
sds_fac = 42. ! Breaking dissipation factor
sds_power = 2.5 ! Saturation spectrum power
mss_fac = 120. ! Mean-square-slope adjustment to Sds
snl_fac = 5. ! Wave energy downshifting factor
sdt_fac = 0.01 ! Dissipation due to turbulence factor
sbf_fac = 0.001 ! Bottom friction coefficient [m/s]
sbp_fac = 0.001 ! Bottom percolation coefficient [m/s]
/
&GRID
gridid = 2 ! Set 1 for Cartesian or 2 for spherical projection
delx = 0.0625 ! Grid spacing in x and y in: meters for gridid=1, MSMデータの格子間隔を記入します。
dely = 0.05 ! degrees for gridid=2
lat_s = 24.85 ! Latitude of bottom left point; used for gridid=2 緯度・経度の開始点を記入します。
lon_w = 132.4375 ! Longitude of bottom left point; used for gridid=2
/
&INPUT_FLAGS ! Set to 1 for input from file, or 0 for homogenous
! fields input from &INPUT list below;
dflag = 0 ! Bathymetry
excom = 0 ! Set to 1 to exclude user chosen basins/lakes
wflag = 1 ! Winds inputファイルを使用するので 1。
pflag = 0 ! Surface pressure
tflag = 0 ! Air temperature
cflag = 0 ! Currents
rflag = 0 ! Water density
/
&INPUT
depth = 3000 ! Water depth [m] 深度を3000m に固定。
u0 = 25 ! Wind magnitude [m/s]
wdir0 = 0. ! Wind direction [rad]
sfcp0 = 101325. ! Surface pressure [Pa]
temp0 = 300. ! Surface air temperature [K]
uc0 = 0. ! x-component ocean current [m/s]
vc0 = 0. ! y-component ocean current [m/s]
rhow0 = 1025. ! Water density [kg/m^3]
/
&OUTPUT
outgrid = .TRUE. ! Gridded output
outspec = .FALSE. ! Spectrum output
xpl = 45 ! Grid cell in x for stdout (screen)
ypl = 6 ! Grid cell in y for stdout (screen)
restart_out = .TRUE. ! Restart output
基本的にはデフォルトの状態で問題ありません。物理量を変更したいときはここの項目を変更することで反映させることができます。
これで、 /umwm-1.0.0ディレクトリで実行コマンド ./umwm を打ち込み実験を開始させられるのですが、ここで実験の計算量についてです。今回はUMWMのコンパイルで触れたようにMPI環境で計算します。しかし、計算量が膨大であるためMPI環境での計算（いわゆるスパコン）を使用してもかなり時間がかかりますので、MPI環境を使用しない場合はかなりの時間がかかってしまいます。
このようなことから main.nmlで記述したようにリスタートオプションをTRUEとすることで、途中から計算を行えるようにしているのです。
また、スパコンでの計算は機械によって仕様が異なるため、確実な説明をできませんが基本的にはノード（CPUやメモリの塊）とコア（CPU）を設定してsubmitコマンドで計算を行います。
私の例を記述しますと、2ノード・24コアで 1時間計算をさせました。1時間と区切った理由はスパコンの計算開始までの時間を短くしたいことと、バグが出た場合にすぐ修正をするためです。
・UMWMの実験結果
計算が問題なく終了すれば /umwm-1.0.0/output ディレクトリに umwmout- というNetCDFファイルが作成されていると思います。今回の実験での個人的な忘備録としては計算開始時刻のoutputデータの値が0になってしまうため、対処が必要なことでした。）それでは、.ctlファイルを作成して図を描いていきます。
.ctlファイルはこちら
今回作成するアニメーションはこちらです。
まず、アニメーションについての説明ですが、色で示しているのが Significant wave height (有義波高)です。有義波高とは有義波の平均で単位は [m]です。有義波とはある一定期間に観測される波の波高と周期を高いものから順に全体の 1/3 までを抽出し、平均したものです。よく天気予報で言われる波高とはこのことを指します。
このアニメーションでは矢印で inputした海上風も示しているので、台風の進行にともなって有義波高が移動しているのがよく分かります。また、風の吹き込む台風の右側の風が強く、有義波高も高いということが分かります。
では最後に、このアニメーションの作成スクリプトについて解説をします。スクリプトはこちら
NetCDFファイルを書き換えた時と同様にループを使用するので「''」で囲ってあります。また、MSMのデータを描写した時のスクリプトもほとんど同じなので、参考にしてみてください。
'reinit'
'open umwmout.ctl'
tt=121 開始時刻の指定をします。
'set display color white'
while (tt<181) 終了時刻の指定をします。
'c'
'set t 'tt
'q dims'
say tt
result2 = sublin(result,5) 作成する図に時刻を表示させるため、組み込み変数を使用します。
say result2
mydate = subwrd(result2,6)
say mydate
'set lon 133 150' 経度の領域を指定します。
'set lat 25 39' 緯度の領域を指定します。
'set grid off'
'set grads off'
'set rgb 19 130 0 220'
'set rgb 20 30 60 255'
'set rgb 21 0 160 255'
'set rgb 22 0 200 200'
'set rgb 23 0 210 140'
'set rgb 24 0 220 0'
'set rgb 25 160 230 50'
'set rgb 26 230 220 50'
'set rgb 27 230 175 40'
'set rgb 28 240 130 40'
'set rgb 29 250 95 50'
'set rgb 30 250 60 60'
'set rgb 31 240 0 130'
'set rgb 32 240 0 200'
'set rgb 33 240 0 240'
'set clevs 0 1 2 3 4 5 6 7 8 9 10 11' 有義波高の色分けの設定です。他にも様々な方法があるので
'set ccols 20 21 22 23 24 25 26 27 28 29 30 31 32' いろいろと調べてみてください。
'set gxout shaded' 色分けの設定です。
'd swh' 有義波高をここで描写しています。
'run cbarn.gs 0.9 1 9.7 4.3' カラーバーの表示設定です。
'draw ylab Latitude'
'draw xlab Longitude'
'set ccolor 1'
'set cthick 10'
'set arrscl 1 40'
'd skip(uwnd,12,12);skip(vwnd,12,12)' ここで、矢印の設定をした後、inputした風速を表示しています。
'draw title Swh (color) and Wind speed (arrows) 'mydate''
'basemap.gs L 0' basemap.gsというスクリプトを使用して陸域をマスクしています。
'printim im-'tt'.png x1000 y800' 画像の保存をしています。
tt=tt+1 ループの設定をしています。
endwhile
'!convert -delay 20 im-* swh.gif' アニメーションの作成コマンドです。
'!rm im-*'
陸域のマスクに basemap.gsを使用していますが、このあたりのスクリプト（カラースケールなど）はネットで公開されていますので調べてみると良いと思います。