SLTM X1 海洋一層数値モデル

2012年 8月 14日掲載，2012年 9月 25日修正

はじめに

”パソコン”で気軽にためして遊べるような海洋数値モデルを作っている．

前作 SLTM（Single Layer Toy Model）は，簡素で改造しやすいのが売りであった．この SLTM をベースにして，実用的な数値モデルに必要な機能を実験的に追加していくのが SLTM Xシリーズである．

Xシリーズ第1弾の SLTM X1 では，おそらくもっとも需要が多いと思われる，陸地をあつかう機能を追加してみた．

SLTM X1 について

SLTM をベースに開発した，陸地の影響をあつかえる 1層（1.5層）モデルである．

開発の方針としては，やはり改造のしやすさを優先することとし，地形に関する部分と力学計算部分の分離を追求した．具体的には，地形をあつかうコード map.c と， SLTM由来の力学計算コード equations.c を分離してある．

今回掲載した SLTM X1（Single Layer Toy Model eXperimental 1）は， SLTM に地形をあつかう機能を追加したかわりに，

移流項
完全解説

をあきらめた．

SLTM X1 では，陸岸境界の向こう側の仮想セルを， 1列しか用意していないので，3次の上流差分が使えず，かといって今さら 1次の上流差分や中央差分を使う気にもならなかったので，省略することにした．そのため，移流の影響が無視できない現象には使えなくなった．

SLTM X1 のコード量は，1170行と SLTM の約2倍になった．使い方と，変数の意味や関数が何をしているかくらいは説明したが，地形をあつかう部分は，かなりオリジナルな設計になっているので，解説するのがむずかしい．

ダウンロード

ソースコードとサンプルの設定ファイルをまとめて圧縮したファイルを用意した．

sltmx1.tgz

ダウンロードして適当なディレクトリで展開してください．

ファイルの改行コードは LF，コードのコメントの文字コードは UTF-8 です．

使用条件・配布条件

以下の免責事項に同意していただくことを，使用条件とします．

免責事項：
SLTM X1 を使うことにより損害が発生しても，作成者である KDO は責任を負いません．

SLTM X1 は，どんどん改造してもらうことや，海洋物理の学習のために使っていただくことを想定しているので，このままでも，改造したものでも，自由に配布していただいて結構です．ただし，コード冒頭の著作権のコメント

/*
 * SLTM X1 - Single Layer Toy Model X1 (C) 2012 KDO
 */

は残していただくよう，お願いします．改造して配布される方は，その旨コメントを追加していってください．

使い方（コンパイル，引数，出力）

コードは，次の 4つのファイルである．

header.h ：ヘッダファイル
main.c ： main関数
map.c ：地形をあつかう関数群
equations.c ： SLTM 由来の力学計算の関数群

4つのファイルを一つのディレクトリに置いて， -lm オプションを付けてコンパイルする．

$ cc -o sltmx1 main.c map.c equations.c -lm

これは，配布する Makefile にも記述されている．

実行時には，2つのファイルを指定する必要がある．必要なのは，

第1引数：計算領域の範囲と格子間隔を記した「range file」
第2引数：層厚分布の設定の「layer file」

「range file」には，以下の情報を記しておく．

1行目：計算領域の西端の経度，単位は度
2行目：計算領域の東端の経度，単位は度
3行目：東西方向の格子間隔，単位は分
4行目：計算領域の南端の緯度，単位は度
5行目：計算領域の北端の緯度，単位は度
6行目：南北方向の格子間隔，単位は分

「layer file」の層厚の設定は，1行に，

経度，単位は度
緯度，単位は度
層厚，単位は m

を，タブで区切って書く．層厚を設定していないセルが，陸地としてあつかわれる．「range file」の範囲外の層厚設定は無視される．

配布するサンプルを使った実行例を以下に示す．

$ ./sltmx1 range.txt layer.txt
0001
0002
0003
 :
 ;

計算結果のフォーマットは，SLTM と同じで，1行に

経度 λ，単位は度
緯度 φ，単位は度
流速の東方成分 u，単位は ms^-1
流速の北方成分 v，単位は ms^-1
境界面変位 h，単位は m

が，タブ区切りでならんでいる．

コードの解説

SLTM のような詳細な解説はできないが，書けるだけのことは書いてみよう．

header.h

header.h には，計算パラメタを設定する #define文，構造体と関数のプロトタイプを記述している．

パラメタは，以下の 13個．

OUTPUT_NUM	出力ファイルの個数
OUTPUT_STEP	何ステップごとに出力するか
OUTPUT_FILENAME_LENGTH	出力ファイル名の長さ
OUTPUT_FILENAME_FORMAT	出力ファイル名のフォーマット
DT_SEC	1計算ステップの時間間隔，単位は s
MATSUNO_STEP	何ステップごとに松野スキームを使うか
PI	円周率 π
GRAVITY	重力加速度 g，単位は ms^-2 （1.5層モデルでは reduced gravity）
EARTH_ROTATION	地球の自転角速度 Ω，単位は s^-1
EARTH_RADIUS	地球の半径 a，単位は m
RAYLEIGH_FRICTION	摩擦係数 R，単位は s^-1
NEWTONIAN_COOLING	鉛直拡散係数 γ，単位は s^-1
HORIZONTAL_VISCOSITY	水平粘性係数 A_h，単位は m²s^-1

構造体は，以下のものがある．

var	物理変数 u，v，h
cell	隣接するセルのセルID と，経緯度，層厚
map	計算領域の範囲と，セルID のマップ
param	力学計算ルーチンのパラメタ
tmp_cell	一時セル，隣接するセルの海陸フラグ
boundary	陸地境界の情報
boundary_list	陸地に隣接するセルと，境界の向こう側の仮想セルのセルID のペア

説明中に出てくる，マップとセルID については，次節で解説する．

main.c

main関数．地形のセッティング，時間発展の計算，メモリの確保・解放を行っている．

map.c

今回追加された地形をあつかう機能．以下の 19個の関数からなる．

create_map	range file を読み込み，map に計算範囲と格子間隔をセット．セルID のマップを作成して初期化する．
delete_map	map の消去．
create_param	力学計算ルーチンのパラメタ param の作成・初期化．
delete_param	param の消去．
create_tmp_cell	一時セルの作成・初期化．
delete_tmp_cell	一時セルの消去．
read_layer_thickness	layer file を読み込み，一時セルに層厚をセット．
set_tmp_cell_flag	一時セルに，隣接するセルの海陸フラグをセット．
set_map_id_and_param	一時セルのフラグを参考にして，u，v，h それぞれの計算セルが連続するように，セルID をマップに割り付ける．u，v，h の計算セルID の始点と終点を param にセット．
create_cell	計算セルと陸地境界の向こうの仮想セルの分の cell を作成・初期化．
delete_cell	cell の消去．
create_boundary	陸地境界処理が必要な u，v，h をカウントし，boundary_list を作成・初期化．
create_boundary_list	boundary_list を作成・初期化する関数．
delete_boundary_list	boundary_list の消去．
delete_boundary	boundary の消去．
set_cell_param	各セルに経緯度，層厚をセット．
set_ewns_id	各計算セルについて，東西南北のセルが海であれば，そのセルID をセット．
set_boundary	計算セルの東西南北のセルが陸地であれば，仮想セルを割り当て，計算セルと仮想セルのセルID のペアのリストを作成する．
set_es_wn_id	南東と北西のセルID をセット．

equations.c

SLTM 由来の力学計算関数群である．ここではくわしく解説しないので，SLTM の解説を読んでほしい． SLTM との主な違いは，

var の配列を作成・消去する create_var，delete_var を追加した．
calc_momentam_equations_add_advection がなくなった．
セルの経緯度情報，隣接セルの情報が，座標からの計算式でなく，cell のメンバになった．
力学計算のループが 1回になった．

SLTM で検証のすんでいる計算式には手をつけていないので安心である．

セルID・マップ

SLTM X1 では，計算セルに「セルID」という番号をふっている．各セルごとに隣接するセルのセルID を，構造体 cell のメンバとして持っている．この仕組みにより，陸岸境界の向こう側にある仮想セルを干渉を気にせず設定できるようになった（図）．

仮想セル

また，u，v，h を計算するセルがそれぞれ連続するように並べることにより，力学計算ループがシンプルになっている（速いかどうかはわからない）．

各セルの経緯度は，構造体 cell のメンバに経緯度がふくまれているので，簡単に取得できる．経緯度からセルID を取得するためには，構造体 map を使う．構造体 map のメンバを使って，2次元の経緯度座標から，対応するセルIDを求めることができる．

struct map *m;

ilon = (long)((lon_deg - m->lon_west_deg) * 60.0 / m->dlon_min + 0.5);
ilat = (long)((lat_deg - m->lat_south_deg) * 60.0 / m->dlat_min + 0.5);
igrid = m->ilon_num * ilat + ilon;
cell_id = m->id[igrid];

SLTM X1 の方程式

運動方程式 2つと連続の式 1つの 3つの式から，u，v，h の時間差分を計算する．

du/dt
dv/dt
dh/dt

ここで，

各変数は，

u：経度方向の流速（東向き正）
v：緯度方向の流速（北向き正）
h：境界面変位（層厚が厚くなる向きが正）
t：時刻
λ：経度
φ：緯度
H₀：層厚
a：地球の半径
Ω：自転角速度
g：重力加速度（1.5層モデルでは，reduced gravity）
R：摩擦係数
A_h：水平粘性係数
γ：鉛直拡散係数

実行例

配布しているファイルを解凍，コンパイルして実行すると，中緯度の5度×5度の範囲にある，円形の海での重力波の伝播のシミュレーションがはじまる．格子間隔は 0.02度（1.2分），層厚は 2000m．

$ cc -o sltmx1 main.c map.c equations.c -lm
$ ./sltmx1 range.txt layer.txt

海面変位分布（赤：正，青：負）．陸岸で反射する重力波（経過時間 120分）

重力波が円形に広がり，陸岸に反射してもどってくる様子を確認できる．

120分の計算をするのに，AMD Phenom X4 マシンで 8分かかる．

改造してみよう！

計算範囲と地形については外部ファイル化してあるが，初期条件はコード内で作っているので，最初にいじるのは equations.c の set_initial_condition になるだろう．

いろいろパラメタを変えて実験してみたいときは， header.h の #define文を書き換えてみることになるのかな．

次にいじるとしたら equations.c の力学計算関数だろう．摩擦項をもうちょっとマシなのに交換するとか，いろいろ遊べると思う．

map.c はちょっと手がつけられないかな．ここの改造のアイディアは，Xシリーズでためしていくつもりです．

更新履歴

2012年 8月 14日

公開．

2012年 9月 25日

このページのスペルミスを修正しました．

equation.c → equations.c
HORIZONTAL_VISCOUSITY → HORIZONTAL_VISCOSITY

[ KDOホームにもどる ]