!--
!----------------------------------------------------------------------
! Copyright(c) 2009-2016 SPMDODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  tee_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/20  竹広真一  Polmagbc ルーチンバグ fix
!
!凡例
!      データ種類と index
!        x : 水平(X)        y : 水平(Y)        z : 鉛直
!        e : フーリエ変換スペクトル
!        l : フーリエ関数スペクトル(X 方向波数)
!        m : フーリエ関数スペクトル(Y 方向波数)
!        t : チェビシェフ関数スペクトル
!        a : 任意の次元
!
!        zyx : 3 次元格子点データ
!        yx  : 水平 2 次元格子点データ
!        zy  : 鉛直 2 次元格子点データ
!        zx  : 鉛直 2 次元格子点データ
!
!        zee : 水平スペクトル鉛直格子点データ
!        tee : スペクトルデータ鉛直チェビシェフデータ
!
!++
module tee_module
  !
  != tee_module
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  !    spml/tee_module モジュールは平行平板間での 3 次元流体運動を
  !    スペクトル法によって数値計算するための Fortran90 関数を提供する
  !    ものである.
  !
  !    水平方向にフーリエ数変換および上下の境界壁を扱うための
  !    チェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな
  !    関数を提供する.
  !
  !    内部で ee_module, at_module を用いている. 最下部ではフーリエ
  !    およびチェビシェフ変換のエンジンとして ISPACK の Fortran77
  !    サブルーチンを用いている.
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (tee_, zyx_, zee_, ee_, yx_, x_, y_, z_, a_) は,
  !   返す値の形を示している.
  !   tee_  :: スペクトルデータ(2 重フーリエ・チェビシェフ変換)
  !   zyx_ :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・)
  !   zee_ :: 水平スペクトル, 鉛直格子点データ
  !   e2a_ :: 1 次元化した水平スペクトル, 任意座標データ
  !   aee_ :: 任意座標データ, 水平スペクトルデータ
  !
  ! * 関数名の間の文字列(Dx, Dy, Dz, Lapla,..)
  !   は, その関数の作用を表している.
  !
  ! * 関数名の最後 (tee_, zyx_, zee_, ee_, yx_, x_, y_, z_, a_) は, 入力変数の
  !   形がスペクトルデータおよび格子点データであることを示している.
  !   _tee :: スペクトルデータ(2 重フーリエ・チェビシェフ変換)
  !   _zyx :: 3 次元格子点データ(水平 2 次元鉛直 1 次元・)
  !   _zee :: 水平スペクトル, 鉛直格子点データ
  !   _e2a :: 1 次元化した水平スペクトル, 任意座標データ
  !   _aee :: 任意座標データ, 水平スペクトルデータ
  !
  !=== 各データの種類の説明
  !
  ! * zyx : 3 次元格子点データ(鉛直, 水平 2 次元)
  !   * 変数の種類と次元は real(8), dimension(0:km,0:jm-1,0:im-1).
  !   * im, jm, km はそれぞれ水平 X, Y, 鉛直 Z 座標の格子点数であり,
  !     サブルーチン tee_Initial にてあらかじめ設定しておく.
  !
  ! * tee : スペクトルデータ
  !   * 変数の種類と次元は real(8), dimension(0:nm,-mm:mm,-lm:lm).
  !   * lm, mm は X,Y 方向最大波数, nm はチェビシェフ多項式の最大次数
  !     であり, サブルーチン tee_Initial にてあらかじめ設定しておく.
  !
  ! * zee : 水平スペクトル, 鉛直格子点データ.
  !   * 変数の種類と次元は real(8), dimension(0:km,-mm:mm,-lm:lm).
  !
  ! * e2a : 1 次元化した水平スペクトル, 任意座標データ
  !   * 変数の種類と次元は real(8), dimension((2*mm+1)*(2*lm+1),:)
  !
  ! * aee :   任意座標データ, 水平スペクトルデータ
  !   * 変数の種類と次元は real(8), dimension(:,-mm:mm,-lm:lm).
  !
  ! * tee_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * zyx_ で始まる関数が返す値は 3 次元格子点データに同じ.
  !
  ! * zee_ で始まる関数が返す値は水平スペクトル, 動径格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! tee_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  !
  !==== 座標変数
  !
  ! x_X, y_Y, z_Z                :: 格子点座標(水平 X,Y, 鉛直 Z 座標)を
  !                                 格納した1 次元配列
  ! x_X_Weight, y_X_Weight, z_Z_Weight :: 重み座標を格納した 1 次元配列
  ! zyx_X, zyx_Y, zyx_Z          :: 格子点データの水平鉛直座標(X,Y,Z)
  !                                 (格子点データ型 3 次元配列)
  ! yx_X, yx_Y                   :: 格子点データの水平座標(X,Y)
  ! zy_Z, zy_Y                   :: 格子点データの鉛直水平座標(Z,Y)
  ! zx_Z, zx_X                   :: 格子点データの鉛直水平座標(Z,X)
  !                                 (格子点データ型 2 次元配列)
  !
  !==== 基本変換
  !
  ! zyx_tee, tee_zyx :: スペクトルデータと 3 次元格子データの間の変換
  !                     (2 重フーリエ, チェビシェフ変換)
  ! zyx_zee, zee_zyx :: 3 次元格子データと水平スペクトル・鉛直格子データとの間
  !                     の変換 (2 重フーリエ変換)
  ! zee_tee, tee_zee :: スペクトルデータと水平スペクトル・鉛直格子データとの間
  !                     の変換 (チェビシェフ変換)
  ! ee_yx, yx_ee     :: スペクトルデータと 2 次元水平格子データの間の変換
  !                     (2 重フーリエ変換)
  ! az_at, at_az     :: 同時に複数個行う (チェビシェフ変換)格子データと
  !                     チェビシェフデータの間の変換を
  ! e2a_aee, aee_e2a :: 水平スペクトル軸を 1 次元化し転置, 転置 ２次元化する.
  !
  !==== 微分
  !
  ! tee_Dx_tee          :: スペクトルデータに動径微分∂/∂x を作用させる
  ! tee_Dy_tee          :: スペクトルデータに動径微分∂/∂y を作用させる
  ! tee_Dz_tee          :: スペクトルデータに動径微分∂/∂z を作用させる
  !
  ! tee_Lapla_tee       :: スペクトルデータにラプラシアンを作用させる
  ! tee_LaplaH_tee      :: スペクトルデータに水平ラプラシアンを作用させる
  ! tee_LaplaHInv_tee   :: スペクトルデータに逆水平ラプラシアンを作用させる
  !
  ! tee_Div_zyx_zyx_zyx :: ベクトル成分である 3 つの格子データに
  !                        発散を作用させる
  !
  !==== トロイダルポロイダル計算用微分
  !
  ! tee_ZRot_zyx_zyx         :: ベクトル v の渦度と動径ベクトル r の内積
  !                             z・(▽×v) を計算する
  ! tee_ZRotRot_zyx_zyx_zyx  :: ベクトルの v の z・(▽×▽×v) を計算する
  ! tee_Potential2Vector     :: トロイダルポロイダルポテンシャルから
  !                             ベクトル場を計算する
  ! tee_Potential2Rotation   :: トロイダルポロイダルポテンシャルで表される
  !                             非発散ベクトル場の回転の各成分を計算する
  !
  !==== ポロイダル/トロイダルモデル用スペクトル解析
  !
  ! zee_ToroidalEnergySpectrum_tee, zk_ToroidalEnergySpectrum_tee   ::
  !     トロイダルポテンシャルからエネルギーのフーリエ各成分を計算する
  ! zee_PoloidalEnergySpectrum_tee, zk_PoloidalEnergySpectrum_tee   ::
  !     ポロイダルポテンシャルからエネルギーのフーリエ各成分を計算する
  !
  !==== 境界値問題
  !
  ! tee_BoundariesTau, tee_BoundariesGrid, tee_Boundaries                ::
  !     ディリクレ, ノイマン境界条件を適用する(タウ法, 選点法)
  !
  ! tee_TorBoundariesTau, tee_TorBoundariesGrid, tee_TorBoundaries       ::
  !     速度トロイダルポテンシャルの境界条件を適用する(タウ法,選点法)
  !
  ! tee_LaplaPol2PolTau_tee, zee_LaplaPol2Pol_zee, tee_LaplaPol2PolGrid_tee ::
  !     速度ポロイダルポテンシャルΦを▽^2Φから求める
  !     (入出力がそれぞれチェビシェフ格子点,チェビシェフ係数)
  !
  ! tee_TorMagBoundariesTau, tee_TorMagBoundariesGrid, tee_TorMagBoundaries ::
  !     磁場トロイダルポテンシャルの境界条件を適用する(タウ法, 選点法)
  !
  ! tee_PolMagBoundariesTau, tee_PolMagBoundariesGrid, tee_PolMagBoundaries ::
  !     磁場トロイダルポテンシャル境界の境界条件を適用する(タウ法, 選点法)
  !
  !==== 積分・平均(3 次元データ)
  !
  ! IntZYX_zyx, AvrZYX_zyx     :: 3 次元格子点データの全領域積分および平均
  ! z_IntYX_zyx, z_AvrYX_zyx   :: 3 次元格子点データの水平積分および平均
  ! y_IntZX_zyx, y_AvrZX_zyx   :: 3 次元格子点データのZX積分および平均
  ! z_IntZY_zyx, z_AvrZY_zyx   :: 3 次元格子点データのZY積分および平均
  ! zy_IntX_zyx, zy_AvrX_zyx   :: 3 次元格子点データの水平X方向積分および平均
  ! zx_IntY_zyx, zx_AvrY_zyx   :: 3 次元格子点データの水平Y方向積分および平均
  ! zx_IntZ_zyx, zx_AvrZ_zyx   :: 3 次元格子点データの鉛直方向積分および平均
  !
  !==== 積分・平均(2 次元データ)
  !
  ! IntYX_yx, AvrYX_yx :: 2 次元格子点データの水平積分および平均
  ! IntZX_zx, AvrZX_zx :: 2 次元(ZX)格子点データのZX積分および平均
  ! IntZY_zy, AvrZY_zy :: 2 次元(ZY)格子点データのZY積分および平均
  ! y_IntX_yx, y_AvrX_yx   :: 水平 2 次元格子点データのX方向積分および平均
  ! x_IntY_yx, x_AvrY_yx   :: 水平2 次元格子点データのY方向積分および平均
  ! z_IntX_zx, z_AvrX_zx   :: 2 次元(ZX)格子点データのX方向積分および平均
  ! x_IntZ_zx, x_AvrZ_zx   :: 2 次元(ZX)格子点データのZ方向積分および平均
  ! z_IntY_zy, z_AvrY_zy   :: 2 次元(ZY)格子点データのY方向積分および平均
  ! y_IntZ_zy, y_AvrZ_zy   :: 2 次元(ZY)格子点データのZ方向積分および平均
  !
  !==== 積分・平均(1 次元データ)
  !
  ! IntX_x, AvrX_x  :: 1 次元(X)格子点データのX方向積分および平均
  ! IntY_y, AvrY_y  :: 1 次元(Y)格子点データのY方向積分および平均
  ! IntZ_z, AvrZ_z  :: 1 次元(Z)格子点データのZ方向積分および平均
  !
  !==== 補間計算
  !
  ! Interpolate_tee :: スペクトルデータから任意の点の値を補間する.
  !
  use dc_message
  use lumatrix
  use ee_module
  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_Initial

  public x_X, x_X_Weight
  public y_Y, y_Y_Weight
  public z_Z, z_Z_Weight
  public yx_X, yx_Y, zy_Z, zy_Y, zx_X, zx_Z
  public zyx_X, zyx_Y, zyx_Z
  public zee_Z
  public tee_VMiss

  public ee_yx, yx_ee
  public at_Dz_at, t_Dz_t, az_at, at_az, z_t, t_z
  public zyx_tee, tee_zyx, zyx_zee, zee_zyx, zee_tee, tee_zee
  public e2z_e2t, e2t_e2z 
  public aee_e2a, e2a_aee, ee_e2, e2_ee
  public tee_Dx_tee, tee_Dy_tee, tee_Dz_tee
  public tee_Lapla_tee, tee_LaplaH_tee, tee_LaplaHInv_tee
  public tee_Div_zyx_zyx_zyx

  public zy_IntX_zyx, zx_IntY_zyx, yx_IntZ_zyx
  public x_IntZY_zyx, y_IntZX_zyx, z_IntYX_zyx
  public IntZYX_zyx

  public x_IntY_yx, y_IntX_yx, IntYX_yx
  public z_IntY_zy, y_IntZ_zy, IntZY_zy
  public z_IntX_zx, x_IntZ_zx, IntZX_zx
  public IntX_x, IntY_y, IntZ_z

  public zy_AvrX_zyx, zx_AvrY_zyx, yx_AvrZ_zyx
  public x_AvrZY_zyx, y_AvrZX_zyx, z_AvrYX_zyx
  public AvrZYX_zyx

  public x_AvrY_yx, y_AvrX_yx, AvrYX_yx
  public z_AvrY_zy, y_AvrZ_zy, AvrZY_zy
  public z_AvrX_zx, x_AvrZ_zx, AvrZX_zx
  public AvrX_x, AvrY_y, AvrZ_z

  public tee_ZRot_zyx_zyx, tee_ZRotRot_zyx_zyx_zyx
  public tee_Potential2vector, tee_Potential2Rotation

  public Interpolate_tee

  public zee_ToroidalEnergySpectrum_tee ! nz_ToroidalEnergySpectrum_wt
  public zee_PoloidalEnergySpectrum_tee ! nz_PoloidalEnergySpectrum_wt

  public tee_Boundaries, zee_LaplaPol2Pol_zee
  public tee_TorBoundaries, tee_LaplaPol2Pol_tee
  public tee_TormagBoundaries, tee_PolmagBoundaries

  public tee_BoundariesTau, tee_LaplaPol2PolTau_tee, tee_TorBoundariesTau
  public tee_TormagBoundariesTau, tee_PolmagBoundariesTau

  public tee_BoundariesGrid, tee_LaplaPol2PolGrid_tee, tee_TorBoundariesGrid
  public tee_TormagBoundariesGrid, tee_PolmagBoundariesGrid

  public at_Boundaries, t_Boundaries

  interface tee_Boundaries
     module procedure tee_BoundariesTau
  end interface

  interface tee_TorBoundaries
     module procedure tee_TorBoundariesTau
  end interface

  interface tee_LaplaPol2Pol_tee
     module procedure tee_LaplaPol2PolTau_tee
  end interface

  interface tee_TorMagBoundaries
     module procedure tee_TorMagBoundariesTau
  end interface

  interface tee_PolMagBoundaries
     module procedure tee_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 :: zyx_X, zyx_Y, zyx_Z    ! 座標
  real(8), dimension(:,:),   allocatable :: zy_Z, zy_Y             ! 座標
  real(8), dimension(:,:),   allocatable :: zx_Z, zx_X             ! 座標
  real(8), dimension(:,:,:), allocatable :: zee_Z                  ! 座標

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

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

  contains
  !--------------- 初期化 -----------------
   subroutine tee_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            :: id

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

     call ee_Initial(im,jm,lm,mm,xmin,xmax,ymin,ymax,id)

     call at_Initial(km,nm,zmin,zmax)

     allocate(zyx_X(0:km,0:jm-1,0:im-1))
     allocate(zyx_Y(0:km,0:jm-1,0:im-1))
     allocate(zyx_Z(0:km,0:jm-1,0:im-1))

     allocate(zy_Z(0:km,0:jm-1))
     allocate(zy_Y(0:km,0:jm-1))
     allocate(zx_Z(0:km,0:im-1))
     allocate(zx_X(0:km,0:im-1))

     allocate(zee_Z(0:km,-mm:mm,-lm:lm))

     zyx_X = spread(yx_X,1,km+1)
     zyx_Y = spread(yx_Y,1,km+1)
     zyx_Z = spread(spread(z_Z,2,jm),3,im)

     zy_Z = spread(z_Z,2,jm)
     zy_Y = spread(y_Y,1,km+1)
     zx_Z = spread(z_Z,2,im)
     zx_X = spread(x_X,1,km+1)

     zee_Z = spread(spread(z_Z,2,2*mm+1),3,2*lm+1)

     call MessageNotify('M','tee_initial', &
          'tee_module (2016/09/20) is initialized')

   end subroutine tee_Initial

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

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

      zyx_tee = zyx_zee(zee_tee(tee))

    end function zyx_tee

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

      tee_zyx = tee_zee(zee_zyx(zyx))

    end function tee_zyx

    function zyx_zee(zee)
      !
      ! 水平スペクトル・鉛直格子点データから 3 次元格子点データへ(逆)変換する.
      !
      real(8), dimension(0:km,-mm:mm,-lm:lm), intent(in) :: zee
      !(in) 2 次元水平スペクトル・鉛直格子点データ
      real(8), dimension(0:km,0:jm-1,0:im-1)             :: zyx_zee
      !(out) 3 次元格子点データ
      real(8), dimension(-mm:mm,-lm:lm)                  :: ee
      real(8), dimension(0:jm-1,0:im-1)                  :: yx

      integer :: k

      do k = 0, km
         !zyx_zee(k,:,:) = yx_ee(zee(k,:,:))
         !
         ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
         !
         ee = zee(k,-mm:mm,-lm:lm)
         yx = yx_ee(ee)
         zyx_zee(k,:,:) = yx
      enddo

    end function zyx_zee

    function zee_zyx(zyx)
      !
      ! 3 次元格子データから水平スペクトル・鉛直格子点データへ(正)変換する.
      !
      real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
      !(in) 3 次元格子点データ
      real(8), dimension(0:km,-mm:mm,-lm:lm)             :: zee_zyx
      !(out) 2 次元スペクトル・鉛直格子点データ
      real(8), dimension(-mm:mm,-lm:lm)                  :: ee
      real(8), dimension(0:jm-1,0:im-1)                  :: yx

      integer :: k

      do k = 0, km
         !zee_zyx(k,:,:) = ee_yx(zyx(k,:,:))
         !
         ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
         !
         yx = zyx(k,0:jm-1,0:im-1)
         ee = ee_yx(yx)
         zee_zyx(k,:,:) = ee
      enddo

    end function zee_zyx

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

      zee_tee = aee_e2a(e2z_e2t(e2a_aee(tee)))

    end function zee_tee

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

      tee_zee = aee_e2a(e2t_e2z(e2a_aee(zee)))

    end function tee_zee

    function e2z_e2t(e2t)
      !
      ! 鉛直スペクトルを格子点に変換する
      !
      real(8), dimension((2*lm+1)*(2*mm+1),0:nm)  :: e2t
      !(in)  1 次元化された水平スペクトル・鉛直スペクトル
      real(8), dimension((2*lm+1)*(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*lm+1)*(2*mm+1),0:km)  :: e2z
      !(in)  1 次元化された水平スペクトル・鉛直格子座標
      real(8), dimension((2*lm+1)*(2*mm+1),0:nm)  :: e2t_e2z
      !(out) 1 次元化された水平スペクトル・鉛直スペクトル

      e2t_e2z = at_az(e2z)

    end function e2t_e2z

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

      integer :: j,k,l,m

      if ( size(aee,2) /= 2*mm+1 ) &
           call MessageNotify('E','e2a_aee',&
                              '2nd dimension of input data invalid')
      if ( size(aee,3) /= 2*lm+1 ) &
           call MessageNotify('E','e2a_aee',&
                              '3rd dimension of input data invalid')

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

      do k=1,size(aee,1)
         do l=-lm,lm
            do m=-mm,mm
               j = (m+mm+1)+ (2*mm+1)*(l+lm)
               e2a_aee(j,k) = aee(k,m,l)
            enddo
         enddo
      enddo

    end function e2a_aee

    function aee_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,-lm:lm)      :: aee_e2a
      !(out) 任意座標・水平 2 次元スペクトルデータ

      integer :: j,k,l,m

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

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

      do k=1,size(e2a,2)
         do l=-lm,lm
            do m=-mm,mm
               j = (m+mm+1) + (2*mm+1)*(l+lm)
               aee_e2a(k,m,l) = e2a(j,k)
            enddo
         enddo
      enddo

    end function aee_e2a

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

      integer :: j,l,m

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

      do l=-lm,lm
         do m=-mm,mm
            j = (m+mm+1)+ (2*mm+1)*(l+lm)
            e2_ee(j) = ee(m,l)
         enddo
      enddo

    end function e2_ee

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

      integer :: j,l,m

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

      do l=-lm,lm
         do m=-mm,mm
            j = (m+mm+1) + (2*mm+1)*(l+lm)
            ee_e2(m,l) = e2(j)
         enddo
      enddo

    end function ee_e2

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

      real(8), dimension(0:nm,-mm:mm,-lm:lm)             :: tee_Dx_tee
      !(in) X微分されたスペクトルデータ

      integer :: n

      do n=0,nm
         tee_Dx_tee(n,:,:) = ee_Dx_ee(tee(n,:,:))
      enddo

    end function tee_Dx_tee

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

      real(8), dimension(0:nm,-mm:mm,-lm:lm)             :: tee_Dy_tee
      !(in) X微分されたスペクトルデータ

      integer :: n

      do n=0,nm
         tee_Dy_tee(n,:,:) = ee_Dy_ee(tee(n,:,:))
      enddo

    end function tee_Dy_tee

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

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

      tee_Dz_tee = aee_e2a(at_Dz_at(e2a_aee(tee)))

    end function tee_Dz_tee

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

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

      tee_Lapla_tee = tee_LaplaH_tee(tee) + tee_Dz_tee(tee_Dz_tee(tee))

    end function tee_Lapla_tee

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

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

      integer :: n

      do n=0,nm
         tee_LaplaH_tee(n,:,:) = ee_Lapla_ee(tee(n,:,:))
      enddo

    end function tee_LaplaH_tee

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

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

      integer :: n

      do n=0,nm
         tee_LaplaHInv_tee(n,:,:) = ee_LaplaInv_ee(tee(n,:,:))
      enddo

    end function tee_LaplaHInv_tee

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

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

      real(8), dimension(0:nm,-mm:mm,-lm:lm)     :: tee_Div_zyx_zyx_zyx
      !(out) ベクトル場の発散

      tee_Div_zyx_zyx_zyx =   tee_Dx_tee(tee_zyx(zyx_VX)) &
                            + tee_Dy_tee(tee_zyx(zyx_VY)) &
                            + tee_Dz_tee(tee_zyx(zyx_VZ))

    end function tee_Div_zyx_zyx_zyx

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

      real(8), dimension(0:km,0:jm-1)  :: zy_IntX_zyx
      !(out) X方向積分された 2 次元ZY格子点データ

      integer :: i

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

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

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

      integer :: j

      zx_IntY_zyx = 0.0d0
      do j=0,jm-1
         zx_IntY_zyx(:,:) = zx_IntY_zyx(:,:) &
                       + zyx(:,j,:) * y_Y_Weight(j)
      enddo
    end function zx_IntY_zyx

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

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

      integer :: k

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

    function x_IntZY_zyx(zyx)
      !
      ! 3 次元格子点データのZY積分
      !
      real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
      !(in) 3 次元ZYX格子点データ

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

      integer :: j, k

      x_IntZY_zyx = 0.0D0
      do j=0,jm-1
         do k=0,km
            x_IntZY_zyx = x_IntZY_zyx &
                 + zyx(k,j,:) * y_Y_Weight(j) * z_Z_Weight(k)
         enddo
      enddo
    end function x_IntZY_zyx

    function y_IntZX_zyx(zyx)
      !
      ! 3 次元格子点データのZX積分.
      !
      real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
      !(in) 3 次元ZYX格子点データ

      real(8), dimension(0:jm-1)       :: y_IntZX_zyx
      !(out) ZX積分された 1 次元Y格子点データ

      integer :: i, k

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

    function z_IntYX_zyx(zyx)  ! YX(水平)積分
      !
      ! 3 次元格子点データのYX(水平)積分
      !
      real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
      !(in) 3 次元ZYX格子点データ

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

      integer :: i, j

      z_IntYX_zyx = 0
      do j=0,jm-1
         do i=0,im-1
            z_IntYX_zyx = z_IntYX_zyx &
                 + zyx(:,j,i) * x_X_Weight(i) * y_Y_Weight(j)
         enddo
      enddo
    end function z_IntYX_zyx

    function IntZYX_zyx(zyx) ! ZYX(全領域)積分
      !
      ! 3 次元格子点データのZYX(全領域)積分
      !
      real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
      !(in) 3 次元格子点データ

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

      integer :: i, j, k

      IntZYX_zyx = 0
      do i=0,im-1
         do j=0,jm-1
            do k=0,km
               IntZYX_zyx = IntZYX_zyx &
                    + zyx(k,j,i) * x_X_Weight(i) &
                         * y_Y_Weight(j) * z_Z_Weight(k)
            enddo
         enddo
      enddo
    end function IntZYX_zyx

    !----(入力データ zy)---
    function z_IntY_zy(zy)  ! Y積分
      !
      ! 2 次元(ZY)格子点データのY方向積分.
      !
      real(8), dimension(0:km,0:jm-1), intent(in) :: zy
      !(in) 2 次元ZY(子午面)格子点データ

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

      integer :: j

      z_IntY_zy = 0.0d0
      do j=0,jm-1
         z_IntY_zy(:) = z_IntY_zy(:) + zy(:,j) * y_Y_Weight(j)
      enddo
    end function z_IntY_zy

    function y_IntZ_zy(zy)  ! Z積分
      !
      ! 2 次元(ZY)格子点データのZ方向域積分.
      !
      real(8), dimension(0:km,0:jm-1), intent(in) :: zy
      !(in) 2 次元ZY格子点データ

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

      integer :: k

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

    function IntZY_zy(zy)
      !
      ! 2 次元(ZY)格子点データのZY積分
      !
      real(8), dimension(0:km,0:jm-1), intent(in) :: zy
      !(in) 2 次元ZY(子午面)格子点データ

      real(8)                   :: IntZY_zy
      !(out) 積分値
      integer :: j, k

      IntZY_zy = 0
      do j=0,jm-1
         do k=0,km
            IntZY_zy = IntZY_zy &
                 + zy(k,j) * y_Y_Weight(j) * z_Z_Weight(k)
         enddo
      enddo
    end function IntZY_zy

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

      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 次元ZY格子点データ

      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 次元ZY格子点データ

      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

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

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

      zy_AvrX_zyx = zy_IntX_zyx(zyx)/sum(x_X_Weight)

    end function zy_AvrX_zyx

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

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

      zx_AvrY_zyx = zx_IntY_zyx(zyx)/sum(y_Y_Weight)

    end function zx_AvrY_zyx

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

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

      yx_AvrZ_zyx = yx_IntZ_zyx(zyx)/sum(z_Z_Weight)

    end function yx_AvrZ_zyx

    function x_AvrZY_zyx(zyx)  ! ZY積分
      !
      ! 3 次元格子点データのZY平均
      !
      real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
      !(in) 3 次元ZYX格子点データ

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

      x_AvrZY_zyx = x_IntZY_zyx(zyx) &
                   /( sum(y_Y_Weight)*sum(z_Z_Weight) )

    end function x_AvrZY_zyx

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

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

      y_AvrZX_zyx = y_IntZX_zyx(zyx) &
                 /(sum(x_X_Weight)*sum(z_Z_Weight))

    end function y_AvrZX_zyx

    function z_AvrYX_zyx(zyx)  ! YX(水平)積分
      !
      ! 3 次元格子点データのYX(水平)積分
      !
      real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
      !(in) 3 次元格子点データ

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

      z_AvrYX_zyx = z_IntYX_zyx(zyx) &
                 /(sum(x_X_Weight)*sum(y_Y_Weight))

    end function z_AvrYX_zyx

    function AvrZYX_zyx(zyx) ! ZYX(全領域)積分
      !
      ! 3 次元格子点データのZYX(全領域)積分
      !
      real(8), dimension(0:km,0:jm-1,0:im-1), intent(in) :: zyx
      !(in) 3 次元ZYX格子点データ

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

      AvrZYX_zyx = IntZYX_zyx(zyx) &
            /(sum(x_X_Weight)*sum(y_Y_Weight) * sum(z_Z_Weight))

    end function AvrZYX_zyx

    !----(入力データ zy)---
    function z_AvrY_zy(zy)
      !
      ! 2 次元(ZY)格子点データのY方向平均.
      !
      real(8), dimension(0:km,0:jm-1), intent(in) :: zy
      !(in) 2 次元ZY格子点データ

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

      z_AvrY_zy = z_IntY_zy(zy)/sum(y_Y_Weight)

    end function z_AvrY_zy

    function y_AvrZ_zy(zy)
      !
      ! 2 次元(ZY)格子点データのZ方向平均.
      !
      real(8), dimension(0:km,0:jm-1), intent(in) :: zy
      !(in) 2 次元ZY格子点データ

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

      y_AvrZ_zy = y_IntZ_zy(zy)/sum(z_Z_Weight)

    end function y_AvrZ_zy

    function AvrZY_zy(zy)  ! ZY平均
      !
      ! 2 次元(ZY)格子点データのZY平均
      !
      real(8), dimension(0:km,0:jm-1), intent(in) :: zy
      !(in) 2 次元ZY(子午面)格子点データ

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

      AvrZY_zy = IntZY_zy(zy)/(sum(y_Y_Weight)*sum(z_Z_Weight))

    end function AvrZY_zy

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

      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 次元ZY格子点データ

      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 zyx_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 次元スペクトルデータ
!!$
!!$      xzy_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 tee_ZRot_zyx_zyx(zyx_VX,zyx_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,0:jm-1,0:im-1), intent(in) :: zyx_VX
      !(in) ベクトルのX成分

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

      real(8), dimension(0:nm,-mm:mm,-lm:lm)     :: tee_ZRot_zyx_zyx
      !(out) ベクトルの渦度とZベクトルの内積

      tee_ZRot_zyx_zyx = tee_DX_tee(tee_zyx(zyx_VY)) &
                       - tee_DY_tee(tee_zyx(zyx_VX))

    end function tee_ZRot_zyx_zyx

    function tee_ZRotRot_zyx_zyx_zyx(zyx_VX,zyx_VY,zyx_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,0:jm-1,0:im-1), intent(in) :: zyx_VX
      !(in) ベクトルのX成分

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

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

      real(8), dimension(0:nm,-mm:mm,-lm:lm)     :: tee_ZRotRot_zyx_zyx_zyx
      !(out) ベクトル v の z・(▽×▽×v)

      tee_ZRotRot_zyx_zyx_zyx = &
               tee_DZ_tee(   tee_DX_tee(tee_zyx(zyx_VX))   &
                           + tee_DY_tee(tee_zyx(zyx_VY)) ) &
             - tee_LaplaH_tee(tee_zyx(zyx_VZ))

    end function tee_ZRotRot_zyx_zyx_zyx

    subroutine tee_Potential2Vector(&
         zyx_VX,zyx_VY,zyx_VZ,tee_Torvel,tee_Polvel)
      !
      ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
      !
      !     v = ▽x(Ψz) + ▽x▽x(Φz)
      !
      ! の各成分を計算する
      !
      real(8), dimension(0:km,0:jm-1,0:im-1)     :: zyx_VX
      !(out) ベクトル場のX成分

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

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

      real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: tee_Torvel
      !(in) トロイダルポテンシャル

      real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: tee_Polvel
      !(in) ポロイダルポテンシャル

      zyx_VX = zyx_tee(   tee_DY_tee(tee_Torvel) &
                        + tee_DX_tee(tee_DZ_tee(tee_Polvel))  )
      zyx_VY = zyx_tee( - tee_DX_tee(tee_Torvel) &
                        + tee_DY_tee(tee_DZ_tee(tee_Polvel))  )
      zyx_VZ = -zyx_tee(tee_LaplaH_tee(tee_Polvel))

    end subroutine tee_Potential2Vector

    subroutine tee_Potential2Rotation(&
       zyx_RotVX,zyx_RotVY,zyx_RotVZ,tee_Torvel,tee_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,0:jm-1,0:im-1), intent(OUT) :: zyx_RotVX
      !(out) 回転のX成分

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

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

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

      real(8), dimension(0:nm,-mm:mm,-lm:lm), intent(in) :: tee_Polvel
      !(in) ポロイダルポテンシャル

      call tee_Potential2Vector( &
           zyx_RotVX,zyx_RotVY,zyx_RotVZ, &
           -tee_Lapla_tee(tee_Polvel), tee_Torvel)

    end subroutine tee_Potential2Rotation

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

      Interpolate_tee = &
           Interpolate_ee(ee_e2(a_Interpolate_at(e2a_aee(tee_data),z)),x,y)

    end function Interpolate_tee

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

    function zee_ToroidalEnergySpectrum_tee(tee_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,-lm:lm), intent(in) :: tee_TORPOT
      !(in) トロイダルポテンシャル

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

      real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Work
      integer :: l, m

      zee_Work = zee_tee(tee_Torpot)
      do l=-lm,lm
         do m=-mm,mm
            zee_ToroidalEnergySpectrum_tee(:,m,l) &
              = 0.5 * ((2*PI*l/xl)**2+(2*PI*m/yl)**2) &
              * ( zee_Work(:,m,l)**2 + zee_Work(:,-m,-l)**2 )
         enddo
      enddo

    end function zee_ToroidalEnergySpectrum_tee

!!$    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 zee_PoloidalEnergySpectrum_tee(tee_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,-lm:lm), intent(in) :: tee_POLPOT
      !(in) ポロイダルポテンシャル

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


      real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_Data   ! 作業領域
      real(8), dimension(0:km,-mm:mm,-lm:lm) :: zee_DData  ! 作業領域
      integer :: l, m

      zee_Data = zee_tee(tee_POLPOT)
      zee_DData = zee_tee(tee_DZ_tee(tee_POLPOT))

      do l=-lm,lm
         do m=-mm,mm
            zee_PoloidalEnergySpectrum_tee(:,m,l) =                   &
                 + 0.5* ((2*pi*l/xl)**2+(2*pi*m/yl)**2)               &
                 *(   zee_DData(:,m,l)**2 + zee_DData(:,-m,-l)**2     &
                    + ((2*pi*l/xl)**2+(2*pi*m/yl)**2)                 &
                         *( zee_Data(:,m,l)**2 + zee_Data(:,-m,-l)**2))
         enddo
      enddo

    end function zee_PoloidalEnergySpectrum_tee

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

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

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

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

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

      e2t = e2a_aee(tee)

      if (present(values)) then
         e22_values = e2a_aee(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','tee_BoundariesTau','B.C. not supported')
      end select

      tee = aee_e2a(e2t)

    end subroutine tee_BoundariesTau

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

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

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

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

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

      e2t = e2a_aee(tee)

      if (present(values)) then
         e22_values = e2a_aee(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','tee_BoundariesGrid','B.C. not supported')
      end select

      tee = aee_e2a(e2t)

    end subroutine tee_BoundariesGrid

    subroutine tee_TorBoundariesTau(tee_TOR,cond,new)

      ! トロイダル速度ポテンシャルに対して境界条件を適用する.
      ! Chebyshev 空間での境界条件適用
      !
      ! チェビシェフ空間において境界条件を満たすべく高次の係数を定める方法を
      ! とっている(タウ法). トロイダルポテンシャルの境界条件は
      !
      !     ψ = 0 at boundaries           (粘着条件)
      !     ∂ψ/∂z = 0 at boundaries    (応力なし条件)
      !
      ! であるから tee_Boundaries で対応可能だが, 将来のため別途作成しておく.
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(8), dimension(0:nm,-mm:mm,-lm:lm),intent(inout)   :: tee_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*lm+1)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension((2*lm+1)*(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','tee_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*lm+1)*(2*mm+1),0:nm,0:nm),kp((2*lm+1)*(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','tee_TorBoundariesTau',&
                            'Matrix to apply  b.c. newly produced.')
      endif

      e2t_work = e2a_aee(tee_Tor)

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

      tee_Tor = aee_e2a(lusolve(alu,kp,e2t_work))

    end subroutine tee_TorBoundariesTau

    subroutine tee_TorBoundariesGrid(tee_TOR,cond,new)
      !
      ! トロイダル速度ポテンシャルに対して境界条件を適用する.
      ! 鉛直実空間での境界条件適用.
      !
      ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
      ! 条件を課している(選点法). このルーチンを用いるためには
      ! tee_Initial にて設定するチェビシェフ切断波数(nm)と鉛直格子点数(km)を
      ! 等しくしておく必要がある.
      !
      ! トロイダル速度ポテンシャルの境界条件は
      !
      !     ψ = 0 at boundaries           (粘着条件)
      !     ∂ψ/∂z = 0 at boundaries    (応力なし条件)
      !
      ! であるので tee_Boundaries で対応可能だが, 将来のため別途作成しておく
      !
      ! 最初に呼ばれるときはオプショナル引数 new に関係なく行列が設定される.
      !
      real(8), dimension(0:nm,-mm:mm,-lm:lm),intent(inout)   :: tee_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*lm+1)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension((2*lm+1)*(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','tee_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','tee_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*lm+1)*(2*mm+1),0:km,0:nm),kp((2*lm+1)*(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_aee(tee_Tor))
      e2z_work(:,0)  = 0.0D0
      e2z_work(:,km) = 0.0D0
      tee_TOR = aee_e2a(lusolve(alu,kp,e2z_work))

    end subroutine tee_TorBoundariesGrid

    function tee_LaplaPol2PolTau_tee(tee,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,-lm:lm),intent(in) :: tee
              !(in) 入力▽^2φ分布

      real(8), dimension(0:nm,-mm:mm,-lm:lm)      :: tee_LaplaPol2PolTau_tee
              !(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*lm+1)*(2*mm+1),0:nm)  :: e2t_work
      real(8), dimension((2*lm+1)*(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','tee_LaplaPol2PolTau_tee','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*lm+1)*(2*mm+1),0:nm,0:nm),kp((2*lm+1)*(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_aee(tee_lapla_tee(aee_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','tee_LaplaPol2PolTau_tee',&
                           'Matrix to apply  b.c. newly produced.')
      endif

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

      tee_LaplaPol2PolTau_tee = aee_e2a(lusolve(alu,kp,e2t_work))

    end function tee_LaplaPol2PolTau_tee

    function zee_LaplaPol2Pol_zee(zee,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,-lm:lm),intent(in)  :: zee
              !(in) 入力▽^2φ分布

      real(8), dimension(0:km,-mm:mm,-lm:lm)             :: zee_LaplaPol2Pol_zee
              !(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*lm+1)*(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','zee_laplapol2pol_zee','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','zee_LaplaPol2Pol_zee', &
             'Chebyshev truncation and number of grid points should be same.')
         endif

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

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

            ! 各水平波数に関して独立の式
            alu(:,:,k) &
                 = e2a_aee(zee_tee(tee_lapla_tee(tee_zee(aee_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','zee_LaplaPol2Pol_zee',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      e2z_work        = e2a_aee(zee)
      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)

      zee_laplapol2pol_zee = aee_e2a(e2z_work)

    end function zee_LaplaPol2Pol_zee

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

      real(8), dimension(0:nm,-mm:mm,-lm:lm)      :: tee_LaplaPol2PolGrid_tee
              !(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*lm+1)*(2*mm+1),0:nm)  :: e2t_work
      real(8), dimension((2*lm+1)*(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','tee_LaplaPol2PolGrid_tee','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','tee_LaplaPol2PolGrid_tee', &
             'Chebyshev truncation and number of grid points should be same.')
         endif

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

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

            ! 各水平波数に関して独立の式
            alu(:,:,n) = e2a_aee(zee_tee(tee_lapla_tee(aee_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','tee_LaplaPol2PolGrid_tee',&
                           'Matrix to apply  b.c. newly produced.')
      endif

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

      tee_LaplaPol2PolGrid_tee = aee_e2a(lusolve(alu,kp,e2z_work))

    end function tee_LaplaPol2PolGrid_tee

    subroutine tee_TormagBoundariesTau(tee_TOR,new)

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

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

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

      real(8), dimension((2*lm+1)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension((2*lm+1)*(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*lm+1)*(2*mm+1),0:nm,0:nm),kp((2*lm+1)*(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','tee_TormagBoundariesTau',&
                            'Matrix to apply  b.c. newly produced.')
      endif

      e2t_work = e2a_aee(tee_Tor)

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


      tee_Tor = aee_e2a(lusolve(alu,kp,e2t_work))

    end subroutine tee_TormagBoundariesTau

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

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

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

      real(8), dimension((2*lm+1)*(2*mm+1),0:nm) :: e2t_work
      real(8), dimension((2*lm+1)*(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','tee_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*lm+1)*(2*mm+1),0:km,0:nm),kp((2*lm+1)*(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_aee(tee_Tor))
      e2z_work(:,0)  = 0.0D0
      e2z_work(:,km) = 0.0D0
      tee_TOR = aee_e2a(lusolve(alu,kp,e2z_work))

    end subroutine tee_TormagBoundariesGrid

    subroutine tee_PolmagBoundariesTau(tee_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,-lm:lm),intent(inout)   :: tee_POL
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

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

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

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

      real(8), dimension((2*lm+1)*(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*lm+1)*(2*mm+1),0:nm,0:nm),kp((2*lm+1)*(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
               e2index = (2*mm+1)*(l+lm) + (m+mm+1)
               rl = 2*PI/xl*l ; rm = 2*PI/yl*m
               e2t_kh(e2index,:) = sqrt(rl**2+rm**2)
            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(:,lm-1,n) = e2z_work(:,0)

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

         call ludecomp(alu,kp)

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

      e2t_work = e2a_aee(tee_POL)
      e2t_work(:,lm-1) = 0.0D0
      e2t_work(:,lm)   = 0.0D0
      tee_POL = aee_e2a(lusolve(alu,kp,e2t_work))

    end subroutine tee_PolmagBoundariesTau

    subroutine tee_PolmagBoundariesGrid(tee_POL,new)
      !
      ! 磁場ポロイダルポテンシャルに対して境界条件を適用する.
      ! 鉛直実空間での境界条件適用.
      !
      ! 鉛直実格子点空間において内部領域の値と境界条件を満たすように
      ! 条件を課している(選点法). このルーチンを用いるためには
      ! tee_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,-lm:lm),intent(inout)   :: tee_POL
              !(inout) 境界条件を適用するデータ. 修正された値を返す.

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

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

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

      real(8), dimension((2*lm+1)*(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','tee_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*lm+1)*(2*mm+1),0:km,0:nm),kp((2*lm+1)*(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
               e2index = (2*mm+1)*(l+lm) + (m+mm+1)
               rl = 2*PI/xl*l ; rm = 2*PI/yl*m
               e2t_kh(e2index,:) = sqrt(rl**2+rm**2)
            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','tee_PolmagBoundariesGrid',&
                           'Matrix to apply  b.c. newly produced.')
      endif

      e2z_work = az_at(e2a_aee(tee_POL))
      e2z_work(:,0)  = 0.0D0
      e2z_work(:,km) = 0.0D0
      tee_POL = aee_e2a(lusolve(alu,kp,e2z_work))

    end subroutine tee_PolmagBoundariesGrid

end module tee_module
