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

熱帯太平洋(1980-2000:3ヶ月移動平均) 熱帯太平洋(2000-2018:3ヶ月移動平均)

Ertel渦位(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を行ってください。

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

dgwspeed-sample.dta
dertelpv.YYYYMMDD.dta

この中身を確認するにはgradsを起動してcgwspeed-sample.ctlとceterlpv-sampleをそれぞれopenしてください。1つ目のファイルは、表層の重力波の速度の分布(海洋の場合3m/s程度で水平一様でもかまいません)を与え、それと同時に陸海の分布(値がマイナスになっているところを陸や海岸線直上とみなします)を与える役割があります。サンプルデータの場合は下図のIS0×JS0の範囲が720×241グリッドになっています。この図では紫塗四角のところで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)は大気海洋力学の世界ではストレッチング項と呼ばれる項で、変形半径の逆数の2乗が係数としてついています。その値代入は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
のようにして実行します。

今後の計画

Level-0用のインバージョンプログラムのデバッグ、共同研究の体制整備(科研費申請など)