!--
!----------------------------------------------------------------------
! Copyright (c) 2017 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  wsc_mpi_module_supack
!
!    spml/wsc_mpi_module_supack モジュールは球面上および球殻内での流体運動を
!    スペクトル法と MPI 並列化によって数値計算するための Fortran90 
!    関数を提供するものである. 
!
!    水平方向に球面調和函数変換および鉛直方向に sin/cos 関数変換を
!    用いる場合のスペクトル計算のためのさまざまな関数を提供する. 
!
!    関数, サブルーチンの名前と機能は wsc_mpi_module のものと同じである.
!    したがって use 文を wsc_mpi_module から wsc_mpi_module_supack に
!    変更するだけで SUPACK-MPI の機能が使えるようになる.
!
!    ただし l_nm, nm_l の使い方には注意されたい. wt_mpi_module の l_nm, nm_l は
!    wt_Initial で初期化しなくとも用いることができる(結果が切断波数に依らない)が,
!    wsc_mpi_module_supack のものは初期化したのちにしか使うことができない.
!
!    内部で wa_mpi_module_supack を用いている. 
!    最下部では球面調和変換およびチェビシェフ変換のエンジンとして 
!    ISPACK2 の Fortran77 サブルーチンを用いている.
!
!
!履歴  2017/05/30  竹広真一  wt_mpi_module_supack より wsc 用に改造
!                             
!凡例
!      データ種類と index
!        x : 経度         y : 緯度    v : 緯度(分散格子)     z : 動径
!        w : 球面調和関数スペクトル
!        n : 球面調和関数スペクトル(水平全波数)
!        m : 球面調和関数スペクトル(帯状波数)
!        s : sin 関数スペクトル
!        c : cos 関数スペクトル
!        a : 任意の次元
!
!        xyz : 3 次元格子点データ
!        xvz : 3 次元分散格子点データ
!        xy  : 水平 2 次元格子点データ
!        yz  : 子午面 2 次元格子点データ
!        xz  : 緯度面 2 次元格子点データ
!
!        wz  : 水平スペクトル動径格子点データ
!        ws  : スペクトルデータ(鉛直 sin)
!        wc  : スペクトルデータ(鉛直 cos)
!
!++
module wsc_mpi_module_supack
  !
  != wsc_mpi_module_supack
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id: wt_mpi_module_supack.f90 598 2013-08-20 03:23:44Z takepiro $
  ! Copyright&License:: See COPYRIGHT[link:../../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/wt_mpi_module_supack モジュールは球面上および球殻内での流体運動を
  ! スペクトル法と MPI 並列化によって数値計算するための Fortran90 
  ! 関数を提供するものである. 
  !
  ! 水平方向に球面調和函数変換および鉛直方向に sin/cos 関数変換を
  ! 用いる場合のスペクトル計算のためのさまざまな関数を提供する. 
  !
  ! 内部で wa_mpi_module_supack を用いている. 
  ! 最下部では球面調和変換およびチェビシェフ変換のエンジンとして 
  ! ISPACK2 の Fortran77 サブルーチンを用いている.
  !
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (ws_, wc_, nmz_,nz_,xyz_,wz_,w_,xy_,x_,y_,s_,c_,a_) は, 
  !   返す値の形を示している.
  !   ws_,wc_ :: スペクトルデータ(球面調和函数・三角関数変換)
  !   nmz_    :: 水平スペクトルデータ(全波数 n, 帯状波数各成分, 動径)
  !   nz_     :: 水平スペクトルデータ(全波数 n, 動径)
  !   xyz_    :: 3 次元格子点データ(経度・緯度・動径)
  !   wz_     :: 水平スペクトル, 動径格子点データ
  !
  ! * 関数名の間の文字列(DLon, GradLat, GradLat, DivLon, DivLat, Lapla,..)
  !   は, その関数の作用を表している.
  !
  ! * 関数名の最後 (ws_,wc_,xyz_,wz_,w_,xy_,x_,y_,z_,a_) は, 入力変数の
  !   形がスペクトルデータおよび格子点データであることを示している.
  !   _ws,_wc  :: スペクトルデータ
  !   _xyz     :: 3 次元格子点データ
  !   _xyz_xyz :: 2 つの3 次元格子点データ, ...
  !
  !=== 各データの種類の説明
  !
  ! * xyz : 3 次元格子点データ(経度・緯度・動径)
  !   * 変数の種類と次元は real(8), dimension(0:im-1,jm,0:km). 
  !   * im, jm, km はそれぞれ経度, 緯度, 動径座標の格子点数であり, 
  !     サブルーチン wt_mpi_Initial にてあらかじめ設定しておく.
  !
  ! * xvz : 3 次元分散格子点データ(経度・緯度・動径)
  !   * 変数の種類と次元は real(8), dimension(0:im-1,jc,0:km). 
  !   * im, km はそれぞれ経度, 動径座標の格子点数であり, 
  !   * jc はこのプロセスで保有する緯度格子点数である. 
  !     サブルーチン wsc_mpi_Initial を呼ぶと jc が設定される. 
  !
  ! * ws : スペクトルデータ(鉛直 sin)
  !   * 変数の種類と次元は real(8), dimension(nc,lm). 
  !   * nm は球面調和函数の最大全波数, lm は sin 変換の最大次数であり,
  !     サブルーチン wsc_mpi_Initial にてあらかじめ設定しておく. 
  !   * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる
  !     ことができる.
  !
  ! * wc : スペクトルデータ(鉛直 cos)
  !   * 変数の種類と次元は real(8), dimension(nc,0:lm). 
  !   * nm は球面調和函数の最大全波数, lm は cos 変換の最大次数であり,
  !     サブルーチン wsc_mpi_Initial にてあらかじめ設定しておく. 
  !   * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる
  !     ことができる.
  !
  ! * nmz : 水平スペクトルデータの並んだ 3 次元配列.
  !   * 変数の種類と次元は real(8), dimension(0:nm,-nm:nm,0:km). 
  !   * 第 1 次元が水平全波数, 第 2 次元が帯状波数, 第 3 次元が動径座標を表す. 
  !   * nm は球面調和函数の最大全波数であり, サブルーチン wt_Initial にて
  !     あらかじめ設定しておく.
  !
  ! * nz : スペクトルデータの並んだ 2 次元配列.
  !   * 変数の種類と次元は real(8), dimension(0:nm,0:km). 
  !     第 1 次元が水平全波数を表す. 
  !   * nm は球面調和函数の最大全波数であり, 
  !     サブルーチン wt_Initial にてあらかじめ設定しておく.
  !
  ! * wz : 水平スペクトル, 動径格子点データ.
  !   * 変数の種類と次元は real(8), dimension(nc,0:km).
  !
  ! * ws_, wc_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * xyz_ で始まる関数が返す値は 3 次元格子点データに同じ.
  !
  ! * xvz_ で始まる関数が返す値は 3 次元分散格子点データに同じ.
  !
  ! * wz_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化 
  !
  ! wsc_mpi_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  ! 
  !==== 座標変数
  !
  ! x_Lon, x_Lon_Weight,  ::  格子点座標(経度)と重みを格納した 1 次元配列
  ! y_Lat, y_Lat_Weight   ::  格子点座標(緯度)と重みを格納した 1 次元配列
  ! v_Lat, v_Lat_Weight   ::  分散格子点座標(緯度)と重みを格納した 1 次元配列
  ! z_Rad, z_Rad_Weight   ::  格子点座標(同形)と重みを格納した 1 次元配列
  ! xyz_Lon, xyz_Lat, xyz_Rad :: 格子点データの経度・緯度・動径座標(X,Y,Z) (格子点データ型 3 次元配列)
  ! xvz_Lon, xvz_Lat, xvz_Rad :: 格子点データの経度・緯度・動径座標(X,Y,Z) (分散格子点データ型 3 次元配列)
  !
  !==== 基本変換
  !
  ! xyz_ws, ws_xyz :: スペクトルデータと 3 次元格子データの間の変換
  !                   (球面調和函数, sin 変換)
  ! xyz_wc, wc_xyz :: スペクトルデータと 3 次元格子データの間の変換
  !                   (球面調和函数, cos 変換)
  ! xvz_ws, ws_xvz :: スペクトルデータと 3 次元分散格子データの間の変換
  !                   (球面調和函数, sin 変換)
  ! xvz_wc, wc_xvz :: スペクトルデータと 3 次元分散格子データの間の変換
  !                   (球面調和函数, cos 変換)
  ! xyz_wz, wz_xyz :: 3 次元格子データと水平スペクトル・動径格子データとの間
  !                   の変換 (球面調和函数)
  ! xvz_wz, wz_xvz :: 3 次元分散格子データと水平スペクトル・動径格子データとの間
  !                   の変換 (球面調和函数)
  ! wz_ws, ws_wz   :: スペクトルデータと水平スペクトル・動径格子データとの間
  !                   の変換 (sin 変換)
  ! wz_wc, wc_wz   :: スペクトルデータと水平スペクトル・動径格子データとの間
  !                   の変換 (cos 変換)
  ! w_xy, xy_w     :: スペクトルデータと 2 次元水平格子データの間の変換
  !                   (球面調和函数変換) 
  ! w_xv, xv_w     :: スペクトルデータと 2 次元水平分散格子データの間の変換
  !                   (球面調和函数変換) 
  ! az_as, as_az   :: 格子データと sin 変換スペクトルデータの間の変換を
  !                   同時に複数個行う
  ! az_as, as_az   :: 格子データと cos 変換スペクトル間の変換を
  !                   同時に複数個行う
  ! l_nm, nm_l     :: スペクトルデータの格納位置と全波数・帯状波数の変換 
  !
  ! ws_DRad_wc, wc_DRad_ws     :: スペクトルデータに動径微分∂/∂r を作用させる
  ! ws_DivRad_wc, wc_DivRad_ws :: スペクトルデータに発散型動径微分
  !                               ∂/∂r 
  ! ws_RotRad_wc, wc_RotRad_ws :: スペクトルデータに回転型動径微分
  !                               ∂/∂r を作用させる
  ! ws_Lapla_ws, wc_Lapla_wc   :: スペクトルデータにラプラシアンを作用させる
  ! ws_LaplaInv_ws, wc_LaplaInv_wc  :: スペクトルデータに逆ラプラシアンを作用させる  
  ! xyz_GradLon_ws, xyz_GradLon_wc :: スペクトルデータに勾配型経度微分
  !                                   1/(a cosφ)・∂/∂λを作用させる
  ! xyz_GradLat_wt, xyz_GradLat_wc :: スペクトルデータに勾配型緯度微分
  !                                   1/a・∂/∂φを作用させる
  ! ws_DivLon_xyz, wc_DivLon_xyz   :: 格子データに発散型経度微分
  !                                   1/(a cosφ)・∂/∂λを作用させる
  ! ws_DivLat_xyz, wc_DivLat_xyz   :: 格子データに発散型緯度微分
  !                                   (1/a cosφ)・∂(g cosφ)/∂φを作用させる
  ! wc_Div_xyz_xyz_xyz  :: ベクトル成分である 3 つの格子データに
  !                        発散を作用させる
  ! ws_Div_xyz_xyz_xyz  :: ベクトル成分である 3 つの格子データに
  !                        発散を作用させる
  ! xyz_Div_xyz_xyz_xyz :: ベクトル成分である 3 つの格子データに
  !                        発散を作用させる
  ! wc_RotRad_xyz_xyz   :: ベクトル場の回転の動径成分を計算する
  !
  !==== トロイダルポロイダル計算用微分
  !
  ! ws_KxRGrad_ws, wc_KxRGrad_wc :: スペクトルデータに
  !                                 経度微分 k×r・▽ = ∂/∂λを作用させる
  ! ws_L2_ws, wc_L2_wc       :: スペクトルデータに 
  !                             L2 演算子 = -水平ラプラシアンを作用させる
  ! ws_L2Inv_ws, ws_L2Inv_ws :: スペクトルデータに 
  !                             L2 演算子の逆 = -逆水平ラプラシアンを作用させる
  ! wc_RadRot_xyz_xyz        :: ベクトル v の渦度と動径ベクトル r の内積
  !                             r・(▽×v) を計算する
  ! ws_RadRotRot_xyz_xyz_xyz :: ベクトルの v の r・(▽×▽×v) を計算する
  ! wsc_Potential2Vector      :: トロイダルポロイダルポテンシャルから
  !                              ベクトル場を計算する
  ! wsc_Potential2Rotation    :: トロイダルポロイダルポテンシャルで表される
  !                              非発散ベクトル場の回転の各成分を計算する
  !
  !==== ポロイダル/トロイダルモデル用スペクトル解析
  !
  ! nmz_ToroidalEnergySpectrum_wt, :: トロイダルポテンシャルからエネルギーの
  ! nz_ToroidalEnergySpectrum_wt   :: 球面調和函数各成分を計算する
  ! 
  ! nmz_PoloidalEnergySpectrum_wt, :: ポロイダルポテンシャルからエネルギーの
  ! nz_PoloidalEnergySpectrum_wt   :: 球面調和函数各成分を計算する
  !
  !==== 積分・平均(3 次元データ)
  !
  ! IntLonLatRad_xyz, AvrLonLatRad_xyz :: 3 次元格子点データの全領域積分および平均
  !
  ! z_IntLonLat_xyz, z_AvrLonLat_xyz :: 3 次元格子点データの緯度経度(水平・球面)積分および平均               
  !
  ! y_IntLonRad_xyz, y_AvrLonRad_xyz :: 3 次元格子点データの緯度動径積分および平均
  !
  ! z_IntLatRad_xyz, z_AvrLatRad_xyz :: 3 次元格子点データの経度動径(子午面)積分および平均              
  !
  ! yz_IntLon_xyz, yz_AvrLon_xyz :: 3 次元格子点データの経度方向積分および平均
  ! xz_IntLat_xyz, xz_AvrLat_xyz :: 3 次元格子点データの緯度方向積分および平均
  ! xz_IntRad_xyz, xz_AvrRad_xyz :: 3 次元格子点データの動径方向積分および平均
  !
  ! IntLonLatRad_xvz, AvrLonLatRad_xvz :: 3 次元格子点データの全領域積分および平均
  !
  ! z_IntLonLat_xvz, z_AvrLonLat_xvz :: 3 次元格子点データの緯度経度(水平・球面)積分および平均               
  !
  ! v_IntLonRad_xvz, v_AvrLonRad_xvz :: 3 次元格子点データの緯度動径積分および平均
  !
  ! z_IntLatRad_xvz, z_AvrLatRad_xvz :: 3 次元格子点データの経度動径(子午面)積分および平均              
  !
  ! vz_IntLon_xvz, vz_AvrLon_xvz :: 3 次元格子点データの経度方向積分および平均
  ! xz_IntLat_xvz, xz_AvrLat_xvz :: 3 次元格子点データの緯度方向積分および平均
  ! xy_IntRad_xvz, xy_AvrRad_xvz :: 3 次元格子点データの動径方向積分および平均
  !
  !==== 積分・平均(2 次元データ)
  !
  ! IntLonLat_xy, AvrLonLat_xy :: 2 次元格子点データの水平(球面)積分および平均
  ! IntLonRad_xz, AvrLonRad_xz :: 2 次元(XZ)格子点データの経度動径積分
  !                            :: および平均
  ! IntLatRad_yz, AvrLatRad_yz :: 2 次元(YZ)格子点データの緯度動径(子午面)
  !                            :: 積分および平均 
  ! y_IntLon_xy, y_AvrLon_xy   :: 水平 2 次元(球面)格子点データの経度方向
  !                            :: 積分および平均
  ! x_IntLat_xy, x_AvrLat_xy   :: 水平2 次元(球面)格子点データの緯度方向積分
  !                            :: および平均
  ! z_IntLon_xz, z_AvrLon_xz   :: 2 次元(XZ)格子点データの経度方向積分および
  !                            :: 平均
  ! x_IntRad_xz, x_AvrRad_xz   :: 2 次元(XZ)格子点データの動径方向積分および
  !                            :: 平均
  ! z_IntLat_yz, z_AvrLat_yz   :: 2 次元(YZ)格子点データの緯度方向積分および
  !                            :: 平均
  ! y_IntRad_yz, y_AvrRad_yz   :: 2 次元(YZ)格子点データの動径方向積分および
  !                            :: 平均                  
  !
  ! IntLonLat_xv, AvrLonLat_xv :: 2 次元格子点データの水平(球面)積分および平均
  ! IntLonRad_xz, AvrLonRad_xz :: 2 次元(XZ)格子点データの経度動径積分
  !                            :: および平均
  ! IntLatRad_vz, AvrLatRad_vz :: 2 次元(YZ)格子点データの緯度動径(子午面)
  !                            :: 積分および平均 
  ! v_IntLon_xv, v_AvrLon_xv   :: 水平 2 次元(球面)格子点データの経度方向
  !                            :: 積分および平均
  ! v_IntLat_xv, x_AvrLat_xv   :: 水平2 次元(球面)格子点データの緯度方向積分
  !                            :: および平均
  ! z_IntLon_xz, z_AvrLon_xz   :: 2 次元(XZ)格子点データの経度方向積分および
  !                            :: 平均
  ! x_IntRad_xz, x_AvrRad_xz   :: 2 次元(XZ)格子点データの動径方向積分および
  !                            :: 平均
  ! z_IntLat_vz, z_AvrLat_vz   :: 2 次元(YZ)格子点データの緯度方向積分および
  !                            :: 平均
  ! v_IntRad_vz, v_AvrRad_vz   :: 2 次元(YZ)格子点データの動径方向積分および
  !                            :: 平均                  
  !
  !==== 積分・平均(1 次元データ)
  !
  ! IntLon_x, AvrLon_x  :: 1 次元(X)格子点データの経度方向積分および平均
  ! IntLat_y, AvrLat_y  :: 1 次元(Y)格子点データの緯度方向積分および平均
  ! IntRad_z, AvrRad_z  :: 1 次元(Z)格子点データの動径方向積分および平均
  !
  ! IntLat_v, AvrLat_v  :: 1 次元(Y)格子点データの緯度方向積分および平均
  !
  !==== 補間計算
  !
  ! Interpolate_ws, Interpolate_wc :: スペクトルデータから任意の点の値を補間する  ! 
  use dc_message
!  use lumatrix
  use mpi
  use wa_mpi_module_supack
  use asc_module, z_RAD => g_X, z_RAD_WEIGHT => g_X_WEIGHT, &
                  as_az => as_ag, az_as => ag_as, &
                  ac_az => ac_ag, az_ac => ag_ac, &
                  s_z => s_g, z_s => g_s, c_z => c_g, z_c => g_c, &
                  s_Dr_c => s_Dx_c, as_Dr_ac => as_Dx_ac, &
                  c_Dr_s => s_Dx_c, ac_Dr_as => ac_Dx_as
  implicit none
  integer :: ierr

  private

  public wsc_mpi_Initial
  public jc
  public nc

  public x_Lon, x_Lon_Weight
  public v_Lat, v_Lat_Weight
  public z_Rad, z_Rad_Weight
  public l_nm, nm_l
  public xv_Lon, xv_Lat
  public xvz_Lon, xvz_Lat, xvz_Rad
  public wz_Rad
  public wsc_VMiss

  public w_xv, xv_w
  public az_ac, ac_az, az_as, as_az
  public wz_wc, wc_wz, wz_ws, ws_wz
  public xvz_wc, wc_xvz, xvz_ws, ws_xvz, xvz_wz, wz_xvz
  public as_Dr_ac, s_Dr_c, ac_Dr_as, c_Dr_s
  public ws_DRad_wc, ws_DivRad_wc, ws_RotRad_wc, wc_Lapla_wc, wc_LaplaInv_wc
  public wc_DRad_ws, wc_DivRad_ws, wc_RotRad_ws, ws_Lapla_ws, ws_LaplaInv_ws
  public xvz_GradLon_wc, xvz_Gradlat_wc, xvz_GradLon_ws, xvz_Gradlat_ws
  public wc_DivLon_xvz, ws_DivLat_xvz
  public wc_Div_xvz_xvz_xvz, ws_Div_xvz_xvz_xvz, xvz_Div_xvz_xvz_xvz
  public wc_RotRad_xvz_xvz

!!$  public xvz_RotLon_wt_wt, xvz_RotLat_wt_wt, wt_RotRad_xvz_xvz

  public vz_IntLon_xvz, xz_IntLat_xvz, xv_IntRad_xvz
  public x_IntLatRad_xvz, v_IntLonRad_xvz, z_IntLonLat_xvz
  public IntLonLatRad_xvz

  public z_IntLon_xz, x_IntRad_xz, IntLonRad_xz

  public x_IntLat_xv, v_IntLon_xv, IntLonLat_xv
  public z_IntLat_vz, v_IntRad_vz, IntLatRad_vz
  public IntLon_x, IntLat_v, IntRad_z

  public vz_AvrLon_xvz, xz_AvrLat_xvz, xv_AvrRad_xvz
  public x_AvrLatRad_xvz, v_AvrLonRad_xvz, z_AvrLonLat_xvz
  public AvrLonLatRad_xvz

  public z_AvrLon_xz, x_AvrRad_xz, AvrLonRad_xz

  public x_AvrLat_xv, v_AvrLon_xv, AvrLonLat_xv
  public z_AvrLat_vz, v_AvrRad_vz, AvrLatRad_vz
  public AvrLon_x, AvrLat_v, AvrRad_z

  public wc_KxRGrad_wc, ws_KxRGrad_ws
  public wc_L2_wc, wc_L2Inv_wc, ws_L2_ws, ws_L2Inv_ws
  public wc_RadRot_xvz_xvz, ws_RadRotRot_xvz_xvz_xvz
  public wsc_Potential2vectorMPI, wsc_Potential2RotationMPI

  public Interpolate_ws, Interpolate_wc

  public nmz_ToroidalEnergySpectrum_wc, nz_ToroidalEnergySpectrum_wc
  public nmz_PoloidalEnergySpectrum_ws, nz_PoloidalEnergySpectrum_ws


  integer            :: im=64, jm=32, km=16  ! 格子点の設定(経度, 緯度, 動径)
  integer            :: nm=21, lm=16         ! 切断波数の設定(水平, 動径)
  real(8)            :: ri=0.0, ro=1.0       ! 球殻内外半径
  real(8), parameter :: pi=3.1415926535897932385D0

  real(8), dimension(:,:,:), allocatable :: xvz_LON, xvz_LAT, xvz_RAD ! 座標
  real(8), dimension(:,:), allocatable   :: wz_RAD                    ! 座標

  real(8)             :: Lat_Weight_Sum

  real(8) :: wsc_VMiss = -999.0        ! 欠損値
  
  save im, jm, km, nm, lm, ri, ro, Lat_Weight_Sum, wsc_VMiss

contains
  !--------------- 初期化 -----------------
  subroutine wsc_mpi_Initial(i,j,k,n,l,r_in,r_out,np,wa_init)
    !
    ! スペクトル変換の格子点数, 波数, 動径座標の範囲を設定する.
    !
    ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
    ! しなければならない. 
    !
    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_in           ! 球殻内半径
    real(8),intent(in) :: r_out          ! 球殻外半径

    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
    ri = r_in ; ro = r_out

    if ( present(wa_init) ) then
       wa_initialize = wa_init
    else
       wa_initialize = .true.
    endif

    if ( present(np) ) then
       call MessageNotify('M','wt_mpi_Initial', &
            'Optional argument "NP_IN" is dummy in wt_mpi_module_supack')
    endif

    if ( wa_initialize ) then
       call wa_mpi_Initial(nm,im,jm,km+1)
    endif

    call asc_Initial(km,lm,r_in,r_out)

    allocate(xvz_Lon(0:im-1,jc,0:km))
    allocate(xvz_Lat(0:im-1,jc,0:km))
    allocate(xvz_Rad(0:im-1,jc,0:km))

    allocate(wz_Rad(nc,0:km))

    wz_Rad = spread(z_Rad,1,nc)

    z_Rad_Weight = z_Rad_Weight * Ro**2           ! ro^2 dr の積分重み

    xvz_Lon = spread(xv_Lon,3,km+1)
    xvz_Lat = spread(xv_Lat,3,km+1)
    xvz_Rad = spread(spread(z_Rad,1,jc),1,im)

    CALL MPI_ALLREDUCE(sum(v_Lat_weight),Lat_Weight_Sum,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_WORLD,IERR)


    call MessageNotify('M','wsc_mpi_initial', &
         'wsc_mpi_module_supack (2017/05/30) is initialized')

  end subroutine wsc_mpi_Initial

  !--------------- 基本変換 -----------------

  function xvz_wc(wc)
    !
    ! スペクトルデータから 3 次元格子点データへ(逆)変換する.
    !
    real(8), dimension(0:im-1,jc,0:km)      :: xvz_wc
    !(out) 3 次元経度緯度動径格子点データ

    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    xvz_wc = xva_wa(wz_wc(wc))

  end function xvz_wc

  function wc_xvz(xvz)
    !
    ! 3 次元格子点データからスペクトルデータへ(正)変換する.
    !
    real(8), dimension(nc,0:lm)                    :: wc_xvz
    !(out) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    wc_xvz = wc_wz(wa_xva(xvz))

  end function wc_xvz

  function xvz_ws(ws)
    !
    ! スペクトルデータから 3 次元格子点データへ(逆)変換する.
    !
    real(8), dimension(0:im-1,jc,0:km)      :: xvz_ws
    !(out) 3 次元経度緯度動径格子点データ

    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ

    xvz_ws = xva_wa(wz_ws(ws))

  end function xvz_ws

  function ws_xvz(xvz)
    !
    ! 3 次元格子点データからスペクトルデータへ(正)変換する.
    !
    real(8), dimension(nc,lm)                    :: ws_xvz
    !(out) 2 次元球面調和函数 sin スペクトルデータ

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    ws_xvz = ws_wz(wa_xva(xvz))

  end function ws_xvz

  function xvz_wz(wz)
    !
    ! 水平スペクトル・動径格子点データから 3 次元格子点データへ(逆)変換する.
    !
    real(8), dimension(0:im-1,jc,0:km)                     :: xvz_wz
    !(out) 3 次元経度緯度動径格子点データ

    real(8), dimension(nc,0:km), intent(in) :: wz
    !(in) 2 次元球面調和函数スペクトル・動径格子点データ

    xvz_wz = xva_wa(wz)

  end function xvz_wz

  function wz_xvz(xvz)
    !
    ! 3 次元格子データから水平スペクトル・動径格子点データへ(正)変換する.
    !
    real(8), dimension(nc,0:km)             :: wz_xvz
    !(out) 2 次元球面調和函数スペクトル・動径格子点データ

    real(8), dimension(0:im-1,jc,0:km), intent(in)         :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    wz_xvz = wa_xva(xvz)

  end function wz_xvz

  function wz_wc(wc)
    !
    ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ
    real(8), dimension(nc,0:km)             :: wz_wc
    !(out) 2 次元球面調和函数スペクトル・動径格子点データ

    wz_wc = az_ac(wc)

  end function wz_wc

  function wc_wz(wz)
    !
    ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
    !
    real(8), dimension(nc,0:km), intent(in) :: wz
    !(in) 2 次元球面調和函数スペクトル・動径格子点データ
    real(8), dimension(nc,0:lm)             :: wc_wz
    !(out) 2 次元球面調和函数 cos スペクトルデータ

    wc_wz = ac_az(wz)

  end function wc_wz

  function wz_ws(ws)
    !
    ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ
    real(8), dimension(nc,0:km)           :: wz_ws
    !(out) 2 次元球面調和函数スペクトル・動径格子点データ

    wz_ws = az_as(ws)

  end function wz_ws

  function ws_wz(wz)
    !
    ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
    !
    real(8), dimension(nc,0:km), intent(in) :: wz
    !(in) 2 次元球面調和函数スペクトル・動径格子点データ
    real(8), dimension(nc,lm)             :: ws_wz
    !(out) 2 次元球面調和函数 sin スペクトルデータ

    ws_wz = as_az(wz)

  end function ws_wz

  !--------------- 微分計算 -----------------
  function ws_DRad_wc(wc)
    !
    ! 入力スペクトルデータに動径微分 ∂/∂r を作用する.
    !
    ! スペクトルデータの動径微分とは, 対応する格子点データに動径微分を
    ! 作用させたデータのスペクトル変換のことである.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(nc,lm)               :: ws_DRad_wc
    !(in) 動径微分された2 次元球面調和函数 sin スペクトルデータ

    ws_DRad_wc = as_Dr_ac(wc)

  end function ws_DRad_wc

  function wc_DRad_ws(ws)
    !
    ! 入力スペクトルデータに動径微分 ∂/∂r を作用する.
    !
    ! スペクトルデータの動径微分とは, 対応する格子点データに動径微分を
    ! 作用させたデータのスペクトル変換のことである.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ

    real(8), dimension(nc,0:lm)           :: wc_DRad_ws
    !(in) 動径微分された2 次元球面調和函数 cos スペクトルデータ

    wc_DRad_ws = ac_Dr_as(ws)

  end function wc_DRad_ws

  function ws_DivRad_wc(wc)
    !
    ! 入力スペクトルデータに近似的発散型動径微分 ∂/∂r を作用する.
    !
    ! スペクトルデータの発散型動径微分とは, 対応する格子点データに
    ! 発散型動径微分を作用させたデータのスペクトル変換のことである.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(nc,lm)               :: ws_DivRad_wc
    !(out) 発散型動径微分を作用された 2 次元スペクトルデータ

    ws_DivRad_wc = ws_Drad_wc(wc)

  end function ws_DivRad_wc

  function wc_DivRad_ws(ws)
    !
    ! 入力スペクトルデータに近似的発散型動径微分 ∂/∂r を作用する.
    !
    ! スペクトルデータの発散型動径微分とは, 対応する格子点データに
    ! 発散型動径微分を作用させたデータのスペクトル変換のことである.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ

    real(8), dimension(nc,0:lm)           :: wc_DivRad_ws
    !(out) 発散型動径微分を作用された 2 次元スペクトルデータ

    wc_DivRad_ws = wc_Drad_ws(ws)

  end function wc_DivRad_ws

  function ws_RotRad_wc(wc)
    !
    ! 入力スペクトルデータに近似的回転型動径微分 ∂/∂r を作用する.
    !
    ! スペクトルデータの回転型動径微分とは, 対応する格子点データに
    ! 回転型動径微分を作用させたデータのスペクトル変換のことである.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(nc,lm)               :: ws_RotRad_wc
    !(out) 回転型動径微分を作用された 2 次元スペクトルデータ

    ws_RotRad_wc = ws_Drad_wc(wc)

  end function ws_RotRad_wc

  function wc_RotRad_ws(ws)
    !
    ! 入力スペクトルデータに近似的回転型動径微分 ∂/∂r を作用する.
    !
    ! スペクトルデータの回転型動径微分とは, 対応する格子点データに
    ! 回転型動径微分を作用させたデータのスペクトル変換のことである.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ

    real(8), dimension(nc,0:lm)           :: wc_RotRad_ws
    !(out) 回転型動径微分を作用された 2 次元スペクトルデータ

    wc_RotRad_ws = wc_Drad_ws(ws)

  end function wc_RotRad_ws

  function wc_Lapla_wc(wc)
    !
    ! 入力スペクトルデータに近似的ラプラシアン
    !
    !     ▽^2 =   1/a^2 cos^2φ・∂^2/∂λ^2 
    !            + 1/a^2 cosφ・∂/∂φ(cosφ∂/∂φ) 
    !            + ∂^2/∂r^2
    !
    ! を作用する.
    !
    ! スペクトルデータのラプラシアンとは, 対応する格子点データに
    ! ラプラシアンを作用させたデータのスペクトル変換のことである.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(nc,0:lm)             :: wc_Lapla_wc
    !(out) ラプラシアンを作用された 2 次元スペクトルデータ

    wc_Lapla_wc = wc_DivRad_ws(ws_Drad_wc(wc)) + wa_Lapla_wa(wc)/Ro**2

  end function wc_Lapla_wc

  function ws_Lapla_ws(ws)
    !
    ! 入力スペクトルデータに近似的ラプラシアン
    !
    !     ▽^2 =   1/a^2 cos^2φ・∂^2/∂λ^2 
    !            + 1/a^2 cosφ・∂/∂φ(cosφ∂/∂φ) 
    !            + ∂^2/∂r^2
    !
    ! を作用する.
    !
    ! スペクトルデータのラプラシアンとは, 対応する格子点データに
    ! ラプラシアンを作用させたデータのスペクトル変換のことである.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ

    real(8), dimension(nc,lm)             :: ws_Lapla_ws
    !(out) ラプラシアンを作用された 2 次元スペクトルデータ

    ws_Lapla_ws = ws_DivRad_wc(wc_Drad_ws(ws)) + wa_Lapla_wa(ws)/Ro**2

  end function ws_Lapla_ws

  function wc_LaplaInv_wc(wc)
    !
    ! スペクトルデータに逆ラプラシアンを作用させる. 
    !
    real(8), dimension(nc,0:lm),intent(in) :: wc
    real(8), dimension(nc,0:lm)            :: wc_LaplaInv_wc

    real(8), allocatable, dimension(:,:)   :: wc_LaplaInvFactor
    logical :: first=.true.
    integer :: n, m, l, nn, nmary(2)

    save first, wc_LaplaInvFactor

    if ( first ) then
       first = .false.
       allocate(wc_LaplaInvFactor(nc,0:lm))
       do nn=1,nc
          nmary=nm_l(nn)
          n = nmary(1) ; m = nmary(2)
          do l=0,lm
             if ( n==0 .and. l==0 ) then
                wc_LaplaInvFactor(nn,l) = 0.0D0
             else
                wc_LaplaInvFactor(nn,l) &
                     = -1.0D0/( n*(n+1)/Ro**2 + (PI/(Ro-Ri)*l)**2 )
             end if
          end do
       end do
    end if

    wc_LaplaInv_wc = wc_LaplaInvFactor * wc

  end function wc_LaplaInv_wc

  function ws_LaplaInv_ws(ws)
    !
    ! スペクトルデータに逆ラプラシアンを作用させる. 
    !
    real(8), dimension(nc,lm),intent(in) :: ws
    real(8), dimension(nc,lm)            :: ws_LaplaInv_ws

    real(8), allocatable, dimension(:,:) :: ws_LaplaInvFactor
    logical :: first=.true.
    integer :: n, m, l, nn, nmary(2)

    save first, ws_LaplaInvFactor

    if ( first ) then
       first = .false.
       allocate(ws_LaplaInvFactor(nc,lm))
       do nn=1,nc
          nmary=nm_l(nn)
          n = nmary(1) ; m = nmary(2)
          do l=1,lm
             ws_LaplaInvFactor(nn,l) &
                  = -1.0D0/( n*(n+1)/Ro**2 + (PI/(Ro-Ri)*l)**2 )
          end do
       end do
    end if

    ws_LaplaInv_ws = ws_LaplaInvFactor * ws

  end function ws_LaplaInv_ws

  function xvz_GradLon_wc(wc)
    !
    ! スペクトルデータに近似的勾配型経度微分 1/(a cosφ)・∂/∂λ
    ! を作用させる.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(0:im-1,jc,0:km)                     :: xvz_GradLon_wc
    !(out) 勾配型経度微分を作用された 2 次元スペクトルデータ

    xvz_GradLon_wc = xva_GradLon_wa(wz_wc(wc))/Ro

  end function xvz_GradLon_wc

  function xvz_GradLon_ws(ws)
    !
    ! スペクトルデータに近似的勾配型経度微分 1/(a cosφ)・∂/∂λ
    ! を作用させる.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(0:im-1,jc,0:km)    :: xvz_GradLon_ws
    !(out) 勾配型経度微分を作用された 2 次元スペクトルデータ

    xvz_GradLon_ws = xva_GradLon_wa(wz_ws(ws))/Ro

  end function xvz_GradLon_ws

  function xvz_GradLat_wc(wc) 
    !
    ! スペクトルデータに近似的勾配型経度微分 1/a ∂/∂φ を作用させる.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(0:im-1,jc,0:km)      :: xvz_GradLat_wc
    !(out) 勾配型緯度微分を作用された 2 次元スペクトルデータ

    xvz_GradLat_wc = xva_GradLat_wa(wz_wc(wc))/Ro
  end function xvz_GradLat_wc

  function xvz_GradLat_ws(ws) 
    !
    ! スペクトルデータに近似的勾配型経度微分 1/a ∂/∂φ を作用させる.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ

    real(8), dimension(0:im-1,jc,0:km)    :: xvz_GradLat_ws
    !(out) 勾配型緯度微分を作用された 2 次元スペクトルデータ

    xvz_GradLat_ws = xva_GradLat_wa(wz_ws(ws))/Ro
  end function xvz_GradLat_ws

  function wc_DivLon_xvz(xvz)
    ! 
    ! 格子点データに近似的発散型経度微分 1/(a cosφ)・∂/∂λ を作用させた
    ! スペクトルデータを返す.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(nc,0:lm)                    :: wc_DivLon_xvz
    !(out) 発散型経度微分を作用された 2 次元スペクトルデータ

    wc_DivLon_xvz = wc_wz(wa_DivLon_xva(xvz/Ro))
  end function wc_DivLon_xvz

  function wc_DivLat_xvz(xvz)
    !
    ! 格子データに近似的発散型緯度微分 1/(a cosφ)・∂(f cosφ)/∂φ を
    ! 作用させたスペクトルデータを返す.
    !
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in)   :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(nc,0:lm)       :: wc_DivLat_xvz
    !(out) 発散型緯度微分を作用された 2 次元スペクトルデータ

    wc_DivLat_xvz = wc_wz(wa_divlat_xva(xvz/Ro))
  end function wc_DivLat_xvz

  function ws_DivLon_xvz(xvz)
    ! 
    ! 格子点データに近似的発散型経度微分 1/(a cosφ)・∂/∂λ を作用させた
    ! スペクトルデータを返す.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(nc,lm)                      :: ws_DivLon_xvz
    !(out) 発散型経度微分を作用された 2 次元スペクトルデータ

    ws_DivLon_xvz = ws_wz(wa_DivLon_xva(xvz/Ro))
  end function ws_DivLon_xvz

  function ws_DivLat_xvz(xvz)
    !
    ! 格子データに近似的発散型緯度微分 1/(a cosφ)・∂(f cosφ)/∂φ を
    ! 作用させたスペクトルデータを返す.
    !
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(nc,lm)                      :: ws_DivLat_xvz
    !(out) 発散型緯度微分を作用された 2 次元スペクトルデータ

    ws_DivLat_xvz = ws_wz(wa_divlat_xva(xvz/Ro))
  end function ws_DivLat_xvz

  function wc_Div_xvz_xvz_xvz(xvz_Vlon,xvz_Vlat,xvz_Vrad)
    !
    ! べクトル成分である 3 つの格子データに近似的発散を作用させた
    ! スペクトルデータを返す.
    !
    ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, 
    ! 動径成分を表し, 近似的発散は 
    !
    !   1/(a cosφ)・∂u/∂λ + 1/(a cosφ)・∂(v cosφ)/∂φ + ∂w/∂r 
    !
    ! と計算される.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vlon
    !(in) ベクトル場の経度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vlat
    !(in) ベクトル場の緯度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vrad
    !(in) ベクトル場の動径成分

    real(8), dimension(nc,0:lm)     :: wc_Div_xvz_xvz_xvz
    !(out) ベクトル場の発散

    wc_Div_xvz_xvz_xvz =   wc_DivLon_xvz(xvz_Vlon) &
         + wc_DivLat_xvz(xvz_Vlat) &
         + wc_DivRad_ws(ws_xvz(xvz_Vrad))

  end function wc_Div_xvz_xvz_xvz

  function ws_Div_xvz_xvz_xvz(xvz_Vlon,xvz_Vlat,xvz_Vrad)
    !
    ! べクトル成分である 3 つの格子データに近似的発散を作用させた
    ! スペクトルデータを返す.
    !
    ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, 
    ! 動径成分を表し, 近似的発散は 
    !
    !   1/(a cosφ)・∂u/∂λ + 1/(a cosφ)・∂(v cosφ)/∂φ + ∂w/∂r 
    !
    ! と計算される.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vlon
    !(in) ベクトル場の経度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vlat
    !(in) ベクトル場の緯度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vrad
    !(in) ベクトル場の動径成分

    real(8), dimension(nc,lm)     :: ws_Div_xvz_xvz_xvz
    !(out) ベクトル場の発散

    ws_Div_xvz_xvz_xvz =   ws_DivLon_xvz(xvz_Vlon) &
         + ws_DivLat_xvz(xvz_Vlat) &
         + ws_DivRad_wc(wc_xvz(xvz_Vrad))

  end function ws_Div_xvz_xvz_xvz

  function xvz_Div_xvz_xvz_xvz(xvz_Vlon,xvz_Vlat,xvz_Vrad)
    !
    ! ベクトル成分である 3 つの格子データに近似的発散を作用させる.
    !
    ! 第 1, 2 ,3 引数(u,v,w)がそれぞれベクトルの経度成分, 緯度成分, 
    ! 動径成分を表す.
    !
    ! 極の特異性を回避するためにベクトル場に cosφ/r の重みをかけて
    ! 計算している. 
    !
    !      div V = (1/cosφ)・div (Vcosφ) + V_φtanφ/a
    ! 
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vlon
    !(in) ベクトル場の経度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vlat
    !(in) ベクトル場の緯度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vrad
    !(in) ベクトル場の動径成分

    real(8), dimension(0:im-1,jc,0:km)             :: xvz_Div_xvz_xvz_xvz
    !(out) ベクトル場の発散

    xvz_Div_xvz_xvz_xvz &
         = 1.0D0/cos(xvz_Lat) &
         * xvz_wc(wc_Div_xvz_xvz_xvz(xvz_VLon*cos(xvz_Lat),  &
         xvz_VLat*cos(xvz_Lat),  &
         xvz_VRad*cos(xvz_Lat))) &
         + xvz_VLat*tan(xvz_Lat)/Ro

  end function xvz_Div_xvz_xvz_xvz

!!$    function xvz_RotLon_wt_wt(wt_Vrad,wt_Vlat) 
!!$      !
!!$      ! ベクトル場の動径成分, 緯度成分である第 1, 2 引数 Vrad, Vlat から
!!$      ! 回転の経度成分 
!!$      !
!!$      !    1/r ∂Vrad/∂φ-1/r ∂(r Vlat)/∂r を計算する.
!!$      !
!!$      ! を計算する
!!$      !
!!$      real(8), dimension(nc,0:lm), intent(in) :: wt_Vrad
!!$      !(in) ベクトル場の動径成分
!!$
!!$      real(8), dimension(nc,0:lm), intent(in) :: wt_Vlat
!!$      !(in) ベクトル場の緯度成分
!!$
!!$      real(8), dimension(0:im-1,jc,0:km)                     :: xvz_RotLon_wt_wt
!!$      !(out) ベクトル場の回転の経度成分
!!$
!!$        xvz_RotLon_wt_wt =   xvz_GradLat_wt(wt_Vrad) &
!!$                           - xvz_wt(wt_RotRad_wt(wt_Vlat))
!!$
!!$    end function xvz_RotLon_wt_wt
!!$
!!$    function xvz_RotLat_wt_wt(wt_Vlon,wt_Vrad) 
!!$      !
!!$      ! ベクトル場の経度成分, 動径成分である第 1, 2 引数 Vlon, Vrad から
!!$      ! 回転の緯度成分 
!!$      !
!!$      !    1/r ∂(r Vlon)/∂r - 1/rcosφ・∂Vrad/∂λ
!!$      !
!!$      ! を計算する.
!!$      !
!!$      real(8), dimension(nc,0:lm), intent(in) :: wt_Vlon
!!$      !(in) ベクトル場の経度成分
!!$
!!$      real(8), dimension(nc,0:lm), intent(in) :: wt_Vrad
!!$      !(in) ベクトル場の動径成分
!!$
!!$      real(8), dimension(0:im-1,jc,0:km)                     :: xvz_RotLat_wt_wt
!!$      !(out) ベクトル場の回転の緯度成分
!!$
!!$        xvz_RotLat_wt_wt =   xvz_wt(wt_RotRad_wt(wt_Vlon)) &
!!$                           - xvz_GradLon_wt(wt_Vrad) 
!!$
!!$    end function xvz_RotLat_wt_wt

  function wc_RotRad_xvz_xvz(xvz_Vlat,xvz_Vlon) 
    !
    ! ベクトルの緯度成分, 経度成分である第 1, 2 引数 Vlat, Vlon に対して
    ! ベクトル場の近似的回転の動径成分 
    !
    !    1/(a cosφ)・∂Vlat/∂λ - 1/(a cosφ)・∂(Vlon cosφ)/∂φ
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vlat
    !(in) ベクトル場の緯度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_Vlon
    !(in) ベクトル場の経度成分

    real(8), dimension(nc,0:lm)     :: wc_RotRad_xvz_xvz
    !(out) ベクトル場の回転の動径成分

    wc_RotRad_xvz_xvz =   wc_DivLon_xvz(xvz_Vlat) &
         - wc_DivLat_xvz(xvz_Vlon)

  end function wc_RotRad_xvz_xvz

  !--------------- ポロイダル/トロイダルモデル用微分 -----------------
  function wc_KxRGrad_wc(wc)
    !
    ! 入力スペクトルデータに経度微分 k×r・▽ = ∂/∂λを作用する.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(nc,0:lm)             :: wc_KxRGrad_wc
    !(out) 経度微分を作用された 2 次元スペクトルデータ

    wc_KxRGrad_wc =  wa_Dlon_wa(wc)

  end function wc_KxRGrad_wc

  function ws_KxRGrad_ws(ws)
    !
    ! 入力スペクトルデータに経度微分 k×r・▽ = ∂/∂λを作用する.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ

    real(8), dimension(nc,lm)             :: ws_KxRGrad_ws
    !(out) 経度微分を作用された 2 次元スペクトルデータ

    ws_KxRGrad_ws =  wa_Dlon_wa(ws)

  end function ws_KxRGrad_ws

  function wc_L2_wc(wc)
    !
    ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する.
    !
    ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる.
    !  入力スペクトルデ ータに対応する格子点データに演算子
    !
    !     L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ)
    !
    ! を作用させたデータのスペクトル変換が返される.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(nc,0:lm)             :: wc_L2_wc
    !(out) L^2 演算子を作用された 2 次元スペクトルデータ

    wc_L2_wc = -wa_Lapla_wa(wc)

  end function wc_L2_wc

  function ws_L2_ws(ws)
    !
    ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する.
    !
    ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる.
    !  入力スペクトルデ ータに対応する格子点データに演算子
    !
    !     L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ)
    !
    ! を作用させたデータのスペクトル変換が返される.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ

    real(8), dimension(nc,lm)             :: ws_L2_ws
    !(out) L^2 演算子を作用された 2 次元スペクトルデータ

    ws_L2_ws = -wa_Lapla_wa(ws)

  end function ws_L2_ws

  function wc_L2Inv_wc(wc)
    !
    ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を
    ! 作用する.
    !
    ! スペクトルデータに L^2 演算子を作用させる関数 wc_L2_wc の逆計算を
    ! 行う関数である.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc
    !(in) 2 次元球面調和函数 cos スペクトルデータ

    real(8), dimension(nc,0:lm)             :: wc_L2Inv_wc
    !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ

    wc_L2Inv_wc = -wa_LaplaInv_wa(wc)

  end function wc_L2Inv_wc

  function ws_L2Inv_ws(ws)
    !
    ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を
    ! 作用する.
    !
    ! スペクトルデータに L^2 演算子を作用させる関数 ws_L2_ws の逆計算を
    ! 行う関数である.
    !
    real(8), dimension(nc,lm), intent(in) :: ws
    !(in) 2 次元球面調和函数 sin スペクトルデータ

    real(8), dimension(nc,lm)             :: ws_L2Inv_ws
    !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ

    ws_L2Inv_ws = -wa_LaplaInv_wa(ws)

  end function ws_L2Inv_ws

  function wc_RadRot_xvz_xvz(xvz_VLON,xvz_VLAT)  ! r・(▽×v)
    !
    ! ベクトルの渦度と動径ベクトルの内積 r・(▽×v) を計算する.
    !
    ! 第 1, 2 引数(v[λ], v[φ])がそれぞれベクトルの経度成分, 緯度成分を表す.
    !
    !    r・(▽×v) = 1/cosφ・∂v[φ]/∂λ - 1/cosφ・∂(v[λ] cosφ)/∂φ
    !
    ! のスペクトル データが返される.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_VLON
    !(in) ベクトルの経度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_VLAT
    !(in) ベクトルの緯度成分

    real(8), dimension(nc,0:lm)     :: wc_RadRot_xvz_xvz
    !(out) ベクトルの渦度と動径ベクトルの内積

    wc_RadRot_xvz_xvz = wc_wz(wa_DivLon_xva(xvz_VLAT) &
         - wa_DivLat_xva(xvz_VLON))

  end function wc_RadRot_xvz_xvz

  function ws_RadRotRot_xvz_xvz_xvz(xvz_VLON,xvz_VLAT,xvz_VRAD) 
    ! 
    ! ベクトル v に対して近似的な r・(▽×▽×v) を計算する.
    !
    ! 第 1, 2, 3 引数(v[λ], v[φ], v[r])がそれぞれベクトルの経度成分, 
    ! 緯度成分, 動径成分を表す. 
    !
    !    r・(▽×▽×v)  = ∂/∂r ((  1/cosφ・∂v[λ]/∂λ 
    !                               + 1/cosφ・∂(v[φ] cosφ)/∂φ ) ) 
    !                     + L^2 v[r]/a 
    !
    ! のスペクトルデータが返される.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_VLON
    !(in) ベクトルの経度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_VLAT
    !(in) ベクトルの緯度成分

    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz_VRAD
    !(in) ベクトルの動径成分

    real(8), dimension(nc,lm)     :: ws_RadRotRot_xvz_xvz_xvz
    !(out) ベクトル v の r・(▽×▽×v) 

    ws_RadRotRot_xvz_xvz_xvz = &
         ws_RotRad_wc(wc_wz( &
         (wa_DivLon_xva(xvz_VLON)+ wa_DivLat_xva(xvz_VLAT)))) &
         + ws_L2_ws(ws_xvz(xvz_VRAD/Ro))

  end function ws_RadRotRot_xvz_xvz_xvz

  subroutine wsc_Potential2VectorMPI(&
       xvz_VLON,xvz_VLAT,xvz_VRAD,wc_TORPOT,ws_POLPOT)
    !
    ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
    !
    !     v = ▽x(Ψr) + ▽x▽x(Φr) 
    !
    ! の各成分を計算する
    !
    real(8), dimension(0:im-1,jc,0:km)     :: xvz_VLON
    !(out) ベクトル場の経度成分

    real(8), dimension(0:im-1,jc,0:km)     :: xvz_VLAT
    !(out) ベクトル場の緯度成分

    real(8), dimension(0:im-1,jc,0:km)     :: xvz_VRAD
    !(out) ベクトル場の動径成分

    real(8), dimension(nc,0:lm), intent(in) :: wc_TORPOT
    !(in) トロイダルポテンシャル

    real(8), dimension(nc,lm), intent(in)   :: ws_POLPOT
    !(in) ポロイダルポテンシャル

    xvz_VLON =   Ro * xvz_GradLat_wc(wc_TORPOT) &
         + xva_GradLon_wa(wz_wc(wc_RotRad_ws(ws_POLPOT)))
    xvz_VLAT = - Ro * xvz_GradLon_wc(wc_TORPOT) &
         + xva_GradLat_wa(wz_wc(wc_RotRad_ws(ws_POLPOT)))
    xvz_VRAD = xvz_ws(ws_L2_ws(ws_POLPOT))/Ro

  end subroutine wsc_Potential2VectorMPI

  subroutine wsc_Potential2RotationMPI(&
       xvz_RotVLON,xvz_RotVLAT,xvz_RotVRAD,wc_TORPOT,ws_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,jc,0:km), intent(OUT) :: xvz_RotVLON
    !(out) 回転の経度成分

    real(8), dimension(0:im-1,jc,0:km), intent(OUT) :: xvz_RotVLAT
    !(out) 回転の緯度成分

    real(8), dimension(0:im-1,jc,0:km), intent(OUT) :: xvz_RotVRAD
    !(out) 回転の動径成分

    ! 入力ベクトル場を表すポテンシャル
    real(8), dimension(nc,0:lm), intent(in) :: wc_TORPOT
    !(in) トロイダルポテンシャル

    real(8), dimension(nc,lm), intent(in)   :: ws_POLPOT
    !(in) ポロイダルポテンシャル

    xvz_RotVLON =   Ro * xvz_GradLat_ws(-ws_Lapla_ws(ws_POLPOT)) &
         + xva_GradLon_wa(wz_ws(ws_RotRad_wc(wc_TORPOT)))
    xvz_RotVLAT = - Ro * xvz_GradLon_ws(-ws_Lapla_ws(ws_POLPOT)) &
         + xva_GradLat_wa(wz_ws(ws_RotRad_wc(wc_TORPOT)))
    xvz_RotVRAD = xvz_wc(wc_L2_wc(wc_TORPOT))/Ro

  end subroutine wsc_Potential2RotationMPI

  !--------------- 積分計算 -----------------
  !----(入力データ xvz)---
  function vz_IntLon_xvz(xvz)  ! 経度(帯状)積分
    !
    ! 3 次元格子点データの経度方向(帯状)積分.
    !
    ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(jc,0:km)  :: vz_IntLon_xvz
    !(out) 経度方向(帯状)積分された 2 次元子午面格子点データ

    integer :: i

    vz_IntLon_xvz = 0.0d0
    do i=0,im-1
       vz_IntLon_xvz(:,:) = vz_IntLon_xvz(:,:) &
            + xvz(i,:,:) * x_Lon_Weight(i)
    enddo
  end function vz_IntLon_xvz

  function xz_IntLat_xvz(xvz)
    !
    ! 3 次元格子点データの緯度方向域積分.
    !
    ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) cosφ dφ を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(0:im-1,0:km)  :: xz_IntLat_xvz        ! 緯度円格子点データ
    !(out) 緯度積分された 2 次元緯度動径格子点データ

    real(8), dimension(0:im-1,0:km)  :: xz_IntLatTMP
    integer :: j

    xz_IntLat_xvz = 0.0d0
    do j=1,jc
       xz_IntLat_xvz(:,:) = xz_IntLat_xvz(:,:) &
            + xvz(:,j,:) * v_Lat_Weight(j)
    enddo
    xz_IntLatTmp=xz_IntLat_xvz
    CALL MPI_ALLREDUCE(xz_IntLatTMP,xz_IntLat_xvz,im*(km+1),MPI_REAL8, &
         MPI_SUM,MPI_COMM_WORLD,IERR)

  end function xz_IntLat_xvz

  function xv_IntRad_xvz(xvz)  ! 動径積分
    !
    ! 3 次元格子点データの動径方向域積分.
    !
    ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dr を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(0:im-1,jc)  :: xv_IntRad_xvz
    !(out) 動径積分された 2 次元経度緯度(水平, 球面)格子点データ

    integer :: k

    xv_IntRad_xvz = 0.0d0
    do k=0,km
       xv_IntRad_xvz(:,:) = xv_IntRad_xvz(:,:) &
            + xvz(:,:,k) * z_Rad_Weight(k) 
    enddo
  end function xv_IntRad_xvz

  function x_IntLatRad_xvz(xvz)
    !
    ! 3 次元格子点データの緯度動径(子午面)積分
    !
    ! 3 次元データ f(λ,φ,r) に対して
    !
    !    ∫f(λ,φ,r) r^2cosφ dφdr 
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(0:im-1)     :: x_IntLatRad_xvz
    !(out) 緯度動径(子午面)積分された 1 次元経度格子点データ

    real(8), dimension(0:im-1)     :: x_IntLatRadTMP
    integer :: j, k

    x_IntLatRad_xvz = 0
    do k=0,km
       do j=1,jc
          x_IntLatRad_xvz = x_IntLatRad_xvz &
               + xvz(:,j,k) * v_Lat_Weight(j) * z_Rad_Weight(k)
       enddo
    enddo

    x_IntLatRadTmp=x_IntLatRad_xvz
    CALL MPI_ALLREDUCE(x_IntLatRadTMP,x_IntLatRad_xvz,im,MPI_REAL8, &
         MPI_SUM,MPI_COMM_WORLD,IERR)

  end function x_IntLatRad_xvz

  function v_IntLonRad_xvz(xvz)
    !
    ! 3 次元格子点データの経度動径(緯度円)積分.
    !
    ! 3 次元データ f(λ,φ,r) に対して∫f(λ,φ,r) r^2dλdr を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(jc)       :: v_IntLonRad_xvz
    !(out) 経度動径(緯度円)積分された 1 次元緯度格子点データ

    integer :: i, k

    v_IntLonRad_xvz = 0
    do k=0,km
       do i=0,im-1
          v_IntLonRad_xvz = v_IntLonRad_xvz &
               + xvz(i,:,k) * x_Lon_Weight(i) * z_Rad_Weight(k)
       enddo
    enddo
  end function v_IntLonRad_xvz

  function z_IntLonLat_xvz(xvz)  ! 緯度経度(水平)積分
    !
    ! 3 次元格子点データの緯度経度(水平, 球面)積分
    ! 
    ! 3 次元データ f(λ,φ,r) に対して
    !
    !    ∫f(λ,φ,r) cosφ dλdφ 
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(0:km)     :: z_IntLonLat_xvz
    !(out) 緯度経度(水平, 球面)積分された 1 次元動径格子点データ

    real(8), dimension(0:km)     :: z_IntLonLatTMP
    integer :: i, j

    z_IntLonLat_xvz = 0
    do j=1,jc
       do i=0,im-1
          z_IntLonLat_xvz = z_IntLonLat_xvz &
               + xvz(i,j,:) * x_Lon_Weight(i) * v_Lat_Weight(j)
       enddo
    enddo

    z_IntLonLatTmp=z_IntLonLat_xvz
    CALL MPI_ALLREDUCE(z_IntLonLatTMP,z_IntLonLat_xvz,km+1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_WORLD,IERR)

  end function z_IntLonLat_xvz

  function IntLonLatRad_xvz(xvz) ! 緯度経度動径(全球)積分
    !
    ! 3 次元格子点データの緯度経度動径(全球)積分
    !
    ! 3 次元データ f(λ,φ,r) に対して
    !
    !     ∫f(λ,φ,r) r^2cosφ dλdφdr 
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz 
    !(in) 3 次元経度緯度動径格子点データ

    real(8)                     :: IntLonLatRad_xvz 
    !(out) 全球積分値

    real(8)                     :: IntLonLatRadTMP
    integer :: i, j, k

    IntLonLatRad_xvz = 0
    do k=0,km
       do j=1,jc
          do i=0,im-1
             IntLonLatRad_xvz = IntLonLatRad_xvz &
                  + xvz(i,j,k) * x_Lon_Weight(i) &
                  * v_Lat_Weight(j) * z_Rad_Weight(k)
          enddo
       enddo
    enddo

    IntLonLatRadTmp=IntLonLatRad_xvz
    CALL MPI_ALLREDUCE(IntLonLatRadTMP,IntLonLatRad_xvz,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_WORLD,IERR)

  end function IntLonLatRad_xvz

  !----(入力データ vz)---
  function z_IntLat_vz(vz)  ! 緯度積分
    !
    ! 2 次元(VZ)格子点データの緯度方向域積分.
    !
    ! 2 次元データ f(φ,r) に対して∫f(φ,r) cosφ dφ を計算する.
    !
    real(8), dimension(jc,0:km), intent(in) :: vz
    !(in) 2 次元緯度動径(子午面)格子点データ

    real(8), dimension(0:km)  :: z_IntLat_vz
    !(out) 緯度積分された 1 次元動径格子点データ

    real(8), dimension(0:km)  :: z_IntLatTMP
    integer :: j

    z_IntLat_vz = 0.0d0
    do j=1,jc
       z_IntLat_vz(:) = z_IntLat_vz(:) + vz(j,:) * v_Lat_Weight(j)
    enddo
    z_IntLatTmp=z_IntLat_vz
    CALL MPI_ALLREDUCE(z_IntLatTMP,z_IntLat_vz,km+1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_WORLD,IERR)

  end function z_IntLat_vz

  function v_IntRad_vz(vz)  ! 動径積分
    !
    ! 2 次元(VZ)格子点データの動径方向域積分.
    !
    ! 2 次元データ f(φ,r) に対して∫f(φ,r) r^2dr を計算する.
    !
    real(8), dimension(jc,0:km), intent(in) :: vz
    !(in) 2 次元緯度動径(子午面)格子点データ

    real(8), dimension(jc)  :: v_IntRad_vz
    !(out) 動径積分された 1 次元緯度格子点データ

    integer :: k

    v_IntRad_vz = 0.0d0
    do k=0,km
       v_IntRad_vz(:) = v_IntRad_vz(:) &
            + vz(:,k) * z_Rad_Weight(k) 
    enddo
  end function v_IntRad_vz

  function IntLatRad_vz(vz)
    !
    ! 2 次元(VZ)格子点データの緯度動径積分(子午面)および平均
    !
    ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2cosφ dφdr を計算する.
    !
    real(8), dimension(jc,0:km), intent(in) :: vz
    !(in) 2 次元緯度動径(子午面)格子点データ

    real(8)                   :: IntLatRad_vz
    !(out) 積分値

    real(8)                   :: IntLatRadTMP
    integer :: j, k

    IntLatRad_vz = 0
    do k=0,km
       do j=1,jc
          IntLatRad_vz = IntLatRad_vz &
               + vz(j,k) * v_Lat_Weight(j) * z_Rad_Weight(k)
       enddo
    enddo

    IntLatRadTmp=IntLatRad_vz
    CALL MPI_ALLREDUCE(IntLatRadTMP,IntLatRad_vz,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_WORLD,IERR)

  end function IntLatRad_vz

  !----(入力データ xz)---
  function z_IntLon_xz(xz)
    !
    ! 2 次元(XZ)格子点データの経度方向積分.
    !
    ! 2 次元データ f(λ,r) に対して ∫f(λ,r)dλ を計算する.
    !
    real(8), dimension(0:im-1,0:km), intent(in) :: xz
    !(in) 2 次元緯度動径格子点データ

    real(8), dimension(0:km)  :: z_IntLon_xz
    !(out) 経度積分された 1 次元動径格子点データ

    integer :: i

    z_IntLon_xz = 0.0d0
    do i=0,im-1
       z_IntLon_xz(:) = z_IntLon_xz(:) + xz(i,:) * x_Lon_Weight(i)
    enddo

  end function z_IntLon_xz

  function x_IntRad_xz(xz)
    !
    ! 2 次元(XZ)格子点データの動径方向域積分.
    !
    ! 2 次元データ f(λ,r) に対して ∫f(λ,r) r^2dr を計算する.
    !
    real(8), dimension(0:im-1,0:km), intent(in) :: xz
    !(in) 2 次元緯度動径格子点データ

    real(8), dimension(0:im-1)  :: x_IntRad_xz
    !(out) 動径積分された 1 次元経度格子点データ

    integer :: k

    x_IntRad_xz = 0.0d0
    do k=0,km
       x_IntRad_xz(:) = x_IntRad_xz(:) &
            + xz(:,k) * z_Rad_Weight(k)
    enddo

  end function x_IntRad_xz

  function IntLonRad_xz(xz)  ! 経度動径(緯度円)積分
    !
    ! 2 次元(XZ)格子点データの経度動径積分
    !
    ! 2 次元データ f(λ,r) に対して∫f(λ,r) r^2dλdr を計算する.
    !
    real(8), dimension(0:im-1,0:km), intent(in) :: xz
    !(in) 2 次元緯度動径格子点データ

    real(8)                                 :: IntLonRad_xz
    !(out) 積分値

    integer :: i, k

    IntLonRad_xz = 0
    do k=0,km
       do i=0,im-1
          IntLonRad_xz = IntLonRad_xz &
               + xz(i,k) * x_Lon_Weight(i) * z_Rad_Weight(k)
       enddo
    enddo
  end function IntLonRad_xz

  !----(入力データ z)---
  function IntRad_z(z)  ! 動径積分
    !
    ! 1 次元(Z)格子点データの動径方向域積分.
    !
    ! 1 次元データ f(r) に対して ∫f(r) r^2dr を計算する.
    !
    real(8), dimension(0:km), intent(in) :: z
    !(in) 1 次元動径格子点データ

    real(8)                              :: IntRad_z
    !(out) 積分値

    integer :: k

    IntRad_z = 0.0d0
    do k=0,km
       IntRad_z = IntRad_z + z(k) * z_Rad_Weight(k)
    enddo
  end function IntRad_z

  !--------------- 平均計算 -----------------
  !----(入力データ xvz)---
  function vz_AvrLon_xvz(xvz)  ! 経度(帯状)積分
    !
    ! 3 次元格子点データの経度方向(帯状)平均.
    !
    ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)dλ/2π を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(jc,0:km)  :: vz_AvrLon_xvz
    !(out) 経度方向(帯状)平均された 2 次元子午面格子点データ

    vz_AvrLon_xvz = vz_IntLon_xvz(xvz)/sum(x_Lon_Weight)

  end function vz_AvrLon_xvz

  function xz_AvrLat_xvz(xvz)  ! 緯度積分
    !
    ! 3 次元格子点データの緯度方向域平均.
    !
    ! 3 次元データ f(λ,φ,r) に対して ∫f(λ,φ,r)cosφ dφ/2 を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(im,0:km)  :: xz_AvrLat_xvz
    !(out) 緯度平均された 2 次元緯度動径格子点データ

    xz_AvrLat_xvz = xz_IntLat_xvz(xvz)/Lat_Weight_Sum

  end function xz_AvrLat_xvz

  function xv_AvrRad_xvz(xvz)
    !
    ! 3 次元格子点データの動径方向域平均.
    !
    ! 3 次元データ f(λ,φ,r) に対して 
    !
    !    ∫f(λ,φ,r) r^2dr/((r[o]^3-r[i]^3)/3) 
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(0:im-1,jc)  :: xv_AvrRad_xvz          ! 水平格子点データ
    !(out) 動径平均された 2 次元経度緯度(水平, 球面)格子点データ

    xv_AvrRad_xvz = xv_IntRad_xvz(xvz)/sum(z_Rad_Weight)

  end function xv_AvrRad_xvz

  function x_AvrLatRad_xvz(xvz)  ! 緯度動径(子午面)積分
    !
    ! 3 次元格子点データの緯度動径(子午面)平均
    !
    ! 3 次元データ f(λ,φ,r) に対して
    !
    !    ∫f(λ,,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) 
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(0:im-1)     :: x_AvrLatRad_xvz
    !(out) 緯度動径(子午面)平均された 1 次元経度格子点データ

    x_AvrLatRad_xvz = x_IntLatRad_xvz(xvz) &
         /( Lat_Weight_Sum*sum(z_Rad_Weight) )

  end function x_AvrLatRad_xvz

  function v_AvrLonRad_xvz(xvz)  ! 経度動径(緯度円)積分
    !
    ! 3 次元格子点データの経度動径(緯度円)平均.
    !
    ! 3 次元データ f(λ,φ,r) に対して
    !
    !     ∫f(λ,φ,r) r^2dλdr /(2π(r[o]^3-r[i]^3)/3) 
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(jc)       :: v_AvrLonRad_xvz
    !(out) 経度動径(緯度円)平均された 1 次元緯度格子点データ

    v_AvrLonRad_xvz = v_IntLonRad_xvz(xvz) &
         /(sum(x_Lon_Weight)*sum(z_Rad_Weight))

  end function v_AvrLonRad_xvz

  function z_AvrLonLat_xvz(xvz)  ! 緯度経度(水平)積分
    !
    ! 3 次元格子点データの緯度経度(水平, 球面)積分
    ! 
    ! 3 次元データ f(λ,φ,r) に対して
    !
    !    ∫f(λ,φ,r) cosφ dλdφ /4π 
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8), dimension(0:km)     :: z_AvrLonLat_xvz
    !(out) 緯度経度(水平, 球面)平均された 1 次元動径格子点データ

    z_AvrLonLat_xvz = z_IntLonLat_xvz(xvz) &
         /(sum(x_Lon_Weight)*Lat_Weight_Sum)

  end function z_AvrLonLat_xvz

  function AvrLonLatRad_xvz(xvz) ! 緯度経度動径(全球)積分
    !
    ! 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,jc,0:km), intent(in) :: xvz
    !(in) 3 次元経度緯度動径格子点データ

    real(8)                     :: AvrLonLatRad_xvz
    !(out) 全球平均値

    AvrLonLatRad_xvz = IntLonLatRad_xvz(xvz) &
         /(sum(x_Lon_Weight)*Lat_Weight_Sum * sum(z_Rad_Weight))

  end function AvrLonLatRad_xvz

  !----(入力データ vz)---
  function z_AvrLat_vz(vz)
    !
    ! 2 次元(VZ)格子点データの緯度方向域平均.
    !
    ! 2 次元データ f(φ,r) に対して ∫f(φ,r) cosφ dφ/2 を計算する.
    !
    real(8), dimension(jc,0:km), intent(in) :: vz
    !(in) 2 次元緯度動径(子午面)格子点データ

    real(8), dimension(0:km)  :: z_AvrLat_vz
    !(out) 緯度平均された 1 次元動径格子点データ

    z_AvrLat_vz = z_IntLat_vz(vz)/Lat_Weight_Sum
  end function z_AvrLat_vz

  function v_AvrRad_vz(vz)
    !
    ! 2 次元(VZ)格子点データの動径方向域平均.
    !
    ! 2 次元データ f(φ,r) に対して ∫f(φ,r) r^2dr /((r[o]^3-r[i]^3)/3) 
    ! を計算する.
    !
    real(8), dimension(jc,0:km), intent(in) :: vz
    !(in) 2 次元緯度動径(子午面)格子点データ

    real(8), dimension(jc)  :: v_AvrRad_vz
    !(out) 動径平均された 1 次元緯度格子点データ

    v_AvrRad_vz = v_IntRad_vz(vz)/sum(z_Rad_Weight)

  end function v_AvrRad_vz

  function AvrLatRad_vz(vz)  ! 緯度動径(子午面)積分
    !
    ! 2 次元(VZ)格子点データの緯度動径(子午面)平均
    !
    ! 2 次元データ f(φ,r) に対して
    !
    !    ∫f(φ,r) r^2cosφ dφdr /(2(r[o]^3-r[i]^3)/3) 
    !
    ! を計算する.
    !
    real(8), dimension(jc,0:km), intent(in) :: vz
    !(in) 2 次元緯度動径(子午面)格子点データ

    real(8)                   :: AvrLatRad_vz
    !(out) 平均値

    AvrLatRad_vz = IntLatRad_vz(vz)/(Lat_Weight_Sum*sum(z_Rad_Weight))

  end function AvrLatRad_vz

  !----(入力データ xz)---
  function z_AvrLon_xz(xz)  ! 経度(帯状)積分
    !
    ! 2 次元(XZ)格子点データの経度方向平均.
    !
    ! 2 次元データ f(λ,r) に対して ∫f(λ,r)dλ/2π を計算する.
    !
    real(8), dimension(0:im-1,0:km), intent(in) :: xz
    !(in) 2 次元緯度動径格子点データ

    real(8), dimension(0:km)  :: z_AvrLon_xz
    !(out) 経度平均された 1 次元動径格子点データ

    z_AvrLon_xz = z_IntLon_xz(xz)/sum(x_Lon_Weight)

  end function z_AvrLon_xz

  function x_AvrRad_xz(xz)  ! 動径積分
    !
    ! 2 次元(XZ)格子点データの動径方向域平均.
    !
    ! 2 次元データ f(λ,r) に対して
    !
    !   ∫f(λ,r) r^2dr /((r[o]^3-r[i]^3)/3)
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,0:km), intent(in) :: xz
    !(in) 2 次元緯度動径格子点データ

    real(8), dimension(0:im-1)  :: x_AvrRad_xz
    !(out) 動径平均された 1 次元経度格子点データ

    x_AvrRad_xz = x_IntRad_xz(xz)/sum(z_Rad_Weight)

  end function x_AvrRad_xz

  function AvrLonRad_xz(xz)  ! 経度動径(緯度円)積分
    !
    ! 2 次元(XZ)格子点データの経度動径平均
    !
    ! 2 次元データ f(λ,r) に対して
    !
    !    ∫f(λ,r) r^2dλdr /(2π(r[o]^3-r[i]^3)/3)
    !
    ! を計算する.
    !
    real(8), dimension(0:im-1,0:km), intent(in) :: xz
    ! (in) 2 次元格子点データ
    real(8)                                 :: AvrLonRad_xz
    ! 積分値

    AvrLonRad_xz = IntLonRad_xz(xz)/(sum(x_Lon_Weight)*sum(z_Rad_Weight))

  end function AvrLonRad_xz

  !----(入力データ z)---
  function AvrRad_z(z)
    !
    ! 1 次元(Z)格子点データの動径方向域平均.
    !
    ! 1 次元データ f(r) に対して ∫f(r) r^2dr /((r[o]^3-r[i]^3)/3) を
    ! 計算する.
    !
    real(8), dimension(0:km), intent(in) :: z
    !(in) 1 次元動径格子点データ
    real(8)                              :: AvrRad_z
    !(out) 平均値

    AvrRad_z = IntRad_z(z)/sum(z_Rad_Weight)

  end function AvrRad_z

  !--------------- 補間計算 -----------------
  function Interpolate_wc(wc_data,alon,alat,arad)
    !
    ! 緯度 alon, 経度 alat 動径 arad における関数値を
    ! その球面調和変換係数 wa_data から補間計算する
    !
    real(8), intent(IN) :: wc_data(nc,0:lm)         ! スペクトルデータ
    real(8), intent(IN) :: alon                     ! 補間する位置(経度)
    real(8), intent(IN) :: alat                     ! 補間する位置(緯度)
    real(8), intent(IN) :: arad                     ! 補間する位置(動径)
    real(8) :: Interpolate_wc                       ! 補間した値

    Interpolate_wc = &
         Interpolate_w(a_Interpolate_ac(wc_data,arad),alon,alat)

  end function Interpolate_wc

  function Interpolate_ws(ws_data,alon,alat,arad)
    !
    ! 緯度 alon, 経度 alat 動径 arad における関数値を
    ! その球面調和変換係数 wa_data から補間計算する
    !
    real(8), intent(IN) :: ws_data(nc,lm)           ! スペクトルデータ
    real(8), intent(IN) :: alon                     ! 補間する位置(経度)
    real(8), intent(IN) :: alat                     ! 補間する位置(緯度)
    real(8), intent(IN) :: arad                     ! 補間する位置(動径)
    real(8) :: Interpolate_ws                       ! 補間した値

    Interpolate_ws = &
         Interpolate_w(a_Interpolate_as(ws_data,arad),alon,alat)

  end function Interpolate_ws

  !--------------- ポロイダル/トロイダルモデル用スペクトル解析 -------------

  function nmz_ToroidalEnergySpectrum_wc(wc_TORPOT)
    !
    ! トロイダルポテンシャルから, トロイダルエネルギーの
    ! 球面調和函数全波数 n, 帯状波数 m の各成分を計算する
    !
    !  * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分
    !    ψ(n,m,r)から全波数 n, 帯状波数 m 成分のトロイダルエネルギー
    !    スペクトルは  (1/2)n(n+1)4πa^2|ψ(n,m,r)|^2  と計算される.
    !
    !  * 全てのエネルギースペクトル成分の和を動径積分したもの(a^2の重み無し)
    !    が球殻内での全エネルギーに等しい.
    !
    !  * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
    !    wsc_VMiss によって設定できる (初期値は -999.0)
    !
    real(8), intent(in)   :: wc_TORPOT(nc,0:lm)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_ToroidalEnergySpectrum_wc
    !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間)

    real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_Work ! 作業領域

    real(8), dimension(nc,0:km) ::wz_DATA      ! 作業領域

    integer :: l, n, m, nm_ary(2)

    integer :: ierr

    wz_DATA = wz_wc(wc_TORPOT)

    nmz_work = 0.0D0
    do l=1,nc
       nm_ary = nm_l(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             nmz_Work(n,0,:) = nmz_Work(n,0,:) &
                  +  0.5 * n*(n+1) * (4*pi) * Ro**2 * wz_Data(l,:)**2
          else
             nmz_Work(n,m,:) = nmz_Work(n,m,:) &
                  +  0.5 * n*(n+1) * (4*pi) * Ro**2 * wz_Data(l,:)**2
             nmz_Work(n,-m,:) = nmz_Work(n,-m,:) &
                  +  0.5 * n*(n+1) * (4*pi) * Ro**2 * wz_Data(l,:)**2
          endif
       endif
    enddo

    CALL MPI_ALLREDUCE(nmz_Work,nmz_ToroidalEnergySpectrum_wc, &
         (nm+1)*(2*nm+1)*(km+1), &
         MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

    do m=1,nm
       do n=0,m-1
          nmz_ToroidalEnergySpectrum_wc(n,m,:)  = wsc_Vmiss
          nmz_ToroidalEnergySpectrum_wc(n,-m,:) = wsc_Vmiss
       enddo
    enddo

  end function nmz_ToroidalEnergySpectrum_wc

  function nz_ToroidalEnergySpectrum_wc(wc_TORPOT)
    !
    ! トロイダルポテンシャルから, トロイダルエネルギーの
    ! 球面調和函数全波数の各成分を計算する.
    !
    !  * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分
    !    ψ(n,m,r)から全波数 n 成分のトロイダルエネルギースペクトルは
    !    Σ[m=-n]^n(1/2)n(n+1)4πa^2|ψ(n,m,r)|^2 と計算される.
    !
    ! * 全てのエネルギースペクトル成分の和を動径積分したもの(a^2の重み無し)
    !    が球殻内での全エネルギーに等しい.
    !
    real(8), dimension(nc,0:lm), intent(in) :: wc_TORPOT
    !(in) トロイダルポテンシャル

    real(8), dimension(0:nm,0:km) :: nz_ToroidalEnergySpectrum_wc
    !(out) エネルギースペクトルトロイダル成分

    real(8), dimension(nc,0:km) ::wz_DATA   ! 作業領域
    real(8), dimension(0:nm,0:km)     ::nz_Work
    integer :: l, n, m, nm_ary(2)
    real(8) :: factor

    nz_Work = 0.0D0
    wz_DATA = wz_wc(wc_TORPOT)
    do l=1,nc
       nm_ary = nm_l(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             factor = 1.0D0
          else
             factor = 2.0D0
          endif
          nz_Work(n,:) = nz_Work(n,:) &
               + factor * 0.5 * n*(n+1)* (4*pi) * Ro**2 * wz_Data(l,:)**2
          !              + 0.5 * n*(n+1)* (4*pi) * z_Rad**2 * wz_Data(l,:)**2
       endif
    enddo

    CALL MPI_ALLREDUCE(nz_Work,nz_ToroidalEnergySpectrum_wc,(nm+1)*(km+1), &
         MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

  end function nz_ToroidalEnergySpectrum_wc

  function nmz_PoloidalEnergySpectrum_ws(ws_POLPOT)
    !
    ! ポロイダルポテンシャルから, ポロイダルエネルギーの
    ! 球面調和函数全波数 n, 帯状波数 m の各成分を計算する.
    !
    !  * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分
    !    φ(n,m,r)から全波数 n, 帯状波数 m 成分のポロイダルエネルギー
    !    スペクトルは
    !
    !      (1/2)n(n+1)4π{|d(aφ(n,m,r))/dr|^2 + n(n+1)|φ(n,m,r)|^2} 
    !
    !    と計算される.
    !
    !  * 全てのエネルギースペクトル成分の和を動径積分したもの(a^2の重み無し)
    !    が球殻内での全エネルギーに等しい.
    !
    !  * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
    !    欠損値の値はモジュール変数 wsc_VMiss によって設定できる
    !    (初期値は -999.0)
    !
    real(8), intent(in)   :: ws_POLPOT(nc,lm)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_PoloidalEnergySpectrum_ws
    !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間)

    real(8), dimension(0:nm,-nm:nm,0:km) :: nmz_Work ! 作業領域

    real(8), dimension(nc,0:km) ::wz_DATA1   ! 作業領域
    real(8), dimension(nc,0:km) ::wz_DATA2   ! 作業領域

    integer :: l, n, m, nm_ary(2)

    integer :: ierr

    wz_Data1 = wz_ws(ws_POLPOT)
    wz_Data2 = Ro*wz_wc(wc_DRad_ws(ws_POLPOT))  ! d(rφ)/dr ~ a dφ/dr

    nmz_work = 0.0D0
    do l=1,nc
       nm_ary = nm_l(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             nmz_Work(n,0,:) = nmz_Work(n,0,:) &
                  +  0.5 * n*(n+1) * (4*pi) &
                  * ( wz_Data2(l,:)**2 +  n*(n+1)*wz_Data1(l,:)**2 )
          else
             nmz_Work(n,m,:) = nmz_Work(n,m,:) &
                  +  0.5 * n*(n+1) * (4*pi) &
                  * ( wz_Data2(l,:)**2 +  n*(n+1)*wz_Data1(l,:)**2 )
             nmz_Work(n,-m,:) = nmz_Work(n,-m,:) &
                  +  0.5 * n*(n+1) * (4*pi) &
                  * ( wz_Data2(l,:)**2 +  n*(n+1)*wz_Data1(l,:)**2 )
          endif
       endif
    enddo

    CALL MPI_ALLREDUCE(nmz_Work,nmz_PoloidalEnergySpectrum_ws, &
         (nm+1)*(2*nm+1)*(km+1), &
         MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

    do m=1,nm
       do n=0,m-1
          nmz_PoloidalEnergySpectrum_ws(n,m,:)  = wsc_Vmiss
          nmz_PoloidalEnergySpectrum_ws(n,-m,:) = wsc_Vmiss
       enddo
    enddo

  end function nmz_PoloidalEnergySpectrum_ws

  function nz_PoloidalEnergySpectrum_ws(ws_POLPOT)
    !
    ! ポロイダルポテンシャルから, ポロイダルエネルギーの
    ! 球面調和函数全波数の各成分を計算する
    !
    !  * 全波数 n, 帯状波数 m のポロイダルポテンシャルのスペクトル成分
    !    φ(n,m,r)から全波数 n 成分のポロイダルエネルギースペクトルは
    !
    !      Σ[m=-n]^n ((1/2)n(n+1)4π{|d(rφ(n,m,r))/dr|^2 
    !                 + n(n+1)|φ(n,m,r)|^2} 
    !
    !    と計算される.
    !
    !  * 全ての全波数に対してのエネルギースペクトル成分の和を動径積分したもの
    !    (a^2の重み無し)が球殻内での全エネルギーに等しい.
    !
    real(8), dimension(nc,lm), intent(in) :: ws_POLPOT
    !(in) ポロイダルポテンシャル

    real(8), dimension(0:nm,0:km) :: nz_PoloidalEnergySpectrum_ws
    !(out) エネルギースペクトルポロイダル成分

    real(8), dimension(nc,0:km) ::wz_DATA1   ! 作業領域
    real(8), dimension(nc,0:km) ::wz_DATA2   ! 作業領域
    real(8), dimension(0:nm,0:km)     ::nz_Work
    integer :: l, n, m, nm_ary(2)
    real(8) :: factor

    nz_Work = 0.0D0
    wz_Data1 = wz_ws(ws_POLPOT)
    wz_Data2 = Ro*wz_wc(wc_DRad_ws(ws_POLPOT))  ! d(rφ)/dr ~ a dφ/dr

    do l=1,nc
       nm_ary = nm_l(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             factor = 1.0D0
          else
             factor = 2.0D0
          endif
          nz_Work(n,:) = nz_Work(n,:) &
               + factor * 0.5* n*(n+1)* (4*pi) &
               *( wz_Data2(l,:)**2 + n*(n+1)*wz_Data1(l,:)**2 )
       endif
    enddo

    CALL MPI_ALLREDUCE(nz_Work,nz_PoloidalEnergySpectrum_ws,(nm+1)*(km+1), &
         MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

  end function nz_PoloidalEnergySpectrum_ws
    
end module wsc_mpi_module_supack
  
