2010年 4月 1日 公開， 2012年 1月 9日 修正
”パソコン”で気軽にためして遊べるような簡易な海洋数値モデルを作ってみた． 非常にシンプルなプログラムなので， ぜひ動かして，改造して，海洋物理を楽しみながら学んでほしい．
SLTM（Single Layer Toy Model）は， 初学者向けのシンプルな海洋一層矩形モデルである． 取っつきやすさと，改造のしやすさを追求して開発した．
まずは，SLTM のプログラム的な特徴をいくつか挙げてみる．
実行には，ソースファイルが 1つあればよい． 適当なソースファイル（たとえば sltm_rossby.c）をダウンロードして，
$ cc sltm_rossby.c -lm $ ./a.out
の，たった 2行の操作で，シミュレーションがはじまる． これ以上簡単な操作があるだろうか． 本当はコンパイルオプションも不要にしたかったのだが， さすがに -lm は取り去れなかった．
コードは素朴な C で書いてあるので，移植性も抜群，というか， おそらく Cコンパイラのインストールしてあるほとんどのコンピュータで， 何もしなくてもコンパイルできると思う．
コードは 600行程度と短い上，高速化をまったく考えずに， 方程式をそのままの形でコーディングしてあるので， 処理の内容がわかりやすく，改造もしやすくなっている．
初期値の作成，パラメータのセットもコード内で行っているので， 変数やパラメータの追加を気軽にできる． 初期値ファイルやパラメータファイルのフォーマットを覚える必要がないのも楽でいい． ファイル入力がないことは，コードのシンプル化にも貢献している．
数値モデルとしての特徴も挙げておこう．
1層の矩形モデルなので，概念的な実験が SLTM の主な用途になるだろう．
前節で紹介したように，初期条件，パラメータが 1つのソースファイルに組みこまれている． 今回は，4つの現象について 4つのソースファイルを用意した．
それぞれの解説は，
で行っている．
実行方法は，ダウンロードしたソースファイル（たとえば sltm_rossby.c）を適当なディレクトリに置き，
$ cc sltm_rossby.c -lm $ ./a.out
これだけである．
出力ファイルはタブ区切りのアスキーファイルで， 1行に，経度λ，緯度φ，経度方向の流速 u，緯度方向の流速 v，境界面の変位 h が順に並んでいる．
以下の免責事項に同意していただくことを，SLTM の使用条件とします．
損害が発生しそうな用途に使う場合は，ソースを良く読んで，納得してから自己責任でお使いください．
SLTM は，どんどん改造してもらうことや， 海洋物理の学習のために使っていたただくことを想定しているので， このままでも，改造したものでも，自由に配布していただいて結構です． ただし，ソース冒頭の作成者情報のコメント
/* * SLTM - Single Layer Toy Model (C) 2010 KDO */
は残していただくよう，お願いします． 改造して配布される方は，その旨コメントを追加していってください．
SLTM に実装している方程式について解説する．
運動方程式 2つと連続の式 1つの 3つの式から，u，v，h の時間差分を計算する．
ここで，
各変数は，
である．
1層（1.5層）モデルの模式図（上面が境界面の場合）
時間積分は，リープフロッグにときどき松野スキームを混ぜている．
リープフロッグ：
松野スキーム：
ここで，u は状態ベクトル， n は時刻ステップ数， f は du/dt を計算する関数．
コードは 600行程度と短いので，全部解説してしまおう．
パラメータの類は #define でセットしている． 順に解説していこう．
| OUTPUT_NUM | 計算結果をいくつ出力するか． |
| OUTPUT_STEP | 計算結果を何ステップごとに出力するか． |
| OUTPUT_FILENAME_LENGTH OUTPUT_FILENAME_FORMAT |
出力ファイル名のフォーマットを指定する．関数 output で使う． |
| DT_SEC | Δt：計算 1ステップの時間．単位は s． |
| MATSUNO_STEP | SLTM の時間積分は，基本的にはリープフロッグだが，MATSUNO_STEP ごとに松野スキームを混ぜている． |
| ILON_NUM | 経度方向のグリッド数． |
| ILAT_NUM | 緯度方向のグリッド数． |
| LON0_DEG | 西端の経度．単位は度． |
| LAT0_DEG | 南端の緯度．単位は度． |
| DLON_DEG | Δλ：経度方向のグリッド間隔．単位は度． |
| DLAT_DEG | Δφ：緯度方向のグリッド間隔．単位は度． |
| LON_DEG | 座標(ilon，ilad) のグリッドの経度の定義式．単位は度． |
| LAT_DEG | 座標(ilon，ilad) のグリッドの緯度の定義式．単位は度． |
| PI | 円周率 π． |
| DLON_RAD | Δλ：経度方向のグリッド間隔 DLON_DEG の単位をラジアンに変換した値．計算にはこれを使う． |
| DLAT_RAD | Δφ：緯度方向のグリッド間隔 DLAT_DEG の単位をラジアンに変換した値．計算にはこれを使う． |
| LON_RAD | 座標(ilon，ilad) のグリッドの経度 LON_DEG の単位をラジアンに変換した値．計算には使われないが一応． |
| LAT_RAD | 座標(ilon，ilad) のグリッドの緯度 LAT_DEG の単位をラジアンに変換した値．計算にはこれを使う． |
| GRAVITY | 重力加速度 g．単位は ms-2．1.5層モデルでは，reduced gravity の値をセットする． |
| DENSITY | 海水密度 ρ0．単位は kgm-3． |
| EARTH_ROTATION | 地球自転角速度 Ω．単位は s-1． |
| EARTH_RADIUS | 地球の半径 a．単位は m． |
| RAYLEIGH_FRICTION | 摩擦係数 R．単位は s-1 |
| NEWTONIAN_COOLING | 鉛直拡散係数 γ．単位は s-1 |
| HORIZONTAL_VISCOSITY | 水平粘性係数 Ah．単位は m2s-1 |
| LAYER_THICKNESS | 層厚 H0 の基本値．単位は m．関数 set_layer_thickness で層厚をセットするとき参考に使う． |
このあとに，計算する座標(ilon, ilat) のグリッドと， それに隣接するグリッドのグリッド番号（変数配列の添字）を示す記号の定義を行っている． ここが，SLTM の売りのひとつである．
#define C ILON_NUM*ilat+ilon
#define W C-1
#define E C+1
#define S C-ILON_NUM
#define N C+ILON_NUM
#define ES C+1-ILON_NUM
#define WN C-1+ILON_NUM
#define WW W-1
#define EE E+1
#define SS S-ILON_NUM
#define NN N+ILON_NUM
/*
NN
|
WN N
\|
WW--W--C--E--EE
|\
S ES
|
SS
*/
座標(ilon, ilat) のグリッドを「C」として， 北隣のグリッドは「N」，東隣のグリッドは「E」という具合に定義していく． これらの定義により，方程式を実装する際に， グリッドの座標やグリッド番号を意識する必要がなくなり， シンプルで読みやすく，間違いの少ないコードを書くことが可能となった．
2次元モデルをプログラムするときは，変数を 2次元配列で書くのがふつうだろうが， そうすると，たとえば 2階微分 d2u/dx2 のコードは，
(v[j][i+1].u - 2.0 * v[j][i].u + v[j][i-1].u) / dx / dx
のように繁雑になってしまう．これが，SLTM では，
(v[E].u - 2.0 * v[C].u + v[W].u) / dx / dx
のようにシンプルに書ける． グリッドの座標で記述するのではなく， 相対的な位置関係で記述するので， コーディング時に間違えにくいし，読みやすい．
変数とグリッド情報を構造体で記述している．
struct var に，経度方向の流速 u，緯度方向の流速 v， 境界面の変位 h をまとめて記述している． u，v の単位は ms-1，h の単位は m． SLTM では，Arakawa C-grid を採用しており， var に属する u，v，h の位置関係は，コードにも書いてある通り，
| | -+----+- |v | | | |h u| -+----+- | |
である． LON_DEG, LAT_DEG, LON_RAD, LAT_RAD は h の位置の値である．
struct cell には，各グリッドの層厚 layer_thicness とエクマン湧昇 ekman_upwelling を記述している． layer_thickness は，関数 set_layer_thickness で， ekman_upwelling は，関数 set_ekman_upwelling でそれぞれセットする．
各関数をざっと解説しよう．まあ，関数名が，そのままやっていることではあるのだが．
| set_layer_thickness | 各グリッドの層厚 layer_thickness をセット． |
| set_ekman_upwelling | 各グリッドのエクマン湧昇 ekman_upwelling をセット． |
| set_initial_condition | 初期条件 u，v，h をセット． |
| make_mount | cos形の山を作るのに使う関数．グリッド座標（ilon，ilat），山の中心座標，山の半径（どちらもグリッド単位），山の高さ（m）を指定すると，座標（ilon，ilat）の標高が返ってくる．set_initial_condition で使う． |
| calc_boundary_condition | 境界条件を計算する関数．かなり簡単な処理で済ませている． |
| calc_momentum_equations_first | u，v の運動方程式の一部を計算．メトリック項，圧力傾度力項，コリオリ項を計算している．時間差分 f の値を上書きするので，calc_momentum_equations_add_advection，calc_momentum_equations_add_viscosity よりも先に書くこと． |
| calc_momentum_equations_add_advection | u，v の運動方程式の移流項を計算．3次精度の上流差分法を使っている．時間差分 f に加算するようになっているので calc_momentum_equations_first の後に書くこと． |
| calc_momentum_equations_add_viscosity | u，v の運動方程式の水平粘性項を計算．時間差分 f に加算するようになっているので calc_momentum_equations_first の後に書くこと． |
| calc_momentum_equations_add_friction | u，v の運動方程式の摩擦項を計算．時間差分 f に加算するようになっているので calc_momentum_equations_first の後に書くこと． |
| calc_continuity_equation | h の連続の式を計算． |
| calc_matsuno | 松野スキームの式 1本分の計算（1計算ステップに 2回呼び出す）． |
| calc_leap_frog | リープフロッグの計算． |
| output | 計算結果を出力する．フォーマットは，タブ区切りで，経度，緯度，u，v，h の順に出力している． |
| main | main関数．メモリの確保と開放，各種セッティング，時間積分ループの実行を行っている． |
set_* と calc_* では，最初の引数に計算結果を代入あるいは加算している．
main の変数のうち，主なものを解説する．
| long itime | 計算ステップ数． |
| struct cell *c | グリッド情報． |
| struct var *v_pre | 1ステップ前の状態（un-1，vn-1，hn-1） |
| struct var *v_now | 現在のステップの状態（un，vn，hn）． |
| struct var *v_next | 次のステップの状態（un+1，vn+1，hn+1）． |
| struct var *v_matsuno | 松野スキームで使う仮の予測値（u*，v*，h*） |
| struct var *f | 時間差分値（du/dt，dv/dt，dh/dt）． |
方程式のコーディングでは，グリッドを相対位置で指定するので， グリッドの座標を意識する必要はないが， 各変数の計算範囲を決めたり，境界条件処理をコーディングする際には， 境界付近のグリッドの座標の値を把握している必要がある．
以下に，(ilon, ilat)座標空間での境界と計算領域を図示する． 水色の実線が境界で，この境界上にある u, v の値は 0 である． 水色の網掛け部分が計算領域になる． 境界に隣接する変数の計算には， 計算領域外の赤色網掛け部分の変数も使う必要があり， この変数の値は，関数 calc_boundary_condition で計算する．
SLTM の出力は，タブ区切りのアスキーファイルで， 1行にλ，φ，u，v，h の値が順に並んでいる．
1層（1.5層）の SLTM では， 境界面の変位 h の等値線が流線とほぼ一致するので， h の分布がわかれば，流れの様子がわかる． ここでは h の分布をプロットする方法を 2つ紹介しよう．
gnuplot を起動し，
> splot "0120.dat" using 1:2:5
とすると，出力ファイル 0120.dat の h の分布を見ることができる．
たとえば，以下のようなシェルスクリプトを書く．
#!/bin/sh
FILE=0120.dat
HEIGHT=100
EPS=${FILE%.dat}.eps
TMP=tmp.$$.dat
GRD=tmp.$$.grd
CPT=tmp.$$.cpt
makecpt -Cpolar -T-${HEIGHT}/${HEIGHT}/${HEIGHT} -Z > ${CPT}
awk '{print $1, $2, $5}' ${FILE} > ${TMP}
xyz2grd `minmax -I1 ${FILE}` -I0.2 -G${GRD} ${TMP}
grdimage ${GRD} -Jm1 -C${CPT} -Ba1:."${FILE}": -P -K > ${EPS}
grdcontour ${GRD} -J -C10 -A50 -P -O >> ${EPS}
rm ${TMP} ${GRD} ${CPT}
このスクリプトを実行すると，出力ファイル 0120.dat から， EPSファイル（GMT の設定によっては PSファイル） 0120.eps ができる．
xyz2grd の -Iオプションの数値は，グリッド間隔の値に応じて， 適宜変更すること（この例では 0.2度）．
FILE，HEIGHT の値は直接書かずに，$1，$2 としておいて， スクリプトへの引数で与えるようにすると便利である．
日本海を模した中緯度の 10x10度の海域で， 重力波のシミュレーションを行う． 主なパラメータの値は，
重力波の位相速度 ( gH0 )1/2 が 140ms-1 と速いので，Δt は 1s にセットする．
set_initial_condition で make_mount を使い， 海域の中心に高さ 10m，半径 0.2度（10グリッド）の初期波形を与えた．
シミュレーションを実行すると， 重力波が広がって，岸で反射する様子を見ることができる．
海面変位分布（赤：正，青：負）．重力波が東西岸で反射したところ（経過時間 60分）．
120分の計算をするのに，AMD Phenom X4 9350e マシンで約90分かかる．
日本海を模した中緯度の 10x10度の海域で， ロスビー波による中規模渦の西進のシミュレーションを行う． 主なパラメータの値は，
表層の 1.5層モデルで，g は reduced gravity，h は内部の境界面の変位（下向き正）になる．
set_initial_condition で make_mount を使い， 海域の中心に変位 100m，半径 2度（20グリッド）の初期波形を与えた．
シミュレーションを実行すると， 地衡流平衡した高気圧性渦が，ゆっくりと西進する様子を見ることができる．
境界面変位分布（m）．海域の中心に与えた渦が，西進しているところ（経過時間 100日）．
120日分の計算をするのに，AMD Phenom X4 9350e マシンでは 8分しかかからない．
太平洋の西半分を模した， 北緯 10〜50度，東経 120〜180度の海域で， 風成循環のシミュレーションを行う． 本当は，太平洋全域で行うつもりだったのだが， 広すぎて計算に時間がかかりすぎるのであきらめた．
主なパラメータの値は，
表層の 1.5層モデルで，g は reduced gravity，h は内部の境界面の変位（下向き正）になる．
set_ekman_upwelling で，エクマン湧昇をセットする．
低緯度と高緯度で西向き，中緯度で東向きの風応力 τ
を与えると，エクマン湧昇 wE は，
となる（で，合ってるよね？）．今回は，τ0 = 0.1Nm-2 とした．
シミュレーションを実行すると， 時計回りの亜熱帯循環と，反時計回りの亜寒帯循環が形成される様子を見ることができる．
境界面変位分布（m）．亜熱帯循環と亜寒帯循環が形成される（経過時間 20年）．
20年分の計算をするのに，AMD Phenom X4 9350e マシンで約25時間かかる．
※ 本当は，もっと格子間隔を小さくした方がよいのだが， この範囲をパソコンで計算するには 0.2度間隔が精一杯である （0.1度間隔にすると計算に 9日かかる）．
太平洋の西半分を模した， 北緯 10〜50度，東経 120〜180度の海域で， 古典的な深層循環のシミュレーションを行う． 本当は，太平洋全域で行うつもりだったのだが， 広すぎて計算に時間がかかりすぎるのであきらめた．
主なパラメータの値は，
深層の 1.5層モデルで，g は reduced gravity，h は内部の境界面の変位（上向き正）になる．
calc_boundary_condition で， 海域の南西端に，20Sv の流入に相当する 0.2ms-1 の北向きの流速を与える．
シミュレーションを実行すると，西岸と北岸に強い流れが生じ， 他のほとんどの領域では，ゆっくりと北東向きに流れている様子を見ることができる．
境界面変位分布（m）．ほとんどの領域でゆるやかに北向きに流れる（経過時間 40年）．
40年分の計算をするのに，AMD Phenom X4 9350e マシンで約8時間かかる．
せっかく短く，シンプルにコーディングしたのだから， どんどん改造して遊んでもらいたい．
最初はパラメータをいじって遊んでみよう．
次は，初期値をいじってみるのかなあ．
方程式もいじってみよう．まずは手始めに，
さらに，
もっと高度な改造は，
こんな簡単なモデルでも，遊び方はいくらでもある．
Enjoy Happy Hacking!
コードにまちがいがみつかりましたので修正しました． 移流項の上流差分の判定に使う流速が，半グリッドずれていました． 修正箇所は，以下の一行です．
/* u dv/dx */
u_at_v = (v[W].u + v[C].u + v[WN].u + v[N].u) * 0.25;
誤 if(ilon == ILON_NUM - 2 || (ilon != 1 && v[C].u > 0.0)){
正 if(ilon == ILON_NUM - 2 || (ilon != 1 && u_at_v > 0.0)){
f[C].v += - u_at_v
* (2.0 * v[E].v + 3.0 * v[C].v - 6.0 * v[W].v + v[WW].v)
/ 6.0 / DLON_RAD / EARTH_RADIUS / cos(LAT_RAD + DLAT_RAD * 0.5);
ご迷惑をおかけして申しわけありません．
コードに大量に混入していた「タブ」を除去しました． 動作への影響はありません．
境界付近の処理を大きく変更しました．
境界の向こう側の変数は境界に隣接するもののみ計算していた （参考：変更前のグリッド配置）． この方法だと，境界に隣接する u，v の移流項の計算については， 陸地側の 2つ隣りの変数が足りないので，3次の上流差分法をそのまま使えない． そこで，陸地側が上流になることはないと仮定して計算していた．
計算領域を東西南北 1セルづつせばめ， 境界の向こう側の値を，2セル分計算するようにし，境界に隣接する u，v の特別扱いをやめた．
実行例 2，3，4 について，水平粘性係数 Ah の値が不適当だった（小さすぎた）のを修正しました． 計算時間を短縮するために格子間隔を大きくしたとき，この係数をなおすのをわすれていました． おはずかしいかぎりです．
先日変更した境界付近の処理に，不自然な部分があったので修正しました．
2つの関数名をわかりやすいものに変更しました．動作の変更はありません．
境界条件の処理を修正しました．
境界条件の処理を修正しました． 境界付近の計算に，メモリの初期化していない部分を使っていました． コードを整理したとき，処理をけずりすぎていました．
摩擦項がなぜか中心差分だったのを，前進差分に修正しました．
あわせて，パラメタ名を 2つ変更しました．
[ KDOホームにもどる ]
(C) 2010, 2011 KDO
