! -*- mode: f90; coding: utf-8 -*-
!----------------------------------------------------------------------
! Copyright(c) 2019 SPMDODEL Development Group. All rights reserved.
!
! 履歴  2019/03/28 竹広真一
!       2019/09/28 佐々木洋平
!       2019/10/09 佐々木洋平
!----------------------------------------------------------------------
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @copyright Copyright (C) SPMODEL Development Group, 2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @brief 1 次元有限領域: チェビシェフ変換
!>
!> @details
!> 1 次元有限領域の下での流体運動をチェビシェフ変換によりスペクトル数値
!> 計算するための Fortran90 関数を提供する.
!>
!> 2 次元データの 1 次元に関して同時にスペクトル計算を実行するための
!> 関数も提供しており, 2, 3 次元領域での計算のベースも提供する.
!>
!> このモジュールは内部で ISPACK3/FXPACK の Fortran77 サブルーチンを
!> 呼んでいる. スペクトルデータおよび格子点データの格納方法については
!> ISPACK3/FXPACK のマニュアルを参照されたい.
!>
!> ## 関数・変数の名前と型について
!>
!> ### 命名法
!>
!> * 関数名の先頭 (`t_`, `g_`, `at_`, `ag_`) は, 返す値の形を示している.
!>   * `t_`  :: チェビシェフデータ
!>   * `g_`  :: 1 次元格子点データ
!>   * `at_` :: 1 次元チェビシェフデータが複数並んだ 2 次元データ
!>   * `ag_` :: 1 次元格子点データが複数並んだ 2 次元データ.
!>
!> * 関数名の間の文字列(`Dx`)は, その関数の作用を表している.
!>
!> * 関数名の最後 (`_e`,`_at`,`_g`,`_ag`) は,
!>   入力変数の形がチェビシェフデータおよび格子点データであることを示している.
!>   * `_t`  :: チェビシェフデータ
!>   * `_g`  :: 1 次元格子点データ
!>   * `_at` :: 1 次元チェビシェフデータが複数並んだ 2 次元データ
!>   * `_ag` :: 1 次元格子点データが複数並んだ 2 次元データ
!>
!> ### 各データの種類の説明
!>
!> * 以下, `DP  = kind(1.0D)` である
!> * `g` : 1 次元格子点データ.
!>   * 変数の種類と次元は `real(DP), dimension(0:im)`.
!>   * `im` は X 座標の格子点数であり, サブルーチン at_initial() にて
!>     あらかじめ設定しておく.
!>
!> * `t` : チェビシェフデータ.
!>   * 変数の種類と次元は `real(DP), dimension(0:km)`.
!>   * `km` は X 方向の最大波数であり, サブルーチン at_initial() にて
!>     あらかじめ設定しておく.
!>
!> * `ag` : 1 次元(X)格子点データの並んだ 2 次元データ.
!>   * 変数の種類と次元は `real(DP), dimension(:,0:im)`.
!>     第 2 次元が X 方向を表す.
!>
!> * `at` : 1 次元チェビシェフデータの並んだ 2 次元データ.
!>   * 変数の種類と次元は `real(DP), dimension(:,0:km)`.
!>     第 2 次元がスペクトルを表す.
!>
!> * `g_` で始まる関数が返す値は 1 次元格子点データに同じ.
!>
!> * `t_` で始まる関数が返す値はチェビシェフデータに同じ.
!>
!> * `ag_` で始まる関数が返す値は 1 次元格子点データの並んだ
!>   2 次元データに同じ.
!>
!> * `at_` で始まる関数が返す値は 1 次元チェビシェフデータの並んだ
!>   2 次元データに同じ.
!>
!> * チェビシェフデータに対する微分等の作用とは, 対応する格子点データに
!>   微分などを作用させたデータをチェビシェフ変換したもののことである.

module at_module

  use dc_types, only: DP
  use dc_message, only: MessageNotify
  use spml_utils, only: pi
  use lumatrix, only: lusolve, ludecomp
  use ae_module, only: ae_initial, ae_ag, ag_ae
  implicit none

  private
  public g_X
  public g_X_Weight
  public at_Initial
  public at_Finalize
  public ag_at
  public at_ag
  public g_t
  public t_g
  public at_Dx_at
  public t_Dx_t
  public Interpolate_t
  public a_Interpolate_at
  public at_Boundaries
  public t_Boundaries
  public at_Boundaries_DD
  public at_Boundaries_DN
  public at_Boundaries_ND
  public at_Boundaries_NN
  public at_BoundariesTau_DD
  public at_BoundariesTau_DN
  public at_BoundariesTau_ND
  public at_BoundariesTau_NN
  public at_BoundariesGrid
  public t_BoundariesGrid
  public at_BoundariesGrid_DD
  public at_BoundariesGrid_DN
  public at_BoundariesGrid_ND
  public at_BoundariesGrid_NN
  public a_Int_ag
  public Int_g
  public a_Avr_ag
  public Avr_g

  !----------------------------------------------------------------------
  !> 両端ディリクレ型境界条件の適用(タウ法). 両境界での値を与える.
  interface at_Boundaries_DD
     module procedure at_BoundariesTau_DD_1d
     module procedure at_BoundariesTau_DD_2d
  end interface at_Boundaries_DD

  !----------------------------------------------------------------------
  !> ディリクレ・ノイマン型境界条件の適用(タウ法).
  !> i=0 で値, i=im で勾配の値を与える.
  interface at_Boundaries_DN
     module procedure at_BoundariesTau_DN_1d
     module procedure at_BoundariesTau_DN_2d
  end interface at_Boundaries_DN

  !----------------------------------------------------------------------
  !> ディリクレ・ノイマン型境界条件の適用(タウ法).
  !> i=0 で勾配の値, i=im で値を与える.
  interface at_Boundaries_ND
     module procedure at_BoundariesTau_ND_1d
     module procedure at_BoundariesTau_ND_2d
  end interface at_Boundaries_ND

  !----------------------------------------------------------------------
  !> 両端ノイマン型境界条件の適用(タウ法).
  !> 両境界で勾配の値を与える.
  interface at_Boundaries_NN
     module procedure at_BoundariesTau_NN_1d
     module procedure at_BoundariesTau_NN_2d
  end interface at_Boundaries_NN

  !----------------------------------------------------------------------
  !> 両端ディリクレ型境界条件の適用(タウ法).
  !> 両境界での値を与える.
  interface at_BoundariesTau_DD
     module procedure at_BoundariesTau_DD_1d
     module procedure at_BoundariesTau_DD_2d
  end interface at_BoundariesTau_DD

  !----------------------------------------------------------------------
  !> ディリクレ・ノイマン型境界条件の適用(タウ法).
  !> i=0 で値, i=im で勾配の値を与える.
  interface at_BoundariesTau_DN
     module procedure at_BoundariesTau_DN_1d
     module procedure at_BoundariesTau_DN_2d
  end interface at_BoundariesTau_DN

  !----------------------------------------------------------------------
  !> ノイマン・ディリクレ型境界条件の適用(タウ法).
  !> i=0 で勾配の値, i=im で値を与える.
  interface at_BoundariesTau_ND
     module procedure at_BoundariesTau_ND_1d
     module procedure at_BoundariesTau_ND_2d
  end interface at_BoundariesTau_ND

  !----------------------------------------------------------------------
  !> 両端ノイマン型境界条件の適用(タウ法).
  !> 両境界で勾配の値を与える.
  interface at_BoundariesTau_NN
     module procedure at_BoundariesTau_NN_1d
     module procedure at_BoundariesTau_NN_2d
  end interface at_BoundariesTau_NN

  !----------------------------------------------------------------------
  !> 両端ディリクレ型境界条件の適用(実空間での評価).
  !> 両境界での値を与える.
  interface at_BoundariesGrid_DD
     module procedure at_BoundariesGrid_DD_1d
     module procedure at_BoundariesGrid_DD_2d
  end interface at_BoundariesGrid_DD

  !----------------------------------------------------------------------
  !> ディリクレ・ノイマン型境界条件の適用(実空間での評価).
  !> i=0 で値, i=im で勾配の値を与える.
  interface at_BoundariesGrid_DN
     module procedure at_BoundariesGrid_DN_1d
     module procedure at_BoundariesGrid_DN_2d
  end interface at_BoundariesGrid_DN

  !----------------------------------------------------------------------
  !> ノイマン・ディリクレ型境界条件の適用(実空間での評価).
  !> i=0 で勾配の値, i=0 で値を与える.
  interface at_BoundariesGrid_ND
     module procedure at_BoundariesGrid_ND_1d
     module procedure at_BoundariesGrid_ND_2d
  end interface at_BoundariesGrid_ND

  !----------------------------------------------------------------------
  !> 両端ノイマン型境界条件の適用(実空間での評価).
  !> 両境界で勾配の値を与える.
  interface at_BoundariesGrid_NN
     module procedure at_BoundariesGrid_NN_1d
     module procedure at_BoundariesGrid_NN_2d
  end interface at_BoundariesGrid_NN

  !----------------------------------------------------------------------
  !> 格子点座標(X)を格納した 1 次元配列
  real(DP), allocatable :: g_X(:)
  !> 重み座標を格納した 1 次元配列
  real(DP), allocatable :: g_X_Weight(:)

  !----------------------------------------------------------------------
  integer(8), dimension(:), allocatable  :: it
  real(DP), dimension(:), allocatable    :: t
  !> 格子点数
  integer(8) :: im
  !> 切断波数
  integer(8) :: km
  !> 領域の大きさ
  real(DP)   :: xl
  !> 初期化スイッチ
  logical :: at_Initialize = .false.
  !> 境界条件ルーチン初期化スイッチ
  logical :: at_BoundariesTau_DD_first = .true.
  !> 境界条件ルーチン初期化スイッチ
  logical :: at_BoundariesTau_DN_first = .true.
  !> 境界条件ルーチン初期化スイッチ
  logical :: at_BoundariesTau_ND_first = .true.
  !> 境界条件ルーチン初期化スイッチ
  logical :: at_BoundariesTau_NN_first = .true.
  !> 境界条件ルーチン初期化スイッチ
  logical :: at_BoundariesGrid_DD_first = .true.
  !> 境界条件ルーチン初期化スイッチ
  logical :: at_BoundariesGrid_DN_first = .true.
  !> 境界条件ルーチン初期化スイッチ
  logical :: at_BoundariesGrid_ND_first = .true.
  !> 境界条件ルーチン初期化スイッチ
  logical :: at_BoundariesGrid_NN_first = .true.

  ! save :: im, km, xl, it, t, g_X, g_X_Weight, at_Initialize
  save :: im, km, xl, g_X, g_X_Weight, at_Initialize
  save :: at_BoundariesTau_DD_first
  save :: at_BoundariesTau_DN_first
  save :: at_BoundariesTau_ND_first
  save :: at_BoundariesTau_NN_first
  save :: at_BoundariesGrid_DD_first
  save :: at_BoundariesGrid_DN_first
  save :: at_BoundariesGrid_ND_first
  save :: at_BoundariesGrid_NN_first

contains
  !-----------------------------------------------------------------------
  !> チェビシェフ変換の格子点数, 波数, 領域の大きさを設定する.
  !>
  !> 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
  !> 初期設定をしなければならない.
  subroutine at_Initial(i,k,xmin,xmax)
    !> 格子点数
    integer,intent(in) :: i
    !> 切断波数
    integer,intent(in) :: k
    !> 座標の範囲
    real(DP),intent(in) :: xmin, xmax
    integer(8) :: ii,kk

    im=i
    km=k
    xl = xmax-xmin

    if ( im <= 0 .or. km <= 0 ) then
       call MessageNotify('E','at_initial', &
            & 'Number of grid points and waves should be positive')
    elseif ( mod(im,2_8) /= 0 ) then
       call MessageNotify('E','at_initial', &
            & 'Number of grid points should be even')
    elseif ( km > im ) then
       call MessageNotify('E','at_initial', &
            & 'KM shoud be less equal IM')
    endif

    allocate(it(im),t(3*im))
    call fxrini(2*im,it,t)

    allocate(g_X(0:im))
    do ii=0,im
       g_X(ii) = (xmax+xmin)/2 + xl/2 * cos(pi*ii/im)
    end do

    allocate(g_X_Weight(0:im))

    do ii=0,im
       g_X_Weight(ii) = 1.0D0
       do kk=2,km,2
          g_X_Weight(ii) = g_X_Weight(ii) &
               + 2.0D0/(1.0D0 - kk**2) * cos(kk*ii*pi/im)
       enddo
       ! 最後の和は factor 1/2.
       if ( (km == im) .and. (mod(im,2_8) == 0) ) then
          g_X_Weight(ii) = g_X_Weight(ii) &
               - 1.0D0/(1.0D0 - km**2)* cos(km*ii*pi/im)
       endif
       g_X_Weight(ii) = 2.0D0/im * g_X_Weight(ii) * (xl/2.0D0)
    enddo
    ! g_X_Weight(0)  = g_X_Weight(0) / 2
    ! g_X_Weight(im) = g_X_Weight(im) / 2
    g_X_Weight(0)  = g_X_Weight(0) * 0.5D0
    g_X_Weight(im) = g_X_Weight(im) * 0.5D0

    at_Initialize = .true.

    call MessageNotify('M','at_initial',&
         & 'at_module (2019/09/28) is initialized')
  end subroutine at_Initial

  !-----------------------------------------------------------------------
  !> チェビシェフデータから格子データへ変換する(2 次元配列用).
  function ag_at(at_data)
    !> チェビシェフデータ
    real(DP), dimension(:,0:), intent(in)     :: at_data
    !> 格子点データ
    real(DP), dimension(size(at_data,1),0:im) :: ag_at

    real(DP), dimension(size(at_data,1),0:2*im-1) :: y
    integer(8) :: nm, k

    nm=size(at_data,1)

    y = 0.0D0
    y(:,0) = at_data(:,0)
    do k=1,min(km,im-1)
       y(:,2*k)   = at_data(:,k)
       y(:,2*k+1) = 0.0D0
    enddo
    if ( km == im ) then
       y(:,1) = at_data(:,km)
    endif

    call fxrtba(nm,2*im,y,it,t)
    ! ag_at = y(:,0:im)/2
    ag_at = y(:,0:im)*0.50D0

  end function ag_at

  !-----------------------------------------------------------------------
  !> チェビシェフデータから格子データへ変換する(1 次元配列用).
  function g_t(t_data)
    !> チェビシェフデータ
    real(DP), dimension(:), intent(in)  :: t_data
    !> 格子点データ
    real(DP), dimension(0:im)           :: g_t

    real(DP), dimension(1,size(t_data)) :: t_work
    real(DP), dimension(1,0:im)         :: g_work

    t_work(1,:) = t_data
    g_work = ag_at(t_work)
    g_t = g_work(1,:)

  end function g_t

  !-----------------------------------------------------------------------
  !> 格子データからチェビシェフデータへ変換する(2 次元配列用).
  function at_ag(ag_data)
    !> 格子点データ
    real(DP), dimension(:,0:), intent(in)     :: ag_data
    !> チェビシェフデータ
    real(DP), dimension(size(ag_data,1),0:km) :: at_ag

    real(DP), dimension(size(ag_data,1),0:2*im-1) :: y
    integer(8) :: nm, k

    nm=size(ag_data,1)

    y(:,0) = ag_data(:,0)
    do k=1,im-1
       y(:,2*k)   = ag_data(:,k)
       y(:,2*k+1) = 0.0D0
    enddo
    y(:,1) = ag_data(:,im)

    call fxrtba(nm,2*im,y,it,t)
    at_ag = y(:,0:km)/im

  end function at_ag

  !-----------------------------------------------------------------------
  !> 格子データからチェビシェフデータへ変換する(1 次元配列用).
  function t_g(g_data)
    !> 格子点データ
    real(DP), dimension(:), intent(in)     :: g_data
    !>  チェビシェフデータ
    real(DP), dimension(0:km)              :: t_g


    real(DP), dimension(1,size(g_data)) :: ag_work
    real(DP), dimension(1,0:km)         :: at_work

    ag_work(1,:) = g_data
    at_work = at_ag(ag_work)
    t_g = at_work(1,:)

  end function t_g

  !-----------------------------------------------------------------------
  !> 入力チェビシェフデータに X 微分を作用する(2 次元配列用).
  !>
  !> チェビシェフデータの X 微分とは, 対応する格子点データに X 微分を
  !> 作用させたデータのチェビシェフ変換のことである.
  function at_Dx_at(at_data)
    !> 入力チェビシェフデータ
    real(DP), dimension(:,0:), intent(in)                    :: at_data
    !> チェビシェフデータの X 微分
    real(DP), dimension(size(at_data,1),0:size(at_data,2)-1) :: at_Dx_at

    integer(8) :: m, k
    integer(8) :: nm, kmax

    nm=size(at_data,1)
    kmax=size(at_data,2)-1

    ! TODO: loop 併合できないか?
    if ( kmax == im ) then
       do m=1,nm
          at_Dx_at(m,kmax)   = 0.0D0
          ! at_Dx_at(m,kmax-1) = 2.0D0 * km * at_data(m,kmax) /2
          at_Dx_at(m,kmax-1) = 2.0D0 * km * at_data(m,kmax)* 0.5D0
       enddo
    else
       do m=1,nm
          at_Dx_at(m,kmax)   = 0.0D0
          ! スタートはグリッド対応最大波数未満. Factor 1/2 不要
          at_Dx_at(m,kmax-1) = 2.0D0 * km * at_data(m,kmax)
       enddo
    endif

    do k=kmax-2,0,-1
       do m=1,nm
          at_Dx_at(m,k) = at_Dx_at(m,k+2) + 2.0D0*(k+1)*at_data(m,k+1)
       enddo
    enddo

    do k=0,kmax
       do m=1,nm
          at_Dx_at(m,k) = 2.0D0/xl * at_Dx_at(m,k)
       enddo
    enddo

  end function at_Dx_at


  !-----------------------------------------------------------------------
  !> 入力チェビシェフデータに X 微分を作用する(1 次元配列用).
  !>
  !> チェビシェフデータの X 微分とは, 対応する格子点データに X 微分を
  !> 作用させたデータのチェビシェフ変換のことである.
  function t_Dx_t(t_data)
    !> 入力チェビシェフデータ
    real(DP), dimension(:), intent(in)   :: t_data
    !> チェビシェフデータの X 微分
    real(DP), dimension(size(t_data))    :: t_Dx_t

    real(DP), dimension(1,size(t_data))  :: at_work

    at_work(1,:) = t_data
    at_work = at_Dx_at(at_work)
    t_Dx_t = at_work(1,:)

  end function t_Dx_t

  !-----------------------------------------------------------------------
  !> 1 次元格子点データが並んだ 2 次元配列の積分
  function a_Int_ag(ag)
    !> 入力格子点データ
    real(DP), dimension(:,0:), intent(in)     :: ag
    !> 積分したデータ
    real(DP), dimension(size(ag,1))           :: a_Int_ag

    integer(8) :: i

    if ( size(ag,2) < im+1 ) then
       call MessageNotify('E','ae_ag', &
            'The Grid points of input data too small.')
    elseif ( size(ag,2) > im+1 ) then
       call MessageNotify('W','ae_ag', &
            'The Grid points of input data too large.')
    endif

    a_Int_ag = 0.0d0
    do i=0,im
       a_Int_ag(:) = a_Int_ag(:) + ag(:,i)*g_X_Weight(i)
    enddo
  end function a_Int_ag

  !-----------------------------------------------------------------------
  !> 1 次元格子点データの積分
  function Int_g(g)
    !> 格子点データ
    real(DP), dimension(0:im), intent(in)   :: g
    !> 積分値
    real(DP)                                :: Int_g
    Int_g = sum(g*g_X_Weight)
  end function Int_g

  !-----------------------------------------------------------------------
  !> 1 次元格子点データが並んだ 2 次元配列の平均
  function a_Avr_ag(ag)
    !> 入力格子点データ
    real(DP), dimension(:,0:), intent(in)   :: ag
    !> 平均したデータ
    real(DP), dimension(size(ag,1))         :: a_Avr_ag
    a_Avr_ag = a_Int_ag(ag)/sum(g_X_Weight)

  end function a_Avr_ag

  !-----------------------------------------------------------------------
  !> 1 次元格子点データの平均
  function Avr_g(g)
    !> 格子点データ
    real(DP), dimension(0:im), intent(in)   :: g
    !> 積分値
    real(DP)                                :: Avr_g

    Avr_g = Int_g(g)/sum(g_X_Weight)
  end function Avr_g

  !-----------------------------------------------------------------------
  !> 1 次元データの補間
  function Interpolate_t(t_data,xval)
    !> 入力チェビシェフデータ
    real(DP), dimension(0:), intent(in)  :: t_data
    !> 補間する点の座標
    real(DP), intent(in)                 :: xval
    !> 補間した結果の値
    real(DP)                             :: Interpolate_t

    integer :: kmax
    ! Crenshow's reccurence formula 計算用変数
    real(DP) :: y2, y1, y0, x
    integer :: k
    ! DO 文変数

    kmax = size(t_data)-1

    ! x =(xval -(g_X(0)+g_X(im))/2 )/(g_X(0)-g_X(im))*2
    x =(xval -(g_X(0)+g_X(im))*0.50D0 )/(g_X(0)-g_X(im))*2.0D0

    y2 = 0.0D0
    y1 = 0.0D0
    do k=kmax,1,-1
       y0 = 2.0D0*x*y1 - y2 + t_data(k)
       y2 = y1
       y1 = y0
    enddo

    ! Interpolate_t = - y2 + x*y1 + t_data(0)/2
    Interpolate_t = - y2 + x*y1 + t_data(0)*0.5D0
    if ( kmax == int(im) ) then
       Interpolate_t = Interpolate_t -t_data(kmax)/2.0D0*cos(kmax*acos(x))
    endif

  end function Interpolate_t

  !-----------------------------------------------------------------------
  !> 1 次元データの並んだ 2次元配列の補間
  function a_Interpolate_at(at_data,xval)
    !> 入力チェビシェフデータ
    real(DP), dimension(:,0:), intent(in) :: at_data
    !> 補間する点の座標
    real(DP), intent(in)                  :: xval
    !> 補間した結果の値
    real(DP), dimension(size(at_data,1))  :: a_Interpolate_at


    integer :: kmax
    ! 入力配列の最大次数
    real(DP), dimension(size(at_data,1))  :: y2, y1, y0
    real(DP)                              :: x
    ! Crenshow's reccurence formula 計算用変数
    integer :: k
    ! DO 文変数

    kmax = size(at_data,2)-1

    ! x =(xval -(g_X(0)+g_X(im))/2 )/(g_X(0)-g_X(im))*2
    x =(xval -(g_X(0)+g_X(im))*0.5D0 )/(g_X(0)-g_X(im))*2.0D0

    y2 = 0.0D0
    y1 = 0.0D0
    do k=kmax,1,-1
       y0 = 2.0D0*x*y1 - y2 + at_data(:,k)
       y2 = y1
       y1 = y0
    enddo

    ! a_Interpolate_at = - y2 + x*y1 + at_data(:,0)/2
    a_Interpolate_at = - y2 + x*y1 + at_data(:,0)*0.5D0
    if ( kmax == int(im) ) then
       a_Interpolate_at = a_Interpolate_at &
            & - at_data(:,kmax)/2.0D0*cos(kmax*acos(x))
    endif
  end function a_Interpolate_at

  !-----------------------------------------------------------------------
  !> 境界条件の適用(タウ法, 2 次元配列用)
  !>
  !> 境界条件と両境界での値をオプションで与える.
  !> デフォルトは DD, values = 0
  subroutine at_Boundaries(at_data,cond,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), intent(INOUT)                  :: at_data(:,:)
    !> 境界条件.
    !> * DD : 両端ディリクレ
    !> * DN : 上端ディリクレ, 下端ノイマン
    !> * ND : 上端ノイマン,   下端ディリクレ
    !> * NN : 両端ノイマン
    character(len=2), intent(IN), optional  :: cond
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values

    real(DP)          :: BCvalues(size(at_data),2)
    character(len=2) :: BC

    if ( present(cond) ) then
       BC = cond
    else
       BC = 'DD'
    endif

    if ( present(values) ) then
       BCvalues = values
    else
       BCvalues = 0.0D0
    endif
    select case(BC)
    case ('DD')
      call at_Boundaries_DD(at_data,BCvalues)
    case ('DN')
      call at_Boundaries_DN(at_data,BCValues)
    case ('ND')
      call at_Boundaries_ND(at_data,BCValues)
    case ('NN')
      call at_Boundaries_NN(at_data,BCValues)
    case default
      call at_Boundaries_DD(at_data,BCvalues)
    end select

  end subroutine at_Boundaries

  !---------------------------------------------------------------------
  !> 境界条件の適用(タウ法, 1 次元配列用)
  !> 境界条件と両境界での値をオプションで与える.
  !> デフォルトは DD, values = 0
  subroutine t_Boundaries(t_data,cond,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), intent(INOUT)                  :: t_data(:)
    !> 境界条件.
    !> * DD : 両端ディリクレ
    !> * DN : 上端ディリクレ, 下端ノイマン
    !> * ND : 上端ノイマン,   下端ディリクレ
    !> * NN : 両端ノイマン
    character(len=2), intent(IN), optional  :: cond
    !> 境界値(2)
    real(DP), dimension(2), intent(in), optional  :: values

    real(DP)          :: BCvalues(2)
    character(len=2) :: BC


    if ( present(cond) ) then
      BC = cond
    else
      BC = 'DD'
    endif

    if ( present(values) ) then
      BCvalues = values
    else
      BCvalues = 0.0D0
    endif

    select case(BC)
    case ('DD' )
      call at_Boundaries_DD(t_data,BCvalues)
    case ('DN' )
      call at_Boundaries_DN(t_data,BCvalues)
    case ('ND' )
      call at_Boundaries_ND(t_data,BCvalues)
    case ('NN' )
      call at_Boundaries_NN(t_data,BCvalues)
    case default
      call at_Boundaries_DD(t_data,BCvalues)
    end select

  end subroutine t_Boundaries

  !---------------------------------------------------------------------
  !> 境界条件の適用(実空間での評価, 2 次元配列用)
  !> 境界条件と両境界での値をオプションで与える.
  !> デフォルトは DD, values = 0
  subroutine at_BoundariesGrid(at_data,cond,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), intent(INOUT)                  :: at_data(:,:)
    !> 境界条件.
    !> * DD : 両端ディリクレ
    !> * DN : 上端ディリクレ, 下端ノイマン
    !> * ND : 上端ノイマン,   下端ディリクレ
    !> * NN : 両端ノイマン
    character(len=2), intent(IN), optional  :: cond
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values

    real(DP)          :: BCvalues(size(at_data),2)
    character(len=2) :: BC

    if ( present(cond) ) then
       BC = cond
    else
       BC = 'DD'
    endif

    if ( present(values) ) then
       BCvalues = values
    else
       BCvalues = 0.0D0
    endif

    select case(BC)
    case ('DD')
      call at_BoundariesGrid_DD(at_data,BCvalues)
    case ('DN')
      call at_BoundariesGrid_DN(at_data,BCValues)
    case ('ND')
      call at_BoundariesGrid_ND(at_data,BCValues)
    case ('NN')
      call at_BoundariesGrid_NN(at_data,BCValues)
    end select

  end subroutine at_BoundariesGrid

  !---------------------------------------------------------------------
  !> 境界条件の適用(タウ法, 1 次元配列用)
  !> 境界条件と両境界での値をオプションで与える.
  !> デフォルトは DD, values = 0
  subroutine t_BoundariesGrid(t_data,cond,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), intent(INOUT)                  :: t_data(:)
    !> 境界条件.
    !> * DD : 両端ディリクレ
    !> * DN : 上端ディリクレ, 下端ノイマン
    !> * ND : 上端ノイマン,   下端ディリクレ
    !> * NN : 両端ノイマン
    character(len=2), intent(IN), optional  :: cond
    !> 境界値(2)
    real(DP), dimension(2), intent(in), optional  :: values

    real(DP)          :: BCvalues(2)
    character(len=2) :: BC

    if ( present(cond) ) then
       BC = cond
    else
       BC = 'DD'
    endif

    if ( present(values) ) then
       BCvalues = values
    else
       BCvalues = 0.0D0
    endif
    select case(BC)
    case ('DD')
      call at_BoundariesGrid_DD(t_data,BCvalues)
    case ('DN')
      call at_BoundariesGrid_DN(t_data,BCvalues)
    case ('ND')
      call at_BoundariesGrid_ND(t_data,BCvalues)
    case ('NN')
      call at_BoundariesGrid_NN(t_data,BCvalues)
    end select

  end subroutine t_BoundariesGrid

  !----------------------------------------------------------------------
  !> モジュールの終了処理(割り付け配列の解放)をおこなう.
  subroutine at_Finalize
    if ( .not. at_initialize ) then
      call MessageNotify('W','at_Finalize',&
        'at_module not initialized yet')
      return
    endif

    deallocate(t)                  ! 変換用配列
    deallocate(it)                 ! 変換用配列
    deallocate(g_X)
    deallocate(g_X_Weight)         ! 格子点座標格納配列

    at_initialize = .false.

    at_BoundariesTau_DD_first = .true.
    at_BoundariesTau_DN_first = .true.
    at_BoundariesTau_ND_first = .true.
    at_BoundariesTau_NN_first = .true.
    at_BoundariesGrid_DD_first = .true.
    at_BoundariesGrid_DN_first = .true.
    at_BoundariesGrid_ND_first = .true.
    at_BoundariesGrid_NN_first = .true.

    call MessageNotify('M','at_Finalize',&
         & 'at_module (2019/09/29) is finalized')

  end subroutine at_Finalize

  !----------------------------------------------------------
  ! Private function
  !----------------------------------------------------------

  !----------------------------------------------------------
  !> Dirichlet 型境界条件の適用(タウ法, 2 次元配列用)
  !> 両境界での値を与える.
  subroutine at_BoundariesTau_DD_2d(at_data,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), dimension(:,0:),intent(inout)         :: at_data
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values

    real(DP), dimension(:,:), allocatable  :: alu
    integer, dimension(:), allocatable     :: kp
    real(DP), dimension(0:km,0:km)         :: tt_data
    real(DP), dimension(0:km,0:im)         :: tg_data
    real(DP), dimension(size(at_data,1))   :: value1, value2  ! 境界値

    integer(8) :: k
    save       :: alu, kp

    if ( size(at_data,2)-1 < km ) then
       call MessageNotify('E','at_BoundariesTau_DD', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(at_data,2)-1 > km ) then
       call MessageNotify('W','at_BoundariesTau_DD', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0.0D0
       value2=0.0D0
    else
       value1 = values(:,1)
       value2 = values(:,2)
    endif

    if ( at_BoundariesTau_DD_first ) then
       at_BoundariesTau_DD_first = .false.

       if (allocated(alu)) deallocate(alu)
       if (allocated(kp))  deallocate(kp)

       allocate(alu(0:km,0:km),kp(0:km))

       tt_data=0.0D0
       do k=0,km
          tt_data(k,k)=1.0D0
       enddo
       alu = tt_data

       tg_data = ag_at(tt_data)
       alu(km-1,:) = tg_data(:,0)
       alu(km,:)   = tg_data(:,im)

       call ludecomp(alu,kp)
    endif

    at_data(:,km-1) = value1
    at_data(:,km)   = value2
    at_data = lusolve(alu,kp,at_data)

  end subroutine at_BoundariesTau_DD_2d

  !----------------------------------------------------------
  !> Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)
  !> 両境界での値を与える.
  subroutine at_BoundariesTau_DD_1d(t_data,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), dimension(0:km),intent(inout)       :: t_data
    !> 境界値
    real(DP), dimension(2), intent(in), optional  :: values

    real(DP), dimension(1,0:km)                   :: at_work
    real(DP), dimension(1,2)                      :: vwork

    if (.not. present(values)) then
       vwork(1,1)=0.0D0
       vwork(1,2)=0.0D0
    else
       vwork(1,:) = values
    endif

    at_work(1,:)=t_data
    call at_BoundariesTau_DD_2d(at_work,vwork)
    t_data=at_work(1,:)

  end subroutine at_BoundariesTau_DD_1d


  !----------------------------------------------------------
  !> Dirichlet/Neumann 型境界条件の適用(タウ法, 2 次元配列用)
  !> i=0 で値, i=im で勾配の値を与える.
  subroutine at_BoundariesTau_DN_2d(at_data,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), dimension(:,0:),intent(inout)         :: at_data
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values

    real(DP), dimension(:,:), allocatable  :: alu
    integer, dimension(:), allocatable     :: kp
    real(DP), dimension(0:km,0:km)         :: tt_data
    real(DP), dimension(0:km,0:im)         :: tg_data
    real(DP), dimension(size(at_data,1))   :: value1, value2

    integer(8) :: k
    save       :: alu, kp

    if ( size(at_data,2)-1 < km ) then
       call MessageNotify('E','at_BoundariesTau_DN', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(at_data,2)-1 > km ) then
       call MessageNotify('W','at_BoundariesTau_DN', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0.0D0
       value2=0.0D0
    else
       value1 = values(:,1)
       value2 = values(:,2)
    endif

    if (  at_BoundariesTau_DN_first ) then
       at_BoundariesTau_DN_first = .false.

       if (allocated(alu)) deallocate(alu)
       if (allocated(kp)) deallocate(kp)

       allocate(alu(0:km,0:km),kp(0:km))

       tt_data = 0.0D0
       do k=0,km
          tt_data(k,k)=1.0D0
       enddo
       alu = tt_data

       tg_data = ag_at(tt_data)
       alu(km-1,:) = tg_data(:,0)
       tg_data = ag_at(at_Dx_at(tt_data))
       alu(km,:)   = tg_data(:,im)

       call ludecomp(alu,kp)
    endif

    at_data(:,km-1) = value1
    at_data(:,km)   = value2
    at_data = lusolve(alu,kp,at_data)

  end subroutine at_BoundariesTau_DN_2d

  !----------------------------------------------------------
  !> Dirichlet/Neumann 型境界条件の適用(タウ法, 1 次元配列用)
  !> i=0 で値, i=im で勾配の値を与える.
  subroutine at_BoundariesTau_DN_1d(t_data,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), dimension(0:km),intent(inout)       :: t_data
    !> 境界値

    real(DP), dimension(2), intent(in), optional  :: values
    real(DP), dimension(1,0:km)                   :: at_work
    real(DP), dimension(1,2)                      :: vwork

    if (.not. present(values)) then
       vwork(1,1)=0.0D0
       vwork(1,2)=0.0D0
    else
       vwork(1,:) = values
    endif

    at_work(1,:)=t_data
    call at_BoundariesTau_DN_2d(at_work,vwork)
    t_data=at_work(1,:)

  end subroutine at_BoundariesTau_DN_1d

  !----------------------------------------------------------
  !> Neumann/Dirichlet 型境界条件の適用(タウ法, 2 次元配列用)
  !> i=0 で勾配の値, i=im で値を与える.
  subroutine at_BoundariesTau_ND_2d(at_data,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), dimension(:,0:),intent(inout)         :: at_data
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values

    real(DP), dimension(:,:), allocatable  :: alu
    integer, dimension(:), allocatable     :: kp
    real(DP), dimension(0:km,0:km)         :: tt_data
    real(DP), dimension(0:km,0:im)         :: tg_data
    real(DP), dimension(size(at_data,1))   :: value1, value2

    integer(8) :: k
    save       :: alu, kp

    if ( size(at_data,2)-1 < km ) then
       call MessageNotify('E','at_BoundariesTau_ND', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(at_data,2)-1 > km ) then
       call MessageNotify('W','at_BoundariesTau_ND', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0.0D0
       value2=0.0D0
    else
       value1 = values(:,1)
       value2 = values(:,2)
    endif

    if ( at_BoundariesTau_ND_first ) then
       at_BoundariesTau_ND_first = .false.

       if (allocated(alu)) deallocate(alu)
       if (allocated(kp)) deallocate(kp)

       allocate(alu(0:km,0:km),kp(0:km))

       tt_data = 0.0D0
       do k=0,km
          tt_data(k,k)=1.0D0
       enddo
       alu = tt_data

       tg_data = ag_at(at_Dx_at(tt_data))
       alu(km-1,:) = tg_data(:,0)
       tg_data = ag_at(tt_data)
       alu(km,:)   = tg_data(:,im)

       call ludecomp(alu,kp)
    endif

    at_data(:,km-1) = value1
    at_data(:,km)   = value2
    at_data = lusolve(alu,kp,at_data)

  end subroutine at_BoundariesTau_ND_2d

  !----------------------------------------------------------
  !> Neumann/Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)
  !> i=0 で勾配の値, i=im で値を与える.
  subroutine at_BoundariesTau_ND_1d(t_data,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), dimension(0:km),intent(inout)       :: t_data
    !> 境界値
    real(DP), dimension(2), intent(in), optional  :: values

    real(DP), dimension(1,0:km)                   :: at_work
    real(DP), dimension(1,2)                      :: vwork

    if (.not. present(values)) then
       vwork(1,1)=0.0D0
       vwork(1,2)=0.0D0
    else
       vwork(1,:) = values
    endif

    at_work(1,:)=t_data
    call at_BoundariesTau_ND_2d(at_work,vwork)
    t_data=at_work(1,:)

  end subroutine at_BoundariesTau_ND_1d

  !----------------------------------------------------------
  !> Neumann/Dirichlet 型境界条件の適用(タウ法, 2 次元配列用)
  !> 両境界で勾配の値を与える.
  subroutine at_BoundariesTau_NN_2d(at_data,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), dimension(:,0:),intent(inout)         :: at_data
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values

    real(DP), dimension(:,:), allocatable  :: alu
    integer, dimension(:), allocatable     :: kp
    real(DP), dimension(0:km,0:km)         :: tt_data
    real(DP), dimension(0:km,0:im)         :: tg_data
    real(DP), dimension(size(at_data,1))   :: value1, value2

    integer(8) :: k
    save       :: alu, kp

    if ( size(at_data,2)-1 < km ) then
       call MessageNotify('E','at_BoundariesTau_NN', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(at_data,2)-1 > km ) then
       call MessageNotify('W','at_BoundariesTau_NN', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1 = 0.0D0
       value2 = 0.0D0
    else
       value1 = values(:,1)
       value2 = values(:,2)
    endif

    if ( at_BoundariesTau_NN_first ) then
       at_BoundariesTau_NN_first = .false.

       if (allocated(alu)) deallocate(alu)
       if (allocated(kp)) deallocate(kp)

       allocate(alu(0:km,0:km),kp(0:km))

       tt_data = 0.0D0
       do k=0,km
          tt_data(k,k)=1.0D0
       enddo
       alu = tt_data

       tg_data = ag_at(at_Dx_at(tt_data))
       alu(km-1,:) = tg_data(:,0)
       alu(km,:)   = tg_data(:,im)

       call ludecomp(alu,kp)
    endif

    at_data(:,km-1) = value1
    at_data(:,km)   = value2
    at_data = lusolve(alu,kp,at_data)

  end subroutine at_BoundariesTau_NN_2d

  !----------------------------------------------------------
  !> Neumann/Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)
  !> 両境界で勾配の値を与える.
  subroutine at_BoundariesTau_NN_1d(t_data,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), dimension(0:km),intent(inout)       :: t_data
    !> 境界値
    real(DP), dimension(2), intent(in), optional  :: values

    real(DP), dimension(1,0:km)                   :: at_work
    real(DP), dimension(1,2)                      :: vwork

    if (.not. present(values)) then
       vwork(1,1)=0.0D0
       vwork(1,2)=0.0D0
    else
       vwork(1,:) = values
    endif

    at_work(1,:)=t_data
    call at_BoundariesTau_NN_2d(at_work,vwork)
    t_data=at_work(1,:)

  end subroutine at_BoundariesTau_NN_1d

  !----------------------------------------------------------
  !> Dirichlet 型境界条件の適用(実空間での評価, 2 次元配列用)
  !> 両境界での値を与える.
  subroutine at_BoundariesGrid_DD_2d(at_data,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), dimension(:,0:),intent(inout)         :: at_data
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values

    real(DP), dimension(:,:), allocatable     :: alu
    integer, dimension(:), allocatable        :: kp
    real(DP), dimension(size(at_data,1),0:im) :: ag_data
    real(DP), dimension(0:km,0:km)            :: tt_data
    real(DP), dimension(0:km,0:im)            :: tg_data
    real(DP), dimension(size(at_data,1))      :: value1, value2  ! 境界値

    integer(8) :: k
    save       :: alu, kp

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

    if ( size(at_data,2)-1 < km ) then
       call MessageNotify('E','at_BoundariesGrid_DD', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(at_data,2)-1 > km ) then
       call MessageNotify('W','at_BoundariesGrid_DD', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0.0D0
       value2=0.0D0
    else
       value1 = values(:,1)
       value2 = values(:,2)
    endif

    if ( at_BoundariesGrid_DD_first ) then
       at_BoundariesGrid_DD_first = .false.

       if (allocated(alu)) deallocate(alu)
       if (allocated(kp)) deallocate(kp)

       allocate(alu(0:im,0:km),kp(0:im))

       tt_data = 0.0D0
       do k=0,km
          tt_data(k,k)=1.0D0
       enddo
       tg_data = ag_at(tt_data)
       alu = transpose(tg_data)
       !       alu(km-1,:) = tg_data(:,0)
       !       alu(km,:)   = tg_data(:,im)

       call ludecomp(alu,kp)
    endif

    ag_data = ag_at(at_data)
    ag_data(:,0)  = value1
    ag_data(:,im) = value2
    at_data = lusolve(alu,kp,ag_data)

  end subroutine at_BoundariesGrid_DD_2d

  !----------------------------------------------------------
  !> Dirichlet 型境界条件の適用(実空間での評価, 1 次元配列用)
  !> 両境界での値を与える.
  subroutine at_BoundariesGrid_DD_1d(t_data,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), dimension(0:km),intent(inout)       :: t_data
    !> 境界値
    real(DP), dimension(2), intent(in), optional  :: values

    real(DP), dimension(1,0:km)                   :: at_work
    real(DP), dimension(1,2)                      :: vwork           ! 境界値

    if (.not. present(values)) then
       vwork(1,1)=0.0D0
       vwork(1,2)=0.0D0
    else
       vwork(1,:) = values
    endif

    at_work(1,:)=t_data
    call at_BoundariesGrid_DD_2d(at_work,vwork)
    t_data=at_work(1,:)

  end subroutine at_BoundariesGrid_DD_1d

  !----------------------------------------------------------
  !> Dirichlet/Neumann 型境界条件の適用(実空間での評価, 2 次元配列用)
  !> i=0 で値, i=im で勾配の値を与える.
  subroutine at_BoundariesGrid_DN_2d(at_data,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), dimension(:,0:),intent(inout)         :: at_data
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values

    real(DP), dimension(:,:), allocatable     :: alu
    integer, dimension(:), allocatable        :: kp
    real(DP), dimension(size(at_data,1),0:im) :: ag_data
    real(DP), dimension(0:km,0:km)            :: tt_data
    real(DP), dimension(0:km,0:im)            :: tg_data
    real(DP), dimension(size(at_data,1))      :: value1, value2  ! 境界値

    integer(8) :: k
    save       :: alu, kp

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

    if ( size(at_data,2)-1 < km ) then
       call MessageNotify('E','at_BoundariesGrid_DN', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(at_data,2)-1 > km ) then
       call MessageNotify('W','at_BoundariesGrid_DN', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0.0D0
       value2=0.0D0
    else
       value1 = values(:,1)
       value2 = values(:,2)
    endif

    if ( at_BoundariesGrid_DN_first ) then
       at_BoundariesGrid_DN_first = .false.

       if (allocated(alu)) deallocate(alu)
       if (allocated(kp)) deallocate(kp)

       allocate(alu(0:im,0:km),kp(0:im))

       tt_data = 0.0D0
       do k=0,km
          tt_data(k,k)=1.0D0
       enddo
       tg_data = ag_at(tt_data)
       alu = transpose(tg_data)

       tg_data = ag_at(at_dx_at(tt_data))
       alu(im,:) = tg_data(:,im)

       call ludecomp(alu,kp)
    endif

    ag_data = ag_at(at_data)
    ag_data(:,0)  = value1
    ag_data(:,im) = value2
    at_data = lusolve(alu,kp,ag_data)

  end subroutine at_BoundariesGrid_DN_2d

  !----------------------------------------------------------------------
  !> Dirichlet/Neumann 型境界条件の適用(実空間での評価, 1 次元配列用)
  !> i=0 で値, i=im で勾配の値を与える.
  subroutine at_BoundariesGrid_DN_1d(t_data,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), dimension(0:km),intent(inout)       :: t_data
    !> 境界値
    real(DP), dimension(2), intent(in), optional  :: values

    real(DP), dimension(1,0:km)                   :: at_work
    real(DP), dimension(1,2)                      :: vwork           ! 境界値

    if (.not. present(values)) then
       vwork(1,1)=0.0D0
       vwork(1,2)=0.0D0
    else
       vwork(1,:) = values
    endif

    at_work(1,:)=t_data
    call at_BoundariesGrid_DN_2d(at_work,vwork)
    t_data=at_work(1,:)

  end subroutine at_BoundariesGrid_DN_1d

  !----------------------------------------------------------------------
  !> Neumann/Dirichlet 型境界条件の適用(実空間での評価, 2 次元配列用)
  !> i=0 で勾配の値, i=im で値を与える.
  !----------------------------------------------------------------------
  subroutine at_BoundariesGrid_ND_2d(at_data,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), dimension(:,0:),intent(inout)         :: at_data
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values

    real(DP), dimension(:,:), allocatable     :: alu
    integer, dimension(:), allocatable        :: kp
    real(DP), dimension(size(at_data,1),0:im) :: ag_data
    real(DP), dimension(0:km,0:km)            :: tt_data
    real(DP), dimension(0:km,0:im)            :: tg_data
    real(DP), dimension(size(at_data,1))      :: value1, value2  ! 境界値

    integer(8) :: k
    save       :: alu, kp

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

    if ( size(at_data,2)-1 < km ) then
       call MessageNotify('E','at_BoundariesGrid_ND', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(at_data,2)-1 > km ) then
       call MessageNotify('W','at_BoundariesGrid_DD', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0.0D0
       value2=0.0D0
    else
       value1 = values(:,1)
       value2 = values(:,2)
    endif

    if ( at_BoundariesGrid_ND_first ) then
       at_BoundariesGrid_ND_first = .false.

       if (allocated(alu)) deallocate(alu)
       if (allocated(kp)) deallocate(kp)

       allocate(alu(0:im,0:km),kp(0:im))

       tt_data = 0.0D0
       do k=0,km
          tt_data(k,k)=1.0D0
       enddo
       tg_data = ag_at(tt_data)
       alu = transpose(tg_data)

       tg_data = ag_at(at_dx_at(tt_data))
       alu(0,:)  = tg_data(:,0)

       call ludecomp(alu,kp)
    endif

    ag_data = ag_at(at_data)
    ag_data(:,0)  = value1
    ag_data(:,im) = value2
    at_data = lusolve(alu,kp,ag_data)

  end subroutine at_BoundariesGrid_ND_2d

  !----------------------------------------------------------------------
  !> Neumann 型境界条件の適用(実空間での評価, 1 次元配列用)
  !> i=0 で勾配の値, i=im で値を与える.
  !----------------------------------------------------------------------
  subroutine at_BoundariesGrid_ND_1d(t_data,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), dimension(0:km),intent(inout)       :: t_data
    !> 境界値
    real(DP), dimension(2), intent(in), optional  :: values

    real(DP), dimension(1,0:km)                   :: at_work
    real(DP), dimension(1,2)                      :: vwork           ! 境界値

    if (.not. present(values)) then
       vwork(1,1)=0.0D0
       vwork(1,2)=0.0D0
    else
       vwork(1,:) = values
    endif

    at_work(1,:)=t_data
    call at_BoundariesGrid_ND_2d(at_work,vwork)
    t_data=at_work(1,:)

  end subroutine at_BoundariesGrid_ND_1d

  !----------------------------------------------------------------------
  !> Neumann 型境界条件の適用(実空間での評価, 2 次元配列用)
  !> 両境界で勾配の値を与える.
  !----------------------------------------------------------------------
  subroutine at_BoundariesGrid_NN_2d(at_data,values)
    !> 境界条件を適用するチェビシェフデータ(m,0:km)
    real(DP), dimension(:,0:),intent(inout)         :: at_data
    !> 境界値(m,2)
    real(DP), dimension(:,:), intent(in), optional  :: values


    real(DP),   dimension(:,:), allocatable     :: alu
    integer,    dimension(:), allocatable       :: kp
    real(DP),   dimension(size(at_data,1),0:im) :: ag_data
    real(DP),   dimension(0:km,0:km)            :: tt_data
    real(DP),   dimension(0:km,0:im)            :: tg_data
    real(DP),   dimension(size(at_data,1))      :: value1, value2  ! 境界値

    integer(8) :: k
    save       :: alu, kp

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

    if ( size(at_data,2)-1 < km ) then
       call MessageNotify('E','at_BoundariesGrid_NN', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(at_data,2)-1 > km ) then
       call MessageNotify('W','at_BoundariesGrid_NN', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0.0D0
       value2=0.0D0
    else
       value1 = values(:,1)
       value2 = values(:,2)
    endif

    if ( at_BoundariesGrid_NN_first ) then
       at_BoundariesGrid_NN_first = .false.

       if (allocated(alu)) deallocate(alu)
       if (allocated(kp)) deallocate(kp)

       allocate(alu(0:im,0:km),kp(0:im))

       tt_data = 0.0D0
       do k=0,km
          tt_data(k,k)=1.0D0
       enddo
       tg_data = ag_at(tt_data)
       alu = transpose(tg_data)

       tg_data = ag_at(at_dx_at(tt_data))
       alu(0,:)  = tg_data(:,0)
       alu(im,:) = tg_data(:,im)

       call ludecomp(alu,kp)
    endif

    ag_data = ag_at(at_data)
    ag_data(:,0)  = value1
    ag_data(:,im) = value2
    at_data = lusolve(alu,kp,ag_data)

  end subroutine at_BoundariesGrid_NN_2d

  !----------------------------------------------------------------------
  !> Neumann 型境界条件の適用(実空間での評価, 1 次元配列用)
  !> 両境界で勾配の値を与える.
  !----------------------------------------------------------------------
  subroutine at_BoundariesGrid_NN_1d(t_data,values)
    !> 境界条件を適用するチェビシェフデータ(0:km)
    real(DP), dimension(0:km),intent(inout)       :: t_data
    !> 境界値
    real(DP), dimension(2), intent(in), optional  :: values
    real(DP), dimension(1,0:km)                   :: at_work
    real(DP), dimension(1,2)                      :: vwork  ! 境界値

    if (.not. present(values)) then
       vwork(1,1)=0.0D0
       vwork(1,2)=0.0D0
    else
       vwork(1,:) = values
    endif

    at_work(1,:)=t_data
    call at_BoundariesGrid_NN_2d(at_work,vwork)
    t_data=at_work(1,:)

  end subroutine at_BoundariesGrid_NN_1d

end module at_module
