著名な海洋物理モデル Princeton Ocean Model (POM) を使っていて、1ヶ月以上にわたり悩まされた問題について、メモ代わりに残しておきたい。予め書いておくと、自分は初期値として水平に密度勾配があるという条件について数値実験していた。そしてこのことが問題となっていた。

この海洋モデルでは座標系のうち鉛直方向を、一般的な直交座標系ではなくシグマ座標系というのを採用している。シグマ座標についてはほかにいくらでも書かれているので、ここでは詳細に記載しない。メリットもあり、デメリットもある。

デメリットの一つに、Pressure Gradient Error という誤差がある。海洋において圧力とは、主に水の密度である。この密度が水平に異なるとき、高密度の方から低密度の方へ流速が起きるような力が生まれる。これを圧力傾度力と呼ぶ。シグマ座標系において圧力傾度力の計算で起きる密度の打ち切り誤差は、圧力傾度力に対して過大に影響する。これが Pressure Gradient Error で、シグマ座標系の問題として長く議論されている。

POM ではこの誤差を軽減するため、圧力傾度力を計算する前に密度から密度の平均値を一旦減算 (rho = rho - rmean) し、圧力傾度力の計算が終わると再び加算 (rho = rho + rmean) する。こうすることで、密度の絶対値ではなく相対値によって圧力傾度力が計算され、打ち切り誤差が軽減される。ここで密度の変数は rho、密度の平均は rmean である。

問題は、この rmean の初期化にある。初期化する部分のデフォルトのコードは以下のようになっている。

見ての通り、rmean = rho であって、平均でも何でも無い。コメント部分を読むと、『密度の初期値が水深のみの関数であるときはこれで大丈夫。水平に密度勾配があるときは平均するように。』というようなことが書かれている。これを見逃していたために、圧力傾度力が正しく計算されていなかった。これが問題のすべてである。

解決するためには、rmean に適切な平均値を入れるか、またはいっそ 0 を代入しておいても良い。Pressure Gradient Error は大きくなるものの、大した問題ではない。

ということで、POM を使って初期値に水平の密度勾配がある条件の数値実験をする際には、rmean を適切に初期化しなければならない、という教訓です。1ヶ月くらい気付かなかったし、すごく辛い思いをした。普通に、平均という名前の変数には、最初から平均を入れておいてほしい。