!-- !---------------------------------------------------------------------- ! Copyright(c) 2008-2012 SPMDODEL Development Group. All rights reserved. !---------------------------------------------------------------------- !表題 wq_mpi_module ! ! spml/wq_mpi_module モジュールは球内での流体運動をスペクトル法と ! MPI 並列化によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向に球面調和函数変換および動径方向に ! Matsushima and Marcus (1994) で提唱された多項式を用いた ! スペクトル計算のためのさまざまな関数を提供する. ! ! 内部で wa_mpi_module, aq_module を用いている. ! 最下部では球面調和変換変換のエンジンとして ! ISPACK の Fortran77 サブルーチンを用いている. ! ! Matsushima and Marcus (1994) の多項式に関する説明は ! doc/spectral_radial.tex を参照のこと. ! ! !履歴 2008/04/03 竹広真一 wu_module より改変 ! 2008/05/02 竹広真一 コメント追加 ! 2008/07/01 佐々木洋平 水平方向の格子点を 0:im-1,1:jm に修正 ! 2008/07/04 佐々木洋平 コメントを RDoc 用に微修正 ! 2008/07/20 竹広真一 wq_Rad2_wq, wq_Rad2Inv_wq, wq_Lapla_wq 追加 ! 2008/07/21 竹広真一 行列設定計算変更 ! 2009/01/09 竹広真一 wu_Initial メッセージに日付を追加 ! 2009/01/29 佐々木洋平 コメントを RDoc 用に微修正 ! 2009/07/31 竹広真一 境界条件計算用配列を threadprivate 指定(OpenMP) ! 2009/12/06 竹広真一 先頭コメントを修正 ! 2010/03/10 佐々木洋平 threadprivate 削除(コンパイラ依存) ! 2010/04/17 竹広真一 wr_DivRad_xyr 追加, wr_Div_xyr_xyr_xyr を変更 ! wq_LaplaPol2Pol_wq の計算法変更 ! 2010/05/08 竹広真一 wq_LaplaPol2Pol_wq の計算法元にもどす. ! 2011/09/08 竹広真一 MPI 並列化 ! 2012/05/25 竹広真一 wq_QoperatorMPI_wq を追加 ! !凡例 ! データ種類と index ! x : 経度 y : 緯度 v : 緯度(分散格子) r : 動径 ! w : 球面調和関数スペクトル ! n : 球面調和関数スペクトル(水平全波数) ! m : 球面調和関数スペクトル(帯状波数) ! q : スペクトル関数スペクトル ! a : 任意の次元 ! ! xyr : 3 次元格子点データ ! xvr : 3 次元分散格子点データ ! xy : 水平 2 次元格子点データ ! xv : 水平 2 次元分散格子点データ ! yr : 子午面 2 次元格子点データ ! vr : 子午面 2 次元分散格子点データ ! xr : 緯度面 2 次元格子点データ ! ! wr : 水平スペクトル動径格子点データ ! wq : スペクトルデータ ! !++ module wq_mpi_module ! != wq_mpi_module ! ! Authors:: Shin-ichi Takehiro, Youhei SASAKI ! Version:: $Id$ ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT] ! !== 概要 ! ! spml/wq_mpi_module モジュールは球内での流体運動をスペクトル法と ! MPI 並列化によって数値計算するための Fortran90 関数を提供する ! ものである. ! ! 水平方向に球面調和函数変換および動径方向に ! Matsushima and Marcus (1994) で提唱された多項式を用いた ! スペクトル計算のためのさまざまな関数を提供する. ! ! 内部で wa_mpi_module, aq_module を用いている. ! 最下部では球面調和変換のエンジンとして ! ISPACK の Fortran77 サブルーチンを用いている. ! !== 関数・変数の名前と型について ! !=== 命名法 ! ! * 関数名の先頭 (wq_, nmr_, nr_, xyr_, wr_, w_, xy_, xv_, yr_, vr_, xr_, ! x_, y_, v_, r_, a_) は, 返す値の形を示している. ! wq_ :: スペクトルデータ(球面調和函数・チェビシェフ変換) ! nmr_ :: 水平スペクトルデータ(全波数 n, 帯状波数各成分, 動径) ! nr_ :: 水平スペクトルデータ(全波数 n, 動径) ! xyr_ :: 3 次元格子点データ(経度・緯度・動径) ! xvr_ :: 3 次元分散格子点データ(経度・緯度・動径) ! wr_ :: 水平スペクトル, 動径格子点データ ! xy_ :: 2 次元格子点データ(経度・緯度) ! xv_ :: 2 次元分散格子点データ(経度・緯度) ! yr_ :: 2 次元格子点データ(緯度・動径) ! vr_ :: 2 次元分散格子点データ(緯度・動径) ! xr_ :: 2 次元格子点データ(経度・動径) ! x_ :: 1 次元格子点データ(経度) ! y_ :: 1 次元格子点データ(緯度) ! v_ :: 1 次元分散格子点データ(緯度) ! r_ :: 1 次元格子点データ(動径) ! a_ :: 任意の次元 ! ! * 関数名の間の文字列(DLon, GradLat, GradLat, DivLon, DivLat, Lapla,..) ! は, その関数の作用を表している. ! ! * 関数名の最後 (_wq, _nmr, _nr, _xyr, _wr, _w, _xy, _xv, _yr, _vr, _xr, ! _x, _y, _v, _r, _a) は, 入力変数の形がスペクトルデータおよび ! 格子点データであることを示している. ! _wq :: スペクトルデータ ! _xyr :: 3 次元格子点データ ! _xyr_xyr :: 2 つの3 次元格子点データ, ... ! _xvr :: 3 次元分散格子点データ ! _xvr_xvr :: 2 つの3 次元分散格子点データ, ... ! !=== 各データの種類の説明 ! ! * xyr : 3 次元格子点データ(経度・緯度・動径) ! * 変数の種類と次元は real(8), dimension(0:im-1,1:jm,km). ! * im, jm, km はそれぞれ経度, 緯度, 動径座標の格子点数であり, ! サブルーチン wq_mpi_Initial にてあらかじめ設定しておく. ! ! * xvr : 3 次元分散格子点データ(経度・緯度・動径) ! * 変数の種類と次元は real(8), dimension(0:im-1,jc,km). ! * im, km はそれぞれ経度, 動径座標の格子点数であり, ! * jc はこのプロセスで保有する緯度格子点数である. ! サブルーチン wq_mpi_Initial を呼ぶと jc が設定される. ! ! * wq : スペクトルデータ ! * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1),0:lm). ! * nm は球面調和函数の最大全波数, lm はチェビシェフ多項式の最大次数 ! であり, サブルーチン wq_mpi_Initial にてあらかじめ設定しておく. ! * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる ! ことができる. ! * 動径スペクトルデータの格納方法については aq_module.f90 を ! 参照のこと. km < 2*im でなければならない. ! ! * nmr : 水平スペクトルデータの並んだ 3 次元配列. ! * 変数の種類と次元は real(8), dimension(0:nm,-nm:nm,km). ! * 第 1 次元が水平全波数, 第 2 次元が帯状波数, 第 3 次元が動径座標を表す. ! * nm は球面調和函数の最大全波数であり, サブルーチン wq_Initial にて ! あらかじめ設定しておく. ! ! * nr : スペクトルデータの並んだ 2 次元配列. ! * 変数の種類と次元は real(8), dimension(0:nm,km). ! * 第 1 次元が水平全波数を表す. nm は球面調和函数の最大全波数であり, ! サブルーチン wq_Initial にてあらかじめ設定しておく. ! ! * wr : 水平スペクトル, 動径格子点データ. ! * 変数の種類と次元は real(8), dimension((nm+1)*(nm+1),km). ! ! * wq_ で始まる関数が返す値はスペクトルデータに同じ. ! ! * xyr_ で始まる関数が返す値は 3 次元格子点データに同じ. ! ! * xvr_ で始まる関数が返す値は 3 次元分散格子点データに同じ. ! ! * wr_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ. ! ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに ! 微分などを作用させたデータをスペクトル変換したものことである. ! ! !== 変数・手続き群の要約 ! !==== 初期化 ! ! wq_mpi_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定 ! !==== 座標変数 ! ! x_Lon, x_Lon_Weight, :: 格子点座標(経度)と重みを格納した 1 次元配列 ! y_Lat, y_Lat_Weight :: 格子点座標(緯度)と重みを格納した 1 次元配列 ! v_Lat, v_Lat_Weight :: 分散格子点座標(緯度)と重みを格納した 1 次元配列 ! r_Rad, r_Rad_Weight :: 格子点座標(同形)と重みを格納した 1 次元配列 ! xyr_Lon, xyr_Lat, xyr_Rad :: 格子点データの経度・緯度・動径座標(X,Y,Z) (格子点データ型 3 次元配列) ! xvr_Lon, xvr_Lat, xvr_Rad :: 格子点データの経度・緯度・動径座標(X,Y,Z) (分散格子点データ型 3 次元配列) ! !==== 基本変換 ! ! xyr_wq, wq_xyr :: スペクトルデータと 3 次元格子データの間の変換 ! (球面調和函数, チェビシェフ変換) ! ! xvr_wq, wq_xvr :: スペクトルデータと 3 次元分散格子データの間の変換 ! (球面調和函数, チェビシェフ変換) ! ! xyr_wr, wr_xyr :: 3 次元格子データと水平スペクトル・動径格子データとの ! 間の変換 (球面調和函数) ! ! xvr_wr, wr_xvr :: 3 次元分散格子データと水平スペクトル・動径格子データとの ! 間の変換 (球面調和函数) ! ! wr_wq, wq_wr :: スペクトルデータと水平スペクトル・動径格子データとの ! 間の変換 (チェビシェフ変換) ! ! w_xy, xy_w :: スペクトルデータと 2 次元水平格子データの ! 間の変換(球面調和函数変換) ! ! w_xv, xv_w :: スペクトルデータと 2 次元水平分散格子データの ! 間の変換(球面調和函数変換) ! ! l_nm, nm_l :: スペクトルデータの格納位置と全波数・帯状波数の変換 ! !==== 微分 ! ! wq_RadDRad_wq :: スペクトルデータに動径微分 r∂/∂r を作用させる ! wr_DivRad_wq :: スペクトルデータに発散型動径微分 ! 1/r^2 ∂/∂r r^2 = ∂/∂r + 2/r を作用させる ! wr_DivRad_xyr :: 格子点データに発散型動径微分 ! 1/r^2 ∂/∂r r^2 = ∂/∂r + 2/r を作用させる ! wr_DivRad_xvr :: 格子点データに発散型動径微分 ! 1/r^2 ∂/∂r r^2 = ∂/∂r + 2/r を作用させる ! wr_RotDRad_wq :: スペクトルデータに回転型動径微分 ! 1/r ∂/∂rr = ∂/∂r + 1/r を作用させる ! wr_RotDRad_wr :: スペクトルデータに回転型動径微分 ! 1/r ∂/∂rr = ∂/∂r + 1/r を作用させる ! wq_RotDRad_wr :: スペクトルデータに回転型動径微分 ! 1/r ∂/∂rr = ∂/∂r + 1/r を作用させる ! wq_Lapla_wq :: スペクトルデータにラプラシアンを作用させる ! xyr_GradLon_wq :: スペクトルデータに勾配型経度微分 ! 1/rcosφ・∂/∂λを作用させる ! xvr_GradLon_wq :: スペクトルデータに勾配型経度微分 ! 1/rcosφ・∂/∂λを作用させる ! xyr_GradLat_wq :: スペクトルデータに勾配型緯度微分 ! 1/r・∂/∂φを作用させる ! xvr_GradLat_wq :: スペクトルデータに勾配型緯度微分 ! 1/r・∂/∂φを作用させる ! wr_DivLon_xyr :: 格子データに発散型経度微分 ! 1/rcosφ・∂/∂λを作用させる ! wr_DivLon_xvr :: 格子データに発散型経度微分 ! 1/rcosφ・∂/∂λを作用させる ! wr_DivLat_xyr :: 格子データに発散型緯度微分 ! 1/rcosφ・∂(g cosφ)/∂φを作用させる ! wr_DivLat_xvr :: 格子データに発散型緯度微分 ! 1/rcosφ・∂(g cosφ)/∂φを作用させる ! wr_Div_xyr_xyr_xyr :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! xyr_Div_xyr_xyr_xyr :: ベクトル成分である 3 つの格子データに ! 発散を作用させる ! xyr_RotLon_wq_wq :: ベクトル場の回転の経度成分を計算する ! xvr_RotLon_wq_wq :: ベクトル場の回転の経度成分を計算する ! xyr_RotLat_wq_wq :: ベクトル場の回転の緯度成分を計算する ! xvr_RotLat_wq_wq :: ベクトル場の回転の緯度成分を計算する ! wr_RotRad_xyr_xyr :: ベクトル場の回転の動径成分を計算する ! wr_RotRad_xvr_xvr :: ベクトル場の回転の動径成分を計算する ! !==== トロイダルポロイダル計算用微分 ! ! wq_KxRGrad_wq :: スペクトルデータに経度微分 ! k×r・▽ = ∂/∂λを作用させる ! xyr_KGrad_wq :: スペクトルデータに軸方向微分 ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r を作用させる ! xvr_KGrad_wq :: スペクトルデータに軸方向微分 ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r を作用させる ! wq_L2_wq :: スペクトルデータに L2 演算子 = -水平ラプラシアンを ! 作用させる ! wq_L2Inv_wq :: スペクトルデータに L2 演算子の逆 = -逆水平ラプラシアン ! を作用させる ! wq_QOperator_wq :: スペクトルデータに演算子 ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2)) を作用させる ! wq_QOperatorMPI_wq:: スペクトルデータに演算子 ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2)) を作用させる ! wr_RadRot_xyr_xyr :: ベクトル v の渦度と動径ベクトル r の内積 r・(▽×v) を ! 計算する ! wr_RadRot_xvr_xvr :: ベクトル v の渦度と動径ベクトル r の内積 r・(▽×v) を ! 計算する ! wr_RadRotRot_xyr_xyr_xyr :: ベクトルの v の r・(▽×▽×v) を計算する ! wr_RadRotRot_xvr_xvr_xvr :: ベクトルの v の r・(▽×▽×v) を計算する ! wq_RadRotRot_xyr_xyr_xyr :: ベクトルの v の r・(▽×▽×v) を計算する ! wq_RadRotRot_xvr_xvr_xvr :: ベクトルの v の r・(▽×▽×v) を計算する ! wq_Potential2Vector :: トロイダルポロイダルポテンシャルから ! ベクトル場を計算する ! wq_Potential2VectorMPI :: トロイダルポロイダルポテンシャルから ! ベクトル場を計算する ! wq_Potential2Rotation :: トロイダルポロイダルポテンシャルで表される ! 非発散ベクトル場の回転の各成分を計算する ! wq_Potential2RotationMPI :: トロイダルポロイダルポテンシャルで表される ! 非発散ベクトル場の回転の各成分を計算する ! !==== ポロイダル/トロイダルモデル用スペクトル解析 ! ! nmr_ToroidalEnergySpectrum_wq, nr_ToroidalEnergySpectrum_wq :: ! トロイダルポテンシャルからエネルギーの球面調和函数各成分を計算する ! nmr_PoloidalEnergySpectrum_wq, nr_PoloidalEnergySpectrum_wq :: ! ポロイダルポテンシャルからエネルギーの球面調和函数各成分を計算する ! !==== 境界値問題 ! ! wq_BoundaryTau, wr_BoundaryGrid, wq_Boundary :: ! ディリクレ, ノイマン境界条件を適用する(タウ法, 選点法) ! wq_TorBoundaryTau, wr_TorBoundaryGrid, wq_TorBoundary :: ! 速度トロイダルポテンシャルの境界条件を適用する(タウ法,選点法) ! wq_LaplaPol2Pol_wq, wq_LaplaPol2PolTau_wq :: ! 速度ポロイダルポテンシャルΦを▽^2Φから求める ! (入出力がそれぞれ格子点およびスペクトル係数) ! wq_TorMagBoundaryTau, wr_TorMagBoundaryGrid, wq_TorMagBoundary :: ! 磁場トロイダルポテンシャルの境界条件を適用する(タウ法, 選点法) ! wq_PolMagBoundaryTau, wr_PolMagBoundaryGrid, wq_PolMagBoundary :: ! 磁場トロイダルポテンシャル境界の境界条件を適用する(タウ法, 選点法) ! !==== 積分・平均(3 次元データ) ! ! IntLonLatRad_xyr, AvrLonLatRad_xyr :: 3 次元格子点データの ! 全領域積分および平均 ! IntLonLatRad_xvr, AvrLonLatRad_xvr :: 3 次元格子点データの ! 全領域積分および平均 ! r_IntLonLat_xyr, r_AvrLonLat_xyr :: 3 次元格子点データの ! 緯度経度(水平・球面)積分および平均 ! r_IntLonLat_xvr, r_AvrLonLat_xvr :: 3 次元格子点データの ! 緯度経度(水平・球面)積分および平均 ! y_IntLonRad_xyr, y_AvrLonRad_xyr :: 3 次元格子点データの ! 緯度動径積分および平均 ! v_IntLonRad_xvr, v_AvrLonRad_xvr :: 3 次元格子点データの ! 緯度動径積分および平均 ! x_IntLatRad_xyr,xr_AvrLatRad_xyr :: 3 次元格子点データの ! 経度動径(子午面)積分および平均 ! x_IntLatRad_xvr,xr_AvrLatRad_xvr :: 3 次元格子点データの ! 経度動径(子午面)積分および平均 ! yr_IntLon_xyr, yr_AvrLon_xyr :: 3 次元格子点データの ! 経度方向積分および平均 ! vr_IntLon_xvr, vr_AvrLon_xvr :: 3 次元格子点データの ! 経度方向積分および平均 ! xr_IntLat_xyr, xr_AvrLat_xyr :: 3 次元格子点データの ! 緯度方向積分および平均 ! xr_IntLat_xvr, xr_AvrLat_xvr :: 3 次元格子点データの ! 緯度方向積分および平均 ! xr_IntRad_xyr, xr_AvrRad_xyr :: 3 次元格子点データの ! 動径方向積分および平均 ! xr_IntRad_xvr, xr_AvrRad_xvr :: 3 次元格子点データの ! 動径方向積分および平均 ! !==== 積分・平均(2 次元データ) ! ! IntLonLat_xy, AvrLonLat_xy :: 2 次元格子点データの水平(球面)積分および平均 ! IntLonLat_xv, AvrLonLat_xv :: 2 次元格子点データの水平(球面)積分および平均 ! IntLonRad_xr, AvrLonRad_xr :: 2 次元(XZ)格子点データの経度動径積分 ! および平均 ! IntLatRad_yr, AvrLatRad_yr :: 2 次元(YZ)格子点データの緯度動径(子午面) ! 積分および平均 ! IntLatRad_vr, AvrLatRad_vr :: 2 次元(YZ)格子点データの緯度動径(子午面) ! 積分および平均 ! y_IntLon_xy, y_AvrLon_xy :: 水平 2 次元(球面)格子点データの経度方向 ! 積分および平均 ! v_IntLon_xv, v_AvrLon_xv :: 水平 2 次元(球面)格子点データの経度方向 ! 積分および平均 ! x_IntLat_xy, x_AvrLat_xy :: 水平2 次元(球面)格子点データの緯度方向積分 ! および平均 ! x_IntLat_xv, x_AvrLat_xv :: 水平2 次元(球面)格子点データの緯度方向積分 ! および平均 ! r_IntLon_xr, r_AvrLon_xr :: 2 次元(XZ)格子点データの経度方向積分および ! 平均 ! x_IntRad_xr, x_AvrRad_xr :: 2 次元(XZ)格子点データの動径方向積分および ! 平均 ! r_IntLat_yr, r_AvrLat_yr :: 2 次元(YZ)格子点データの緯度方向積分および ! 平均 ! r_IntLat_vr, r_AvrLat_vr :: 2 次元(YZ)格子点データの緯度方向積分および ! 平均 ! y_IntRad_yr, y_AvrRad_yr :: 2 次元(YZ)格子点データの動径方向積分および ! 平均 ! v_IntRad_vr, v_AvrRad_vr :: 2 次元(YZ)格子点データの動径方向積分および ! 平均 ! !==== 積分・平均(1 次元データ) ! ! IntLon_x, AvrLon_x :: 1 次元(X)格子点データの経度方向積分および平均 ! IntLat_y, AvrLat_y :: 1 次元(Y)格子点データの緯度方向積分および平均 ! IntLat_v, AvrLat_v :: 1 次元(Y)格子点データの緯度方向積分および平均 ! IntRad_r, AvrRad_r :: 1 次元(Z)格子点データの動径方向積分および平均 ! ! use dc_message use lumatrix use w_base_mpi_module use w_deriv_mpi_module use w_integral_mpi_module use wa_base_mpi_module use wa_deriv_mpi_module use wq_module use wa_module, only : wa_Lapla_wa use mpi implicit none integer :: ierr private public wq_mpi_Initial public jc public x_Lon, x_Lon_Weight public y_Lat, y_Lat_Weight public v_Lat, v_Lat_Weight public r_Rad, r_Rad_Weight public l_nm, nm_l public xy_Lon, xy_Lat public xv_Lon, xv_Lat public xyr_Lon, xyr_Lat, xyr_Rad public xvr_Lon, xvr_Lat, xvr_Rad public wr_Rad public wq_VMiss public w_xy, xy_w public w_xv, xv_w public wq_RadDRad_wq, q_RadDRad_q, wr_wq, wq_wr public wq_Rad2_wq, q_Rad2_q, wq_Rad2Inv_wq, q_Rad2Inv_q public xyr_wq, wq_xyr, xyr_wr, wr_xyr public xvr_wq, wq_xvr, xvr_wr, wr_xvr public wr_DivRad_wq, wr_RotDRad_wq, wr_RotDRad_wr, wq_Lapla_wq public wq_RotDRad_wr public xyr_GradLon_wq, xyr_GradLat_wq public xvr_GradLon_wq, xvr_GradLat_wq public wr_DivLon_xyr, wr_DivLat_xyr, wr_DivRad_xyr public wr_DivLon_xvr, wr_DivLat_xvr, wr_DivRad_xvr public wr_Div_xyr_xyr_xyr, xyr_Div_xyr_xyr_xyr public wr_Div_xvr_xvr_xvr, xvr_Div_xvr_xvr_xvr public xyr_RotLon_wq_wq, xyr_RotLat_wq_wq, wr_RotRad_xyr_xyr public xvr_RotLon_wq_wq, xvr_RotLat_wq_wq, wr_RotRad_xvr_xvr public yr_IntLon_xyr, xr_IntLat_xyr, xy_IntRad_xyr public vr_IntLon_xvr, xr_IntLat_xvr, xv_IntRad_xvr public x_IntLatRad_xyr, y_IntLonRad_xyr, r_IntLonLat_xyr public x_IntLatRad_xvr, v_IntLonRad_xvr, r_IntLonLat_xvr public IntLonLatRad_xyr public IntLonLatRad_xvr public x_IntLat_xy, y_IntLon_xy, IntLonLat_xy public x_IntLat_xv, v_IntLon_xv, IntLonLat_xv public r_IntLat_yr, y_IntRad_yr, IntLatRad_yr public r_IntLat_vr, v_IntRad_vr, IntLatRad_vr public r_IntLon_xr, x_IntRad_xr, IntLonRad_xr public IntLon_x, IntLat_y, IntRad_r public IntLat_v public yr_AvrLon_xyr, xr_AvrLat_xyr, xy_AvrRad_xyr public vr_AvrLon_xvr, xr_AvrLat_xvr, xv_AvrRad_xvr public x_AvrLatRad_xyr, y_AvrLonRad_xyr, r_AvrLonLat_xyr public x_AvrLatRad_xvr, v_AvrLonRad_xvr, r_AvrLonLat_xvr public AvrLonLatRad_xyr public AvrLonLatRad_xvr public x_AvrLat_xy, y_AvrLon_xy, AvrLonLat_xy public x_AvrLat_xv, v_AvrLon_xv, AvrLonLat_xv public r_AvrLat_yr, y_AvrRad_yr, AvrLatRad_yr public r_AvrLat_vr, v_AvrRad_vr, AvrLatRad_vr public r_AvrLon_xr, x_AvrRad_xr, AvrLonRad_xr public AvrLon_x, AvrLat_y, AvrRad_r public AvrLat_v public wq_KxRGrad_wq, xyr_KGrad_wq, wq_L2_wq, wq_L2Inv_wq, wq_QOperator_wq public xvr_KGrad_wq, wq_QOperatorMPI_wq public wr_RadRot_xyr_xyr, wr_RadRotRot_xyr_xyr_xyr public wr_RadRot_xvr_xvr, wr_RadRotRot_xvr_xvr_xvr public wq_RadRotRot_xyr_xyr_xyr public wq_RadRotRot_xvr_xvr_xvr public wq_Potential2vector, wq_Potential2Rotation public wq_Potential2vectorMPI, wq_Potential2RotationMPI public nmr_ToroidalEnergySpectrum_wq, nr_ToroidalEnergySpectrum_wq public nmr_PoloidalEnergySpectrum_wq, nr_PoloidalEnergySpectrum_wq public wq_Boundary, wq_TorBoundary, wq_LaplaPol2Pol_wq ! wr_LaplaPol2Pol_wr public wq_TormagBoundary, wq_PolmagBoundary public wq_BoundaryTau, wq_TorBoundaryTau, wq_LaplaPol2PolTau_wq public wq_TormagBoundaryTau, wq_PolmagBoundaryTau public wr_BoundaryGrid, wr_TorBoundaryGrid public wr_TormagBoundaryGrid, wr_PolmagBoundaryGrid integer :: im=64, jm=32, km=16 ! 格子点の設定(経度, 緯度, 動径) integer :: nm=21, lm=31 ! 切断波数の設定(水平, 動径) real(8) :: ra=1.0 ! 球半径 real(8), parameter :: pi=3.1415926535897932385D0 real(8), dimension(:,:,:), allocatable :: xvr_LON, xvr_LAT, xvr_RAD ! 座標 save im, jm, km, nm, lm, ra, xvr_Lon, xvr_Lat, xvr_Rad contains !--------------- 初期化 ----------------- subroutine wq_mpi_Initial(i,j,k,n,l,r,np,wa_init) ! ! スペクトル変換の格子点数, 波数, 動径座標の範囲を設定する. ! ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を ! しなければならない. ! ! np に 1 より大きな値を指定すれば ISPACK の球面調和函数変換 ! OPENMP 並列計算ルーチンが用いられる. 並列計算を実行するには, ! 実行時に環境変数 OMP_NUM_THREADS を np 以下の数字に設定する等の ! システムに応じた準備が必要となる. ! ! np に 1 より大きな値を指定しなければ並列計算ルーチンは呼ばれない. ! ! integer,intent(in) :: i ! 格子点数(経度λ) integer,intent(in) :: j ! 格子点数(緯度φ) integer,intent(in) :: k ! 格子点数(動径 r) integer,intent(in) :: n ! 切断波数(水平全波数) integer,intent(in) :: l ! 切断波数(動径波数) real(8),intent(in) :: r ! 球半径 integer,intent(in), optional :: np ! OPENMP での最大スレッド数 logical,intent(in), optional :: wa_init ! wa_initial スイッチ logical :: wa_initialize=.true. ! wa_initial スイッチ im = i ; jm = j ; km = k nm = n ; lm = l ra = r if ( present(wa_init) ) then wa_initialize = wa_init else wa_initialize = .true. endif if ( present(np) ) then call wq_Initial(im,jm,km,nm,lm,ra,np,wa_initialize) else call wq_Initial(im,jm,km,nm,lm,ra,wa_init=wa_initialize) endif if ( wa_initialize ) then call w_base_mpi_Initial call w_deriv_mpi_Initial call wa_base_mpi_Initial endif allocate(xvr_Lon(0:im-1,jc,km)) allocate(xvr_Lat(0:im-1,jc,km)) allocate(xvr_Rad(0:im-1,jc,km)) xvr_Lon = spread(xv_Lon,3,km) xvr_Lat = spread(xv_Lat,3,km) xvr_Rad = spread(spread(r_Rad,1,jc),1,im) call MessageNotify('M','wq_mpi_initial','wq_mpi_module (2011/09/08) is initialized') end subroutine wq_mpi_Initial !--------------- 基本変換 ----------------- function xvr_wq(wq) ! ! スペクトルデータから 3 次元分散格子点データへ(逆)変換する. ! real(8), dimension(0:im-1,1:jc,km) :: xvr_wq !(out) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq !(in) 2 次元球面調和函数チェビシェフスペクトルデータ xvr_wq = xva_wa(wr_wq(wq)) end function xvr_wq function wq_xvr(xvr) ! ! 3 次元格子点データからスペクトルデータへ(正)変換する. ! real(8), dimension((nm+1)*(nm+1),0:lm) :: wq_xvr !(out) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ wq_xvr = wq_wr(wa_xva(xvr)) end function wq_xvr function xvr_wr(wr) ! ! 水平スペクトル・動径格子点データから 3 次元分散格子点データへ(逆)変換する. ! real(8), dimension(0:im-1,1:jc,km) :: xvr_wr !(out) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),km), intent(in) :: wr !(in) 2 次元球面調和函数スペクトル・動径格子点データ xvr_wr = xva_wa(wr) end function xvr_wr function wr_xvr(xvr) ! ! 3 次元格子データから水平スペクトル・動径格子点データへ(正)変換する. ! real(8), dimension((nm+1)*(nm+1),km) :: wr_xvr !(out) 2 次元球面調和函数スペクトル・動径格子点データ real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ wr_xvr = wa_xva(xvr) end function wr_xvr !--------------- 微分計算 ----------------- function wr_DivRad_xvr(xvr) ! ! 格子点データに発散型動径微分 ! ! 1/r^2 ∂/∂r (r^2 = ∂/∂r + 2/= 1/r∂/∂(r.) + 1/r ! ! を作用する. ! ! この関数の入力は動径次数と水平全波数の偶奇性が異なっていることを ! 仮定している. 偶奇性が一致している場合には wr_DigRad_wq を使用すること. ! real(8), dimension(0:im-1,jc,km), intent(in) :: xvr !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension((nm+1)*(nm+1),km) :: wr_DivRad_xvr !(out) 発散型動径微分を作用された水平スペクトル動径格子点データ wr_DivRad_xvr = wr_wq(wq_RadDRad_wq(wq_xvr(xvr_Rad*xvr)))/wr_Rad**2 & + 1/wr_Rad * wr_xvr(xvr) end function wr_DivRad_xvr function xvr_GradLon_wq(wq) ! ! スペクトルデータに勾配型経度微分 1/rcosφ・∂/∂λ ! を作用させる. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:im-1,1:jc,km) :: xvr_GradLon_wq !(out) 勾配型経度微分を作用された 2 次元スペクトルデータ xvr_GradLon_wq = xva_GradLon_wa(wr_wq(wq))/xvr_Rad end function xvr_GradLon_wq function xvr_GradLat_wq(wq) ! ! スペクトルデータに勾配型経度微分 1/r ∂/∂φ を作用させる. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:im-1,1:jc,km) :: xvr_GradLat_wq !(out) 勾配型緯度微分を作用された 2 次元スペクトルデータ xvr_GradLat_wq = xva_GradLat_wa(wr_wq(wq))/xvr_Rad end function xvr_GradLat_wq function wr_DivLon_xvr(xvr) ! ! 格子点データに発散型経度微分 1/rcosφ・∂/∂λ を作用させた ! スペクトルデータを返す. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),km) :: wr_DivLon_xvr !(out) 発散型経度微分を作用された水平スペクトル動径格子点データ wr_DivLon_xvr = wa_DivLon_xva(xvr/xvr_Rad) end function wr_DivLon_xvr function wr_DivLat_xvr(xvr) ! ! 格子データに発散型緯度微分 1/rcosφ・∂(f cosφ)/∂φ を ! 作用させたスペクトルデータを返す. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension((nm+1)*(nm+1),km) :: wr_DivLat_xvr !(out) 発散型緯度微分を作用された水平スペクトル動径格子点データ wr_DivLat_xvr = wa_DivLat_xva(xvr/xvr_Rad) end function wr_DivLat_xvr function wr_Div_xvr_xvr_xvr(xvr_Vlon,xvr_Vlat,xvr_Vrad) ! ! べクトル成分である 3 つの格子データに発散を作用させた ! スペクトルデータを返す. ! ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, ! 動径成分を表し, 発散は ! ! 1/rcosφ・∂u/∂λ + 1/rcosφ・∂(v cosφ)/∂φ ! + 1/r^2 ∂/∂r (r^2 w) ! ! と計算される. ! ! 第3成分の動径次数と水平全波数の偶奇性がずれていることを仮定している. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_Vlon !(in) ベクトル場の経度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_Vrad !(in) ベクトル場の動径成分 real(8), dimension((nm+1)*(nm+1),km) :: wr_Div_xvr_xvr_xvr !(out) ベクトル場の発散 wr_Div_xvr_xvr_xvr = wr_DivLon_xvr(xvr_Vlon) & + wr_DivLat_xvr(xvr_Vlat) & + wr_DivRad_xvr(xvr_Vrad) end function wr_Div_xvr_xvr_xvr function xvr_Div_xvr_xvr_xvr(xvr_Vlon,xvr_Vlat,xvr_Vrad) ! ! ベクトル成分である 3 つの格子データに発散を作用させる. ! ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, ! 動径成分を表す. ! ! 極の特異性を回避するためにベクトル場に cosφ/r の重みをかけて ! 計算している. ! ! div V = (r/cosφ)・div (Vcosφ/r) + V_φtanφ/r + V_r/r ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_Vlon !(in) ベクトル場の経度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_Vrad !(in) ベクトル場の動径成分 real(8), dimension(0:im-1,1:jc,km) :: xvr_Div_xvr_xvr_xvr !(out) ベクトル場の発散 xvr_Div_xvr_xvr_xvr & = xvr_Rad/cos(xvr_Lat) & * xvr_wr(wr_Div_xvr_xvr_xvr(xvr_VLon*cos(xvr_Lat)/xvr_Rad, & xvr_VLat*cos(xvr_Lat)/xvr_Rad, & xvr_VRad*cos(xvr_Lat)/xvr_Rad ))& + xvr_VLat*tan(xvr_Lat)/xvr_Rad & + xvr_VRad/xvr_Rad end function xvr_Div_xvr_xvr_xvr function xvr_RotLon_wq_wq(wq_Vrad,wq_Vlat) ! ! ベクトル場の動径成分, 緯度成分である第 1, 2 引数 Vrad, Vlat から ! 回転の経度成分 ! ! 1/r ∂Vrad/∂φ-1/r ∂(r Vlat)/∂r を計算する. ! ! を計算する ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq_Vrad !(in) ベクトル場の動径成分 real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jc,km) :: xvr_RotLon_wq_wq !(out) ベクトル場の回転の経度成分 xvr_RotLon_wq_wq = xvr_GradLat_wq(wq_Vrad) & - xvr_wr(wr_RotDRad_wq(wq_Vlat)) end function xvr_RotLon_wq_wq function xvr_RotLat_wq_wq(wq_Vlon,wq_Vrad) ! ! ベクトル場の経度成分, 動径成分である第 1, 2 引数 Vlon, Vrad から ! 回転の緯度成分 ! ! 1/r ∂(r Vlon)/∂r - 1/rcosφ・∂Vrad/∂λ ! ! を計算する. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq_Vlon !(in) ベクトル場の経度成分 real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq_Vrad !(in) ベクトル場の動径成分 real(8), dimension(0:im-1,1:jc,km) :: xvr_RotLat_wq_wq !(out) ベクトル場の回転の緯度成分 xvr_RotLat_wq_wq = xvr_wr(wr_RotDRad_wq(wq_Vlon)) & - xvr_GradLon_wq(wq_Vrad) end function xvr_RotLat_wq_wq function wr_RotRad_xvr_xvr(xvr_Vlat,xvr_Vlon) ! ! ベクトルの緯度成分, 経度成分である第 1, 2 引数 Vlat, Vlon に対して ! ベクトル場の回転の動径成分 ! ! 1/rcosφ・∂Vlat/∂λ - 1/rcosφ・∂(Vlon cosφ)/∂φ ! ! を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_Vlat !(in) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_Vlon !(in) ベクトル場の経度成分 real(8), dimension((nm+1)*(nm+1),km) :: wr_RotRad_xvr_xvr !(out) ベクトル場の回転の動径成分 wr_RotRad_xvr_xvr = wr_DivLon_xvr(xvr_Vlat) & - wr_DivLat_xvr(xvr_Vlon) end function wr_RotRad_xvr_xvr !--------------- 積分計算 ----------------- !----(入力データ xvr)--- function vr_IntLon_xvr(xvr) ! 経度(帯状)積分 ! ! 3 次元格子点データの経度方向(帯状)積分. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jc,km) :: vr_IntLon_xvr !(out) 経度方向(帯状)積分された 2 次元子午面格子点データ integer :: i vr_IntLon_xvr = 0.0d0 do i=0,im-1 vr_IntLon_xvr(:,:) = vr_IntLon_xvr(:,:) & + xvr(i,:,:) * x_Lon_Weight(i) enddo end function vr_IntLon_xvr function xr_IntLat_xvr(xvr) ! ! 3 次元格子点データの緯度方向域積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) cosφ dφ を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,km) :: xr_IntLat_xvr !(out) 緯度積分された 2 次元緯度動径格子点データ ! 緯度円格子点データ real(8), dimension(0:im-1,km) :: xr_IntLatTMP integer :: j xr_IntLat_xvr = 0.0d0 do j=1,jc xr_IntLat_xvr(:,:) = xr_IntLat_xvr(:,:) & + xvr(:,j,:) * v_Lat_Weight(j) enddo xr_IntLatTmp=xr_IntLat_xvr CALL MPI_ALLREDUCE(xr_IntLatTMP,xr_IntLat_xvr,im*km,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function xr_IntLat_xvr function xv_IntRad_xvr(xvr) ! 動径積分 ! ! 3 次元格子点データの動径方向域積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dr を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,1:jc) :: xv_IntRad_xvr !(out) 動径積分された 2 次元経度緯度(水平, 球面)格子点データ integer :: k xv_IntRad_xvr = 0.0d0 do k=1,km xv_IntRad_xvr(:,:) = xv_IntRad_xvr(:,:) & + xvr(:,:,k) * r_Rad_Weight(k) enddo end function xv_IntRad_xvr function x_IntLatRad_xvr(xvr) ! ! 3 次元格子点データの緯度動径(子午面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dφdr ! ! を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1) :: x_IntLatRad_xvr !(out) 緯度動径(子午面)積分された 1 次元経度格子点データ real(8), dimension(0:im-1) :: x_IntLatRadTMP integer :: j, k x_IntLatRad_xvr = 0.0D0 do k=1,km do j=1,jc x_IntLatRad_xvr = x_IntLatRad_xvr & + xvr(:,j,k) * v_Lat_Weight(j) * r_Rad_Weight(k) enddo enddo x_IntLatRadTmp=x_IntLatRad_xvr CALL MPI_ALLREDUCE(x_IntLatRadTMP,x_IntLatRad_xvr,im,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function x_IntLatRad_xvr function v_IntLonRad_xvr(xvr) ! ! 3 次元格子点データの経度動径(緯度円)積分. ! ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dλdr を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jc) :: v_IntLonRad_xvr !(out) 経度動径(緯度円)積分された 1 次元緯度格子点データ integer :: i, k v_IntLonRad_xvr = 0.0D0 do k=1,km do i=0,im-1 v_IntLonRad_xvr = v_IntLonRad_xvr & + xvr(i,:,k) * x_Lon_Weight(i) * r_Rad_Weight(k) enddo enddo end function v_IntLonRad_xvr function r_IntLonLat_xvr(xvr) ! 緯度経度(水平)積分 ! ! 3 次元格子点データの緯度経度(水平, 球面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) cosφ dλdφ ! ! を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(km) :: r_IntLonLat_xvr !(out) 緯度経度(水平, 球面)積分された 1 次元動径格子点データ real(8), dimension(km) :: r_IntLonLatTMP integer :: i, j r_IntLonLat_xvr = 0.0D0 do j=1,jc do i=0,im-1 r_IntLonLat_xvr = r_IntLonLat_xvr & + xvr(i,j,:) * x_Lon_Weight(i) * v_Lat_Weight(j) enddo enddo r_IntLonLatTmp=r_IntLonLat_xvr CALL MPI_ALLREDUCE(r_IntLonLatTMP,r_IntLonLat_xvr,km,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function r_IntLonLat_xvr function IntLonLatRad_xvr(xvr) ! 緯度経度動径(全球)積分 ! ! 3 次元格子点データの緯度経度動径(全球)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dλdφdr ! ! を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8) :: IntLonLatRad_xvr !(out) 全球積分値 real(8) :: IntLonLatRadTMP integer :: i, j, k IntLonLatRad_xvr = 0.0D0 do k=1,km do j=1,jc do i=0,im-1 IntLonLatRad_xvr = IntLonLatRad_xvr & + xvr(i,j,k) * x_Lon_Weight(i) & * v_Lat_Weight(j) * r_Rad_Weight(k) enddo enddo enddo IntLonLatRadTmp=IntLonLatRad_xvr CALL MPI_ALLREDUCE(IntLonLatRadTMP,IntLonLatRad_xvr,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntLonLatRad_xvr !----(入力データ vr)--- function r_IntLat_vr(vr) ! 緯度積分 ! ! 2 次元(VR)格子点データの緯度方向域積分. ! ! 2 次元データ f(φ,r) に対して∫f(φ,r) cosφ dφ を計算する. ! real(8), dimension(1:jc,km), intent(in) :: vr !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(km) :: r_IntLat_vr !(out) 緯度積分された 1 次元動径格子点データ real(8), dimension(km) :: r_IntLatTMP integer :: j r_IntLat_vr = 0.0d0 do j=1,jc r_IntLat_vr(:) = r_IntLat_vr(:) + vr(j,:) * v_Lat_Weight(j) enddo r_IntLatTmp=r_IntLat_vr CALL MPI_ALLREDUCE(r_IntLatTMP,r_IntLat_vr,km,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function r_IntLat_vr function v_IntRad_vr(vr) ! 動径積分 ! ! 2 次元(VR)格子点データの動径方向域積分. ! ! 2 次元データ f(φ,r) に対して∫f(φ,r) r^2dr を計算する. ! real(8), dimension(1:jc,km), intent(in) :: vr !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(1:jc) :: v_IntRad_vr !(out) 動径積分された 1 次元緯度格子点データ integer :: k v_IntRad_vr = 0.0d0 do k=1,km v_IntRad_vr(:) = v_IntRad_vr(:) & + vr(:,k) * r_Rad_Weight(k) enddo end function v_IntRad_vr function IntLatRad_vr(vr) ! ! 2 次元(VR)格子点データの緯度動径積分(子午面)および平均 ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2cosφ dφdr を計算する. ! real(8), dimension(1:jc,km), intent(in) :: vr !(in) 2 次元緯度動径(子午面)格子点データ real(8) :: IntLatRad_vr !(out) 積分値 real(8) :: IntLatRadTMP integer :: j, k IntLatRad_vr = 0.0D0 do k=1,km do j=1,jc IntLatRad_vr = IntLatRad_vr & + vr(j,k) * v_Lat_Weight(j) * r_Rad_Weight(k) enddo enddo IntLatRadTmp=IntLatRad_vr CALL MPI_ALLREDUCE(IntLatRadTMP,IntLatRad_vr,1,MPI_REAL8, & MPI_SUM,MPI_COMM_WORLD,IERR) end function IntLatRad_vr !--------------- 平均計算 ----------------- !----(入力データ xvr)--- function vr_AvrLon_xvr(xvr) ! 経度(帯状)積分 ! ! 3 次元格子点データの経度方向(帯状)平均. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ/2π を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jc,km) :: vr_AvrLon_xvr !(out) 経度方向(帯状)平均された 2 次元子午面格子点データ vr_AvrLon_xvr = vr_IntLon_xvr(xvr)/sum(x_Lon_Weight) end function vr_AvrLon_xvr function xr_AvrLat_xvr(xvr) ! 緯度積分 ! ! 3 次元格子点データの緯度方向域平均. ! ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)cosφ dφ/2 を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,km) :: xr_AvrLat_xvr !(out) 緯度平均された 2 次元緯度動径格子点データ xr_AvrLat_xvr = xr_IntLat_xvr(xvr)/sum(y_Lat_Weight) end function xr_AvrLat_xvr function xv_AvrRad_xvr(xvr) ! ! 3 次元格子点データの動径方向域平均. ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2dr/((r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1,1:jc) :: xv_AvrRad_xvr !(out) 動径平均された 2 次元経度緯度(水平, 球面)格子点データ ! 水平格子点データ xv_AvrRad_xvr = xv_IntRad_xvr(xvr)/sum(r_Rad_Weight) end function xv_AvrRad_xvr function x_AvrLatRad_xvr(xvr) ! 緯度動径(子午面)積分 ! ! 3 次元格子点データの緯度動径(子午面)平均 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(0:im-1) :: x_AvrLatRad_xvr !(out) 緯度動径(子午面)平均された 1 次元経度格子点データ x_AvrLatRad_xvr = x_IntLatRad_xvr(xvr) & /( sum(y_Lat_Weight)*sum(r_Rad_Weight) ) end function x_AvrLatRad_xvr function v_AvrLonRad_xvr(xvr) ! 経度動径(緯度円)積分 ! ! 3 次元格子点データの経度動径(緯度円)平均. ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2dλdr /(2π(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(1:jc) :: v_AvrLonRad_xvr !(out) 経度動径(緯度円)平均された 1 次元緯度格子点データ v_AvrLonRad_xvr = v_IntLonRad_xvr(xvr) & /(sum(x_Lon_Weight)*sum(r_Rad_Weight)) end function v_AvrLonRad_xvr function r_AvrLonLat_xvr(xvr) ! 緯度経度(水平)積分 ! ! 3 次元格子点データの緯度経度(水平, 球面)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) cosφ dλdφ /4π ! ! を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8), dimension(km) :: r_AvrLonLat_xvr !(out) 緯度経度(水平, 球面)平均された 1 次元動径格子点データ r_AvrLonLat_xvr = r_IntLonLat_xvr(xvr) & /(sum(x_Lon_Weight)*sum(y_Lat_Weight)) end function r_AvrLonLat_xvr function AvrLonLatRad_xvr(xvr) ! 緯度経度動径(全球)積分 ! ! 3 次元格子点データの緯度経度動径(全球)積分 ! ! 3 次元データ f(λ,φ,r) に対して ! ! ∫f(λ,φ,r) r^2cosφ dλdφdr /(4π(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr !(in) 3 次元経度緯度動径格子点データ real(8) :: AvrLonLatRad_xvr !(out) 全球平均値 AvrLonLatRad_xvr = IntLonLatRad_xvr(xvr) & /(sum(x_Lon_Weight)*sum(y_Lat_Weight) * sum(r_Rad_Weight)) end function AvrLonLatRad_xvr !----(入力データ vr)--- function r_AvrLat_vr(vr) ! ! 2 次元(VR)格子点データの緯度方向域平均. ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) cosφ dφ/2 を計算する. ! real(8), dimension(1:jc,km), intent(in) :: vr !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(km) :: r_AvrLat_vr !(out) 緯度平均された 1 次元動径格子点データ r_AvrLat_vr = r_IntLat_vr(vr)/sum(y_Lat_Weight) end function r_AvrLat_vr function v_AvrRad_vr(vr) ! ! 2 次元(VR)格子点データの動径方向域平均. ! ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2dr /((r[o]^3-r[i]^3)/3) ! を計算する. ! real(8), dimension(1:jc,km), intent(in) :: vr !(in) 2 次元緯度動径(子午面)格子点データ real(8), dimension(1:jc) :: v_AvrRad_vr !(out) 動径平均された 1 次元緯度格子点データ v_AvrRad_vr = v_IntRad_vr(vr)/sum(r_Rad_Weight) end function v_AvrRad_vr function AvrLatRad_vr(vr) ! 緯度動径(子午面)積分 ! ! 2 次元(VR)格子点データの緯度動径(子午面)平均 ! ! 2 次元データ f(φ,r) に対して ! ! ∫f(φ,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) ! ! を計算する. ! real(8), dimension(1:jc,km), intent(in) :: vr !(in) 2 次元緯度動径(子午面)格子点データ real(8) :: AvrLatRad_vr !(out) 平均値 AvrLatRad_vr = IntLatRad_vr(vr)/(sum(y_Lat_Weight)*sum(r_Rad_Weight)) end function AvrLatRad_vr !--------------- ポロイダル/トロイダルモデル用微分 ----------------- function xvr_KGrad_wq(wq) ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r ! ! 入力スペクトルデータに対応する格子データに軸方向微分 ! ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r ! ! を作用させた格子データが返される. ! ここでベクトル k は球の中心から北極向きの単位ベクトルである. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension(0:im-1,1:jc,km) :: xvr_KGrad_wq !(out) 軸方向微分を作用された 2 次元スペクトルデータ xvr_KGrad_wq = cos(xvr_Lat)*xvr_GradLat_wq(wq) & + sin(xvr_Lat)*xvr_wq(wq_RadDRad_wq(wq))/xvr_Rad end function xvr_KGrad_wq function wq_QOperatorMPI_wq(wq) ! ! 入力スペクトルデータに対応する格子点データに演算子 ! ! Q=(k・▽-1/2(L2 k・▽+ k・▽L2)) ! ! を作用させたデータのスペクトル変換が返される. ! real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq !(in) 2 次元球面調和函数チェビシェフスペクトルデータ real(8), dimension((nm+1)*(nm+1),0:lm) :: wq_QOperatorMPI_wq !(out) Q 演算子を作用された 2 次元スペクトルデータ wq_QOperatorMPI_wq = & wq_xvr(xvr_KGrad_wq(wq) - xvr_KGrad_wq(wq_L2_wq(wq))/2) & - wq_L2_wq(wq_xvr(xvr_KGrad_wq(wq)))/2 end function wq_QOperatorMPI_wq function wr_RadRot_xvr_xvr(xvr_VLON,xvr_VLAT) ! r・(▽×v) ! ! ベクトルの渦度と動径ベクトルの内積 r・(▽×v) を計算する. ! ! 第 1, 2 引数(v[λ], v[φ])がそれぞれベクトルの経度成分, 緯度成分を表す. ! ! r・(▽×v) = 1/cosφ・∂v[φ]/∂λ - 1/cosφ・∂(v[λ] cosφ)/∂φ ! ! のスペクトル データが返される. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_VLON !(in) ベクトルの経度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_VLAT !(in) ベクトルの緯度成分 real(8), dimension((nm+1)*(nm+1),km) :: wr_RadRot_xvr_xvr !(out) ベクトルの渦度と動径ベクトルの内積 wr_RadRot_xvr_xvr = wa_DivLon_xva(xvr_VLAT) - wa_DivLat_xva(xvr_VLON) end function wr_RadRot_xvr_xvr function wr_RadRotRot_xvr_xvr_xvr(xvr_VLON,xvr_VLAT,xvr_VRAD) ! ! ベクトル v に対して r・(▽×▽×v) を計算する. ! ! 第 1, 2, 3 引数(v[λ], v[φ], v[r])がそれぞれベクトルの経度成分, ! 緯度成分, 動径成分を表す. ! ! r・(▽×▽×v) = 1/r ∂/∂r (r・( 1/cosφ・∂v[λ]/∂λ ! + 1/cosφ・∂(v[φ] cosφ)/∂φ ) ) ! + L^2 v[r]/r ! ! のスペクトルデータが返される. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_VLON !(in) ベクトルの経度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_VLAT !(in) ベクトルの緯度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_VRAD !(in) ベクトルの動径成分 real(8), dimension((nm+1)*(nm+1),km) :: wr_RadRotRot_xvr_xvr_xvr !(out) ベクトル v の r・(▽×▽×v) wr_RadRotRot_xvr_xvr_xvr = & wr_RotDRad_wr( & wa_DivLon_xva(xvr_VLON)+ wa_DivLat_xva(xvr_VLAT)) & - wa_Lapla_wa(wr_xvr(xvr_VRAD/xvr_RAD)) end function wr_RadRotRot_xvr_xvr_xvr function wq_RadRotRot_xvr_xvr_xvr(xvr_VLON,xvr_VLAT,xvr_VRAD) ! ! ベクトル v に対して r・(▽×▽×v) を計算する. ! ! 第 1, 2, 3 引数(v[λ], v[φ], v[r])がそれぞれベクトルの経度成分, ! 緯度成分, 動径成分を表す. ! ! r・(▽×▽×v) = 1/r ∂/∂r (r・( 1/cosφ・∂v[λ]/∂λ ! + 1/cosφ・∂(v[φ] cosφ)/∂φ ) ) ! + L^2 v[r]/r ! ! のスペクトルデータが返される. ! real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_VLON !(in) ベクトルの経度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_VLAT !(in) ベクトルの緯度成分 real(8), dimension(0:im-1,1:jc,km), intent(in) :: xvr_VRAD !(in) ベクトルの動径成分 real(8), dimension((nm+1)*(nm+1),0:lm) :: wq_RadRotRot_xvr_xvr_xvr !(out) ベクトル v の r・(▽×▽×v) wq_RadRotRot_xvr_xvr_xvr = & wq_RotDRad_wr( & wa_DivLon_xva(xvr_VLON)+ wa_DivLat_xva(xvr_VLAT)) & - wa_Lapla_wa(wq_xvr(xvr_VRAD/xvr_RAD)) end function wq_RadRotRot_xvr_xvr_xvr subroutine wq_Potential2VectorMPI(& xvr_VLON,xvr_VLAT,xvr_VRAD,wq_TORPOT,wq_POLPOT) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψr) + ▽x▽x(Φr) ! ! の各成分を計算する ! real(8), dimension(0:im-1,1:jc,km) :: xvr_VLON !(out) ベクトル場の経度成分 real(8), dimension(0:im-1,1:jc,km) :: xvr_VLAT !(out) ベクトル場の緯度成分 real(8), dimension(0:im-1,1:jc,km) :: xvr_VRAD !(out) ベクトル場の動径成分 real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq_TORPOT !(in) トロイダルポテンシャル real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq_POLPOT !(in) ポロイダルポテンシャル xvr_VLON = xvr_RAD * xvr_GradLat_wq(wq_TORPOT) & + xva_GradLon_wa(wr_RotDRad_wq(wq_POLPOT)) xvr_VLAT = - xvr_RAD * xvr_GradLon_wq(wq_TORPOT) & + xva_GradLat_wa(wr_RotDRad_wq(wq_POLPOT)) xvr_VRAD = xvr_wq(wq_L2_wq(wq_POLPOT))/xvr_RAD end subroutine wq_Potential2VectorMPI subroutine wq_Potential2RotationMPI(& xvr_RotVLON,xvr_RotVLAT,xvr_RotVRAD,wq_TORPOT,wq_POLPOT) ! ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場 ! ! v = ▽x(Ψr) + ▽x▽x(Φr) ! ! に対して, その回転 ! ! ▽xv = ▽x▽x(Ψr) + ▽x▽x▽x(Φr) = ▽x▽x(Ψr) - ▽x((▽^2Φ)r) ! ! を計算する. ! ベクトル場の回転 real(8), dimension(0:im-1,1:jc,km), intent(OUT) :: xvr_RotVLON !(out) 回転の経度成分 real(8), dimension(0:im-1,1:jc,km), intent(OUT) :: xvr_RotVLAT !(out) 回転の緯度成分 real(8), dimension(0:im-1,1:jc,km), intent(OUT) :: xvr_RotVRAD !(out) 回転の動径成分 ! 入力ベクトル場を表すポテンシャル real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq_TORPOT !(in) トロイダルポテンシャル real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wq_POLPOT !(in) ポロイダルポテンシャル call wq_Potential2VectorMPI( & xvr_RotVLON,xvr_RotVLAT,xvr_RotVRAD, & -wq_Lapla_wq(wq_POLPOT), wq_TORPOT) end subroutine wq_Potential2RotationMPI end module wq_mpi_module