Group-velocity-based energy flux at all latitudes in the ocean and atmosphere

海洋表層のエネルギーフラックス解析結果

Ertel渦位（EPV）のインバージョンプログラム

大気や海洋中の大規模な波動によるエネルギーの伝達経路を気候学的に同定する事によって、熱帯―亜熱帯相互作用（エルニーニョ現象やインド洋ダイポールモード現象が日本付近の気候に及ぼす影響など）を診断することを目的としています。　
Aiki et al. (2017 PEPS)は理論研究によって、EPVのインバージョンから得られる流線関数を使用すれば、浅水方程式で表される全ての種類の波についての、群速度ベクトル分布を全ての緯度帯において連続して計算することができる事を示しました。これは（全ての種類の波についての）オートフォーカス機能と（全ての移動帯についての）シームレス機能の両方を備えた新しい解析手法です。つまり赤道域と中緯度域の力学を融合し、統一的に赤道波・中緯度惑星波・慣性重力波を取り扱うことができます。ポスター
本ホームページでは第一弾としてLevel-2エネルギーフラックス計算用のインバージョンプログラムのソースコード（v20170719）を配布して、AGCM/OGCM結果の試験的な診断を開始してもらっています。
これによって、(i) どのくらいの時間スケールの平均が適切なのか？ (ii) 診断用の１層を定義するにはどの高度帯（AGCMの場合）および深さ帯（OGCMの場合）を選べばよいのか？ (iii) Level-0エネルギーフラックス用のEPVインバージョンプログラムの開発はどのくらい求められているのか？についての情報交換を行い、新しい研究成果を創出していきたいと思います。

クイックスタート

まずインバージョンプログラムをコンパイルするにはunixのプロンプトから
make clean all
と打ちます。次にサンプルデータ（EPVの水平分布・重力波速度の水平分布）を作成します。このためにはgradsを起動して
exec gmksample.gex
とします。そしてgradsを終了してunixのプロンプトに戻ります。最後にインバージョンプログラムを起動するには
./invc2d dertelpv.19990101.dta
のようにします。結果を見るにはgradsを起動して
open cstrmfc-sample.ctl
d stf
とします。

全体の説明

グリッド数や解析範囲などの情報はhclassic2d_mysetup.hにて設定します。配布しているソースファイルではとりあえず東西方向720グリッド南北方向241グリッド、東西境界条件は周期的（PERIODIC）になっています。太平洋などに限定した解析では東西境界条件を閉鎖的（CLOSED）にすることもできます。hclassic2d_mysetup.hを書き換えた後は必ずmake allを行ってください。

入力データは以下の２つのバイナリファイルです。

dgwspeed-sample.dta
dertelpv.YYYYMMDD.dta

この中身を確認するにはgradsを起動してcgwspeed-sample.ctlとceterlpv-sampleをそれぞれopenしてください。１つ目のファイルは、表層の重力波の速度の分布（海洋の場合3m/s程度で水平一様でもかまいません）を与え、それと同時に陸海の分布（値がマイナスになっているところを陸や海岸線直上とみなします）を与える役割があります。サンプルデータの場合は下図のIS0×JS0の範囲が720×241グリッドになっています。この図では紫塗四角のところでEPVが定義されていることを表しています。例えば

海洋大循環モデルMRI.COMの出力（Arakawa's B-grid）を使ってEPVファイルを作る時には、温度と塩分点が下図の紫四角に対応し、水平速度点が下図の青丸に対応するという風に解釈します。
海洋大循環モデルPOMの出力（Arakawa's C-grid）を使ってEPVファイルを作る時には、温度と塩分点が下図の青丸に対応し、相対渦度点が下図の紫四角に対応するという風に解釈します。

いずれにせよ入力ファイル（dertelpv.YYYYMMDD.dta）には、EPVアノマリーの値を単位が相対渦度と同じになるように規格化して、下図のIS0×JS0の範囲で入れておきます。

共役勾配法の説明

インバージョンの式の左辺は楕円型の微分方程式になっているので、BICGSTAB2法を使って解いています。楕円型の微分方程式は差分化してfclassic2d_cgmatrix.F90の

do j = JJS0
do i = IIS0
s(i,j) = ((flx(i,j)-flx(i-1,j)+fly(i,j)-fly(i,j-1))*cidxy_o(i,j)+elc(i,j )*w(i,j))*lm_o(i,j)
end do
end do

の部分に記述してあります。elc(i,j)は大気海洋力学の世界ではストレッチング項と呼ばれる項で、変形半径の逆数の２乗が係数としてついています。その値代入はfclassic2d_main.F90の

do j = JJS0
do i = IIS0
if (gwspd(i,j) > 0.0d0) then
lm_o(i,j) = 1.0d0
elc_cof(i,j) = -bgf_o(i,j)*bgf_o(i,j)/(gwspd(i,j)*gwspd(i,j))
else
lm_o(i,j) = 0.0d0
elc_cof(i,j) = 0.0d0
endif
enddo
enddo

の部分でしています。BICGSTAB法は反復解法ですので、元々の外力と比較してインバージョンの式の誤差がEPS_2CG倍以下に収束したら、反復を打ち切るようになっています。その設定は、fclassic2d_cgmain.F90の

#define EPS_2CG (1.0d-10) /* relative evaluation */

の部分でしています。

Level-2用の流線関数はfclassic2d_*.F90、Level-0用の流線関数はfexact2d_*F90

Level-2のエネルギーフラック計算用の流線関数は古典的（classic）なEPVインバージョンの式を解いて求めます。その実行プログラムはinvc2dで並列計算機でなくても動きます。一方、Level-0のエネルギーフラックス計算用の流線関数は新しい（exact）EPVインバージョンの式を解いて求められています。この式には2階の時間微分項が含まれます。この式を解く実行プログラムはinve2dで、水平方向に2次元、時間方向に1次元、合計3次元の時空間でメモリーを確保し、流線関数の分布をBICGSTAB2法で求めます。前処理には時間方向のLU分解に基づく陰解法を使用しています。inve2dは使用メモリー量負荷が高いので並列計算機を使って動かすように設計されています。実行するには、epvファイル群のリストを先に作っておきます。そのためにはunixのプロンプトから
ls dertel*.dta >! list
と打ちます。その後
mpirun -n 40 ./inve2d list
のようにして実行します。

今後の計画

３次元大気用のインバージョンプログラムの配布、共同研究の体制整備（科研費申請など）
今後の定式化の方針