!--
!----------------------------------------------------------------------
! Copyright(c) 2008-2012 SPMDODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  wq_mpi_module_sjpack
!
!    spml/wq_mpi_module_sjpack モジュールは球内での流体運動をスペクトル法と 
!    MPI 並列化によって数値計算するための Fortran90 関数を提供する
!    ものである. 
!
!    水平方向に球面調和函数変換および動径方向に
!    Matsushima and Marcus (1994) で提唱された多項式を用いた
!    スペクトル計算のためのさまざまな関数を提供する. 
!
!    内部で wa_mpi_module_sjpack, aq_module を用いている. 
!    最下部では球面調和変換変換のエンジンとして 
!    ISPACK の Fortran77 サブルーチンを用いている.
!
!    関数, サブルーチンの名前と機能は wq_mpi_module のものと同じである.
!    したがって use 文を wq_mpi_module から wq_mpi_module_sjpack に
!    変更するだけで SJPACK-MPI の機能が使えるようになる.
!
!    ただし l_nm, nm_l の使い方には注意されたい. wq_module の l_nm,
!    nm_l はwq_Initial で初期化しなくとも用いることができる(結果が切断
!    波数に依らない)が, wq_mpi_module_sjpack のものは初期化したのちにしか使
!    うことができない.
!
!    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/04/02  竹広真一  sjpack 化
!      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_sjpack
  !
  != wq_mpi_module_sjpack
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/wq_mpi_module_sjpack モジュールは球内での流体運動をスペクトル法と 
  ! MPI 並列化によって数値計算するための Fortran90 関数を提供する
  ! ものである. 
  !
  ! 水平方向に球面調和函数変換および動径方向に
  ! Matsushima and Marcus (1994) で提唱された多項式を用いた
  ! スペクトル計算のためのさまざまな関数を提供する. 
  !
  ! 内部で wa_mpi_module_sjpack, 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_sjpack
  use w_deriv_mpi_module_sjpack
  use w_integral_mpi_module_sjpack
  use wa_base_mpi_module_sjpack
  use wa_deriv_mpi_module_sjpack
  use wq_module_sjpack
  use wa_module_sjpack, 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_sjpack (2012/04/02) 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_sjpack
