tee_module

    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

    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

    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

    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 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

    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

    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

    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

    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 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 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 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 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 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

    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_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

    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

    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

   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 (2013/08/20) is initialized')

   end subroutine tee_Initial

    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_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_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

    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 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

    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.
      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)
               e2t_kh(e2index,:) = sqrt(dble(l**2+m**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

    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.
      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)
               e2t_kh(e2index,:) = sqrt(dble(l**2+m**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_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

    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_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

    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_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_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

    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

    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_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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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