!--
!----------------------------------------------------------------------
! Copyright (c) 2019 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  ut_mpi_module
!
!    spml/ut_mpi_module モジュールは球面上および球殻内での流体運動を
!    スペクトル法と MPI 並列化によって数値計算するための Fortran90
!    関数を提供するものである.
!
!    水平方向に球面調和函数変換および上下の境界壁を扱うための
!    チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな
!    関数を提供する.
!
!    内部で ua_mpi_module_syapck を用いている. 最下部では球面調和変換
!    として SYPACK/ISPACK3 を, チェビシェフ変換のエンジンとして ISPACK の
!    サブルーチンを用いている.
!
!
!履歴  2002/05/19  竹広真一  yt_module よりモジュール名変更.
!                            格子点次元を区別する必要があるため.
!      2002/06/10  竹広真一  ポロイダル磁場境界条件ルーチン追加
!      2002/11/10  竹広真一  実空間での境界条件ルーチン追加
!      2002/11/24  竹広真一  チェビシェフ係数入出力の
!                            ラプラシアン逆解きルーチンを追加
!      2002/11/28  竹広真一  VGradV の計算変更.
!      2005/01/09  竹広真一  msgdmp -> MessageNotify に変更
!      2005/01/09  竹広真一  ベクトル場の回転の各成分を計算する関数を追加.
!      2005/02/19  竹広真一  トロイダル粘着境界条件にて非零速度場を
!                             与えられるようにオプションを追加
!      2005/02/19  竹広真一  変数 xy_Lon, xy_Lat を追加
!      2005/03/18  竹広真一  積分・平均関数を追加
!                            z_IntLon_xz, x_IntRad_xz, IntLonRad_xz
!                            z_InRad_xz, y_IntRad_yz, IntLatRad_yz
!                            IntRad_z
!                            z_IntLon_xz, x_IntRad_xz, IntLonRad_xz
!                            z_InRad_xz, y_IntRad_yz, IntLatRad_yz
!                            IntRad_z
!      2005/04/24  竹広真一  スペクトル解析計算ルーチンを追加
!                            nmz_ToroidalEnergySpectrum_wt
!                            nz_ToroidalEnergySpectrum_wt
!                            nmz_PoloidalEnergySpectrum_wt
!                            nz_PoloidalEnergySpectrum_wt
!                            内部変数欠損値 wt_VMiss を追加
!      2005/07/09  竹広真一  OPENMP 版変換ルーチンに対応
!      2006/03/03  竹広真一  配列サイズミス修正
!                            wt_wz, nmz_ToroidalEnergySpectrum_wt,
!                            nmz_PoloidalEnergySpectrum_wt
!      2006/03/08  竹広真一  コメントを RDoc 用に修正
!      2006/03/19  竹広真一  変数・手続き群の要約をコメントに追加
!      2007/08/11  竹広真一  境界値ルーチンのコメントを増強
!      2007/09/15  竹広真一  微分ルーチンのコメントを修正
!      2007/11/02  竹広真一  補間計算ルーチン追加
!      2007/11/11  竹広真一  境界条件関数サブルーチンに行列設定スイッチ導入
!      2007/11/21  竹広真一  初期化サブルーチンメッセージ出力
!      2008/01/07  竹広真一  xyz_RotLon_wt_wt, xyz_RotLat_wt_wt バグ修正
!      2008/01/13  竹広真一  wa_initial スイッチ導入
!      2008/05/29  竹広真一  MPI 並列化
!      2008/06/11  竹広真一  OPENMP スイッチ導入
!      2008/06/11  竹広真一  配列添え字変更
!      2009/02/27  佐々木洋平 RDoc 用のコメントを修正
!      2010/01/07  佐々木洋平  RDoc 用のドキュメント修正,
!                              include 'mpif.h' -> use mpi
!      2011/09/08  竹広真一   コメント修正
!      2011/09/12  竹広真一   wt_mpi_Initial に wa_init スイッチ導入
!      2012/05/25  竹広真一   wt_QoperatorMPI_wt を追加
!      2012/08/24  竹広真一   wt_LaplaPol2Pol_wt を追加
!      2013/08/20  竹広真一  gnu fortran 対応
!      2014/03/10  竹広真一  wt_mpi_module_nm に改造
!      2014/03/17  竹広真一  mint パラメターに対応
!      2014/03/21  竹広真一  OpenMP 版
!      2015/02/15  竹広真一  Interpolate 関数追加
!      2016/03/02  竹広真一  nz_ToroidalEnergySpectrum_wt,
!                            nz_PoloidalEnergySpectrum_wt 追加
!      2016/11/04  竹広真一  nmz_ToroidalEnergySpectrum_wt,
!                            nmz_PoloidalEnergySpectrum_wt 追加
!      2017/02/07  竹広真一  wt_mpi_Finalize 追加
!      2017/10/26  竹広真一  ut_mpi_module_nmomp に改造
!      2017/11/07  竹広真一  ut_mpi_module_nmomp_lp に改造
!      2017/11/11  竹広真一  ut_wh, wh_ut を追加
!      2018/01/28  竹広真一  ut_mpi_module_nmomp に改名
!      2018/06/14  竹広真一  ut_mpi_module_extomp に改名, 変数定義変更
!      2019/05/28  竹広真一  ut_mpi_module に改造
!      2020/01/25  竹広真一  u_xv, xv_u を ua_mpi_module に移動
!      2021/11/15  竹広真一  Potentail2VectorMPI を ua_Potentail2VectorMPI で計算
!                             ut_Vector2VorDivMPI, ut_RadRot_RadRotRotRot を作成
!
!凡例
!      データ種類と index
!        x : 経度   y : 緯度    v : 緯度(分散格子) z : 動径  h : 動径(分散格子)
!        w : 球面調和関数スペクトル
!        n : 球面調和関数スペクトル(水平全波数)
!        m : 球面調和関数スペクトル(帯状波数)
!        t : チェビシェフ関数スペクトル
!        a : 任意の次元
!
!        xvh : 3 次元分散格子点データ
!        xv  : 水平 2 次元格子点データ
!        vh  : 子午面 2 次元格子点データ
!        xh  : 緯度面 2 次元格子点データ
!
!        wz  : 水平スペクトル動径格子点データ
!        wt  : スペクトルデータ
!
!++
module ut_mpi_module
  !
  != ut_mpi_module
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id: wt_mpi_module.f90 598 2013-08-20 03:23:44Z takepiro $
  ! Copyright&License:: See COPYRIGHT[link:../../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/wt_mpi_module モジュールは球面上および球殻内での流体運動を
  ! スペクトル法と MPI 並列化によって数値計算するための Fortran90
  ! 関数を提供するものである.
  !
  ! 水平方向に球面調和函数変換および上下の境界壁を扱うための
  ! チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな
  ! 関数を提供する.
  !
  ! 内部で ua_mpi_module_syapck を用いている. 最下部では球面調和変換
  ! として SYPACK/ISPACK3 を, チェビシェフ変換のエンジンとして ISPACK の
  ! サブルーチンを用いている.
  !
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭(ut_, nmz_, nz_, xyz_, xvh_ uz_, w_, xv_, x_, y_, z_, a_)は,
  !   返す値の形を示している.
  !      ut_ :: スペクトルデータ(球面調和函数・チェビシェフ変換)
  !     xvh_ :: 3 次元分散格子点データ(経度・緯度・動径)
  !      uz_ :: 水平スペクトル, 動径格子点データ
  !
  ! * 関数名の間の文字列(DLon, GradLat, GradLat, DivLon, DivLat, Lapla,..)
  !   は, その関数の作用を表している.
  !
  ! * 関数名の最後 (_ut, _xyz, _xvh, _uz_, _w, _xy, _xv, _x, _y, _z, _a) は,
  !   入力変数の形がスペクトルデータおよび格子点データであることを示している.
  !         _ut :: スペクトルデータ
  !        _xvh :: 3 次元分散格子点データ
  !    _xvh_xvh :: 2 つの3 次元分散格子点データ, ...
  !
  !=== 各データの種類の説明
  !
  ! * xvh : 3 次元分散格子点データ(経度・緯度・動径)
  !   * 変数の種類と次元は real(8), dimension(0:im-1,jc,kc).
  !   * im, は経度の格子点数である.
  !   * jc, kc はこのプロセスで保有する緯度・動径格子点数である.
  !     サブルーチン ut_mpi_Initial を呼ぶと jc が設定される.
  !
  ! * wh : スペクトルデータ
  !   * 変数の種類と次元は real(8), dimension(nc,kc).
  !   * nc はこのプロセスで保有する球面調和函数成分の数,
  !     lm はチェビシェフ多項式の最大次数であり,
  !     サブルーチン ut_mpi_Initial を呼ぶと nc が設定される.
  !   * 水平スペクトルデータの格納のされ方は関数 l_nm, nm_l によって調べる
  !     ことができる.
  !
  ! * ut : スペクトルデータ
  !   * 変数の種類と次元は real(8), dimension(nmc,0:lm).
  !   * nmc はこのプロセスで保有する球面調和函数成分の数,
  !     lm はチェビシェフ多項式の最大次数であり,
  !     サブルーチン ut_mpi_Initial を呼ぶと nc が設定される.
  !   * 水平スペクトルデータの格納のされ方は関数 lu_nm, nm_lu によって調べる
  !     ことができる.
  !
  ! * uz : 水平スペクトル, 動径格子点データ.
  !   * 変数の種類と次元は real(8), dimension(nmc,0:km).
  !
  ! * ut_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * xvh_ で始まる関数が返す値は 3 次元分散格子点データに同じ.
  !
  ! * uz_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! ut_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 次元配列
  ! xvh_Lon, xvh_Lat, xvh_Rad :: 格子点データの経度・緯度・動径座標(X,Y,Z) (分散格子点データ型 3 次元配列)
  !
  !==== 基本変換
  !
  ! xvh_ut, ut_xvh :: スペクトルデータと 3 次元分散格子データの間の変換 (球面調和函数, チェビシェフ変換)
  !
  ! xvh_uz, uz_xvh :: 3 次元分散格子データと水平スペクトル・動径格子データとの間の変換 (球面調和函数)
  !
  ! uz_ut, ut_uz   :: スペクトルデータと水平スペクトル・動径格子データとの間の変換 (チェビシェフ変換)
  !
  ! w_xv, xv_w     :: スペクトルデータと 2 次元水平分散格子データの間の変換(球面調和函数変換)
  !
  ! az_at, at_az   :: 同時に複数個行う (チェビシェフ変換)格子データとチェビシェフデータの間の変換を
  !
  ! l_nm, nm_l     :: スペクトルデータの格納位置と全波数・帯状波数の変換
  !
  !==== 微分
  !
  ! ut_DRad_ut     :: スペクトルデータに動径微分∂/∂r を作用させる
  ! ut_DivRad_ut   :: スペクトルデータに発散型動径微分 1 /r^2 ∂/∂r r^2 = ∂/∂r + 2/r を作用させる
  !
  ! ut_RotRad_ut   :: スペクトルデータに回転型動径微分 1/r ∂/∂rr = ∂/∂r + 1/r を作用させる
  !
  ! ut_Lapla_ut    :: スペクトルデータにラプラシアンを作用させる
  !
  ! xvh_GradLon_ut :: スペクトルデータに勾配型経度微分 1/rcosφ・∂/∂λを作用させる
  !
  ! xvh_GradLat_ut :: スペクトルデータに勾配型緯度微分 1/r・∂/∂φを作用させる(分散格子)
  ! ut_DivLon_xvh  :: 分散格子データに発散型経度微分 1/rcosφ・∂/∂λを作用させる
  ! ut_DivLat_xvh  :: 分散格子データに発散型緯度微分 1/rcosφ・∂(g cosφ)/∂φを作用させる(分散格子)
  !
  ! ut_Div_xvh_xvh_xvh  :: ベクトル成分である 3 つの分散格子データに発散を作用させる
  ! xvh_Div_xvh_xvh_xvh :: ベクトル成分である 3 つの分散格子データに発散を作用させる
  ! xvh_RotLon_ut_ut  :: ベクトル場の回転の経度成分を計算する(分散格子)
  ! xvh_RotLat_ut_ut  :: ベクトル場の回転の緯度成分を計算する(分散格子)
  ! ut_RotRad_xyz_xyz :: ベクトル場の回転の動径成分を計算する
  ! ut_RotRad_xvh_xvh :: ベクトル場の回転の動径成分を計算する(分散格子)
  !
  !==== トロイダルポロイダル計算用微分
  !
  ! ut_L2_ut        :: スペクトルデータに L2 演算子 = -水平ラプラシアンを作用させる
  !
  ! ut_L2Inv_ut     :: スペクトルデータに L2 演算子の逆 = -逆水平ラプラシアンを作用させる
  !
  ! ut_RadRot_xvh_xvh :: ベクトル v の渦度と動径ベクトル r の内積 r・(▽×v) を計算する(分散格子)
  ! ut_RadRotRot_xvh_xvh_xvh :: ベクトルの v の r・(▽×▽×v) を計算する(分散格子)
  ! ut_Potential2VectorMPI   :: トロイダルポロイダルポテンシャルからベクトル場を計算する(分散格子)
  !
  ! ut_Potential2RotationMPI :: トロイダルポロイダルポテンシャルで表される非発散ベクトル場の回転の各成分を計算する(分散格子)
  !
  !==== ポロイダル/トロイダルモデル用スペクトル解析
  !
  ! nmz_ToroidalEnergySpectrum_ut, :: トロイダルポテンシャルからエネルギーの
  ! nz_ToroidalEnergySpectrum_ut   :: 球面調和函数各成分を計算する
  !
  ! nmz_PoloidalEnergySpectrum_ut, :: ポロイダルポテンシャルからエネルギーの
  ! nz_PoloidalEnergySpectrum_ut   :: 球面調和函数各成分を計算する
  !
  !==== 境界値問題
  !
  ! ut_BoundariesTau,    :: ディリクレ, ノイマン境界条件を適用する
  ! ut_BoundariesGrid,   ::(タウ法, 選点法)
  ! ut_Boundaries        ::
  !
  ! ut_TorBoundariesTau,  :: 速度トロイダルポテンシャルの境界条件を
  ! ut_TorBoundariesGrid, :: 適用する(タウ法,選点法)
  ! ut_TorBoundaries      ::
  !
  ! uz_LaplaPol2Pol_uz,   :: 速度ポロイダルポテンシャルΦを▽^2Φから
  ! ut_LaplaPol2Pol_ut    :: 求める (入出力がそれぞれチェビシェフ格子点,
  !                       :: チェビシェフ係数)
  !
  ! ut_TorMagBoundariesTau,  :: 磁場トロイダルポテンシャルの境界条件を
  ! ut_TorMagBoundariesGrid, :: 適用する(タウ法, 選点法)
  ! ut_TorMagBoundaries      ::
  !
  ! ut_PolMagBoundariesTau,  :: 磁場トロイダルポテンシャル境界の境界条件を
  ! ut_PolMagBoundariesGrid, :: 適用する(タウ法, 選点法)
  ! ut_PolMagBoundaries      ::
  !
  !==== 積分・平均(3 次元データ)
  !
  ! IntLonLatRad_xvh, AvrLonLatRad_xvh :: 3 次元格子点データの全領域積分および平均
  !
  ! h_IntLonLat_xvh, h_AvrLonLat_xvh :: 3 次元格子点データの緯度経度(水平・球面)積分および平均
  !
  ! v_IntLonRad_xvh, v_AvrLonRad_xvh :: 3 次元格子点データの緯度動径積分および平均
  !
  ! h_IntLatRad_xvh, h_AvrLatRad_xvh :: 3 次元格子点データの経度動径(子午面)積分および平均
  !
  ! vh_IntLon_xvh, vz_AvrLon_xvh :: 3 次元格子点データの経度方向積分および平均
  ! xh_IntLat_xvh, xz_AvrLat_xvh :: 3 次元格子点データの緯度方向積分および平均
  !
  !==== 積分・平均(2 次元データ)
  !
  ! IntLonLat_xv, AvrLonLat_xv :: 2 次元格子点データの水平(球面)積分および平均
  ! IntLonRad_xh, AvrLonRad_xh :: 2 次元(XZ)格子点データの経度動径積分
  !                            :: および平均
  ! IntLatRad_vh, AvrLatRad_vh :: 2 次元(YZ)格子点データの緯度動径(子午面)
  !                            :: 積分および平均
  ! v_IntLon_xv, v_AvrLon_xv   :: 水平 2 次元(球面)格子点データの経度方向
  !                            :: 積分および平均
  ! v_IntLat_xv, x_AvrLat_xv   :: 水平2 次元(球面)格子点データの緯度方向積分
  !                            :: および平均
  ! h_IntLon_xh, h_AvrLon_xh   :: 2 次元(XZ)格子点データの経度方向積分および
  !                            :: 平均
  ! h_IntRad_xh, x_AvrRad_xh   :: 2 次元(XZ)格子点データの動径方向積分および
  !                            :: 平均
  ! h_IntLat_vh, z_AvrLat_vh   :: 2 次元(YZ)格子点データの緯度方向積分および
  !                            :: 平均
  ! v_IntRad_vh, v_AvrRad_vh   :: 2 次元(YZ)格子点データの動径方向積分および
  !                            :: 平均
  !
  !==== 積分・平均(1 次元データ)
  !
  ! IntLon_x, AvrLon_x  :: 1 次元(X)格子点データの経度方向積分および平均
  ! IntLat_v, AvrLat_v  :: 1 次元(Y)格子点データの緯度方向積分および平均
  ! IntRad_h, AvrRad_h  :: 1 次元(Z)格子点データの動径方向積分および平均
  !
  !==== 補間計算
  !
  ! Interpolate_ut :: スペクトルデータから任意の点の値を補間する.
  !
  use dc_types, only: DP
  use dc_message
  use lumatrix
  use mpi
  use ua_mpi_module, ncc=>nc
  use ua_mpi_module_base, only: &
    & MPI_COMM_H, MPI_COMM_V, iproc, iproc_h, iproc_v
  use w_mpi_module !only : w_xv =>u_xv, xv_w => xv_u,

  use at_module, &
    & z_RAD => g_X, z_RAD_WEIGHT => g_X_WEIGHT, &
    & at_az => at_ag, az_at => ag_at, &
    & t_z => t_g, z_t => g_t, &
    & t_Dr_t => t_Dx_t, at_Dr_at => at_Dx_at

  implicit none

  private

  public ut_mpi_Initial, ut_mpi_Finalize

  public nproc_h, nproc_v
  public iproc_h, iproc_v

  public jc                                   ! 緯度方向分散格子点数
  public ks,ke,kc
  public nc                                   ! MPI 分割水平波数(w_)
  public nms,nme,nmc                          ! MPI 分割水平波数(u_)

  public l_nm, nm_l
  public lu_nm, nm_lu
  public kp_kl

  public x_Lon, x_Lon_Weight
  public y_Lat, y_Lat_Weight
  public v_Lat, v_Lat_Weight
  public z_Rad, z_Rad_Weight
  public h_Rad, h_Rad_Weight
  public xv_Lon, xv_Lat
  public xvh_Lon, xvh_Lat, xvh_Rad
  public wh_Rad
  public uz_Rad
  public ut_VMiss

  public u_xv, xv_u
  public at_Dr_at, t_Dr_t, az_at, at_az
  !public xvh_wt, wt_xvh, xvh_wh, wh_xvh
  public xvh_ut, ut_xvh
  public ut_wh, wh_ut, ua_wa, wa_ua
  public wz_uz, uz_wz, wz_wh, wh_wz
  !public wz_wh, wh_wz, wz_uz, uz_wz, az_ah, ah_az
  public wz_wt, wt_wz, wt_ut, ut_wt, uz_ut, ut_uz
  public ut_DRad_ut, ut_DivRad_ut, ut_RotRad_ut
  public ut_Lapla_ut, ua_LaplaH_ua, ua_LaplaHInv_ua
  public xvh_GradLon_ut, xvh_Gradlat_ut
  public ut_DivLon_xvh, ut_DivLat_xvh
  public ut_Div_xvh_xvh_xvh, xvh_Div_xvh_xvh_xvh
  public xvh_RotLon_ut_ut, xvh_RotLat_ut_ut, ut_RotRad_xvh_xvh

  public vh_IntLon_xvh, xh_IntLat_xvh, xv_IntRad_xvh
  public x_IntLatRad_xvh, v_IntLonRad_xvh, h_IntLonLat_xvh
  public IntLonLatRad_xvh

  public h_IntLon_xh, x_IntRad_xh, IntLonRad_xh
  public IntLon_x, IntRad_h, IntRad_z

  public x_IntLat_xv, v_IntLon_xv, IntLonLat_xv
  public h_IntLat_vh, v_IntRad_vh, IntLatRad_vh
  public IntLat_v

  public vh_AvrLon_xvh, xh_AvrLat_xvh, xv_AvrRad_xvh
  public x_AvrLatRad_xvh, v_AvrLonRad_xvh, h_AvrLonLat_xvh
  public AvrLonLatRad_xvh

  public h_AvrLon_xh, x_AvrRad_xh, AvrLonRad_xh
  public AvrLon_x, AvrRad_h, AvrRad_z

  public x_AvrLat_xv, v_AvrLon_xv, AvrLonLat_xv
  public h_AvrLat_vh, v_AvrRad_vh, AvrLatRad_vh
  public AvrLat_v

  public ut_L2_ut, ut_L2Inv_ut
  public ut_RadRot_xvh_xvh, ut_RadRotRot_xvh_xvh_xvh
  public ut_RadRot_RadRotRot
  public ut_Potential2vectorMPI, ut_Potential2RotationMPI
  public ut_Vector2VorDivMPI

  public Interpolate_ut

  public nmz_ToroidalEnergySpectrum_ut, nz_ToroidalEnergySpectrum_ut
  public nmz_PoloidalEnergySpectrum_ut, nz_PoloidalEnergySpectrum_ut

  public ut_Boundaries, ut_TorBoundaries
  public ut_LaplaPol2Pol_ut, uz_LaplaPol2Pol_uz
  public ut_TormagBoundaries, ut_PolmagBoundaries

  public ut_BoundariesTau, ut_TorBoundariesTau, ut_LaplaPol2PolTau_ut
  public ut_TormagBoundariesTau, ut_PolmagBoundariesTau

  public ut_BoundariesGrid, ut_TorBoundariesGrid, ut_LaplaPol2PolGrid_ut
  public ut_TormagBoundariesGrid, ut_PolmagBoundariesGrid

  interface Interpolate_ut
     !
     ! 緯度 alon, 経度 alat 動径 arad における関数値を
     ! そのスペクトルデータ wt_data から補間計算する
     !
     ! 入力する緯度経度座標は,
     !       １点, 経度複数点緯度動径1点
     ! の 2 種類
     !
     module procedure Interpolate_array000_ut
     module procedure a_Interpolate_array100_ut
  end interface

  interface ut_Boundaries
     module procedure ut_BoundariesTau
  end interface

  interface ut_TorBoundaries
     module procedure ut_TorBoundariesTau
  end interface

  interface ut_LaplaPol2Pol_ut
     module procedure ut_LaplaPol2PolTau_ut
  end interface

  interface ut_TorMagBoundaries
     module procedure ut_TorMagBoundariesTau
  end interface

  interface ut_PolMagBoundaries
     module procedure ut_PolMagBoundariesTau
  end interface

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

  real(DP), dimension(:,:,:), allocatable :: xvh_LON, xvh_LAT, xvh_RAD ! 座標
  real(DP), dimension(:,:), allocatable   :: wh_RAD                    ! 座標
  real(DP), dimension(:,:), allocatable   :: uz_RAD                    ! 座標
  real(DP), dimension(:), allocatable     :: h_rad                     ! 座標
  real(DP), dimension(:), allocatable     :: h_rad_weight              ! 座標

  real(DP) :: ut_VMiss = -999.0        ! 欠損値

!!$  integer              :: ks, ke, kc
!!$  integer, allocatable :: ksp(:), kep(:), kcp(:)
!!$
!!$  integer              :: nc
!!$  integer              :: nms, nme, nmc
!!$  integer, allocatable :: nmsp(:), nmep(:), nmcp(:)

  save im, jm, km, nm, lm, ri, ro, ut_VMiss
!!$  save ks, ke, kc, ksp, kep, kcp, nc, nms, nme, nmc, nmsp, nmep, nmcp

contains
  !--------------- 初期化 -----------------
  subroutine ut_mpi_Initial(i,j,k,n,l,r_in,r_out,np,npv)
     !
     ! スペクトル変換の格子点数, 波数, 動径座標の範囲を設定する.
     !
     ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
     ! しなければならない.
     !
     integer,intent(in) :: i              ! 格子点数(経度λ)
     integer,intent(in) :: j              ! 格子点数(緯度φ)
     integer,intent(in) :: k              ! 格子点数(動径 r)
     integer,intent(in) :: n              ! 切断波数(水平全波数)
     integer,intent(in) :: l              ! 切断波数(動径波数)

     real(DP),intent(in) :: r_in           ! 球殻内半径
     real(DP),intent(in) :: r_out          ! 球殻外半径

     integer,intent(in), optional :: np   ! OPENMP での最大スレッド数(dummy)
     integer,intent(in), optional :: npv  ! 鉛直並列プロセス数

     logical :: wa_initialize=.true.   ! wa_initial スイッチ
     integer :: npv_in                 ! 鉛直並列プロセス数

     im = i  ; jm = j ; km = k
     nm = n  ; lm = l
     ri = r_in ; ro = r_out

     if (present(npv)) then
        npv_in = npv
     else
        npv_in = 1
     endif

     if ( wa_initialize ) then
        if ( present(np) ) then
           call ua_mpi_Initial(nm,im,jm,km+1,np=np,npv=npv_in)
        else
           call ua_mpi_Initial(nm,im,jm,km+1,npv=npv_in)
        endif
     endif
     nc = ncc

     call at_Initial(km,lm,r_in,r_out)

     allocate(h_Rad(kc))
     allocate(h_Rad_Weight(kc))

     allocate(xvh_Lon(0:im-1,jc,kc))
     allocate(xvh_Lat(0:im-1,jc,kc))
     allocate(xvh_Rad(0:im-1,jc,kc))

     allocate(wh_Rad(nc,kc))
     allocate(uz_Rad(nmc,0:km))

     z_Rad_Weight = z_Rad_Weight * z_Rad**2     ! r^2 dr の積分重み
     h_Rad = z_Rad(ks-1:ke-1)
     h_Rad_Weight = z_Rad_Weight(ks-1:ke-1)

     xvh_Lon = spread(xv_Lon,3,kc)
     xvh_Lat = spread(xv_Lat,3,kc)
     xvh_Rad = spread(spread(h_Rad,1,jc),1,im)

     wh_Rad = spread(h_Rad,1,nc)
     uz_Rad = spread(z_Rad,1,nmc)

     call MessageNotify('M','ut_mpi_initial', &
          'ut_mpi_module (2021/11/15) is initialized')

   end subroutine ut_mpi_Initial

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

   function ut_xvh(xvh)
     !
     ! 3 次元格子点データからスペクトルデータへ(正)変換する.
     !
     real(DP), dimension(nmc,0:lm)                   :: ut_xvh
     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ

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

     !ut_xvh = ut_uz(uz_wz(wz_wh(wh_xvh(xvh))))
     ut_xvh = ut_uz(ua_xvb(xvh))

   end function ut_xvh

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

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

     xvh_ut = xvb_ua(uz_ut(ut))

   end function xvh_ut

!!$   function wt_xvh(xvh)
!!$     !
!!$     ! 3 次元格子点データからスペクトルデータへ(正)変換する.
!!$     !
!!$     real(DP), dimension(nc,0:lm)             :: wt_xvh
!!$     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$     real(DP), dimension(0:im-1,jc,kc), intent(in)  :: xvh
!!$     !(in) 3 次元経度緯度動径格子点データ
!!$
!!$     wt_xvh = wt_uz(ua_xvb(xvh))
!!$
!!$   end function wt_xvh
!!$
!!$   function xvh_wt(wt)
!!$     !
!!$     ! 水平スペクトル・動径格子点データから 3 次元格子点データへ(逆)変換する.
!!$     !
!!$     real(DP), dimension(0:im-1,jc,kc)              :: xvh_wt
!!$     !(out) 3 次元経度緯度動径格子点データ
!!$
!!$     real(DP), dimension(nc,0:lm), intent(in) :: wt
!!$     !(in) 2 次元球面調和函数スペクトル・動径格子点データ
!!$
!!$     xvh_wt = xvb_ua(uz_wt(wt))
!!$
!!$   end function xvh_wt
!!$
!!$   function wh_xvh(xvh)
!!$     !
!!$     ! 3 次元格子データから水平スペクトル・動径格子点データへ(正)変換する.
!!$     !
!!$     real(DP), dimension(nc,kc)              :: wh_xvh
!!$     !(out) 2 次元球面調和函数スペクトル・動径格子点データ
!!$
!!$     real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh
!!$     !(in) 3 次元経度緯度動径格子点データ
!!$
!!$     wh_xvh = wb_xvb(xvh)
!!$
!!$   end function wh_xvh
!!$
!!$   function xvh_wh(wh)
!!$     !
!!$     ! 水平スペクトル・動径格子点データから 3 次元格子データへ(逆)変換する.
!!$     !
!!$     real(DP), dimension(0:im-1,jc,kc)            :: xvh_wh
!!$     !(out) 2 次元球面調和函数スペクトル・動径格子点データ
!!$
!!$     real(DP), dimension(nc,kc), intent(in) :: wh
!!$     !(in) 3 次元経度緯度動径格子点データ
!!$
!!$     xvh_wh = xvb_wb(wh)
!!$
!!$   end function xvh_wh
!!$
!!$   function u_xv(xv)
!!$     !
!!$     ! 2 次元格子データから水平スペクトル(正)変換する.
!!$     !
!!$     real(DP), dimension(nmc)                    :: u_xv
!!$     !(out) 2 次元球面調和函数スペクトル・動径格子点データ
!!$
!!$     real(DP), dimension(0:im-1,jc), intent(in) :: xv
!!$     !(in) 3 次元経度緯度動径格子点データ
!!$
!!$     u_xv = reshape(ua_wa(reshape(w_xv(xv),(/nc,1/))),(/nmc/))
!!$
!!$   end function u_xv
!!$
!!$   function xv_u(u)
!!$     !
!!$     ! 水平スペクトルデータから 2 次元格子データへ(逆)変換する.
!!$     !
!!$     real(DP), dimension(0:im-1,jc)      :: xv_u
!!$     !(out) 2 次元球面調和函数スペクトル・動径格子点データ
!!$
!!$     real(DP), dimension(nmc), intent(in) :: u
!!$     !(in) 3 次元経度緯度動径格子点データ
!!$
!!$     xv_u = xv_w(reshape(wa_ua(reshape(u,(/nmc,1/))),(/nc/)))
!!$
!!$   end function xv_u

   function ut_wh(wh)
     !
     ! 水平スペクトル・動径格子点データを集める.
     !
     real(DP), dimension(nc,kc), intent(in) :: wh
     !(in)  2 次元球面調和函数スペクトル・動径格子点データ(分割)
     real(DP), dimension(nmc,0:lm)                 :: ut_wh
     !(out) 2 次元球面調和函数スペクトル・動径格子点データ

     !ut_wh = ut_uz(uz_wz(wz_wh(wh)))
     ut_wh = ut_uz(uz_wh(wh))

   end function ut_wh

   function wh_ut(ut)
     !
     ! 水平スペクトル・動径格子点データを分割する.
     !
     real(DP), dimension(nmc,0:lm), intent(in)       :: ut
     !(in) 2 次元球面調和函数スペクトル・動径格子点データ
     real(DP), dimension(nc,kc)               :: wh_ut
     !(out) 2 次元球面調和函数スペクトル・動径格子点データ(分割)

     wh_ut = wh_uz(uz_ut(ut))

   end function wh_ut

   function wz_wh(wh)
     !
     ! 水平スペクトル・動径格子点データを集める.
     !
     real(DP), dimension(nc,kc), intent(in) :: wh
     !(in)  2 次元球面調和函数スペクトル・動径格子点データ(分割)
     real(DP), dimension(nc,0:km)           :: wz_wh
     !(out) 2 次元球面調和函数スペクトル・動径格子点データ

     wz_wh = aa_ab(wh)

   end function wz_wh

   function wh_wz(wz)
     !
     ! 水平スペクトル・動径格子点データを分割する.
     !
     real(DP), dimension(nc,0:km), intent(in) :: wz
     !(in) 2 次元球面調和函数スペクトル・動径格子点データ
     real(DP), dimension(nc,kc)               :: wh_wz
     !(out) 2 次元球面調和函数スペクトル・動径格子点データ(分割)

     wh_wz = ab_aa(wz)

   end function wh_wz

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

     wz_wt = az_at(wt)

   end function wz_wt

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

     wt_wz = at_az(wz)

   end function wt_wz

   function ut_wt(wt)
     !
     ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する.
     !
     real(DP), dimension(nc,0:lm), intent(in) :: wt
     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
     real(DP), dimension(nmc,0:lm)                   :: ut_wt
     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ(水平再分割)

     ut_wt = ua_wa(wt)

   end function ut_wt

   function wt_ut(ut)
     !
     ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
     !
     real(DP), dimension(nmc,0:lm), intent(in)       :: ut
     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ(水平再分割)
     real(DP), dimension(nc,0:lm)             :: wt_ut
     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ

     wt_ut = wa_ua(ut)

   end function wt_ut

   function uz_wz(wz)
     !
     ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する.
     !
     real(DP), dimension(nc,0:km), intent(in) :: wz
     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
     real(DP), dimension(nmc,0:km)                   :: uz_wz
     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ(水平再分割)

     uz_wz = ua_wa(wz)

   end function uz_wz

   function wz_uz(uz)
     !
     ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
     !
     real(DP), dimension(nmc,0:km), intent(in)       :: uz
     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ(水平再分割)
     real(DP), dimension(nc,0:km)             :: wz_uz
     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ

     wz_uz = wa_ua(uz)

   end function wz_uz

   function uz_ut(ut)
     !
     ! スペクトルデータから水平スペクトル・動径格子点データへ(正)変換する.
     !
     real(DP), dimension(nmc,0:lm), intent(in) :: ut
     !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
     real(DP), dimension(nmc,0:km)             :: uz_ut
     !(out) 2 次元球面調和函数スペクトル・動径格子点データ

     uz_ut = az_at(ut)

   end function uz_ut

   function ut_uz(uz)
     !
     ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
     !
     real(DP), dimension(nmc,0:km), intent(in) :: uz
     !(in) 2 次元球面調和函数スペクトル・動径格子点データ
     real(DP), dimension(nmc,0:lm)             :: ut_uz
     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ

     ut_uz = at_az(uz)

   end function ut_uz

!!$   function wa_ua(ua)
!!$     !
!!$     ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
!!$     !
!!$     real(DP), dimension(:,:), intent(in)       :: ua
!!$     !(in) 2 次元球面調和函数(水平再分割) (nc,:)
!!$     real(DP), dimension(nc,size(ua,2))   :: wa_ua
!!$     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$     real(DP), dimension(size(ua,2),nmc)  :: au
!!$     real(DP), dimension(size(ua,2),nc)   :: aw
!!$
!!$     integer :: a_size
!!$     integer :: IRQSTS
!!$     integer :: IST( MPI_STATUS_SIZE )
!!$
!!$     integer :: ierr, ip_v
!!$     integer, parameter :: itag=110
!!$     ! integer, allocatable ::snd_rank(:)
!!$     integer, allocatable ::snd_rank_v(:)
!!$     ! integer, allocatable ::rcv_rank(:)
!!$     integer, allocatable ::rcv_rank_v(:)
!!$     integer, allocatable ::rcv_idx(:)
!!$     logical :: first=.true.
!!$
!!$     save first, snd_rank_v, rcv_rank_v, rcv_idx !, snd_rank, rcv_rank
!!$
!!$     if ( first ) then
!!$        first = .false.
!!$        ! allocate(snd_rank(1:nproc_v-1))
!!$        allocate(snd_rank_v(1:nproc_v-1))
!!$        ! allocate(rcv_rank(1:nproc_v-1))
!!$        allocate(rcv_rank_v(1:nproc_v-1))
!!$        allocate(rcv_idx(0:nproc_v-1))
!!$
!!$        do ip_v=1,nproc_v-1
!!$           rcv_rank_v(ip_v) = mod(iproc_v+ip_v,nproc_v)
!!$           ! rcv_rank(ip_v) = myid_h + nproc_h*rcv_rank_v(ip_v)
!!$           snd_rank_v(ip_v) = mod(iproc_v-ip_v+nproc_v,nproc_v)
!!$           ! snd_rank(ip_v) = myid_h + nproc_h*snd_rank_v(ip_v)
!!$        enddo
!!$
!!$        rcv_idx(0) = 1
!!$        do ip_v=1,nproc_v-1
!!$           rcv_idx(ip_v) = rcv_idx(ip_v-1)+nmcp(ip_v-1)
!!$        enddo
!!$     endif
!!$
!!$     au = transpose(ua)
!!$     a_size = size(ua,2)
!!$
!!$     call MPI_Barrier(MPI_COMM_V,IERR)
!!$     do ip_v=1,nproc_v-1
!!$        call MPI_ISEND(au, a_size*nmc, MPI_REAL8, rcv_rank_v(ip_v),&
!!$             itag,MPI_COMM_V,IRQSTS,IERR)
!!$        call MPI_RECV(aw(1,rcv_idx(snd_rank_v(ip_v))),  &
!!$             a_size*nmcp(snd_rank_v(ip_v)), MPI_REAL8,   &
!!$             snd_rank_v(ip_v), itag, MPI_COMM_V, IST, IERR)
!!$        CALL MPI_WAIT(IRQSTS,IST,IERR)
!!$     enddo
!!$
!!$     aw(:,rcv_idx(iproc_v):rcv_idx(iproc_v)+nmcp(iproc_v)-1) = au
!!$
!!$     wa_ua = transpose(aw)
!!$
!!$   end function wa_ua
!!$
!!$   function ua_wa(wa)
!!$     !
!!$     ! 水平スペクトルデータを再分割する
!!$     !
!!$     real(DP), dimension(:,:), intent(in) :: wa
!!$     !(in) 2 次元球面調和函数データ (nc,:)
!!$     real(DP), dimension(nmc,size(wa,2))  :: ua_wa
!!$     !(out) 2 次元球面調和函数データ(水平再分割)
!!$
!!$     ua_wa = wa(nms:nme,:)
!!$
!!$   end function ua_wa
!!$

   function uz_wh(wh)
     !
     ! 水平スペクトル・動径格子点データを集める.
     !
     real(DP), dimension(nc,kc), intent(in) :: wh
     !(in)  動径格子点データ(分割)
     real(DP), dimension(nmc,0:km) :: uz_wh
     !(out) 動径格子点データ

     uz_wh = ua_wb(wh)

   end function uz_wh

   function wh_uz(uz)
     !
     ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
     !
     real(DP), dimension(nmc,0:km), intent(in)       :: uz
     !(in) 2 次元球面調和函数(水平再分割) (nc,:)
     real(DP), dimension(nc,kc)               :: wh_uz
     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ

     wh_uz = wb_ua(uz)

   end function wh_uz

   !--------------- 微分計算 -----------------
   function xvh_GradLon_ut(ut)
      !
      ! スペクトルデータに勾配型経度微分 1/rcosφ・∂/∂λ
      ! を作用させる.
      !
      real(DP), dimension(nmc,0:lm), intent(in) :: ut
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

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

      xvh_GradLon_ut = xvb_GradLon_ua(uz_ut(ut))/xvh_Rad

    end function xvh_GradLon_ut

    function xvh_GradLat_ut(ut)
      !
      ! スペクトルデータに勾配型経度微分 1/r ∂/∂φ を作用させる.
      !
      real(DP), dimension(nmc,0:lm), intent(in) :: ut
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

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

      xvh_GradLat_ut = xvb_GradLat_ua(uz_ut(ut))/xvh_Rad
    end function xvh_GradLat_ut

    function ut_DivLon_xvh(xvh)
      !
      ! 格子点データに発散型経度微分 1/rcosφ・∂/∂λ を作用させた
      ! スペクトルデータを返す.
      !
      real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh
      !(in) 3 次元経度緯度動径格子点データ

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

      ut_DivLon_xvh = ut_uz(ua_DivLon_xvb(xvh/xvh_Rad))

    end function ut_DivLon_xvh

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

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

      ut_DivLat_xvh = ut_uz(ua_divlat_xvb(xvh/xvh_Rad))
    end function ut_DivLat_xvh

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

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

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

      real(DP), dimension(nmc,0:lm)                    :: ut_Div_xvh_xvh_xvh
      !(out) ベクトル場の発散

      ut_Div_xvh_xvh_xvh =   ut_DivLon_xvh(xvh_Vlon) &
                           + ut_DivLat_xvh(xvh_Vlat) &
                           + ut_DivRad_ut(ut_xvh(xvh_Vrad))

    end function ut_Div_xvh_xvh_xvh

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

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

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

      real(DP), dimension(0:im-1,jc,kc)             :: xvh_Div_xvh_xvh_xvh
      !(out) ベクトル場の発散

      xvh_Div_xvh_xvh_xvh &
           = xvh_Rad/cos(xvh_Lat) &
                * xvh_ut(ut_Div_xvh_xvh_xvh(xvh_VLon*cos(xvh_Lat)/xvh_Rad,  &
                                            xvh_VLat*cos(xvh_Lat)/xvh_Rad,  &
                                            xvh_VRad*cos(xvh_Lat)/xvh_Rad ))&
             + xvh_VLat*tan(xvh_Lat)/xvh_Rad &
             + xvh_VRad/xvh_Rad

    end function xvh_Div_xvh_xvh_xvh

    function xvh_RotLon_ut_ut(ut_Vrad,ut_Vlat)
      !
      ! ベクトル場の動径成分, 緯度成分である第 1, 2 引数 Vrad, Vlat から
      ! 回転の経度成分
      !
      !    1/r ∂Vrad/∂φ-1/r ∂(r Vlat)/∂r を計算する.
      !
      ! を計算する
      !
      real(DP), dimension(nmc,0:lm), intent(in) :: ut_Vrad
      !(in) ベクトル場の動径成分

      real(DP), dimension(nmc,0:lm), intent(in) :: ut_Vlat
      !(in) ベクトル場の緯度成分

      real(DP), dimension(0:im-1,jc,kc)      :: xvh_RotLon_ut_ut
      !(out) ベクトル場の回転の経度成分

        xvh_RotLon_ut_ut =   xvh_GradLat_ut(ut_Vrad) &
                           - xvh_ut(ut_RotRad_ut(ut_Vlat))

    end function xvh_RotLon_ut_ut

    function xvh_RotLat_ut_ut(ut_Vlon,ut_Vrad)
      !
      ! ベクトル場の経度成分, 動径成分である第 1, 2 引数 Vlon, Vrad から
      ! 回転の緯度成分
      !
      !    1/r ∂(r Vlon)/∂r - 1/rcosφ・∂Vrad/∂λ
      !
      ! を計算する.
      !
      real(DP), dimension(nmc,0:lm), intent(in) :: ut_Vlon
      !(in) ベクトル場の経度成分

      real(DP), dimension(nmc,0:lm), intent(in) :: ut_Vrad
      !(in) ベクトル場の動径成分

      real(DP), dimension(0:im-1,jc,kc)      :: xvh_RotLat_ut_ut
      !(out) ベクトル場の回転の緯度成分

        xvh_RotLat_ut_ut =   xvh_ut(ut_RotRad_ut(ut_Vlon)) &
                           - xvh_GradLon_ut(ut_Vrad)

    end function xvh_RotLat_ut_ut

    function ut_RotRad_xvh_xvh(xvh_Vlat,xvh_Vlon)
      !
      ! ベクトルの緯度成分, 経度成分である第 1, 2 引数 Vlat, Vlon に対して
      ! ベクトル場の回転の動径成分
      !
      !    1/rcosφ・∂Vlat/∂λ - 1/rcosφ・∂(Vlon cosφ)/∂φ
      !
      ! を計算する.
      !
      real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_Vlat
      !(in) ベクトル場の緯度成分

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

      real(DP), dimension(nmc,0:lm)                    :: ut_RotRad_xvh_xvh
      !(out) ベクトル場の回転の動径成分

        ut_RotRad_xvh_xvh =   ut_DivLon_xvh(xvh_Vlat) &
                            - ut_DivLat_xvh(xvh_Vlon)

    end function ut_RotRad_xvh_xvh

  !--------------- ポロイダル/トロイダルモデル用微分 -----------------
    function ut_DRad_ut(ut)
      !
      ! 入力スペクトルデータに動径微分 ∂/∂r を作用する.
      !
      ! スペクトルデータの動径微分とは, 対応する格子点データに動径微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      real(DP), dimension(nmc,0:lm), intent(in) :: ut
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(DP), dimension(nmc,0:lm)             :: ut_DRad_ut
      !(in) 動径微分された2 次元球面調和函数チェビシェフスペクトルデータ

      ut_DRad_ut = at_Dr_at(ut)

    end function ut_DRad_ut

    function ut_DivRad_ut(ut)
      !
      ! 入力スペクトルデータに発散型動径微分
      !
      !       1/r^2 ∂/∂r (r^2 .)= ∂/∂r + 2/r
      !
      ! を作用する.
      !
      ! スペクトルデータの発散型動径微分とは, 対応する格子点データに
      ! 発散型動径微分を作用させたデータのスペクトル変換のことである.
      !
      real(DP), dimension(nmc,0:lm), intent(in) :: ut
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

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

!      ut_DivRad_ut = ut_uz(uz_ut(ut_Drad_ut(ut_uz(uz_Rad**2*uz_ut(ut))))/uz_Rad**2)

      ut_DivRad_ut = ut_Drad_ut(ut) + ut_uz(2/uz_rad*uz_ut(ut))

    end function ut_DivRad_ut

    function ut_RotRad_ut(ut)
      !
      ! 入力スペクトルデータに回転型動径微分
      !
      !      1/r ∂(r.)/∂r = ∂(.)/∂r + (.)/r
      !
      ! を作用する.
      !
      ! スペクトルデータの回転型動径微分とは, 対応する格子点データに
      ! 回転型動径微分を作用させたデータのスペクトル変換のことである.
      !
      real(DP), dimension(nmc,0:lm), intent(in) :: ut
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

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

!      ut_RotRad_ut = ut_uz(uz_ut(ut_Drad_ut(ut_uz(uz_Rad*uz_ut(ut))))/uz_Rad)
      ut_RotRad_ut = ut_Drad_ut(ut) + ut_uz(1/uz_Rad*uz_ut(ut))

    end function ut_RotRad_ut

    function ua_LaplaH_ua(ua_data)
      !
      ! 入力スペクトルデータにラプラシアン
      !
      !    ▽^2 = 1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ)
      !
      ! を作用する(多層用).
      !
      ! スペクトルデータのラプラシアンとは, 対応する格子点データに
      ! ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      real(DP), intent(in)  :: ua_data(:,:)
      !(in) 入力スペクトルデータ
      real(DP)              :: ua_LaplaH_ua(nmc,size(ua_data,2))
      !(out) 入力スペクトルデータのラプラシアン

      ua_LaplaH_ua  = ua_Lapla_ua(ua_data)
    end function ua_LaplaH_ua

    function ua_LaplaHInv_ua(ua_data)
      !
      ! 入力スペクトルデータにラプラシアン
      !
      !    ▽^2 = 1/cos^2φ・∂^2/∂λ^2 + 1/cosφ・∂/∂φ(cosφ∂/∂φ)
      !
      ! を作用する(多層用).
      !
      ! スペクトルデータのラプラシアンとは, 対応する格子点データに
      ! ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      real(DP), intent(in)  :: ua_data(:,:)
      !(in) 入力スペクトルデータ
      real(DP)              :: ua_LaplaHInv_ua(nmc,size(ua_data,2))
      !(out) 入力スペクトルデータのラプラシアン

      ua_LaplaHInv_ua  = ua_LaplaInv_ua(ua_data)

    end function ua_LaplaHInv_ua

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

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

      ut_Lapla_ut = ut_DivRad_ut(ut_Drad_ut(ut)) &
                   + ut_uz(uz_ut(ua_LaplaH_ua(ut))/uz_Rad**2)

    end function ut_Lapla_ut

!!$  !--------------- ポロイダル/トロイダルモデル用微分 -----------------
!!$
!!$    function wt_KxRGrad_wt(wt)
!!$      !
!!$      ! 入力スペクトルデータに経度微分 k×r・▽ = ∂/∂λを作用する.
!!$      !
!!$      real(DP), dimension(nc,0:lm), intent(in) :: wt
!!$      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$      real(DP), dimension(nc,0:lm)             :: wt_KxRGrad_wt
!!$      !(out) 経度微分を作用された 2 次元スペクトルデータ
!!$
!!$      wt_KxRGrad_wt =  wa_Dlon_wa(wt)
!!$
!!$    end function wt_KxRGrad_wt
!!$
!!$    function xvz_KGrad_wt(wt)    ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r
!!$      !
!!$      ! 入力スペクトルデータに対応する格子データに軸方向微分
!!$      !
!!$      !    k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r
!!$      !
!!$      ! を作用させた格子データが返される.
!!$      ! ここでベクトル k は球の中心から北極向きの単位ベクトルである.
!!$      !
!!$      real(DP), dimension(nc,0:lm), intent(in) :: wt
!!$      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$      real(DP), dimension(0:im-1,jc,0:km)                     :: xvz_KGrad_wt
!!$      !(out) 軸方向微分を作用された 2 次元スペクトルデータ
!!$
!!$      xvz_KGrad_wt =  cos(xvz_Lat)*xvz_GradLat_wt(wt) &
!!$                    + sin(xvz_Lat)*xvz_wt(wt_Drad_wt(wt))
!!$    end function xvz_KGrad_wt

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

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

      ut_L2_ut = -ua_LaplaH_ua(ut)

    end function ut_L2_ut

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

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

      ut_L2Inv_ut = -ua_LaplaHInv_ua(ut)

    end function ut_L2Inv_ut

!!$    function wt_QOperatorMPI_wt(wt)
!!$      !
!!$      ! 入力スペクトルデータに対応する格子点データに演算子
!!$      !
!!$      !    Q=(k・▽-1/2(L2 k・▽+ k・▽L2))
!!$      !
!!$      ! を作用させたデータのスペクトル変換が返される.
!!$      !
!!$      real(DP), dimension(nc,0:lm), intent(in) :: wt
!!$      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$      real(DP), dimension(nc,0:lm)             :: wt_QOperatorMPI_wt
!!$      !(out) Q 演算子を作用された 2 次元スペクトルデータ
!!$
!!$      wt_QOperatorMPI_wt = &
!!$             wt_xvz(xvz_KGrad_wt(wt) - xvz_KGrad_wt(wt_L2_wt(wt))/2) &
!!$           - wt_L2_wt(wt_xvz(xvz_KGrad_wt(wt)))/2
!!$
!!$    end function wt_QOperatorMPI_wt

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

      real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLAT
      !(in) ベクトルの緯度成分

      real(DP), dimension(nmc,0:lm)     :: ut_RadRot_xvh_xvh
      !(out) ベクトルの渦度と動径ベクトルの内積

      ut_RadRot_xvh_xvh = ut_uz(   ua_DivLon_xvb(xvh_VLAT) &
                                 - ua_DivLat_xvb(xvh_VLON) )
!!$      ut_RadRot_xvh_xvh = ut_uz(uz_wz(wz_wh(wb_DivLon_xvb(xvh_VLAT) &
!!$                                - wb_DivLat_xvb(xvh_VLON))))

    end function ut_RadRot_xvh_xvh

    function ut_RadRotRot_xvh_xvh_xvh(xvh_VLON,xvh_VLAT,xvh_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(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLON
      !(in) ベクトルの経度成分

      real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLAT
      !(in) ベクトルの緯度成分

      real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VRAD
      !(in) ベクトルの動径成分

      real(DP), dimension(nmc,0:lm)     :: ut_RadRotRot_xvh_xvh_xvh
      !(out) ベクトル v の r・(▽×▽×v)

      ut_RadRotRot_xvh_xvh_xvh = &
               ut_RotRad_ut(ut_uz( &
                   (ua_DivLon_xvb(xvh_VLON)+ ua_DivLat_xvb(xvh_VLAT)) )) &
             + ut_L2_ut(ut_xvh(xvh_VRAD/xvh_RAD))

    end function ut_RadRotRot_xvh_xvh_xvh


    subroutine ut_RadRot_RadRotRot( &
         xvh_VLON,xvh_VLAT,xvh_VRAD, ut_Tor, ut_LaplaPol )
      
      !
      ! ベクトル v に対して r・(▽×v) と r・(▽×▽×v) を計算する.
      !
      ! 第 1, 2, 3 引数(v[λ], v[φ], v[r])がそれぞれベクトルの経度成分,
      ! 緯度成分, 動径成分を表す.
      !
      !    r・(▽×v) = 1/cosφ・∂v[φ]/∂λ - 1/cosφ・∂(v[λ] cosφ)/∂φ
      !
      !    r・(▽×▽×v)  = 1/r ∂/∂r (r・( 1/cosφ・∂v[λ]/∂λ
      !                                  + 1/cosφ・∂(v[φ] cosφ)/∂φ ) )
      !                     + L^2 v[r]/r
      !
      ! のスペクトルデータが返される.
      !
      ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回少なくてすむ.
      !
      real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLON
      !(in) ベクトルの経度成分

      real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VLAT
      !(in) ベクトルの緯度成分

      real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh_VRAD
      !(in) ベクトルの動径成分

      real(DP), dimension(nmc,0:lm)     :: ut_Tor
      !(out) ベクトル v の r・(▽×v)
      
      real(DP), dimension(nmc,0:lm)     :: ut_LaplaPol
      !(out) ベクトル v の r・(▽×▽×v)

      call ut_Vector2VorDivMPI(xvh_VLon, xvh_VLat, ut_Tor, ut_LaplaPol)

      ut_LaplaPol = ut_RotRad_ut(ut_LaplaPol) &
                  + ut_L2_ut(ut_xvh(xvh_VRAD/xvh_RAD))

    end subroutine ut_RadRot_RadRotRot

    subroutine ut_Potential2VectorMPI(&
         xvh_VLON,xvh_VLAT,xvh_VRAD,ut_TORPOT,ut_POLPOT)
      !
      ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
      !
      !     v = ▽x(Ψr) + ▽x▽x(Φr)
      !
      ! の各成分を計算する
      !
      real(DP), dimension(0:im-1,jc,kc)     :: xvh_VLON
      !(out) ベクトル場の経度成分

      real(DP), dimension(0:im-1,jc,kc)     :: xvh_VLAT
      !(out) ベクトル場の緯度成分

      real(DP), dimension(0:im-1,jc,kc)     :: xvh_VRAD
      !(out) ベクトル場の動径成分

      real(DP), dimension(nmc,0:lm), intent(in) :: ut_TORPOT
      !(in) トロイダルポテンシャル

      real(DP), dimension(nmc,0:lm), intent(in) :: ut_POLPOT
      !(in) ポロイダルポテンシャル


      real(DP), dimension(nmc,0:km) :: uz_Tor
      !(in) トロイダルポテンシャル

      real(DP), dimension(nmc,0:km) :: uz_DPolDrad
      !(in) ポロイダルポテンシャル

      uz_Tor = uz_ut(ut_TORPOT)
      uz_DPolDRad = uz_ut(ut_RotRad_ut(ut_POLPOT))

      call ua_StreamPotential2VectorMPI(-uz_Tor, uz_DPolDRad, xvh_VLON, xvh_VLAT)
      xvh_VRAD = xvh_ut(ut_L2_ut(ut_POLPOT))/xvh_RAD


!!$      xvh_VLON =   xvh_RAD * xvh_GradLat_ut(ut_TORPOT) &
!!$                 + xvb_GradLon_ua(uz_ut(ut_RotRad_ut(ut_POLPOT)))
!!$      xvh_VLAT = - xvh_RAD * xvh_GradLon_ut(ut_TORPOT) &
!!$                 + xvb_GradLat_ua(uz_ut(ut_RotRad_ut(ut_POLPOT)))
!!$      xvh_VRAD = xvh_ut(ut_L2_ut(ut_POLPOT))/xvh_RAD
!!$      xvh_VLON =   xvh_RAD * xvh_GradLat_ut(ut_TORPOT) &
!!$                 + xvb_GradLon_wb(wh_wz(wz_uz(uz_ut(ut_RotRad_ut(ut_POLPOT)))))
!!$      xvh_VLAT = - xvh_RAD * xvh_GradLon_ut(ut_TORPOT) &
!!$                 + xvb_GradLat_wb(wh_wz(wz_uz(uz_ut(ut_RotRad_ut(ut_POLPOT)))))
!!$      xvh_VRAD = xvh_ut(ut_L2_ut(ut_POLPOT))/xvh_RAD

    end subroutine ut_Potential2VectorMPI

    subroutine ut_Potential2RotationMPI(&
       xvh_RotVLON,xvh_RotVLAT,xvh_RotVRAD,ut_TORPOT,ut_POLPOT)
      !
      ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
      !
      !     v = ▽x(Ψr) + ▽x▽x(Φr)
      !
      ! に対して, その回転
      !
      !     ▽xv = ▽x▽x(Ψr) + ▽x▽x▽x(Φr) = ▽x▽x(Ψr) - ▽x((▽^2Φ)r)
      !
      ! を計算する.

      ! ベクトル場の回転
      real(DP), dimension(0:im-1,jc,kc), intent(OUT) :: xvh_RotVLON
      !(out) 回転の経度成分

      real(DP), dimension(0:im-1,jc,kc), intent(OUT) :: xvh_RotVLAT
      !(out) 回転の緯度成分

      real(DP), dimension(0:im-1,jc,kc), intent(OUT) :: xvh_RotVRAD
      !(out) 回転の動径成分

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

      real(DP), dimension(nmc,0:lm), intent(in) :: ut_POLPOT
      !(in) ポロイダルポテンシャル

      call ut_Potential2VectorMPI( &
           xvh_RotVLON,xvh_RotVLAT,xvh_RotVRAD, &
           -ut_Lapla_ut(ut_POLPOT), ut_TORPOT)

    end subroutine ut_Potential2RotationMPI

    subroutine ut_Vector2VorDivMPI(xvh_VLon, xvh_VLat, ut_Vor, ut_Div)
      !
      ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に
      ! (正)変換する(多層用)
      !
      ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
      !
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      real(8), intent(in)   :: xvh_VLon(0:im-1,jc,kc)
      !(in) 速度経度成分

      real(8), intent(in)   :: xvh_VLat(0:im-1,jc,kc)
      !(in) 速度緯度成分

      real(8), intent(out)   :: ut_Vor(nmc,0:lm)
      !(out) 渦度
      real(8), intent(out)   :: ut_Div(nmc,0:lm)
      !(out) 発散

      real(8) :: wh_Vor(nc,kc)
      !速度経度成分

      real(8) :: wh_Div(nc,kc)
      !速度緯度成分

      integer :: k

      do k=1,kc
         call w_Vector2VorDivMPI &
              (xvh_VLon(:,:,k), xvh_VLat(:,:,k), wh_Vor(:,k), wh_Div(:,k))
      enddo

      ut_Vor = ut_wh(wh_Vor)
      ut_Div = ut_wh(wh_Div)

    end subroutine ut_Vector2VorDivMPI

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

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

      integer :: i

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

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

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

      xh_IntLat_xvh = xb_IntLat_xvb(xvh)

    end function xh_IntLat_xvh

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

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

      real(DP), dimension(0:im-1,jc)  :: xv_IntRadTmp

      integer :: i,j,k,ierr

      xv_IntRadTmp = 0.0d0
      do k=1,kc
         do j=1,jc
            do i=0,im-1
               xv_IntRadTmp(i,j) = xv_IntRadTmp(i,j) &
                    + xvh(i,j,k) * h_Rad_Weight(k)
            enddo
         enddo
      end do

      CALL MPI_ALLREDUCE(xv_IntRadTmp,xv_IntRad_xvh,im*jc,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_V,IERR)

    end function xv_IntRad_xvh

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

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

      x_IntLatRad_xvh = x_IntRad_xh(xb_IntLat_xvb(xvh))

    end function x_IntLatRad_xvh

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

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

      v_IntLonRad_xvh = v_IntLon_xv(xv_IntRad_xvh(xvh))

    end function v_IntLonRad_xvh

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

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

      h_IntLonLat_xvh = b_IntLonLat_xvb(xvh)

    end function h_IntLonLat_xvh

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

      real(DP)                     :: IntLonLatRad_xvh
      !(out) 全球積分値

      IntLonLatRad_xvh = IntLonLat_xv(xv_IntRad_xvh(xvh))

    end function IntLonLatRad_xvh

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

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

      h_IntLat_vh = b_IntLat_vb(vh)

    end function h_IntLat_vh

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

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

      integer :: j

      do j=1,jc
         v_IntRad_vh(j) = IntRad_h(vh(j,:))
      enddo

    end function v_IntRad_vh

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

      real(DP)                   :: IntLatRad_vh
      !(out) 積分値

      !IntLatRad_vh = IntRad_h(b_IntLat_vb(vh))
      IntLatRad_vh = IntLat_v(v_IntRad_vh(vh))

    end function IntLatRad_vh

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

      real(DP), dimension(kc)  :: h_IntLon_xh
      !(out) 経度積分された 1 次元動径格子点データ

      integer :: i

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

    end function h_IntLon_xh

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

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

      integer :: i

      do i=0,im-1
         x_IntRad_xh(i) = IntRad_h(xh(i,:))
      enddo

    end function x_IntRad_xh

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

      real(DP)                                 :: IntLonRad_xh
      !(out) 積分値

      IntLonRad_xh = IntLon_x(x_IntRad_xh(xh))

    end function IntLonRad_xh

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

      real(DP)                            :: IntRad_h
      !(out) 積分値

      real(DP) :: IntRadTmp
      integer :: k, ierr

      IntRadTmp = 0.0d0
      do k=1,kc
         IntRadTmp = IntRadTmp + h(k) * h_Rad_Weight(k)
      enddo

      CALL MPI_ALLREDUCE(IntRadTmp,IntRad_h,1,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_V,IERR)


    end function IntRad_h

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

      real(DP)                            :: 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

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

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

      vh_AvrLon_xvh = vh_IntLon_xvh(xvh)/sum(x_Lon_Weight)

    end function vh_AvrLon_xvh

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

      real(DP), dimension(im,kc)  :: xh_AvrLat_xvh
      !(out) 緯度平均された 2 次元緯度動径格子点データ

      xh_AvrLat_xvh = xh_IntLat_xvh(xvh)/sum(y_Lat_Weight)

    end function xh_AvrLat_xvh

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

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

      xv_AvrRad_xvh = xv_IntRad_xvh(xvh)/sum(z_Rad_Weight)

    end function xv_AvrRad_xvh

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

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

      x_AvrLatRad_xvh = x_IntLatRad_xvh(xvh) &
                   /( sum(y_Lat_Weight)*sum(z_Rad_Weight) )

    end function x_AvrLatRad_xvh

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

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

      v_AvrLonRad_xvh = v_IntLonRad_xvh(xvh) &
                 /(sum(x_Lon_Weight)*sum(z_Rad_Weight))

    end function v_AvrLonRad_xvh

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

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

      h_AvrLonLat_xvh = h_IntLonLat_xvh(xvh) &
                 /(sum(x_Lon_Weight)*sum(y_Lat_Weight))

    end function h_AvrLonLat_xvh

    function AvrLonLatRad_xvh(xvh) ! 緯度経度動径(全球)積分
      !
      ! 3 次元格子点データの緯度経度動径(全球)積分
      !
      ! 3 次元データ f(λ,φ,r) に対して
      !
      !    ∫f(λ,φ,r) r^2cosφ dλdφdr /(4π(r[o]^3-r[i]^3)/3)
      !
      ! を計算する.
      !
      real(DP), dimension(0:im-1,jc,kc), intent(in) :: xvh
      !(in) 3 次元経度緯度動径格子点データ

      real(DP)                     :: AvrLonLatRad_xvh
      !(out) 全球平均値

      AvrLonLatRad_xvh = IntLonLatRad_xvh(xvh) &
            /(sum(x_Lon_Weight)*sum(y_Lat_Weight) * sum(z_Rad_Weight))

    end function AvrLonLatRad_xvh

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

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

      h_AvrLat_vh = h_IntLat_vh(vh)/sum(y_Lat_Weight)
    end function h_AvrLat_vh

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

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

      v_AvrRad_vh = v_IntRad_vh(vh)/sum(z_Rad_Weight)

    end function v_AvrRad_vh

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

      real(DP)                   :: AvrLatRad_vh
      !(out) 平均値

      AvrLatRad_vh = IntLatRad_vh(vh)/(sum(y_Lat_Weight)*sum(z_Rad_Weight))

    end function AvrLatRad_vh

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

      real(DP), dimension(kc)  :: h_AvrLon_xh
      !(out) 経度平均された 1 次元動径格子点データ

      h_AvrLon_xh = h_IntLon_xh(xh)/sum(x_Lon_Weight)

    end function h_AvrLon_xh

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

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

      x_AvrRad_xh = x_IntRad_xh(xh)/sum(z_Rad_Weight)

    end function x_AvrRad_xh

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

      AvrLonRad_xh = IntLonRad_xh(xh)/(sum(x_Lon_Weight)*sum(z_Rad_Weight))

    end function AvrLonRad_xh

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

      AvrRad_h = IntRad_h(h)/sum(z_Rad_Weight)

    end function AvrRad_h

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

      AvrRad_z = IntRad_z(z)/sum(z_Rad_Weight)

    end function AvrRad_z

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

      real(DP) :: wa_data(nc,1)                  ! 補間した値

      wa_data = wa_ua(reshape(a_Interpolate_at(ut_data,arad),(/nmc,1/)))
      Interpolate_array000_ut = Interpolate_w(wa_data(:,1),alon,alat)

    end function Interpolate_array000_ut

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

      real(DP) :: wa_data(nc,1)                  ! 補間した値

      wa_data = wa_ua(reshape(a_Interpolate_at(ut_data,arad),(/nmc,1/)))
      a_Interpolate_array100_ut = Interpolate_w(wa_data(:,1),alon,alat)

    end function a_Interpolate_array100_ut

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

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

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

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

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

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

      wz_DATA = wz_wt(wt_ut(ut_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) * z_Rad**2 * wz_Data(l,:)**2
                  !  + 0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_Data(l,:)**2
            else
               nmz_Work(n,m,:) = nmz_Work(n,m,:) &
                    +  0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_Data(l,:)**2
                  !  +  0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_Data(l,:)**2
               nmz_Work(n,-m,:) = nmz_Work(n,-m,:) &
                    +  0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_Data(l,:)**2
                  !  +  0.5 * n*(n+1) * (4*pi) * z_Rad**2 * wz_Data(l,:)**2
            endif
         endif
      enddo

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

      do m=1,nm
         do n=0,m-1
            nmz_ToroidalEnergySpectrum_ut(n,m,:)  = ut_Vmiss
            nmz_ToroidalEnergySpectrum_ut(n,-m,:) = ut_Vmiss
         enddo
      enddo

    end function nmz_ToroidalEnergySpectrum_ut

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

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

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

      nz_Work = 0.0D0
      wz_DATA = wz_wt(wt_ut(ut_TORPOT))
      do l=1,nc
         nm_ary = nm_l(l)
         n = nm_ary(1) ; m = abs(nm_ary(2))
         if ( m == 0 ) then
            factor = 1.0D0
         else
            factor = 2.0D0
         endif
         if ( n .ge. 0 ) then
            nz_Work(n,:) = nz_Work(n,:) &
              + factor * 0.5 * n*(n+1)* (4*pi) * z_Rad**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_ut,(nm+1)*(km+1), &
                         MPI_REAL8,MPI_SUM,MPI_COMM_H,IERR)

    end function nz_ToroidalEnergySpectrum_ut

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

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

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

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

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

      integer :: ierr

      wz_Data1 = wz_wt(wt_ut(ut_POLPOT))
      wz_Data2 = wz_uz(uz_Rad*uz_ut(ut_DRad_ut(ut_POLPOT))) &    ! d(rφ)/dr
               + wz_wt(wt_ut(ut_POLPOT))                         ! = rdφ/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 )
                  !  +  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 )
                  !  +  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 )
                  ! +  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_ut, &
                         (nm+1)*(2*nm+1)*(km+1), &
                         MPI_REAL8,MPI_SUM,MPI_COMM_H,IERR)

      do m=1,nm
         do n=0,m-1
            nmz_PoloidalEnergySpectrum_ut(n,m,:)  = ut_Vmiss
            nmz_PoloidalEnergySpectrum_ut(n,-m,:) = ut_Vmiss
         enddo
      enddo

    end function nmz_PoloidalEnergySpectrum_ut

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

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

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

      nz_Work = 0.0D0
      wz_Data1 = wz_wt(wt_ut(ut_POLPOT))
      wz_Data2 = wz_uz(uz_Rad*uz_ut(ut_DRad_ut(ut_POLPOT))) &    ! d(rφ)/dr
               + wz_wt(wt_ut(ut_POLPOT))                         ! = rdφ/dr+φ

      do l=1,nc
         nm_ary = nm_l(l)
         n = nm_ary(1) ; m = abs(nm_ary(2))
         if ( m == 0 ) then
            factor = 1.0D0
         else
            factor = 2.0D0
         endif
         if ( n .ge. 0 ) then
            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_ut,(nm+1)*(km+1), &
                         MPI_REAL8,MPI_SUM,MPI_COMM_H,IERR)

    end function nz_PoloidalEnergySpectrum_ut

    !--------------- 境界値問題 -----------------

    subroutine ut_BoundariesTau(ut,values,cond)
      !
      ! スペクトルデータにディリクレ・ノイマン境界条件を適用する
      ! Chebyshev 空間での境界条件適用(タウ法)
      !
      ! チェビシェフ空間において境界条件を満たすべく高次の係数を
      ! 定める方法をとっている(タウ法).
      !
      real(DP), dimension(nmc,0:lm),intent(inout)      :: ut
      !(inout) 境界条件を適用するデータ. 修正された値を返す.

      real(DP), dimension(nmc,2), intent(in), optional :: values
      !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える.
      !     省略時は値/勾配 0 となる.

      character(len=2), intent(in), optional         :: cond
      !(in) 境界条件. 省略時は 'DD'
      !        DD : 両端ディリクレ条件
      !        DN : 上端ディリクレ, 下端ノイマン条件
      !        ND : 上端ノイマン, 下端ディリクレ条件
      !        NN : 両端ノイマン条件

      if (.not. present(cond)) then
         if (present(values)) then
            call at_BoundariesTau_DD(ut,values)
         else
            call at_BoundariesTau_DD(ut)
         endif
         return
      endif

      select case(cond)
      case ('NN')
         if (present(values)) then
            call at_BoundariesTau_NN(ut,values)
         else
            call at_BoundariesTau_NN(ut)
         endif
      case ('DN')
         if (present(values)) then
            call at_BoundariesTau_DN(ut,values)
         else
            call at_BoundariesTau_DN(ut)
         endif
      case ('ND')
         if (present(values)) then
            call at_BoundariesTau_ND(ut,values)
         else
            call at_BoundariesTau_ND(ut)
         endif
      case ('DD')
         if (present(values)) then
            call at_BoundariesTau_DD(ut,values)
         else
            call at_BoundariesTau_DD(ut)
         endif
      case default
         call MessageNotify('E','ut_BoundariesTau','B.C. not supported')
      end select

    end subroutine ut_BoundariesTau

    subroutine ut_BoundariesGrid(ut,values,cond)
      !
      ! スペクトルデータにディリクレ・ノイマン境界条件を適用する
      ! 実空間での境界条件適用
      !
      ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
      ! 条件を課している(選点法). このルーチンを用いるためには
      ! ut_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を
      ! 等しくしておく必要がある.
      !
      real(DP), dimension(nmc,0:lm),intent(inout)      :: ut
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

      real(DP), dimension(nmc,2), intent(in), optional :: values
              !(in) 境界での 値/勾配 分布を水平スペクトル変換したものを与える.
              !    省略時は値/勾配 0 となる.

      character(len=2), intent(in), optional         :: cond
              !(in) 境界条件. 省略時は 'DD'
              !        DD : 両端ディリクレ条件
              !        DN : 上端ディリクレ, 下端ノイマン条件
              !        ND : 上端ノイマン, 下端ディリクレ条件
              !        NN : 両端ノイマン条件

      if (.not. present(cond)) then
         if (present(values)) then
            call at_boundariesGrid_DD(ut,values)
         else
            call at_boundariesGrid_DD(ut)
         endif
         return
      endif

      select case(cond)
      case ('NN')
         if (present(values)) then
            call at_BoundariesGrid_NN(ut,values)
         else
            call at_BoundariesGrid_NN(ut)
         endif
      case ('DN')
         if (present(values)) then
            call at_BoundariesGrid_DN(ut,values)
         else
            call at_BoundariesGrid_DN(ut)
         endif
      case ('ND')
         if (present(values)) then
            call at_BoundariesGrid_ND(ut,values)
         else
            call at_BoundariesGrid_ND(ut)
         endif
      case ('DD')
         if (present(values)) then
            call at_BoundariesGrid_DD(ut,values)
         else
            call at_BoundariesGrid_DD(ut)
         endif
      case default
         call MessageNotify('E','ut_BoundariesGrid','B.C. not supported')
      end select

    end subroutine ut_BoundariesGrid

    subroutine ut_TorBoundariesTau(ut_TORPOT,values,cond,new)
      !
      ! 速度トロイダルポテンシャルに対して境界条件を適用する.
      ! Chebyshev 空間での境界条件適用.
      !
      ! 速度トロイダルポテンシャルΨに対して与えられる境界条件は
      !
      !   * 粘着条件 : Ψ = Ψb(lon,lat). Ψb は境界球面での速度分布.
      !                                   default は 0(静止状態).
      !
      !   * 応力なし条件 : ∂(Ψ/r)/∂r = 0.
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(DP), dimension(nmc,0:lm),intent(inout)      :: ut_TORPOT
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

      real(DP), dimension(nmc,2), intent(in), optional :: values
              !(in) 両端境界でのトロイダルポテンシャル
              !     粘着条件の時のみ有効

      character(len=2), intent(in), optional  :: cond
              !(in) 境界条件スイッチ. 省略時は 'RR'
              !     RR    : 両端粘着条件
              !     RF    : 上端粘着, 下端応力なし条件
              !     FR    : 上端応力なし, 下端粘着条件
              !     FF    : 両端応力なし条件

      logical, intent(IN), optional :: new
              !(in) true だと境界条件計算用行列を強制的に新たに作る.
              !     default は false.

      real(DP), dimension(:,:), allocatable  :: alu
      integer, dimension(:), allocatable    :: kp
      real(DP), dimension(0:lm)              :: t_data
      real(DP), dimension(0:km)              :: z_data
      logical                               :: rigid1=.true.   ! 境界条件
      logical                               :: rigid2=.true.   ! 境界条件

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer :: l
      save    :: alu, kp, first

      if (.not. present(cond)) then
         rigid1=.TRUE. ; rigid2=.TRUE.
      else
         select case (cond)
         case ('RR')
            rigid1 = .TRUE.  ; rigid2 = .TRUE.
         case ('RF')
            rigid1 = .TRUE.  ; rigid2 = .FALSE.
         case ('FR')
            rigid1 = .FALSE. ; rigid2 = .TRUE.
         case ('FF')
            rigid1 = .FALSE. ; rigid2 = .FALSE.
         case default
            call MessageNotify('E','ut_TorBoundariesTau','B.C. not supported')
         end select
      endif

      if (.not. present(new)) then
         new_matrix=.false.
      else
         new_matrix=new
      endif

      if ( first .OR. new_matrix ) then
         first = .false.

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         allocate(alu(0:lm,0:lm),kp(0:lm))

         do l=0,lm
            t_data = 0.0D0 ; t_data(l) = 1.0D0
            alu(:,l) = t_data

            ! 力学的条件粘着条件
            if ( rigid1 ) then
               z_data = z_t(t_data)
            else
               z_data = z_t(t_dr_t(t_z(z_t(t_data)/z_rad)))
            endif
            alu(lm-1,l) = z_data(0)       ! 境界 k=0 での条件式代入

            if ( rigid2 ) then
               z_data = z_t(t_data)
            else
               z_data = z_t(t_dr_t(t_z(z_t(t_data)/z_rad)))
            endif
            alu(lm,l)   = z_data(km)      ! 境界 k=km での条件式代入
         enddo

         call ludecomp(alu,kp)

         if ( rigid1 .AND. present(values) ) then
            call MessageNotify('M','ut_TorBoundariesTau',&
                 'Toroidal potential at k=0 was given by the optional variable.')
         else if ( rigid1 .AND. (.NOT.present(values)) ) then
            call MessageNotify('M','ut_TorBoundariesTau',&
                 'Toroidal potential at k=0 was set to zero.')
         else if ( (.NOT. rigid1) .AND. present(values) ) then
            call MessageNotify('W','ut_TorBoundariesTau',&
                 'Boundary value k=0 cannot be set under stress-free condition.')
         endif

         if ( rigid2 .AND. present(values) ) then
            call MessageNotify('M','ut_TorBoundariesTau',&
                 'Toroidal potential at k=0 was given by the optional variable.')
         else if ( rigid2 .AND. (.NOT.present(values)) ) then
            call MessageNotify('M','ut_TorBoundariesTau',&
                 'Toroidal potential at k=0 was set to zero.')
         else if ( (.NOT. rigid2) .AND. present(values) ) then
            call MessageNotify('W','ut_TorBoundariesTau',&
                 'Boundary value k=0 cannot be set under stress-free condition.')
         endif

         call MessageNotify('M','ut_TorBoundariesTau',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      if ( rigid1 .AND. present(values) ) then
         ut_torpot(:,lm-1) = values(:,1)
      else
         ut_torpot(:,lm-1) = 0.0D0
      endif
      if ( rigid2 .AND. present(values) ) then
         ut_torpot(:,lm)   = values(:,2)
      else
         ut_torpot(:,lm) = 0.0D0
      endif

      ut_torpot = lusolve(alu,kp,ut_TORPOT)

    end subroutine ut_TorBoundariesTau

    subroutine ut_TorBoundariesGrid(ut_TORPOT,values,cond,new)
      !
      ! 速度トロイダルポテンシャルに対して境界条件を適用する.
      ! 実空間での境界条件適用
      !
      ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
      ! 条件を課している(選点法). このルーチンを用いるためには
      ! ut_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を
      ! 等しくしておく必要がある.
      !
      ! 速度トロイダルポテンシャルΨに対して与えられる境界条件は
      !
      !   * 粘着条件 : Ψ = Ψb(lon,lat). Ψb は境界球面での速度分布.
      !                                   default は 0 (静止状態).
      !
      !   * 応力なし条件 : ∂(Ψ/r)/∂r = 0.
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(DP), dimension(nmc,0:lm),intent(inout)      :: ut_TORPOT
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

      real(DP), dimension(nmc,2), intent(in), optional :: values
              !(in) 両端境界でのトロイダルポテンシャル
              !     粘着条件の時のみ有効

      character(len=2), intent(in), optional  :: cond
              !(in) 境界条件スイッチ. 省略時は 'RR'
              !     RR    : 両端粘着条件
              !     RF    : 上端粘着, 下端応力なし条件
              !     FR    : 上端応力なし, 下端粘着条件
              !     FF    : 両端応力なし条件

      logical, intent(IN), optional :: new
              !(in) true だと境界条件計算用行列を強制的に新たに作る.
              !     default は false.

      real(DP), dimension(nmc,0:km)           :: uz_TORPOT
      real(DP), dimension(:,:), allocatable  :: alu
      integer, dimension(:), allocatable    :: kp
      real(DP), dimension(0:lm)              :: t_data
      real(DP), dimension(0:km)              :: z_data
      logical                               :: rigid1=.true.   ! 境界条件
      logical                               :: rigid2=.true.   ! 境界条件

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer :: l
      save    :: alu, kp, first

      if (.not. present(cond)) then
         rigid1=.TRUE. ; rigid2=.TRUE.
      else
         select case (cond)
         case ('RR')
            rigid1 = .TRUE.  ; rigid2 = .TRUE.
         case ('RF')
            rigid1 = .TRUE.  ; rigid2 = .FALSE.
         case ('FR')
            rigid1 = .FALSE. ; rigid2 = .TRUE.
         case ('FF')
            rigid1 = .FALSE. ; rigid2 = .FALSE.
         case default
            call MessageNotify('E','ut_TorBoundariesGrid','B.C. not supported')
         end select
      endif

      if (.not. present(new)) then
         new_matrix=.false.
      else
         new_matrix=new
      endif

      if ( first .OR. new_matrix ) then
         first = .false.

         if ( lm /= km ) then
            call MessageNotify('E','TorBoundariesGrid', &
             'Chebyshev truncation and number of grid points should be same.')
         endif

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         allocate(alu(0:km,0:lm),kp(0:lm))

         do l=0,lm
            t_data = 0.0D0 ; t_data(l)=1.0D0
            alu(:,l) = z_t(t_data)

            if ( rigid1 ) then
               z_data = z_t(t_data)
            else
               z_data = z_t(t_dr_t(t_z(z_t(t_data)/z_rad)))
            endif
            alu(0,l) = z_data(0)        ! 境界 k=0 での条件式代入

            if ( rigid2 ) then
               z_data = z_t(t_data)
            else
               z_data = z_t(t_dr_t(t_z(z_t(t_data)/z_rad)))
            endif
            alu(km,l)   = z_data(km)    ! 境界 k=km での条件式代入
         end do

         call ludecomp(alu,kp)

         if ( rigid1 .AND. present(values) ) then
            call MessageNotify('M','ut_TorBoundariesGrid',&
                 'Toroidal potential at k=0 was given by the optional variable.')
         else if ( rigid1 .AND. (.NOT.present(values)) ) then
            call MessageNotify('M','ut_TorBoundariesGrid',&
                 'Toroidal potential at k=0 was set to zero.')
         else if ( (.NOT. rigid1) .AND. present(values) ) then
            call MessageNotify('W','ut_TorBoundariesGrid',&
                 'Boundary value at k=0 cannot be set under stress-free condition.')
         endif

         if ( rigid2 .AND. present(values) ) then
            call MessageNotify('M','ut_TorBoundariesGrid',&
                 'Toroidal potential at k=km was given by the optional variable.')
         else if ( rigid2 .AND. (.NOT.present(values)) ) then
            call MessageNotify('M','ut_TorBoundariesGrid',&
                 'Toroidal potential at k=km was set to zero.')
         else if ( (.NOT. rigid2) .AND. present(values) ) then
            call MessageNotify('W','ut_TorBoundariesGrid',&
                 'Boundary value at k=km cannot be set under stress-free condition.')
         endif

         call MessageNotify('M','ut_TorBoundariesGrid',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      uz_TorPot       = uz_ut(ut_TorPot)

      if ( rigid1 .AND. present(values) ) then
         uz_TorPot(:,0)  = values(:,1)
      else
         uz_TorPot(:,0)  = 0.0D0
      endif

      if ( rigid2 .AND. present(values) ) then
         uz_TorPot(:,km) = values(:,2)
      else
         uz_TorPot(:,km) = 0.0D0
      endif

      ut_torpot = lusolve(alu,kp,uz_TorPot)

    end subroutine ut_TorBoundariesGrid

    function uz_LaplaPol2Pol_uz(uz,cond,new)
      !
      ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する.
      !
      ! チェビシェフ格子点空間で境界条件を適用している.
      ! この関数を用いるためには wt_Initial にて設定する
      ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく
      ! しておく必要がある.
      !
      ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は
      !
      !   ▽^2Φ = f
      !     Φ = const. at boundaries.
      !     ∂Φ/∂r = 0 at boundaries           (粘着条件)
      !     or ∂^2Φ/∂r^2 = 0 at boundaries    (応力なし条件)
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(DP), dimension(nmc,0:km),intent(in)  :: uz
      !(in) 入力▽^2φ分布

      real(DP), dimension(nmc,0:km)             :: uz_LaplaPol2Pol_uz
      !(out) 出力ポロイダルポテンシャル分布

      character(len=2), intent(in), optional  :: cond
      !(in) 境界条件スイッチ. 省略時は 'RR'
      !     RR    : 両端粘着条件
      !     RF    : 上端粘着, 下端応力なし条件
      !     FR    : 上端応力なし, 下端粘着条件
      !     FF    : 両端応力なし条件

      logical, intent(IN), optional :: new
      !(in) true だと境界条件計算用行列を強制的に新たに作る.
      !     default は false.

      real(DP), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(DP), dimension(nmc,0:km)            :: uz_work
      real(DP), dimension(0:km,0:km)           :: gg
      real(DP), dimension(0:km,0:km)           :: gg_work
      logical                                 :: rigid1=.true.   ! 境界条件
      logical                                 :: rigid2=.true.   ! 境界条件

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer :: k,n
      save    :: alu, kp, first

      if (.not. present(cond)) then
         rigid1=.TRUE. ; rigid2=.TRUE.
      else
         select case (cond)
         case ('RR')
            rigid1 = .TRUE.  ; rigid2 = .TRUE.
         case ('RF')
            rigid1 = .TRUE.  ; rigid2 = .FALSE.
         case ('FR')
            rigid1 = .FALSE. ; rigid2 = .TRUE.
         case ('FF')
            rigid1 = .FALSE. ; rigid2 = .FALSE.
         case default
            call MessageNotify('E','uz_LaplaPol2Pol_uz','B.C. not supported')
         end select
      endif

      if (.not. present(new)) then
         new_matrix=.false.
      else
         new_matrix=new
      endif

      if ( first .OR. new_matrix ) then
         first = .false.

         if ( lm /= km ) then
            call MessageNotify('E','uz_LaplaPol2Pol_uz', &
                 'Chebyshev truncation and number of grid points should be same.')
         endif

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         allocate(alu(nmc,0:km,0:km),kp(nmc,0:km))

         do k=0,km
            uz_work = 0.0D0 ; uz_work(:,k) = 1.0D0

            ! 各水平波数に関して独立の式
            alu(:,:,k) = uz_ut(ut_lapla_ut(ut_uz(uz_work)))
         enddo

         ! 運動学的条件. 流線は境界で一定
         gg = 0.0D0
         do k=0,km
            gg(k,k)=1.0D0
         enddo
         do n=1,nmc
            alu(n,0,:)   = gg(:,0)
            alu(n,km,:)  = gg(:,km)
         enddo

         ! 力学的条件粘着条件
         if ( rigid1 ) then
            gg_work=az_at(at_dr_at(at_az(gg)))
         else
            gg_work=az_at(at_dr_at(at_dr_at(at_az(gg))))
         endif
         do n=1,nmc
            alu(n,1,:) = gg_work(:,0)
         enddo

         ! 力学的条件粘着条件
         if ( rigid2 ) then
            gg_work=az_at(at_dr_at(at_az(gg)))
         else
            gg_work=az_at(at_dr_at(at_dr_at(at_az(gg))))
         endif
         do n=1,nmc
            alu(n,km-1,:) = gg_work(:,km)
         enddo

         call ludecomp(alu,kp)

         call MessageNotify('M','uz_LaplaPol2Pol_uz',&
              'Matrix to apply  b.c. newly produced.')
      endif

      uz_work         = uz
      uz_work(:,1)    = 0.0D0               ! 力学的条件
      uz_work(:,km-1) = 0.0D0               ! 力学的条件
      uz_work(:,0)    = 0.0D0               ! 運動学的条件
      uz_work(:,km)   = 0.0D0               ! 運動学的条件

      uz_laplapol2pol_uz = lusolve(alu,kp,uz_work)

    end function uz_LaplaPol2Pol_uz

    function ut_LaplaPol2PolTau_ut(ut,cond,new)
      !
      ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する.
      ! チェビシェフタウ法で境界条件を適用している.
      !
      ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は
      !
      !    ▽^2Φ = f
      !      Φ = const. at boundaries.
      !      ∂Φ/∂r = 0 at boundaries          (粘着条件)
      !      or ∂^2Φ/∂r^2 = 0 at boundaries   (応力なし条件)
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(DP), dimension(nmc,0:lm),intent(in)  :: ut
              !(in) 入力▽^2φ分布

      real(DP), dimension(nmc,0:lm)             :: ut_LaplaPol2PolTau_ut
              !(out) 出力ポロイダルポテンシャル分布

      character(len=2), intent(in), optional  :: cond
              !(in) 境界条件スイッチ. 省略時は 'RR'
              !     RR    : 両端粘着条件
              !     RF    : 上端粘着, 下端応力なし条件
              !     FR    : 上端応力なし, 下端粘着条件
              !     FF    : 両端応力なし条件

      logical, intent(IN), optional :: new
              !(in) true だと境界条件計算用行列を強制的に新たに作る.
              !     default は false.

      real(DP), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(DP), dimension(nmc,0:km)  :: uz_work
      real(DP), dimension(nmc,0:lm)  :: ut_work
      logical                       :: rigid1=.true.   ! 境界条件
      logical                       :: rigid2=.true.   ! 境界条件

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer :: l
      save    :: alu, kp, first

      if (.not. present(cond)) then
         rigid1=.TRUE. ; rigid2=.TRUE.
      else
         select case (cond)
         case ('RR')
            rigid1 = .TRUE.  ; rigid2 = .TRUE.
         case ('RF')
            rigid1 = .TRUE.  ; rigid2 = .FALSE.
         case ('FR')
            rigid1 = .FALSE. ; rigid2 = .TRUE.
         case ('FF')
            rigid1 = .FALSE. ; rigid2 = .FALSE.
         case default
            call MessageNotify('E','ut_LaplaPol2PolTau_ut','B.C. not supported')
         end select
      endif

      if (.not. present(new)) then
         new_matrix=.false.
      else
         new_matrix=new
      endif

      if ( first .OR. new_matrix ) then
         first = .false.

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         allocate(alu(nmc,0:lm,0:lm),kp(nmc,0:lm))

         alu = 0.0D0
         do l=0,lm-2
            ut_work = 0.0D0 ; ut_work(:,l) = 1.0D0

            ! 各水平波数に関して独立の式
            alu(:,:,l) = ut_Lapla_ut(ut_work)
         enddo

         ! 運動学的条件. 流線は境界で一定
         do l=0,lm
            ut_work=0.0D0
            ut_work(:,l) = 1.0D0
            uz_work = uz_ut(ut_work)

            alu(:,lm-1,l)  = uz_work(:,0)
            alu(:,lm,l)    = uz_work(:,km)

            ! 力学的条件粘着条件
            if ( rigid1 ) then
               uz_work=uz_ut(ut_DRad_ut(ut_work))
            else
               uz_work=uz_ut(ut_DRad_ut(ut_DRad_ut(ut_work)))
            endif
            alu(:,lm-2,l) = uz_work(:,0)

            ! 力学的条件粘着条件
            if ( rigid2 ) then
               uz_work=uz_ut(ut_DRad_ut(ut_work))
            else
               uz_work=uz_ut(ut_DRad_ut(ut_DRad_ut(ut_work)))
            endif
            alu(:,lm-3,l) = uz_work(:,km)
         end do

         call ludecomp(alu,kp)

         call MessageNotify('M','ut_LaplaPol2PolTau_ut',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      ut_work         = ut
      ut_work(:,lm)   = 0.0D0               ! 運動学的条件
      ut_work(:,lm-1) = 0.0D0               ! 運動学的条件
      ut_work(:,lm-2) = 0.0D0               ! 力学的条件
      ut_work(:,lm-3) = 0.0D0               ! 力学的条件

      ut_LaplaPol2PolTau_ut = lusolve(alu,kp,ut_work)

    end function ut_LaplaPol2PolTau_ut

    function ut_LaplaPol2PolGrid_ut(ut,cond,new)
      !
      ! 速度ポロイダルポテンシャルΦを▽^2Φから計算する.
      ! チェビシェフ格子点空間で境界条件を適用している.
      !
      ! この関数を用いるためには ut_Initial にて設定する
      ! チェビシェフ切断波数(lm)と鉛直格子点数(km)を等しく
      ! しておく必要がある.
      !
      ! 速度ポロイダルポテンシャルΦを f = ▽^2Φから定める式は
      !
      !    ▽^2Φ = f
      !      Φ = const. at boundaries.
      !      ∂Φ/∂r = 0 at boundaries          (粘着条件)
      !      or ∂^2Φ/∂r^2 = 0 at boundaries   (応力なし条件)
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      ! 最終的にチェビシェフ係数の解が欲しい場合には, uz_LaplaPol2Pol_uz に
      ! 比べてチェビシェフ -- 格子点変換が 1 回分少なくて済む.
      !
      real(DP), dimension(nmc,0:lm),intent(in)  :: ut
              !(in) 入力▽^2φ分布

      real(DP), dimension(nmc,0:lm)             :: ut_LaplaPol2PolGrid_ut
              !(out) 出力ポロイダルポテンシャル分布

      character(len=2), intent(in), optional  :: cond
              !(in) 境界条件スイッチ. 省略時は 'RR'
              !     RR    : 両端粘着条件
              !     RF    : 上端粘着, 下端応力なし条件
              !     FR    : 上端応力なし, 下端粘着条件
              !     FF    : 両端応力なし条件

      logical, intent(IN), optional :: new
              !(in) true だと境界条件計算用行列を強制的に新たに作る.
              !     default は false.

      real(DP), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(DP), dimension(nmc,0:km)  :: uz_work
      real(DP), dimension(nmc,0:lm)  :: ut_work
      real(DP), dimension(0:lm,0:lm)           :: tt_I
      real(DP), dimension(0:lm,0:km)           :: tz_work
      logical                                 :: rigid1=.true.   ! 境界条件
      logical                                 :: rigid2=.true.   ! 境界条件

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer :: l,n
      save    :: alu, kp, first

      if (.not. present(cond)) then
         rigid1=.TRUE. ; rigid2=.TRUE.
      else
         select case (cond)
         case ('RR')
            rigid1 = .TRUE.  ; rigid2 = .TRUE.
         case ('RF')
            rigid1 = .TRUE.  ; rigid2 = .FALSE.
         case ('FR')
            rigid1 = .FALSE. ; rigid2 = .TRUE.
         case ('FF')
            rigid1 = .FALSE. ; rigid2 = .FALSE.
         case default
            call MessageNotify('E','ut_LaplaPol2PolGrid_ut','B.C. not supported')
         end select
      endif

      if (.not. present(new)) then
         new_matrix=.false.
      else
         new_matrix=new
      endif

      if ( first .OR. new_matrix ) then
         first = .false.

         if ( lm /= km ) then
            call MessageNotify('E','ut_LaplaPol2PolGrid_ut', &
             'Chebyshev truncation and number of grid points should be same.')
         endif

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         allocate(alu(nmc,0:km,0:lm),kp(nmc,0:lm))

         do l=0,lm
            ut_work = 0.0D0 ; ut_work(:,l) = 1.0D0

            ! 各水平波数に関して独立の式
            alu(:,:,l) = uz_ut(ut_Lapla_ut(ut_work))
         enddo

         ! 運動学的条件. 流線は境界で一定
         tt_I = 0.0D0
         do l=0,lm
            tt_I(l,l)=1.0D0
         enddo

         ! 非電気伝導体
         tz_work = az_at(tt_I)

         do n=1,nmc
            alu(n,0,:)  = tz_work(:,0)
            alu(n,km,:) = tz_work(:,km)
         enddo

         ! 力学的条件粘着条件
         if ( rigid1 ) then
            tz_work=az_at(at_Dr_at(tt_I))
         else
            tz_work=az_at(at_Dr_at(at_Dr_at(tt_I)))
         endif
         do n=1,nmc
            alu(n,1,:) = tz_work(:,0)
         enddo

         ! 力学的条件粘着条件
         if ( rigid2 ) then
            tz_work=az_at(at_Dr_at(tt_I))
         else
            tz_work=az_at(at_Dr_at(at_Dr_at(tt_I)))
         endif
         do n=1,nmc
            alu(n,km-1,:) = tz_work(:,km)
         enddo

         call ludecomp(alu,kp)

         call MessageNotify('M','ut_LaplaPol2PolGrid_ut',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      uz_work         = uz_ut(ut)
      uz_work(:,1)    = 0.0D0               ! 力学的条件
      uz_work(:,km-1) = 0.0D0               ! 力学的条件
      uz_work(:,0)    = 0.0D0               ! 運動学的条件
      uz_work(:,km)   = 0.0D0               ! 運動学的条件

      ut_LaplaPol2PolGrid_ut = lusolve(alu,kp,uz_work)

    end function ut_LaplaPol2PolGrid_ut

    subroutine ut_TormagBoundariesTau(ut_TOR,new)

      ! 磁場トロイダルポテンシャルに対して境界条件を適用する.
      ! Chebyshev 空間での境界条件適用
      !
      ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を
      ! とっている(タウ法). 現在のところ境界物質が非電気伝導体の場合のみ
      ! 対応している. その場合, 磁場トロイダルポテンシャルの境界条件は
      !
      ! 外側
      !    ut_psi = 0   at the outer boundary
      ! 内側
      !    ut_psi = 0       at the inner boundary
      !
      ! であるから ut_Boundaries で対応可能だが, 将来のため別途作成しておく.
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(DP), dimension(nmc,0:lm),intent(inout)   :: ut_TOR
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

      logical, intent(IN), optional :: new
              !(in) true だと境界条件計算用行列を強制的に新たに作る.
              !     default は false.

      real(DP), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(DP), dimension(:,:), allocatable    :: ut_I
      real(DP), dimension(:,:), allocatable    :: uz_PSI

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer  :: l
      save     :: alu, kp, first

      if (.not. present(new)) then
         new_matrix=.false.
      else
         new_matrix=new
      endif

      if ( first .OR. new_matrix ) then
         first = .false.

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         if ( allocated(ut_I) ) deallocate(ut_I)
         if ( allocated(uz_PSI) ) deallocate(uz_PSI)

         allocate(alu(nmc,0:lm,0:lm),kp(nmc,0:lm))
         allocate(ut_I(nmc,0:lm),uz_PSI(nmc,0:km))

         do l=0,lm
            ut_I = 0.0D0 ; ut_I(:,l)=1.0D0
            alu(:,:,l) = ut_I

            ! 非電気伝導体
            uz_PSI = uz_ut(ut_I)

            alu(:,lm-1,l) = uz_PSI(:,0)
            alu(:,lm,l)   = uz_PSI(:,km)
         enddo

         call ludecomp(alu,kp)

         deallocate(ut_I,uz_PSI)

         call MessageNotify('M','TormagBoundariesTau',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      ut_TOR(:,lm-1) = 0.0D0
      ut_TOR(:,lm)   = 0.0D0
      ut_TOR = lusolve(alu,kp,ut_TOR)

    end subroutine ut_TormagBoundariesTau

    subroutine ut_TormagBoundariesGrid(ut_TOR,new)
      !
      ! 磁場トロイダルポテンシャルに対して境界条件を適用する.
      ! 鉛直実空間での境界条件適用.
      !
      ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
      ! 条件を課している(選点法). このルーチンを用いるためには
      ! ut_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を
      ! 等しくしておく必要がある.
      !
      ! 現在のところ境界物質が非電気伝導体の場合のみ対応している.
      ! その場合, 磁場トロイダルポテンシャルの境界条件は
      !
      ! 外側
      !    ut_psi = 0   at the outer boundary
      ! 内側
      !    ut_psi = 0       at the inner boundary
      !
      ! であるので ut_Boundaries で対応可能だが, 将来のため別途作成しておく
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(DP), dimension(nmc,0:lm),intent(inout)   :: ut_TOR
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

      logical, intent(IN), optional :: new
              !(in) true だと境界条件計算用行列を強制的に新たに作る.
              !     default は false.

      real(DP), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(DP), dimension(:,:), allocatable    :: ut_I
      real(DP), dimension(:,:), allocatable    :: uz_PSI
      real(DP), dimension(nmc,0:km)  :: uz_TOR

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer  :: l
      save     :: alu, kp, first

      if (.not. present(new)) then
         new_matrix=.false.
      else
         new_matrix=new
      endif

      if ( first .OR. new_matrix ) then
         first = .false.

         if ( lm /= km ) then
            call MessageNotify('E','TorMagBoundariesGrid', &
             'Chebyshev truncation and number of grid points should be same.')
         endif

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         if ( allocated(ut_I) ) deallocate(ut_I)
         if ( allocated(uz_PSI) ) deallocate(uz_PSI)

         allocate(alu(nmc,0:km,0:lm),kp(nmc,0:lm))
         allocate(ut_I(nmc,0:lm),uz_PSI(nmc,0:km))

         do l=0,lm
            ut_I = 0.0D0 ; ut_I(:,l)=1.0D0
            uz_PSI = uz_ut(ut_I)
            alu(:,:,l) = uz_PSI

            ! 非電気伝導体
            alu(:,0,l)  = uz_PSI(:,0)
            alu(:,km,l) = uz_PSI(:,km)
         enddo

         call ludecomp(alu,kp)

         deallocate(ut_I,uz_PSI)

         call MessageNotify('M','TormagBoundariesGrid',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      uz_TOR       = uz_ut(ut_TOR)
      uz_TOR(:,0)  = 0.0D0
      uz_TOR(:,km) = 0.0D0
      ut_TOR = lusolve(alu,kp,uz_TOR)

    end subroutine ut_TormagBoundariesGrid

    subroutine ut_PolmagBoundariesTau(ut_POL,new)
      !
      ! 磁場ポロイダルポテンシャルに対して境界条件を適用する.
      ! Chebyshev 空間での境界条件適用
      !
      ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を
      ! とっている(タウ法). 現在のところ境界物質が非電気伝導体の場合のみ
      ! 対応している. その場合, 磁場ポロイダルポテンシャルの各水平スペクトル
      ! 成分 h にたいして境界条件が与えられ,
      !
      !  * 外側境界 : dh/dr + (n+1)h/r = 0
      !  * 内側境界 : dh/dr - nh/r = 0
      !
      ! である. ここで n は h の水平全波数である.
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(DP), dimension(nmc,0:lm),intent(inout)   :: ut_POL
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

      logical, intent(IN), optional :: new
              !(in) true だと境界条件計算用行列を強制的に新たに作る.
              !     default は false.

      real(DP), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(DP), dimension(:,:), allocatable    :: ut_I
      real(DP), dimension(:,:), allocatable    :: uz_PSI
      real(DP), dimension(:,:), allocatable    :: uz_DPSIDR
      real(DP), dimension(:), allocatable      :: u_n

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer  :: l, n, nn(2)
      save     :: alu, kp, first

      if (.not. present(new)) then
         new_matrix=.false.
      else
         new_matrix=new
      endif

      if ( first .OR. new_matrix ) then
         first = .false.

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         if ( allocated(ut_I) ) deallocate(ut_I)
         if ( allocated(uz_PSI) ) deallocate(uz_PSI)
         if ( allocated(uz_DPSIDR) ) deallocate(uz_DPSIDR)
         if ( allocated(u_n) ) deallocate(u_n)

         allocate(alu(nmc,0:lm,0:lm),kp(nmc,0:lm))
         allocate(ut_I(nmc,0:lm),u_n(nmc))
         allocate(uz_PSI(nmc,0:km),uz_DPSIDR(nmc,0:km))

         do n=1,nmc
            nn=nm_l(n)
            u_n(n) = nn(1)
         enddo

         do l=0,lm
            ut_I = 0.0D0 ; ut_I(:,l) = 1.0D0
            alu(:,:,l) = ut_I

            uz_PSI = uz_ut(ut_I)
            uz_DPSIDR = uz_ut(ut_DRad_ut(ut_I))
            alu(:,lm-1,l) = uz_DPSIDR(:,0) + (u_n+1) * uz_PSI(:,0)/z_RAD(0)
            alu(:,lm,l)   = uz_DPSIDR(:,km) - u_n * uz_PSI(:,km)/z_RAD(km)
         enddo

         call ludecomp(alu,kp)

         deallocate(ut_I,uz_PSI,uz_DPSIDR,u_n)

         call MessageNotify('M','PolmagBoundariesTau',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      ut_POL(:,lm-1) = 0.0D0
      ut_POL(:,lm)   = 0.0D0
      ut_POL = lusolve(alu,kp,ut_POL)

    end subroutine ut_PolmagBoundariesTau

    subroutine ut_PolmagBoundariesGrid(ut_POL,new)
      !
      ! 磁場ポロイダルポテンシャルに対して境界条件を適用する.
      ! 鉛直実空間での境界条件適用.
      !
      ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
      ! 条件を課している(選点法). このルーチンを用いるためには
      ! ut_Initial にて設定するチェビシェフ切断波数(lm)と鉛直格子点数(km)を
      ! 等しくしておく必要がある.
      !
      ! 現在のところ境界物質が非電気伝導体の場合のみ対応している.
      ! その場合, 磁場ポロイダルポテンシャルの各水平スペクトル成分 h に
      ! たいして境界条件が与えられ,
      !
      !  * 外側境界 : dh/dr + (n+1)h/r = 0
      !  * 内側境界 : dh/dr - nh/r = 0
      !
      ! である. ここで n は h の水平全波数である.
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(DP), dimension(nmc,0:lm),intent(inout)   :: ut_POL
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

      logical, intent(IN), optional :: new
              !(in) true だと境界条件計算用行列を強制的に新たに作る.
              !     default は false.

      real(DP), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(DP), dimension(:,:), allocatable    :: ut_I
      real(DP), dimension(:,:), allocatable    :: uz_PSI
      real(DP), dimension(:,:), allocatable    :: uz_DPSIDR
      real(DP), dimension(:), allocatable      :: u_n

      real(DP), dimension(nmc,0:km)  :: uz_POL

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer  :: l, n, nn(2)
      save     :: alu, kp, first

      if (.not. present(new)) then
         new_matrix=.false.
      else
         new_matrix=new
      endif

      if ( first .OR. new_matrix ) then
         first = .false.

         if ( lm /= km ) then
            call MessageNotify('E','PolMagBoundariesGrid', &
             'Chebyshev truncation and number of grid points should be same.')
         endif

         if ( allocated(alu) )  deallocate(alu)
         if ( allocated(kp) )   deallocate(kp)
         if ( allocated(ut_I) ) deallocate(ut_I)
         if ( allocated(uz_PSI) ) deallocate(uz_PSI)
         if ( allocated(uz_DPSIDR) ) deallocate(uz_DPSIDR)
         if ( allocated(u_n) ) deallocate(u_n)

         allocate(alu(nmc,0:lm,0:lm),kp(nmc,0:lm))
         allocate(ut_I(nmc,0:lm),u_n(nmc))
         allocate(uz_PSI(nmc,0:km),uz_DPSIDR(nmc,0:km))

         do n=1,nmc
            nn=nm_l(n)
            u_n(n) = nn(1)
         enddo

         do l=0,lm
            ut_I = 0.0D0 ; ut_I(:,l) = 1.0D0
            uz_PSI = uz_ut(ut_I)
            alu(:,:,l) = uz_PSI

            uz_DPSIDR = uz_ut(ut_DRad_ut(ut_I))
            alu(:,0,l)  = uz_DPSIDR(:,0) + (u_n+1) * uz_PSI(:,0)/z_RAD(0)
            alu(:,km,l) = uz_DPSIDR(:,km) - u_n * uz_PSI(:,km)/z_RAD(km)
         enddo

         call ludecomp(alu,kp)

         deallocate(ut_I,uz_PSI,uz_DPSIDR,u_n)

         call MessageNotify('M','PolmagBoundariesGrid',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      uz_POL       = uz_ut(ut_POL)
      uz_POL(:,0)  = 0.0D0
      uz_POL(:,km) = 0.0D0
      ut_POL = lusolve(alu,kp,uz_POL)

    end subroutine ut_PolmagBoundariesGrid


  !--------------- 終了処理 -----------------
    subroutine ut_mpi_Finalize
      !
      ! モジュールの終了処理(割り付け配列の解放)をおこなう.
      !
      ! このサブルーチンを単独で用いるのでなく,
      ! 上位サブルーチン wa_mpi_Finalize を使用すること.
      !
      deallocate(h_Rad)
      deallocate(h_Rad_Weight)

      deallocate(xvh_Lon)
      deallocate(xvh_Lat)
      deallocate(xvh_Rad)

      deallocate(wh_Rad)
      deallocate(uz_Rad)

      call ua_mpi_finalize
      call at_Finalize

      call MessageNotify('M','ut_mpi_finalize',&
           'ut_mpi_module (2021/11/15) is finalized')

    end subroutine ut_mpi_Finalize

end module ut_mpi_module
