!--
!----------------------------------------------------------------------
! Copyright(c) 2009-2016 SPMDODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  tee_mpi_module
!
!    spml/tee_module モジュールは平行平板間での 3 次元流体運動を
!    スペクトル法によって数値計算するための Fortran90 関数を提供する
!    ものである.
!
!    水平方向にフーリエ数変換および上下の境界壁を扱うための
!    チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな
!    関数を提供する.
!
!    内部で ee_module, at_module を用いている. 最下部ではフーリエ
!    およびチェビシェフ変換のエンジンとして ISPACK の Fortran77
!    サブルーチンを用いている.
!
!
!履歴  2009/12/19  竹広真一  新規作成
!      2010/03/10  佐々木洋平  threadprivate 削除(コンパイラ依存)
!      2010/03/27  竹広真一  goto 文を削除
!      2011/12/05  竹広真一  tee_Div_zyx_zyx_zyx を public 変数に
!      2012/07/08  竹広真一  tee_LaplaPol2PolTau_tee を追加
!      2012/08/24  竹広真一  tee_LaplaPol2Pol_tee public に追加
!      2013/08/20  竹広真一  gnu fortran 対応
!      2016/09/22  竹広真一  tee_module を mpi 並列化
!
!凡例
!      データ種類と index
!        x : 水平(X)   y : 水平(Y)   v : 水平(Y,MPI 分割)  z : 鉛直
!        e : フーリエ変換スペクトル
!        f : フーリエ変換スペクトル(MPI 分割
!        l : フーリエ関数スペクトル(X 方向波数)
!        m : フーリエ関数スペクトル(Y 方向波数)
!        t : チェビシェフ関数スペクトル
!        a : 任意の次元
!
!        zvx : 3 次元格子点データ(MPI 分割)
!        vx  : 水平 2 次元格子点データ(MPI 分割)
!        zv  : 鉛直 2 次元格子点データ(MPI 分割)
!        zx  : 鉛直 2 次元格子点データ
!
!        zef : 水平スペクトル鉛直格子点データ(MPI 分割)
!        tef : スペクトルデータ鉛直チェビシェフデータ(MPI 分割)
!
!++
module tee_mpi_module
  !
  != tee_mpi_module
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  !    spml/tee_mpi_module モジュールは平行平板間での 3 次元流体運動を
  !    スペクトル法によって数値計算するための Fortran90 関数を提供する
  !    ものである.
  !
  !    水平方向にフーリエ数変換および上下の境界壁を扱うための
  !    チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな
  !    関数を提供する.
  !
  !    内部で ee_mpi_module, at_module を用いている. 最下部ではフーリエ
  !    およびチェビシェフ変換のエンジンとして ISPACK の Fortran77
  !    サブルーチンを用いている.
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (tef_, zvx_, zef_, ef_, vx_, x_, v_, z_, a_) は,
  !   返す値の形を示している.
  !   tef_  :: スペクトルデータ(2 重フーリエ・チェビシェフ変換)
  !   zvx_ :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・)
  !   zef_ :: 水平スペクトル, 鉛直格子点データ
  !   e2a_ :: 1 次元化した水平スペクトル, 任意座標データ
  !   aef_ :: 任意座標データ, 水平スペクトルデータ
  !
  ! * 関数名の間の文字列(Dx, Dy, Dz, Lapla,..)
  !   は, その関数の作用を表している.
  !
  ! * 関数名の最後 (tef_, zvx_, zef_, ef_, vx_, x_, v_, z_, a_) は, 入力変数の
  !   形がスペクトルデータおよび格子点データであることを示している.
  !   _tef :: スペクトルデータ(2 重フーリエ・チェビシェフ変換)
  !   _zvx :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・)
  !   _zef :: 水平スペクトル, 鉛直格子点データ
  !   _e2a :: 1 次元化した水平スペクトル, 任意座標データ
  !   _aef :: 任意座標データ, 水平スペクトルデータ
  !
  !=== 各データの種類の説明
  !
  ! * zvx : 3 次元格子点データ(鉛直, 水平 2 次元, MPI 分割)
  !   * 変数の種類と次元は real(8), dimension(0:km,jc(ip),0:im-1).
  !   * im, km はそれぞれ水平 X, 鉛直 Z 座標の格子点数である. 
  !     サブルーチン tee_mpi_Initial にてあらかじめ設定しておく.
  !   * jc(ip) は tee_mpi_initial にて全 Y 格子点数 jm を設定後に使用できる.
  !
  ! * tef : スペクトルデータ
  !   * 変数の種類と次元は real(8), dimension(0:nm,-mm:mm,2*lc(ip)).
  !   * mm は Y 方向最大波数, nm はチェビシェフ多項式の最大次数
  !     であり, サブルーチン tee_Initial にてあらかじめ設定しておく.
  !   * lc(ip) は tee_mpi_initial にて全 Y 方向最大波数 lm を設定後に使用できる.
  !
  ! * zef : 水平スペクトル, 鉛直格子点データ.
  !   * 変数の種類と次元は real(8), dimension(0:km,-mm:mm,2*lc(ip)).
  !
  ! * e2a : 1 次元化した水平スペクトル, 任意座標データ
  !   * 変数の種類と次元は real(8), dimension((2*mm+1)*2*lc(ip),:)
  !
  ! * aef :   任意座標データ, 水平スペクトルデータ
  !   * 変数の種類と次元は real(8), dimension(:,-mm:mm,2*lc(ip)).
  !
  ! * tef_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * zvx_ で始まる関数が返す値は 3 次元格子点データに同じ.
  !
  ! * zef_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! tee_mpi_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  !
  !==== 座標変数
  !
  ! x_X, v_Y, z_Z                :: 格子点座標(水平 X,Y, 鉛直 Z 座標)を
  !                                 格納した1 次元配列
  ! x_X_Weight, v_X_Weight, z_Z_Weight :: 重み座標を格納した 1 次元配列
  ! zvx_X, zvx_Y, zyx_Z          :: 格子点データの水平鉛直座標(X,Y,Z)
  !                                 (格子点データ型 3 次元配列)
  ! vx_X, vx_Y                   :: 格子点データの水平座標(X,Y)
  ! zv_Z, zv_Y                   :: 格子点データの鉛直水平座標(Z,Y)
  ! zx_Z, zx_X                   :: 格子点データの鉛直水平座標(Z,X)
  !                                 (格子点データ型 2 次元配列)
  !
  !==== 基本変換
  !
  ! zvx_tee, tee_zvx :: スペクトルデータと 3 次元格子データの間の変換
  !                     (2 重フーリエ, チェビシェフ変換)
  ! zvx_zee, zee_zvx :: 3 次元格子データと水平スペクトル・鉛直格子データとの間
  !                     の変換 (2 重フーリエ変換)
  ! zef_tee, tee_zef :: スペクトルデータと水平スペクトル・鉛直格子データとの間
  !                     の変換 (チェビシェフ変換)
  ! ee_vx, vx_ee     :: スペクトルデータと 2 次元水平格子データの間の変換
  !                     (2 重フーリエ変換)
  ! az_at, at_az     :: 同時に複数個行う (チェビシェフ変換)格子データと
  !                     チェビシェフデータの間の変換を
  ! e2a_aef, aef_e2a :: 水平スペクトル軸を 1 次元化し転置, 転置 ２次元化する.
  !
  !==== 微分
  !
  ! tef_Dx_tef          :: スペクトルデータに動径微分∂/∂x を作用させる
  ! tef_Dy_tef          :: スペクトルデータに動径微分∂/∂y を作用させる
  ! tef_Dz_tef          :: スペクトルデータに動径微分∂/∂z を作用させる
  !
  ! tef_Lapla_tef       :: スペクトルデータにラプラシアンを作用させる
  ! tef_LaplaH_tef      :: スペクトルデータに水平ラプラシアンを作用させる
  ! tef_LaplaHInv_tef   :: スペクトルデータに逆水平ラプラシアンを作用させる
  !
  ! tef_Div_zvx_zvx_zvx :: ベクトル成分である 3 つの格子データに
  !                        発散を作用させる
  !
  !==== トロイダルポロイダル計算用微分
  !
  ! tef_ZRot_zvx_zvx         :: ベクトル v の渦度と動径ベクトル r の内積
  !                             z・(▽×v) を計算する
  ! tef_ZRotRot_zvx_zvx_zvx  :: ベクトルの v の z・(▽×▽×v) を計算する
  ! tef_Potential2Vector     :: トロイダルポロイダルポテンシャルから
  !                             ベクトル場を計算する
  ! tef_Potential2Rotation   :: トロイダルポロイダルポテンシャルで表される
  !                             非発散ベクトル場の回転の各成分を計算する
  !
  !==== ポロイダル/トロイダルモデル用スペクトル解析
  !
  ! zef_ToroidalEnergySpectrum_tef, zk_ToroidalEnergySpectrum_tef   ::
  !     トロイダルポテンシャルからエネルギーのフーリエ各成分を計算する
  ! zef_PoloidalEnergySpectrum_tef, zk_PoloidalEnergySpectrum_tef   ::
  !     ポロイダルポテンシャルからエネルギーのフーリエ各成分を計算する
  !
  !==== 境界値問題
  !
  ! tef_BoundariesTau, tef_BoundariesGrid, tef_Boundaries                ::
  !     ディリクレ, ノイマン境界条件を適用する(タウ法, 選点法)
  !
  ! tef_TorBoundariesTau, tef_TorBoundariesGrid, tef_TorBoundaries       ::
  !     速度トロイダルポテンシャルの境界条件を適用する(タウ法,選点法)
  !
  ! tef_LaplaPol2PolTau_tef, zee_LaplaPol2Pol_zee, tef_LaplaPol2PolGrid_tef ::
  !     速度ポロイダルポテンシャルΦを▽^2Φから求める
  !     (入出力がそれぞれチェビシェフ格子点,チェビシェフ係数)
  !
  ! tef_TorMagBoundariesTau, tef_TorMagBoundariesGrid, tef_TorMagBoundaries ::
  !     磁場トロイダルポテンシャルの境界条件を適用する(タウ法, 選点法)
  !
  ! tef_PolMagBoundariesTau, tef_PolMagBoundariesGrid, tef_PolMagBoundaries ::
  !     磁場トロイダルポテンシャル境界の境界条件を適用する(タウ法, 選点法)
  !
  !==== 積分・平均(3 次元データ)
  !
  ! IntZYX_zvx, AvrZYX_zvx     :: 3 次元格子点データの全領域積分および平均
  ! z_IntYX_zvx, z_AvrYX_zvx   :: 3 次元格子点データの水平積分および平均
  ! v_IntZX_zvx, v_AvrZX_zvx   :: 3 次元格子点データのZX積分および平均
  ! z_IntZY_zvx, z_AvrZY_zvx   :: 3 次元格子点データのZY積分および平均
  ! zv_IntX_zvx, zv_AvrX_zvx   :: 3 次元格子点データの水平X方向積分および平均
  ! zx_IntY_zvx, zx_AvrY_zvx   :: 3 次元格子点データの水平Y方向積分および平均
  ! zx_IntZ_zvx, zx_AvrZ_zvx   :: 3 次元格子点データの鉛直方向積分および平均
  !
  !==== 積分・平均(2 次元データ)
  !
  ! IntYX_vx, AvrYX_vx :: 2 次元格子点データの水平積分および平均
  ! IntZX_zx, AvrZX_zx :: 2 次元(ZX)格子点データのZX積分および平均
  ! IntZY_zv, AvrZY_zv :: 2 次元(ZY)格子点データのZY積分および平均
  ! y_IntX_vx, y_AvrX_vx   :: 水平 2 次元格子点データのX方向積分および平均
  ! x_IntY_vx, x_AvrY_vx   :: 水平2 次元格子点データのY方向積分および平均
  ! z_IntX_zx, z_AvrX_zx   :: 2 次元(ZX)格子点データのX方向積分および平均
  ! x_IntZ_zx, x_AvrZ_zx   :: 2 次元(ZX)格子点データのZ方向積分および平均
  ! z_IntY_zv, z_AvrY_zv   :: 2 次元(ZY)格子点データのY方向積分および平均
  ! y_IntZ_zv, y_AvrZ_zv   :: 2 次元(ZY)格子点データのZ方向積分および平均
  !
  !==== 積分・平均(1 次元データ)
  !
  ! IntX_x, AvrX_x  :: 1 次元(X)格子点データのX方向積分および平均
  ! IntY_v, AvrY_v  :: 1 次元(Y)格子点データのY方向積分および平均
  ! IntZ_z, AvrZ_z  :: 1 次元(Z)格子点データのZ方向積分および平均
  !
  !==== 補間計算
  !
  ! Interpolate_tee :: スペクトルデータから任意の点の値を補間する.
  !
  use dc_message
  use mpi
  use lumatrix
  use ee_mpi_module, ls=>ks, le=>ke, lc=>kc, lf_l=>kf_k
  use at_module, only : z_Z => g_X, z_Z_WEIGHT => g_X_WEIGHT, &
                 at_Initial, at_az => at_ag, az_at => ag_at, &
                 t_z => t_g, z_t => g_t, &
                 t_Dz_t => t_Dx_t, at_Dz_at => at_Dx_at, &
                 a_Interpolate_at, &
                 at_Boundaries, t_Boundaries, &
                 at_BoundariesTau_DD, at_BoundariesTau_DN, &
                 at_BoundariesTau_ND, at_BoundariesTau_NN, &
                 at_BoundariesGrid_DD, at_BoundariesGrid_DN, &
                 at_BoundariesGrid_ND, at_BoundariesGrid_NN

  implicit none
  private

  public tee_mpi_Initial

  public ls, le, lc
  public js, je, jc
  public lf_l
  
  public x_X, x_X_Weight
  public v_Y, v_Y_Weight
  public z_Z, z_Z_Weight

  public vx_X, vx_Y, zv_Z, zv_Y, zx_X, zx_Z
  public zvx_X, zvx_Y, zvx_Z
  public zef_Z

  public tee_VMiss

  public ef_vx, vx_ef
  public at_Dz_at, t_Dz_t, az_at, at_az, z_t, t_z
  public zvx_tef, tef_zvx, zvx_zef, zef_zvx, zef_tef, tef_zef
  public e2z_e2t, e2t_e2z 
  public aef_e2a, e2a_aef, ef_e2, e2_ef
  public tef_Dx_tef, tef_Dy_tef, tef_Dz_tef
  public tef_Lapla_tef, tef_LaplaH_tef, tef_LaplaHInv_tef
  public tef_Div_zvx_zvx_zvx

  public zv_IntX_zvx, zx_IntY_zvx, vx_IntZ_zvx
  public x_IntZY_zvx, v_IntZX_zvx, z_IntYX_zvx
  public IntZYX_zvx

  public x_IntY_vx, v_IntX_vx, IntYX_vx
  public z_IntY_zv, v_IntZ_zv, IntZY_zv
  public z_IntX_zx, x_IntZ_zx, IntZX_zx
  public IntX_x, IntY_v, IntZ_z

  public zv_AvrX_zvx, zx_AvrY_zvx, vx_AvrZ_zvx
  public x_AvrZY_zvx, v_AvrZX_zvx, z_AvrYX_zvx
  public AvrZYX_zvx

  public x_AvrY_vx, v_AvrX_vx, AvrYX_vx
  public z_AvrY_zv, v_AvrZ_zv, AvrZY_zv
  public z_AvrX_zx, x_AvrZ_zx, AvrZX_zx
  public AvrX_x, AvrY_v, AvrZ_z

  public tef_ZRot_zvx_zvx, tef_ZRotRot_zvx_zvx_zvx
  public tef_Potential2vector, tef_Potential2Rotation

  public Interpolate_tef

  public zef_ToroidalEnergySpectrum_tef ! nz_ToroidalEnergySpectrum_wt
  public zef_PoloidalEnergySpectrum_tef ! nz_PoloidalEnergySpectrum_wt

  public tef_Boundaries, zef_LaplaPol2Pol_zef
  public tef_TorBoundaries, tef_LaplaPol2Pol_tef
  public tef_TormagBoundaries, tef_PolmagBoundaries

  public tef_BoundariesTau, tef_LaplaPol2PolTau_tef, tef_TorBoundariesTau
  public tef_TormagBoundariesTau, tef_PolmagBoundariesTau

  public tef_BoundariesGrid, tef_LaplaPol2PolGrid_tef, tef_TorBoundariesGrid
  public tef_TormagBoundariesGrid, tef_PolmagBoundariesGrid

  public at_Boundaries, t_Boundaries
  
  interface tef_Boundaries
     module procedure tef_BoundariesTau
  end interface

  interface tef_TorBoundaries
     module procedure tef_TorBoundariesTau
  end interface

  interface tef_LaplaPol2Pol_tef
     module procedure tef_LaplaPol2PolTau_tef
  end interface

  interface tef_TorMagBoundaries
     module procedure tef_TorMagBoundariesTau
  end interface

  interface tef_PolMagBoundaries
     module procedure tef_PolMagBoundariesTau
  end interface

  integer            :: im=32, jm=32, km=16  ! 格子点の設定(X,Y,Z)
  integer            :: lm=10, mm=10, nm=16  ! 切断波数の設定(水平X,Y, 鉛直Z)
  real(8)            :: xl, yl, zl           ! 領域の大きさ(水平X,Y, 鉛直Z)

  real(8), parameter :: pi=3.1415926535897932385D0

  real(8), dimension(:,:,:), allocatable :: zvx_X, zvx_Y, zvx_Z    ! 座標
  real(8), dimension(:,:),   allocatable :: zv_Z, zv_Y             ! 座標
  real(8), dimension(:,:),   allocatable :: zx_Z, zx_X             ! 座標
  real(8), dimension(:,:,:), allocatable :: zef_Z                  ! 座標

  real(8) :: tee_VMiss = -999.0        ! 欠損値

!  integer :: np                                         ! MPI プロセス数
  integer :: ip                                         ! MPI プロセス ID

  save im, jm, km, lm, mm, nm, xl, yl, zl

  contains
  !--------------- 初期化 -----------------
   subroutine tee_mpi_Initial(i,j,k,l,m,n,xmin,xmax,ymin,ymax,zmin,zmax)
     !
     ! スペクトル変換の格子点数, 波数, 各座標の範囲を設定する.
     !
     ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
     ! しなければならない.
     !
     integer,intent(in) :: i              ! 格子点数(水平X)
     integer,intent(in) :: j              ! 格子点数(水平Y)
     integer,intent(in) :: k              ! 格子点数(鉛直Z)
     integer,intent(in) :: l              ! 切断波数(水平X波数)
     integer,intent(in) :: m              ! 切断波数(水平Y波数)
     integer,intent(in) :: n              ! 切断波数(鉛直Z波数)

     real(8),intent(in) :: xmin, xmax     ! X 方向領域
     real(8),intent(in) :: ymin, ymax     ! Y 方向領域
     real(8),intent(in) :: zmin, zmax     ! Z 方向領域

     integer :: ierr     

     im = i  ; jm = j ; km = k
     lm = l  ; mm = m ; nm = n
     xl = xmax - xmin
     yl = ymax - ymin
     zl = zmax - zmin

     CALL MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)
 
     call ee_mpi_Initial(im,jm,lm,mm,xmin,xmax,ymin,ymax)

     call at_Initial(km,nm,zmin,zmax)

     allocate(zvx_X(0:km,jc(ip),0:im-1))
     allocate(zvx_Y(0:km,jc(ip),0:im-1))
     allocate(zvx_Z(0:km,jc(ip),0:im-1))

     allocate(zv_Z(0:km,jc(ip)))
     allocate(zv_Y(0:km,jc(ip)))
     allocate(zx_Z(0:km,0:im-1))
     allocate(zx_X(0:km,0:im-1))

     allocate(zef_Z(0:km,-mm:mm,2*lc(ip)))

     zvx_X = spread(vx_X,1,km+1)
     zvx_Y = spread(vx_Y,1,km+1)
     zvx_Z = spread(spread(z_Z,2,jc(ip)),3,im)

     zv_Z = spread(z_Z,2,jc(ip))
     zv_Y = spread(v_Y,1,km+1)
     zx_Z = spread(z_Z,2,im)
     zx_X = spread(x_X,1,km+1)

     zef_Z = spread(spread(z_Z,2,2*mm+1),3,2*lc(ip))

     call MessageNotify('M','tee_mpi_initial', &
          'tee_mpi_module (2016/09/22) is initialized')

   end subroutine tee_mpi_Initial

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

    function zvx_tef(tef)
      !
      ! スペクトルデータから 3 次元格子点データへ(逆)変換する.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
      !(in) 2 重フーリエチェビシェフスペクトルデータ
      real(8), dimension(0:km,jc(ip),0:im-1)             :: zvx_tef
      !(out) 3 次元格子点データ

      zvx_tef = zvx_zef(zef_tef(tef))

    end function zvx_tef

    function tef_zvx(zvx)
      !
      ! 3 次元格子点データからスペクトルデータへ(正)変換する.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元格子点データ
      real(8), dimension(0:nm,-mm:mm,2*lc(ip))             :: tef_zvx
      !(out) 2 重フーリエチェビシェフスペクトルデータ

      tef_zvx = tef_zef(zef_zvx(zvx))

    end function tef_zvx

    function zvx_zef(zef)
      !
      ! 水平スペクトル・鉛直格子点データから 3 次元格子点データへ(逆)変換する.
      !
      real(8), dimension(0:km,-mm:mm,2*lc(ip)), intent(in) :: zef
      !(in) 2 次元水平スペクトル・鉛直格子点データ
      real(8), dimension(0:km,jc(ip),0:im-1)               :: zvx_zef
      !(out) 3 次元格子点データ
      real(8), dimension(-mm:mm,2*lc(ip))                :: ef
      real(8), dimension(jc(ip),0:im-1)                  :: vx

      integer :: k

      do k = 0, km
         !zyx_zef(k,:,:) = yx_ee(zef(k,:,:))
         !
         ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
         !
         ef = zef(k,-mm:mm,1:2*lc(ip))
         vx = vx_ef(ef)
         zvx_zef(k,:,:) = vx
      enddo

    end function zvx_zef

    function zef_zvx(zvx)
      !
      ! 3 次元格子データから水平スペクトル・鉛直格子点データへ(正)変換する.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in)   :: zvx
      !(in) 3 次元格子点データ
      real(8), dimension(0:km,-mm:mm,2*lc(ip))             :: zef_zvx
      !(out) 2 次元スペクトル・鉛直格子点データ
      real(8), dimension(-mm:mm,2*lc(ip))                :: ef
      real(8), dimension(jc(ip),0:im-1)                  :: vx

      integer :: k

      do k = 0, km
         !zef_zyx(k,:,:) = ee_yx(zyx(k,:,:))
         !
         ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
         !
         vx = zvx(k,1:jc(ip),0:im-1)
         ef = ef_vx(vx)
         zef_zvx(k,:,:) = ef
      enddo

    end function zef_zvx

    function zef_tef(tef)
      !
      ! スペクトルデータから水平スペクトル・鉛直格子点データへ(正)変換する.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
      !(in) 2 次元水平スペクトル鉛直チェビシェフスペクトルデータ
      real(8), dimension(0:km,-mm:mm,2*lc(ip))             :: zef_tef
      !(out) 2 次元水平スペクトル・鉛直格子点データ

      zef_tef = aef_e2a(e2z_e2t(e2a_aef(tef)))

    end function zef_tef

    function tef_zef(zef)
      !
      ! 水平スペクトル・鉛直格子点データからスペクトルデータへ(正)変換する.
      !
      real(8), dimension(0:km,-mm:mm,2*lc(ip)), intent(in) :: zef
      !(in) 2 次元水平スペクトル・鉛直格子点データ
      real(8), dimension(0:nm,-mm:mm,2*lc(ip))             :: tef_zef
      !(out) 2 次元水平鉛直チェビシェフスペクトルデータ

      tef_zef = aef_e2a(e2t_e2z(e2a_aef(zef)))

    end function tef_zef

    function e2z_e2t(e2t)
      !
      ! 鉛直スペクトルを格子点に変換する
      !
      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm)  :: e2t
      !(in)  1 次元化された水平スペクトル・鉛直スペクトル
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km)  :: e2z_e2t
      !(out) 1 次元化された水平スペクトル・鉛直格子座標

      e2z_e2t = az_at(e2t)

    end function e2z_e2t

    function e2t_e2z(e2z)
      !
      ! 鉛直格子点をスペクトルに変換する
      !
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km)  :: e2z
      !(in)  1 次元化された水平スペクトル・鉛直格子座標
      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm)  :: e2t_e2z
      !(out) 1 次元化された水平スペクトル・鉛直スペクトル

      e2t_e2z = at_az(e2z)

    end function e2t_e2z

    function e2a_aef(aef)
      !
      ! 水平スペクトルを 1 次元化転置する.
      !
      real(8), dimension(:,-mm:,:), intent(in)        :: aef
      !(in) 任意座標・水平 2 次元スペクトルデータ dimension(:,-mm:mm,2*lc(ip))
      real(8), dimension(2*lc(ip)*(2*mm+1),size(aef,1))  :: e2a_aef
      !(out) 1 次元化された水平スペクトル・任意座標データ

      integer :: j,k,l,m

      if ( size(aef,2) /= 2*mm+1 ) &
           call MessageNotify('E','e2a_aef',&
                              '2nd dimension of input data invalid')
      if ( size(aef,3) /= 2*lc(ip) ) &
           call MessageNotify('E','e2a_aef',&
                              '3rd dimension of input data invalid')

      !e2a_aee = transpose(reshape(aee,(/size(aee,1),2*lc(ip)*(2*mm+1)/)))
      !
      ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
      !

      do k=1,size(aef,1)
         do l=1,2*lc(ip)
            do m=-mm,mm
               j = (m+mm+1) + (2*mm+1)*(l-1)
               e2a_aef(j,k) = aef(k,m,l)
            enddo
         enddo
      enddo

    end function e2a_aef

    function aef_e2a(e2a)
      !
      ! 水平スペクトルを転置展開する.
      !
      real(8), dimension(:,:),intent(in)                 :: e2a
      !(in) 1 次元化された水平スペクトル・任意座標データ
      !     dimmension((2*mm*1)*(2*lm*1),:)
      real(8), dimension(size(e2a,2),-mm:mm,2*lc(ip))      :: aef_e2a
      !(out) 任意座標・水平 2 次元スペクトルデータ

      integer :: j,k,l,m

      if ( size(e2a,1) /= (2*mm+1)*2*lc(ip) ) &
           call MessageNotify('E','aef_e2a',&
                              '1st dimension of input data invalid')

      !aee_e2a = reshape(transpose(e2a),(/size(e2a,2),2*mm+1,2*lc(ip)/))
      !
      ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
      !

      do k=1,size(e2a,2)
         do l=1,2*lc(ip)
            do m=-mm,mm
               j = (m+mm+1) + (2*mm+1)*(l-1)
               aef_e2a(k,m,l) = e2a(j,k)
            enddo
         enddo
      enddo

    end function aef_e2a

    function e2_ef(ef)
      !
      ! 水平スペクトルを 1 次元化転置する.
      !
      real(8), dimension(-mm:,:), intent(in)        :: ef
      !(in) 任意座標・水平 2 次元スペクトルデータ dimension(-mm:mm,2*lc(ip))
      real(8), dimension(2*lc(ip)*(2*mm+1))            :: e2_ef
      !(out) 1 次元化された水平スペクトル・任意座標データ

      integer :: j,l,m

      !e2_ee = reshape(e2a_aee(reshape(ee,(/1,2*mm+1,2*lc(ip)/))), &
      !                (/lm*(2*mm+1)/))
      !
      ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
      !

      do l=1,2*lc(ip)
         do m=-mm,mm
            j = (m+mm+1)+ (2*mm+1)*(l-1)
            e2_ef(j) = ef(m,l)
         enddo
      enddo

    end function e2_ef

    function ef_e2(e2)
      !
      ! 水平スペクトルを転置展開する.
      !
      real(8), dimension(:),intent(in)                   :: e2
      !(in) 1 次元化された水平スペクトル・任意座標データ
      !     dimmension((2*mm*1)*2*lc(ip),:)
      real(8), dimension(-mm:mm,2*lc(ip))                  :: ef_e2
      !(out) 任意座標・水平 2 次元スペクトルデータ

      integer :: j,l,m

      !ee_e2 = reshape(aee_e2a(reshape(e2,(/(2*mm+1)*2*lc(ip),1/))), &
      !                (/2*mm+1,2*lc(ip)/))
      !
      ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
      !

      do l=1,2*lc(ip)
         do m=-mm,mm
            j = (m+mm+1) + (2*mm+1)*(l-1)
            ef_e2(m,l) = e2(j)
         enddo
      enddo

    end function ef_e2

  !--------------- 微分計算 -----------------
    function tef_Dx_tef(tef)
      !
      ! 入力スペクトルデータに水平微分 ∂/∂x を作用する.
      !
      ! スペクトルデータのX微分とは, 対応する格子点データにX微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
      !(in) 入力スペクトルデータ

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))             :: tef_Dx_tef
      !(in) X微分されたスペクトルデータ

      integer :: n

      do n=0,nm
         tef_Dx_tef(n,:,:) = ef_Dx_ef(tef(n,:,:))
      enddo

    end function tef_Dx_tef

    function tef_Dy_tef(tef)
      !
      ! 入力スペクトルデータに水平微分 ∂/∂y を作用する.
      !
      ! スペクトルデータのY微分とは, 対応する格子点データにY微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
      !(in) 入力スペクトルデータ

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))             :: tef_Dy_tef
      !(in) X微分されたスペクトルデータ

      integer :: n

      do n=0,nm
         tef_Dy_tef(n,:,:) = ef_Dy_ef(tef(n,:,:))
      enddo

    end function tef_Dy_tef

    function tef_Dz_tef(tef)
      !
      ! 入力スペクトルデータに鉛直微分 ∂/∂z を作用する.
      !
      ! スペクトルデータの鉛直微分とは, 対応する格子点データに鉛直微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
      !(in) 入力スペクトルデータ

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))             :: tef_Dz_tef
      !(in) 鉛直微分されたスペクトルデータ

      tef_Dz_tef = aef_e2a(at_Dz_at(e2a_aef(tef)))

    end function tef_Dz_tef

    function tef_Lapla_tef(tef)
      !
      ! 入力スペクトルデータにラプラシアン
      !
      !     ▽^2 = ∂^2/∂X^2 + ∂^2/∂Y^2 + ∂^2/∂Z^2
      !
      ! を作用する.
      !
      ! スペクトルデータのラプラシアンとは, 対応する格子点データに
      ! ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))             :: tef_Lapla_tef
      !(out) ラプラシアンを作用された 2 次元スペクトルデータ

      tef_Lapla_tef = tef_LaplaH_tef(tef) + tef_Dz_tef(tef_Dz_tef(tef))

    end function tef_Lapla_tef

    function tef_LaplaH_tef(tef)
      !
      ! 入力スペクトルデータに水平ラプラシアン
      !
      !     ▽^2_H = ∂^2/∂X^2 + ∂^2/∂Y^2
      !
      ! を作用する.
      !
      ! スペクトルデータの水平ラプラシアンとは, 対応する格子点データに
      ! 水平ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
      !(in) 入力スペクトルデータ

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))             :: tef_LaplaH_tef
      !(out) 水平ラプラシアンを作用された 2 次元スペクトルデータ

      integer :: n

      do n=0,nm
         tef_LaplaH_tef(n,:,:) = ef_Lapla_ef(tef(n,:,:))
      enddo

    end function tef_LaplaH_tef

    function tef_LaplaHInv_tef(tef)
      !
      ! 入力スペクトルデータに逆水平ラプラシアン
      !
      !     ▽^-2_H = (∂^2/∂X^2 + ∂^2/∂Y^2)^-1
      !
      ! を作用する.
      !
      ! スペクトルデータの逆水平ラプラシアンとは, 対応する格子点データに
      ! 逆水平ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef
      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))             :: tef_LaplaHInv_tef
      !(out) 逆水平ラプラシアンを作用された 2 次元スペクトルデータ

      integer :: n

      do n=0,nm
         tef_LaplaHInv_tef(n,:,:) = ef_LaplaInv_ef(tef(n,:,:))
      enddo

    end function tef_LaplaHInv_tef

    function tef_Div_zvx_zvx_zvx(zvx_VX,zvx_VY,zvx_VZ)
      !
      ! べクトル成分である 3 つの格子データに発散を作用させた
      ! スペクトルデータを返す.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VX
      !(in) ベクトル場のX成分
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VY
      !(in) ベクトル場のY成分

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

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))     :: tef_Div_zvx_zvx_zvx
      !(out) ベクトル場の発散

      tef_Div_zvx_zvx_zvx =   tef_Dx_tef(tef_zvx(zvx_VX)) &
                            + tef_Dy_tef(tef_zvx(zvx_VY)) &
                            + tef_Dz_tef(tef_zvx(zvx_VZ))

    end function tef_Div_zvx_zvx_zvx

  !--------------- 積分計算 -----------------
    !----(入力データ zvx)---
    function zv_IntX_zvx(zvx)  ! X積分
      !
      ! 3 次元格子点データのX方向積分.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元格子点データ

      real(8), dimension(0:km,jc(ip))  :: zv_IntX_zvx
      !(out) X方向積分された 2 次元ZV格子点データ

      integer :: i

      zv_IntX_zvx = 0.0d0
      do i=0,im-1
         zv_IntX_zvx(:,:) = zv_IntX_zvx(:,:) &
                       + zvx(:,:,i) * x_X_Weight(i)
      enddo
    end function zv_IntX_zvx

    function zx_IntY_zvx(zvx)
      !
      ! 3 次元格子点データのY方向域積分.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元格子点データ

      real(8), dimension(0:km,0:im-1)  :: zx_IntY_zvx
      !(out) Y積分された 2 次元ZV格子点データ.

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

      zx_IntYTmp = 0.0d0
      do j=1,jc(ip)
         zx_IntYTmp(:,:) = zx_IntYTmp(:,:) &
                       + zvx(:,j,:) * v_Y_Weight(j)
      enddo

      CALL MPI_ALLREDUCE(zx_IntYTmp,zx_IntY_zvx,im*(km+1),MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)

    end function zx_IntY_zvx

    function vx_IntZ_zvx(zvx)  ! Z積分
      !
      ! 3 次元格子点データのZ方向積分.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元格子点データ

      real(8), dimension(jc(ip),0:im-1)  :: vx_IntZ_zvx
      !(out) Z積分された 2 次元VX(水平)格子点データ

      integer :: k

      vx_IntZ_zvx = 0.0d0
      do k=0,km
         vx_IntZ_zvx(:,:) = vx_IntZ_zvx(:,:) &
                       + zvx(k,:,:) * z_Z_Weight(k)
      enddo
    end function vx_IntZ_zvx

    function x_IntZY_zvx(zvx)
      !
      ! 3 次元格子点データのZV積分
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元ZVX格子点データ

      real(8), dimension(0:im-1)     :: x_IntZY_zvx
      !(out) ZV(子午面)積分された 1 次元X格子点データ

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

      x_IntZYtmp = 0.0D0
      do j=1,jc(ip)
         do k=0,km
            x_IntZYTmp = x_IntZYTmp &
                 + zvx(k,j,:) * v_Y_Weight(j) * z_Z_Weight(k)
         enddo
      enddo

      CALL MPI_ALLREDUCE(x_IntZYTmp,x_IntZY_zvx,im,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)
      
    end function x_IntZY_zvx

    function v_IntZX_zvx(zvx)
      !
      ! 3 次元格子点データのZX積分.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元ZVX格子点データ

      real(8), dimension(jc(ip))       :: v_IntZX_zvx
      !(out) ZX積分された 1 次元Y格子点データ

      integer :: i, k

      v_IntZX_zvx = 0
      do i=0,im-1
         do k=0,km
            v_IntZX_zvx = v_IntZX_zvx &
                 + zvx(k,:,i) * x_X_Weight(i) * z_Z_Weight(k)
         enddo
      enddo
    end function v_IntZX_zvx

    function z_IntYX_zvx(zvx)  ! VX(水平)積分
      !
      ! 3 次元格子点データのVX(水平)積分
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元ZVX格子点データ

      real(8), dimension(0:km)     :: z_IntYX_zvx
      !(out) VX(水平)積分された 1 次元Z格子点データ

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

      z_IntYXTmp = 0
      do j=1,jc(ip)
         do i=0,im-1
            z_IntYXTmp = z_IntYXTmp &
                 + zvx(:,j,i) * x_X_Weight(i) * v_Y_Weight(j)
         enddo
      enddo
      
      CALL MPI_ALLREDUCE(z_IntYXTmp,z_IntYX_zvx,km+1,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)
      
    end function z_IntYX_zvx

    function IntZYX_zvx(zvx) ! ZVX(全領域)積分
      !
      ! 3 次元格子点データのZVX(全領域)積分
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元格子点データ

      real(8)                     :: IntZYX_zvx
      !(out) 全領域積分値

      real(8)                     :: IntZYXTmp
      integer :: i, j, k, ierr

      IntZYXTmp = 0
      do i=0,im-1
         do j=1,jc(ip)
            do k=0,km
               IntZYXTmp = IntZYXTmp &
                    + zvx(k,j,i) * x_X_Weight(i) &
                         * v_Y_Weight(j) * z_Z_Weight(k)
            enddo
         enddo
      enddo
      
      CALL MPI_ALLREDUCE(IntZYXTmp,IntZYX_zvx,1,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)
      
    end function IntZYX_zvx

    !----(入力データ zv)---
    function z_IntY_zv(zv)  ! Y積分
      !
      ! 2 次元(ZV)格子点データのY方向積分.
      !
      real(8), dimension(0:km,jc(ip)), intent(in) :: zv
      !(in) 2 次元ZV(子午面)格子点データ

      real(8), dimension(0:km)  :: z_IntY_zv
      !(out) Y積分された 1 次元Z格子点データ

      real(8), dimension(0:km)  :: z_IntYTmp
      integer :: j, ierr

      z_IntYTmp = 0.0d0
      do j=1,jc(ip)
         z_IntYTmp(:) = z_IntYTmp(:) + zv(:,j) * v_Y_Weight(j)
      enddo

      CALL MPI_ALLREDUCE(z_IntYTmp,z_IntY_zv,km+1,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)

    end function z_IntY_zv

    function v_IntZ_zv(zv)  ! Z積分
      !
      ! 2 次元(ZV)格子点データのZ方向域積分.
      !
      real(8), dimension(0:km,jc(ip)), intent(in) :: zv
      !(in) 2 次元ZV格子点データ

      real(8), dimension(jc(ip))  :: v_IntZ_zv
      !(out) Z積分された 1 次元Y格子点データ

      integer :: k

      v_IntZ_zv = 0.0d0
      do k=0,km
         v_IntZ_zv(:) = v_IntZ_zv(:) &
                       + zv(k,:) * z_Z_Weight(k)
      enddo
    end function v_IntZ_zv

    function IntZY_zv(zv)
      !
      ! 2 次元(ZV)格子点データのZV積分
      !
      real(8), dimension(0:km,jc(ip)), intent(in) :: zv
      !(in) 2 次元ZV(子午面)格子点データ

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

      real(8) :: IntZYTmp
      integer :: j, k, ierr

      IntZYTmp = 0
      do j=1,jc(ip)
         do k=0,km
            IntZYTmp = IntZYTmp &
                 + zv(k,j) * v_Y_Weight(j) * z_Z_Weight(k)
         enddo
      enddo

      CALL MPI_ALLREDUCE(IntZYTmp,IntZY_zv,1,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)
      
    end function IntZY_zv

    !----(入力データ zx)---
    function z_IntX_zx(zx)
      !
      ! 2 次元(ZX)格子点データのX方向積分.
      !
      real(8), dimension(0:km,0:im-1), intent(in) :: zx
      !(in) 2 次元ZV格子点データ

      real(8), dimension(0:km)  :: z_IntX_zx
      !(out) X積分された 1 次元Z格子点データ

      integer :: i

      z_IntX_zx = 0.0d0
      do i=0,im-1
         z_IntX_zx(:) = z_IntX_zx(:) + zx(:,i) * x_X_Weight(i)
      enddo

    end function z_IntX_zx

    function x_IntZ_zx(zx)
      !
      ! 2 次元(ZX)格子点データのZ方向積分.
      !
      real(8), dimension(0:km,0:im-1), intent(in) :: zx
      !(in) 2 次元ZV格子点データ

      real(8), dimension(0:im-1)  :: x_IntZ_zx
      !(out) Z積分された 1 次元X格子点データ

      integer :: k

      x_IntZ_zx = 0.0d0
      do k=0,km
         x_IntZ_zx(:) = x_IntZ_zx(:) &
                       + zx(k,:) * z_Z_Weight(k)
      enddo

    end function x_IntZ_zx

    function IntZX_zx(zx)  ! ZX積分
      !
      ! 2 次元(ZX)格子点データのZX積分
      !
      real(8), dimension(0:km,0:im-1), intent(in) :: zx
      !(in) 2 次元ZV格子点データ

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

      integer :: i, k

      IntZX_zx = 0
      do i=0,im-1
         do k=0,km
            IntZX_zx = IntZX_zx &
                 + zx(k,i) * x_X_Weight(i) * z_Z_Weight(k)
         enddo
      enddo
    end function IntZX_zx

    !----(入力データ z)---
    function IntZ_z(z)  ! Z積分
      !
      ! 1 次元(Z)格子点データのZ方向積分.
      !
      real(8), dimension(0:km), intent(in) :: z
      !(in) 1 次元Z格子点データ

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

      integer :: k

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

  !--------------- 平均計算 -----------------
    !----(入力データ zvx)---
    function zv_AvrX_zvx(zvx)  ! X積分
      !
      ! 3 次元格子点データのX方向平均.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元ZVX格子点データ

      real(8), dimension(0:km,jc(ip))  :: zv_AvrX_zvx
      !(out) X方向平均された 2 次元子午面格子点データ

      zv_AvrX_zvx = zv_IntX_zvx(zvx)/sum(x_X_Weight)

    end function zv_AvrX_zvx

    function zx_AvrY_zvx(zvx)  ! Y平均
      !
      ! 3 次元格子点データのY方向域平均.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元格子点データ

      real(8), dimension(0:km,0:im-1)  :: zx_AvrY_zvx
      !(out) Y平均された 2 次元ZX格子点データ

      zx_AvrY_zvx = zx_IntY_zvx(zvx)/yl

    end function zx_AvrY_zvx

    function vx_AvrZ_zvx(zvx)
      !
      ! 3 次元格子点データのZ方向平均.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元格子点データ

      real(8), dimension(jc(ip),0:im-1)  :: vx_AvrZ_zvx
      !(out) Z平均された 2 次元VX(水平)格子点データ

      vx_AvrZ_zvx = vx_IntZ_zvx(zvx)/sum(z_Z_Weight)

    end function vx_AvrZ_zvx

    function x_AvrZY_zvx(zvx)  ! ZV積分
      !
      ! 3 次元格子点データのZV平均
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元ZVX格子点データ

      real(8), dimension(0:im-1)     :: x_AvrZY_zvx
      !(out) ZV平均された 1 次元X格子点データ

      x_AvrZY_zvx = x_IntZY_zvx(zvx)/(yl*zl)

    end function x_AvrZY_zvx

    function v_AvrZX_zvx(zvx)  ! ZX積分
      !
      ! 3 次元格子点データのZX平均.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元格子点データ

      real(8), dimension(jc(ip))       :: v_AvrZX_zvx
      !(out) ZX平均された 1 次元Y格子点データ

      v_AvrZX_zvx = v_IntZX_zvx(zvx) &
                 /(sum(x_X_Weight)*sum(z_Z_Weight))

    end function v_AvrZX_zvx

    function z_AvrYX_zvx(zvx)  ! VX(水平)積分
      !
      ! 3 次元格子点データのVX(水平)積分
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元格子点データ

      real(8), dimension(0:km)     :: z_AvrYX_zvx
      !(out) VX(水平)平均された 1 次元Z格子点データ

      z_AvrYX_zvx = z_IntYX_zvx(zvx)/(xl*yl)

    end function z_AvrYX_zvx

    function AvrZYX_zvx(zvx) ! ZVX(全領域)積分
      !
      ! 3 次元格子点データのZVX(全領域)積分
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx
      !(in) 3 次元ZVX格子点データ

      real(8)                     :: AvrZYX_zvx
      !(out) 全領域平均値

      AvrZYX_zvx = IntZYX_zvx(zvx)/(xl*yl*zl)

    end function AvrZYX_zvx

    !----(入力データ zv)---
    function z_AvrY_zv(zv)
      !
      ! 2 次元(ZV)格子点データのY方向平均.
      !
      real(8), dimension(0:km,jc(ip)), intent(in) :: zv
      !(in) 2 次元ZV格子点データ

      real(8), dimension(0:km)  :: z_AvrY_zv
      !(out) Y平均された 1 次元Z格子点データ

      z_AvrY_zv = z_IntY_zv(zv)/yl

    end function z_AvrY_zv

    function v_AvrZ_zv(zv)
      !
      ! 2 次元(ZV)格子点データのZ方向平均.
      !
      real(8), dimension(0:km,jc(ip)), intent(in) :: zv
      !(in) 2 次元ZV格子点データ

      real(8), dimension(jc(ip))  :: v_AvrZ_zv
      !(out) Z平均された 1 次元Y格子点データ

      v_AvrZ_zv = v_IntZ_zv(zv)/sum(z_Z_Weight)

    end function v_AvrZ_zv

    function AvrZY_zv(zv)  ! ZV平均
      !
      ! 2 次元(ZV)格子点データのZV平均
      !
      real(8), dimension(0:km,jc(ip)), intent(in) :: zv
      !(in) 2 次元ZV(子午面)格子点データ

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

      AvrZY_zv = IntZY_zv(zv)/(yl*zl)

    end function AvrZY_zv

    !----(入力データ zx)---
    function z_AvrX_zx(zx)  ! X積分
      !
      ! 2 次元(ZX)格子点データのX方向平均.
      !
      real(8), dimension(0:km,0:im-1), intent(in) :: zx
      !(in) 2 次元ZV格子点データ

      real(8), dimension(0:km)  :: z_AvrX_zx
      !(out) X平均された 1 次元Z格子点データ

      z_AvrX_zx = z_IntX_zx(zx)/sum(x_X_Weight)

    end function z_AvrX_zx

    function x_AvrZ_zx(zx)  ! Z平均
      !
      ! 2 次元(ZX)格子点データのZ方向平均.
      !
      real(8), dimension(0:km,0:im-1), intent(in) :: zx
      !(in) 2 次元ZV格子点データ

      real(8), dimension(0:im-1)  :: x_AvrZ_zx
      !(out) Z平均された 1 次元X格子点データ

      x_AvrZ_zx = x_IntZ_zx(zx)/sum(z_Z_Weight)

    end function x_AvrZ_zx

    function AvrZX_zx(zx)  ! ZX積分
      !
      ! 2 次元(ZX)格子点データのZX平均
      !
      real(8), dimension(0:km,0:im-1), intent(in) :: zx
      ! (in) 2 次元格子点データ
      real(8)                                 :: AvrZX_zx
      ! 平均値

      AvrZX_zx = IntZX_zx(zx)/(sum(x_X_Weight)*sum(z_Z_Weight))

    end function AvrZX_zx

    !----(入力データ z)---
    function AvrZ_z(z)
      !
      ! 1 次元(Z)格子点データのZ方向平均.
      !
      real(8), dimension(0:km), intent(in) :: z
      !(in) 1 次元Z格子点データ
      real(8)                              :: AvrZ_z
      !(out) 平均値

      AvrZ_z = IntZ_z(z)/sum(z_Z_Weight)

    end function AvrZ_z

  !--------------- ポロイダル/トロイダルモデル用微分 -----------------

!!$    function wt_KxRGrad_wt(wt)
!!$      !
!!$      ! 入力スペクトルデータにX微分 k×r・▽ = ∂/∂λを作用する.
!!$      !
!!$      real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$      real(8), dimension(-lm:lm,-mm,mm,0:nm)             :: wt_KxRGrad_wt
!!$      !(out) X微分を作用された 2 次元スペクトルデータ
!!$
!!$      wt_KxRGrad_wt =  wa_Dlon_wa(wt)
!!$
!!$    end function wt_KxRGrad_wt
!!$
!!$    function zvx_KGrad_wt(wt)    ! k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r
!!$      !
!!$      ! 入力スペクトルデータに対応する格子データに軸方向微分
!!$      !
!!$      !    k・▽ = cosφ/r ∂/∂φ + sinφ∂/∂r
!!$      !
!!$      ! を作用させた格子データが返される.
!!$      ! ここでベクトル k は球の中心から北極向きの単位ベクトルである.
!!$      !
!!$      real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$      real(8), dimension(0:km,0:jm-1,0:im-1)                     :: yzx_KGrad_wt
!!$      !(out) 軸方向微分を作用された 2 次元スペクトルデータ
!!$
!!$      xzv_KGrad_wt =  cos(xyz_Y)*xyz_GradY_wt(wt) &
!!$                    + sin(xyz_Y)*xyz_wt(wt_Drad_wt(wt))
!!$
!!$    end function xyz_KGrad_wt
!!$
!!$    function wt_L2_wt(wt)
!!$      !
!!$      ! 入力スペクトルデータに L^2 演算子(=-水平ラプラシアン)を作用する.
!!$      !
!!$      ! L^2 演算子は単位球面上の水平ラプラシアンの逆符号にあたる.
!!$      !  入力スペクトルデ ータに対応する格子点データに演算子
!!$      !
!!$      !     L^2 = -1/cos^2φ・∂^2/∂λ^2 - 1/cosφ・∂/∂φ(cosφ∂/∂φ)
!!$      !
!!$      ! を作用させたデータのスペクトル変換が返される.
!!$      !
!!$      real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$      real(8), dimension(-lm:lm,-mm,mm,0:nm)             :: wt_L2_wt
!!$      !(out) L^2 演算子を作用された 2 次元スペクトルデータ
!!$
!!$      wt_L2_wt = -wa_Lapla_wa(wt)
!!$
!!$    end function wt_L2_wt
!!$
!!$    function wt_L2Inv_wt(wt)
!!$      !
!!$      ! 入力スペクトルデータに L^2 演算子の逆演算(-逆水平ラプラシアン)を
!!$      ! 作用する.
!!$      !
!!$      ! スペクトルデータに L^2 演算子を作用させる関数 wt_L2_wt の逆計算を
!!$      ! 行う関数である.
!!$      !
!!$      real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$      real(8), dimension(-lm:lm,-mm,mm,0:nm)             :: wt_L2Inv_wt
!!$      !(out) L^2 演算子の逆演算を作用された 2 次元スペクトルデータ
!!$
!!$      wt_L2Inv_wt = -wa_LaplaInv_wa(wt)
!!$
!!$    end function wt_L2Inv_wt
!!$
!!$    function wt_QOperator_wt(wt)
!!$      !
!!$      ! 入力スペクトルデータに対応する格子点データに演算子
!!$      !
!!$      !    Q=(k・▽-1/2(L2 k・▽+ k・▽L2))
!!$      !
!!$      ! を作用させたデータのスペクトル変換が返される.
!!$      !
!!$      real(8), dimension(-lm:lm,-mm,mm,0:nm), intent(in) :: wt
!!$      !(in) 2 次元球面調和函数チェビシェフスペクトルデータ
!!$
!!$      real(8), dimension(-lm:lm,-mm,mm,0:nm)             :: wt_QOperator_wt
!!$      !(out) Q 演算子を作用された 2 次元スペクトルデータ
!!$
!!$      wt_QOperator_wt = &
!!$             wt_xyz(xyz_KGrad_wt(wt) - xyz_KGrad_wt(wt_L2_wt(wt))/2) &
!!$           - wt_L2_wt(wt_xyz(xyz_KGrad_wt(wt)))/2
!!$
!!$    end function wt_QOperator_wt

    function tef_ZRot_zvx_zvx(zvx_VX,zvx_VY)  ! z・(▽×v)
      !
      ! ベクトルの渦度とZベクトルの内積 z・(▽×v) を計算する.
      !
      ! 第 1, 2 引数(v[x], v[y])がそれぞれベクトルのX成分, Y成分を表す.
      !
      !    z・(▽×v) = ∂v[y]/∂x - ∂(v[x])/∂y
      !
      ! のスペクトル データが返される.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VX
      !(in) ベクトルのX成分

      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VY
      !(in) ベクトルのY成分

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))     :: tef_ZRot_zvx_zvx
      !(out) ベクトルの渦度とZベクトルの内積

      tef_ZRot_zvx_zvx = tef_DX_tef(tef_zvx(zvx_VY)) &
                       - tef_DY_tef(tef_zvx(zvx_VX))

    end function tef_ZRot_zvx_zvx

    function tef_ZRotRot_zvx_zvx_zvx(zvx_VX,zvx_VY,zvx_VZ)
      !
      ! ベクトル v に対して z・(▽×▽×v) を計算する.
      !
      ! 第 1, 2, 3 引数(v[x], v[y], v[z])がそれぞれベクトルのX成分,
      ! Y成分, Z成分を表す.
      !
      !    z・(▽×▽×v)  = ∂/∂z (∂v[x]/∂x + ∂(v[y])/∂y ) )
      !                    - ▽_H^2 v[z]
      !
      ! のスペクトルデータが返される.
      !
      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VX
      !(in) ベクトルのX成分

      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VY
      !(in) ベクトルのY成分

      real(8), dimension(0:km,jc(ip),0:im-1), intent(in) :: zvx_VZ
      !(in) ベクトルのZ成分

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))     :: tef_ZRotRot_zvx_zvx_zvx
      !(out) ベクトル v の z・(▽×▽×v)

      tef_ZRotRot_zvx_zvx_zvx = &
               tef_DZ_tef(   tef_DX_tef(tef_zvx(zvx_VX))   &
                           + tef_DY_tef(tef_zvx(zvx_VY)) ) &
             - tef_LaplaH_tef(tef_zvx(zvx_VZ))

    end function tef_ZRotRot_zvx_zvx_zvx

    subroutine tef_Potential2Vector(&
         zvx_VX,zvx_VY,zvx_VZ,tef_Torvel,tef_Polvel)
      !
      ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
      !
      !     v = ▽x(Ψz) + ▽x▽x(Φz)
      !
      ! の各成分を計算する
      !
      real(8), dimension(0:km,jc(ip),0:im-1)     :: zvx_VX
      !(out) ベクトル場のX成分

      real(8), dimension(0:km,jc(ip),0:im-1)     :: zvx_VY
      !(out) ベクトル場のY成分

      real(8), dimension(0:km,jc(ip),0:im-1)     :: zvx_VZ
      !(out) ベクトル場のZ成分

      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Torvel
      !(in) トロイダルポテンシャル

      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Polvel
      !(in) ポロイダルポテンシャル

      zvx_VX = zvx_tef(   tef_DY_tef(tef_Torvel) &
                        + tef_DX_tef(tef_DZ_tef(tef_Polvel))  )
      zvx_VY = zvx_tef( - tef_DX_tef(tef_Torvel) &
                        + tef_DY_tef(tef_DZ_tef(tef_Polvel))  )
      zvx_VZ = -zvx_tef(tef_LaplaH_tef(tef_Polvel))

    end subroutine tef_Potential2Vector

    subroutine tef_Potential2Rotation(&
       zvx_RotVX,zvx_RotVY,zvx_RotVZ,tef_Torvel,tef_Polvel)
      !
      ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
      !
      !     v = ▽x(Ψz) + ▽x▽x(Φz)
      !
      ! に対して, その回転
      !
      !     ▽xv = ▽x▽x(Ψz) + ▽x▽x▽x(Φz) = ▽x▽x(Ψz) - ▽x((▽^2Φ)z)
      !
      ! を計算する.

      ! ベクトル場の回転
      real(8), dimension(0:km,jc(ip),0:im-1), intent(OUT) :: zvx_RotVX
      !(out) 回転のX成分

      real(8), dimension(0:km,jc(ip),0:im-1), intent(OUT) :: zvx_RotVY
      !(out) 回転のY成分

      real(8), dimension(0:km,jc(ip),0:im-1), intent(OUT) :: zvx_RotVZ
      !(out) 回転のZ成分

      ! 入力ベクトル場を表すポテンシャル
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Torvel
      !(in) トロイダルポテンシャル

      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_Polvel
      !(in) ポロイダルポテンシャル

      call tef_Potential2Vector( &
           zvx_RotVX,zvx_RotVY,zvx_RotVZ, &
           -tef_Lapla_tef(tef_Polvel), tef_Torvel)

    end subroutine tef_Potential2Rotation

  !--------------- 補間計算 -----------------
    function Interpolate_tef(tef_data,x,y,z)
      !
      ! (x,y,z) における関数値を
      ! そのスペクトル係数 tef_data から補間計算する
      !
      real(8), intent(IN) :: tef_data(0:nm,-mm:mm,2*lc(ip))! スペクトルデータ
      real(8), intent(IN) :: x                             ! 補間する位置(X)
      real(8), intent(IN) :: y                             ! 補間する位置(Y)
      real(8), intent(IN) :: z                             ! 補間する位置(Z)
      real(8) :: Interpolate_tef                           ! 補間した値

      Interpolate_tef = &
           Interpolate_ef(ef_e2(a_Interpolate_at(e2a_aef(tef_data),z)),x,y)

    end function Interpolate_tef

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

    function zef_ToroidalEnergySpectrum_tef(tef_TORPOT)
      !
      ! トロイダルポテンシャルから, トロイダルエネルギーの
      ! 水平 X 波数 l, Y  波数 m の各成分を計算する
      !
      !  * X 波数 l, Y 波数 m のトロイダルポテンシャルのスペクトル成分
      !    ψ(l,m,z)から水平 X 波数 l, Y 波数 m 成分のトロイダルエネルギー
      !    スペクトルは  (1/2)(l**2+m**2)|ψ(l,m,z)|^2  と計算される.
      !
      !  * 全てのエネルギースペクトル成分の和をZ積分したものが領域内での
      !    水平平均エネルギーに等しい.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_TORPOT
      !(in) トロイダルポテンシャル

      real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_ToroidalEnergySpectrum_tef
      !(out) エネルギースペクトルトロイダル成分

      real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_Work
      integer :: l, m

      zef_ToroidalEnergySpectrum_tef = 0.0D0
      zef_Work = zef_tef(tef_Torpot)
      do l=-lm,lm
         do m=-mm,mm
            if ( lf_l(l) > 0 ) then
               zef_ToroidalEnergySpectrum_tef(:,m,lf_l(l)) &
                    = 0.5 * ((2*PI*l/xl)**2+(2*PI*m/yl)**2) &
                    * ( zef_Work(:,m,lf_l(l))**2 + zef_Work(:,-m,lf_l(-l))**2 )
            endif
         enddo
      enddo

    end function zef_ToroidalEnergySpectrum_tef

!!$    function nz_ToroidalEnergySpectrum_wt(wt_TORPOT)
!!$      !
!!$      ! トロイダルポテンシャルから, トロイダルエネルギーの
!!$      ! 球面調和函数全波数の各成分を計算する.
!!$      !
!!$      !  * 全波数 n, 帯状波数 m のトロイダルポテンシャルのスペクトル成分
!!$      !    ψ(n,m,r)から全波数 n 成分のトロイダルエネルギースペクトルは
!!$      !    Σ[m=-n]^n(1/2)n(n+1)4πr^2ψ(n,m,r)^2 と計算される.
!!$      !
!!$      ! * 全てのエネルギースペクトル成分の和をZ積分したもの(r^2の重み無し)
!!$      !    が球殻内での全エネルギーに等しい.
!!$      !
!!$      real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: wt_TORPOT
!!$      !(in) トロイダルポテンシャル
!!$
!!$      real(8), dimension(0:nm,0:km) :: nz_ToroidalEnergySpectrum_wt
!!$      !(out) エネルギースペクトルトロイダル成分
!!$
!!$      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA   ! 作業領域
!!$      integer :: n, m
!!$
!!$      wz_DATA = wz_wt(wt_TORPOT)
!!$      do n=0,nm
!!$         nz_ToroidalEnergySpectrum_wt(n,:) &
!!$              = 0.5 * n*(n+1)* (4*pi) * z_Z**2 &
!!$                * sum(wz_Data(l_nm(n,(/(m,m=-n,n)/)),:)**2,1)
!!$      enddo
!!$
!!$    end function nz_ToroidalEnergySpectrum_wt

    function zef_PoloidalEnergySpectrum_tef(tef_POLPOT)
      !
      ! ポロイダルポテンシャルから, ポロイダルエネルギーの
      ! 水平 X 波数 l, Y 波数 m の各成分を計算する.
      !
      !  * X 波数 l, Y 波数 m のポロイダルポテンシャルのスペクトル成分
      !    φ(l,m,z)から X 波数 l, Y 波数 m 成分のポロイダルエネルギー
      !    スペクトルは
      !
      !      (1/2)(l**2+m**2){[dφ(n,m,z)/dz]^2 + (l**2+m**2)φ(n,m,z)^2}
      !
      !    と計算される.
      !
      !  * 全てのエネルギースペクトル成分の和をZ積分したもの
      !    が水平平均エネルギーに等しい.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)), intent(in) :: tef_POLPOT
      !(in) ポロイダルポテンシャル

      real(8), dimension(0:km,-nm:nm,2*lc(ip)) :: zef_PoloidalEnergySpectrum_tef
      !(out) エネルギースペクトルポロイダル成分


      real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_Data   ! 作業領域
      real(8), dimension(0:km,-mm:mm,2*lc(ip)) :: zef_DData  ! 作業領域
      integer :: l, m

      zef_PoloidalEnergySpectrum_tef = 0.0D0
      zef_Data = zef_tef(tef_POLPOT)
      zef_DData = zef_tef(tef_DZ_tef(tef_POLPOT))

      do l=-lm,lm
         do m=-mm,mm
            if ( lf_l(l) > 0 ) then
               zef_PoloidalEnergySpectrum_tef(:,m,lf_l(l)) =             &
                    + 0.5* ((2*pi*l/xl)**2+(2*pi*m/yl)**2)               &
                    *(   zef_DData(:,m,lf_l(l))**2 + zef_DData(:,-m,lf_l(-l))**2     &
                       + ((2*pi*l/xl)**2+(2*pi*m/yl)**2)                 &
                           *( zef_Data(:,m,lf_l(l))**2 + zef_Data(:,-m,lf_l(-l))**2))
            endif
         enddo
      enddo

    end function zef_PoloidalEnergySpectrum_tef

!!$    function nz_PoloidalEnergySpectrum_wt(wt_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}
!!$      !
!!$      !    と計算される.
!!$      !
!!$      !  * 全ての全波数に対してのエネルギースペクトル成分の和をZ積分したもの
!!$      !    (r^2の重み無し)が球殻内での全エネルギーに等しい.
!!$      !
!!$      real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: wt_POLPOT
!!$      !(in) ポロイダルポテンシャル
!!$
!!$      real(8), dimension(0:nm,0:km) :: nz_PoloidalEnergySpectrum_wt
!!$      !(out) エネルギースペクトルポロイダル成分
!!$
!!$      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA1   ! 作業領域
!!$      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA2   ! 作業領域
!!$      integer :: n, m
!!$
!!$      wz_Data1 = wz_wt(wt_POLPOT)
!!$      wz_Data2 = wz_Z*wz_wt(wt_DZ_wt(wt_POLPOT)) &    ! d(rφ)/dr
!!$               + wz_wt(wt_POLPOT)                         ! = rdφ/dr+φ
!!$
!!$      do n=0,nm
!!$         nz_PoloidalEnergySpectrum_wt(n,:) = &
!!$              + 0.5* n*(n+1)* (4*pi) &
!!$              *( sum(wz_Data2(l_nm(n,(/(m,m=-n,n)/)),:)**2,1)  &
!!$                + n*(n+1)*sum(wz_Data1(l_nm(n,(/(m,m=-n,n)/)),:)**2,1) )
!!$      enddo
!!$
!!$    end function nz_PoloidalEnergySpectrum_wt

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

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm)         :: e2t

      real(8), dimension(2*lc(ip)*(2*mm+1),2)            :: e22_values
      character(len=2)                                   :: Bcond

      e2t = e2a_aef(tef)

      if (present(values)) then
         e22_values = e2a_aef(values)
      endif

      if (.not. present(cond)) then
         Bcond='DD'
      else
         Bcond = cond
      endif

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

      tef = aef_e2a(e2t)

    end subroutine tef_BoundariesTau

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm)         :: e2t

      real(8), dimension(2*lc(ip)*(2*mm+1),2)            :: e22_values
      character(len=2)                                   :: Bcond

      e2t = e2a_aef(tef)

      if (present(values)) then
         e22_values = e2a_aef(values)
      endif

      if (.not. present(cond)) then
         BCond = 'DD'
      else
         BCond = cond
      endif

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

      tef = aef_e2a(e2t)

    end subroutine tef_BoundariesGrid

    subroutine tef_TorBoundariesTau(tef_TOR,cond,new)

      ! トロイダル速度ポテンシャルに対して境界条件を適用する.
      ! Chebyshev 空間での境界条件適用
      !
      ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を
      ! とっている(タウ法). トロイダルポテンシャルの境界条件は
      !
      !     ψ = 0 at boundaries           (粘着条件)
      !     ∂ψ/∂z = 0 at boundaries    (応力なし条件)
      !
      ! であるから tef_Boundaries で対応可能だが, 将来のため別途作成しておく.
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout)   :: tef_TOR
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work

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

      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','tef_TorBoundariesTau','B.C. not supported')
      end select

      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(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))

         do n=0,nm
            e2t_work= 0.0D0 ; e2t_work(:,n)= 1.0D0
            alu(:,:,n) = e2t_work

            ! 粘着条件
            ! 力学的条件粘着条件
            if ( rigid1 ) then
               e2z_work = az_at(e2t_work)
            else
               e2z_work = az_at(at_Dz_at(e2t_work))
            endif
            alu(:,nm-1,n) = e2z_work(:,0)

            ! 力学的条件粘着条件
            if ( rigid2 ) then
               e2z_work = az_at(e2t_work)
            else
               e2z_work = az_at(at_Dz_at(e2t_work))
            endif
            alu(:,nm,n)   = e2z_work(:,km)
         enddo

         call ludecomp(alu,kp)
         call MessageNotify('M','tef_TorBoundariesTau',&
                            'Matrix to apply  b.c. newly produced.')
      endif

      e2t_work = e2a_aef(tef_Tor)

      e2t_work(:,nm-1) = 0.0D0
      e2t_work(:,nm)   = 0.0D0

      tef_Tor = aef_e2a(lusolve(alu,kp,e2t_work))

    end subroutine tef_TorBoundariesTau

    subroutine tef_TorBoundariesGrid(tef_TOR,cond,new)
      !
      ! トロイダル速度ポテンシャルに対して境界条件を適用する.
      ! 鉛直実空間での境界条件適用.
      !
      ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
      ! 条件を課している(選点法). このルーチンを用いるためには
      ! tef_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を
      ! 等しくしておく必要がある.
      !
      ! トロイダル速度ポテンシャルの境界条件は
      !
      !     ψ = 0 at boundaries           (粘着条件)
      !     ∂ψ/∂z = 0 at boundaries    (応力なし条件)
      !
      ! であるので tef_Boundaries で対応可能だが, 将来のため別途作成しておく
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(8), dimension(0:nm,-mm:mm,2*lc(ip)),intent(inout)   :: tef_TOR
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work

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

      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','tef_TorBoundariesGrid','B.C. not supported')
      end select

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

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

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

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))

         do n=0,nm
            e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0
            e2z_work = az_at(e2t_work)

            alu(:,:,n) = e2z_work          ! 内部領域は値そのまま.


            ! 粘着条件
            ! 力学的条件粘着条件
            if ( rigid1 ) then
               e2z_work = az_at(e2t_work)
            else
               e2z_work=az_at(at_Dz_at(e2t_work))
            endif
            alu(:,0,n) = e2z_work(:,0)

            ! 力学的条件粘着条件
            if ( rigid2 ) then
               e2z_work = az_at(e2t_work)
            else
               e2z_work=az_at(at_Dz_at(e2t_work))
            endif
            alu(:,km,n)   = e2z_work(:,km)

         enddo
         call ludecomp(alu,kp)
         call MessageNotify('M','TorBoundariesGrid',&
                            'Matrix to apply  b.c. newly produced.')
      endif

      e2z_work = az_at(e2a_aef(tef_Tor))
      e2z_work(:,0)  = 0.0D0
      e2z_work(:,km) = 0.0D0
      tef_TOR = aef_e2a(lusolve(alu,kp,e2z_work))

    end subroutine tef_TorBoundariesGrid

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

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))      :: tef_LaplaPol2PolTau_tef
              !(out) 出力ポロイダルポテンシャル分布

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm)  :: e2t_work
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km)  :: e2z_work
      logical                                     :: rigid1 = .true.
      logical                                     :: rigid2 = .true.

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

      if ( present(cond) ) then
         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','tef_LaplaPol2PolTau_tef','B.C. not supported')
         end select
      else
         rigid1 = .TRUE.  ; rigid2 = .TRUE.
      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(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))

         alu = 0.0d0
         do n=0,nm-2
            e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0

            ! 各水平波数に関して独立の式
            alu(:,:,n) = e2a_aef(tef_lapla_tef(aef_e2a(e2t_work)))
         enddo

         do n=0,nm
            e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0
            e2z_work = az_at(e2t_work)

            ! 運動学的条件. 流線は境界で一定
            alu(:,nm,n)   = e2z_work(:,0)
            alu(:,nm-1,n) = e2z_work(:,km)

            ! 力学的条件粘着条件
            e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0
            if ( rigid1 ) then
               e2z_work=az_at(at_Dz_at(e2t_work))
            else
               e2z_work=az_at(at_Dz_at(at_Dz_at(e2t_work)))
            endif
            alu(:,nm-2,n) = e2z_work(:,0)

            ! 力学的条件粘着条件
            e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0
            if ( rigid2 ) then
               e2z_work=az_at(at_Dz_at(e2t_work))
            else
               e2z_work=az_at(at_Dz_at(at_Dz_at(e2t_work)))
            endif
            alu(:,nm-3,n) = e2z_work(:,km)
         enddo

         call ludecomp(alu,kp)

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

      e2t_work         = e2a_aef(tef)
      e2t_work(:,nm)   = 0.0D0               ! 運動学的条件
      e2t_work(:,nm-1) = 0.0D0               ! 運動学的条件
      e2t_work(:,nm-2) = 0.0D0               ! 力学的条件
      e2t_work(:,nm-3) = 0.0D0               ! 力学的条件

      tef_LaplaPol2PolTau_tef = aef_e2a(lusolve(alu,kp,e2t_work))

    end function tef_LaplaPol2PolTau_tef

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

      real(8), dimension(0:km,-mm:mm,2*lc(ip))             :: zef_LaplaPol2Pol_zef
              !(out) 出力ポロイダルポテンシャル分布

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:km)  :: e2z_work
      logical                                     :: rigid1 = .true.
      logical                                     :: rigid2 = .true.

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

      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','zef_laplapol2pol_zef','B.C. not supported')
      end select

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

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

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

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:km),kp(2*lc(ip)*(2*mm+1),0:km))

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

            ! 各水平波数に関して独立の式
            alu(:,:,k) &
                 = e2a_aef(zef_tef(tef_lapla_tef(tef_zef(aef_e2a(e2z_work)))))
         enddo

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

            ! 運動学的条件. 流線は境界で一定
            alu(:,0,k)   = e2z_work(:,0)
            alu(:,km,k)  = e2z_work(:,km)

            ! 力学的条件粘着条件
            e2z_work = 0.0D0 ; e2z_work(:,k) = 1.0D0
            if ( rigid1 ) then
               e2z_work=az_at(at_Dz_at(at_az(e2z_work)))
            else
               e2z_work=az_at(at_Dz_at(at_Dz_at(at_az(e2z_work))))
            endif
            alu(:,1,k) = e2z_work(:,0)

            ! 力学的条件粘着条件
            e2z_work = 0.0D0 ; e2z_work(:,k) = 1.0D0
            if ( rigid2 ) then
               e2z_work=az_at(at_Dz_at(at_az(e2z_work)))
            else
               e2z_work=az_at(at_Dz_at(at_Dz_at(at_az(e2z_work))))
            endif
            alu(:,km-1,k) = e2z_work(:,km)
         enddo

         call ludecomp(alu,kp)

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

      e2z_work        = e2a_aef(zef)
      e2z_work(:,1)    = 0.0D0               ! 力学的条件
      e2z_work(:,km-1) = 0.0D0               ! 力学的条件
      e2z_work(:,0)    = 0.0D0               ! 運動学的条件
      e2z_work(:,km)   = 0.0D0               ! 運動学的条件

      e2z_work = lusolve(alu,kp,e2z_work)

      zef_laplapol2pol_zef = aef_e2a(e2z_work)

    end function zef_LaplaPol2Pol_zef

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

      real(8), dimension(0:nm,-mm:mm,2*lc(ip))      :: tef_LaplaPol2PolGrid_tef
              !(out) 出力ポロイダルポテンシャル分布

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm)  :: e2t_work
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km)  :: e2z_work
      logical                                     :: rigid1 = .true.
      logical                                     :: rigid2 = .true.

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

      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','tef_LaplaPol2PolGrid_tef','B.C. not supported')
      end select

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

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

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

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))

         do n=0,nm
            e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0

            ! 各水平波数に関して独立の式
            alu(:,:,n) = e2a_aef(zef_tef(tef_lapla_tef(aef_e2a(e2t_work))))
         enddo

         do n=0,nm
            e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0
            e2z_work = az_at(e2t_work)

            ! 運動学的条件. 流線は境界で一定
            alu(:,0,n)   = e2z_work(:,0)
            alu(:,km,n)  = e2z_work(:,km)

            ! 力学的条件粘着条件
            e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0
            if ( rigid1 ) then
               e2z_work=az_at(at_Dz_at(e2t_work))
            else
               e2z_work=az_at(at_Dz_at(at_Dz_at(e2t_work)))
            endif
            alu(:,1,n) = e2z_work(:,0)

            ! 力学的条件粘着条件
            e2t_work = 0.0D0 ; e2t_work(:,n) = 1.0D0
            if ( rigid2 ) then
               e2z_work=az_at(at_Dz_at(e2t_work))
            else
               e2z_work=az_at(at_Dz_at(at_Dz_at(e2t_work)))
            endif
            alu(:,km-1,n) = e2z_work(:,km)
         enddo

         call ludecomp(alu,kp)

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

      e2z_work         = az_at(e2a_aef(tef))
      e2z_work(:,1)    = 0.0D0               ! 力学的条件
      e2z_work(:,km-1) = 0.0D0               ! 力学的条件
      e2z_work(:,0)    = 0.0D0               ! 運動学的条件
      e2z_work(:,km)   = 0.0D0               ! 運動学的条件

      tef_LaplaPol2PolGrid_tef = aef_e2a(lusolve(alu,kp,e2z_work))

    end function tef_LaplaPol2PolGrid_tef

    subroutine tef_TormagBoundariesTau(tef_TOR,new)

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer  :: n
      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)
         allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))

         do n=0,nm
            e2t_work= 0.0D0 ; e2t_work(:,n)= 1.0D0
            alu(:,:,n) = e2t_work

            ! 非電気伝導体
            e2z_work = az_at(e2t_work)
            alu(:,nm-1,n) = e2z_work(:,0)
            alu(:,nm,n)   = e2z_work(:,km)

         enddo

         call ludecomp(alu,kp)
         call MessageNotify('M','tef_TormagBoundariesTau',&
                            'Matrix to apply  b.c. newly produced.')
      endif

      e2t_work = e2a_aef(tef_Tor)

      e2t_work(:,nm-1) = 0.0D0
      e2t_work(:,nm)   = 0.0D0


      tef_Tor = aef_e2a(lusolve(alu,kp,e2t_work))

    end subroutine tef_TormagBoundariesTau

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work

      logical :: first = .true.
      logical :: new_matrix = .false.
      integer  :: n
      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 ( nm /= km ) then
            call MessageNotify('E','tef_TorMagBoundariesGrid', &
             'Chebyshev truncation and number of grid points should be same.')
         endif

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) ) deallocate(kp)
         allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))

         do n=0,nm
            e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0
            e2z_work = az_at(e2t_work)

            alu(:,:,n) = e2z_work          ! 内部領域は値そのまま.

            ! 非電気伝導体
            alu(:,0,n)  = e2z_work(:,0)
            alu(:,km,n) = e2z_work(:,km)

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

      e2z_work = az_at(e2a_aef(tef_Tor))
      e2z_work(:,0)  = 0.0D0
      e2z_work(:,km) = 0.0D0
      tef_TOR = aef_e2a(lusolve(alu,kp,e2z_work))

    end subroutine tef_TormagBoundariesGrid

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_kh

      logical :: first = .true.
      logical :: new_matrix = .false.
      real(8) :: rl, rm
      integer :: l, m, n, e2index
      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)

         allocate(alu(2*lc(ip)*(2*mm+1),0:nm,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))

         do n=0,nm
            e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0

            alu(:,:,n) = e2t_work(:,:)    ! 内部領域は同じ
         enddo

         ! 非電気伝導体
         do m=-mm,mm
            do l=-lm,lm
               if ( lf_l(l) > 0 ) then
                  e2index = (2*mm+1)*(lf_l(l)-1) + (m+mm+1)
                  rl = 2*PI/xl*l ; rm = 2*PI/yl*m
                  e2t_kh(e2index,:) = sqrt(rl**2+rm**2)
               endif
            enddo
         enddo
         do n=0,nm
            e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0

            e2z_work = az_at(at_Dz_at(e2t_work) + e2t_kh * e2t_work)
            alu(:,nm-1,n) = e2z_work(:,0)

            e2z_work = az_at(at_Dz_at(e2t_work) - e2t_kh * e2t_work)
            alu(:,nm,n)   = e2z_work(:,km)
         enddo

         call ludecomp(alu,kp)

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

      e2t_work = e2a_aef(tef_POL)
      e2t_work(:,nm-1) = 0.0D0
      e2t_work(:,nm)   = 0.0D0
      tef_POL = aef_e2a(lusolve(alu,kp,e2t_work))

    end subroutine tef_PolmagBoundariesTau

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

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

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

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension(2*lc(ip)*(2*mm+1),0:km) :: e2z_work

      real(8), dimension(2*lc(ip)*(2*mm+1),0:nm) :: e2t_kh

      logical :: first = .true.
      logical :: new_matrix = .false.
      real(8) :: rl, rm
      integer  :: l, m, n, e2index
      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 ( nm /= km ) then
            call MessageNotify('E','tef_PolMagBoundariesGrid', &
             'Chebyshev truncation and number of grid points should be same.')
         endif

         if ( allocated(alu) ) deallocate(alu)
         if ( allocated(kp) )  deallocate(kp)
         allocate(alu(2*lc(ip)*(2*mm+1),0:km,0:nm),kp(2*lc(ip)*(2*mm+1),0:nm))

         do n=0,nm
            e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0
            e2z_work = az_at(e2t_work)

            alu(:,:,n) = e2z_work    ! 内部領域は同じ
         enddo

         ! 非電気伝導体
         do m=-mm,mm
            do l=-lm,lm
               if ( lf_l(l) > 0 ) then
                  e2index = (2*mm+1)*(lf_l(l)-1) + (m+mm+1)
                  rl = 2*PI/xl*l ; rm = 2*PI/yl*m
                  e2t_kh(e2index,:) = sqrt(rl**2+rm**2)
               endif
            enddo
         enddo
         do n=0,nm
            e2t_work = 0.0D0 ; e2t_work(:,n)=1.0D0

            e2z_work = az_at(at_Dz_at(e2t_work) + e2t_kh * e2t_work)
            alu(:,0,n) = e2z_work(:,0)

            e2z_work = az_at(at_Dz_at(e2t_work) - e2t_kh * e2t_work)
            alu(:,km,n)   = e2z_work(:,km)
         enddo

         call ludecomp(alu,kp)

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

      e2z_work = az_at(e2a_aef(tef_POL))
      e2z_work(:,0)  = 0.0D0
      e2z_work(:,km) = 0.0D0
      tef_POL = aef_e2a(lusolve(alu,kp,e2z_work))

    end subroutine tef_PolmagBoundariesGrid

end module tee_mpi_module
