at_module.f90

Go to the documentation of this file.

! -*- mode: f90; coding: utf-8 -*-
!----------------------------------------------------------------------
! Copyright(c) 2019 SPMDODEL Development Group. All rights reserved.
!
! 履歴  2019/03/28 竹広真一
!       2019/09/28 佐々木洋平
!       2019/10/09 佐々木洋平
!----------------------------------------------------------------------

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

  !----------------------------------------------------------------------
  interface at_boundaries_dn
     module procedure at_boundariestau_dn_1d
     module procedure at_boundariestau_dn_2d
  end interface at_boundaries_dn

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

  !----------------------------------------------------------------------
  interface at_boundariestau_dn
     module procedure at_boundariestau_dn_1d
     module procedure at_boundariestau_dn_2d
  end interface at_boundariestau_dn

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

  !----------------------------------------------------------------------
  interface at_boundariesgrid_dn
     module procedure at_boundariesgrid_dn_1d
     module procedure at_boundariesgrid_dn_2d
  end interface at_boundariesgrid_dn

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

  !----------------------------------------------------------------------
  real(DP), allocatable :: g_x(:)
  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

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

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

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

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

  !-----------------------------------------------------------------------
  function at_dx_at(at_data)

    real(DP), dimension(:,0:), intent(in)                    :: at_data
    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


  !-----------------------------------------------------------------------
  function t_dx_t(t_data)

    real(DP), dimension(:), intent(in)   :: t_data
    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

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

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

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

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

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

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

  !-----------------------------------------------------------------------
  subroutine at_boundaries(at_data,cond,values)

    real(DP), intent(INOUT)                  :: at_data(:,:)
    character(len=2), intent(IN), optional  :: cond
    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

  !---------------------------------------------------------------------
  subroutine t_boundaries(t_data,cond,values)

    real(DP), intent(INOUT)                  :: t_data(:)
    character(len=2), intent(IN), optional  :: cond
    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

  !---------------------------------------------------------------------
  subroutine at_boundariesgrid(at_data,cond,values)

    real(DP), intent(INOUT)                  :: at_data(:,:)
    character(len=2), intent(IN), optional  :: cond
    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

  !---------------------------------------------------------------------
  subroutine t_boundariesgrid(t_data,cond,values)

    real(DP), intent(INOUT)                  :: t_data(:)
    character(len=2), intent(IN), optional  :: cond
    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
  !----------------------------------------------------------

  !----------------------------------------------------------
  subroutine at_boundariestau_dd_2d(at_data,values)

    real(DP), dimension(:,0:),intent(inout)         :: at_data
    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

  !----------------------------------------------------------
  subroutine at_boundariestau_dd_1d(t_data,values)

    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


  !----------------------------------------------------------
  subroutine at_boundariestau_dn_2d(at_data,values)

    real(DP), dimension(:,0:),intent(inout)         :: at_data
    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

  !----------------------------------------------------------
  subroutine at_boundariestau_dn_1d(t_data,values)

    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

  !----------------------------------------------------------
  subroutine at_boundariestau_nd_2d(at_data,values)

    real(DP), dimension(:,0:),intent(inout)         :: at_data
    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

  !----------------------------------------------------------
  subroutine at_boundariestau_nd_1d(t_data,values)

    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

  !----------------------------------------------------------
  subroutine at_boundariestau_nn_2d(at_data,values)

    real(DP), dimension(:,0:),intent(inout)         :: at_data
    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

  !----------------------------------------------------------
  subroutine at_boundariestau_nn_1d(t_data,values)

    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

  !----------------------------------------------------------
  subroutine at_boundariesgrid_dd_2d(at_data,values)

    real(DP), dimension(:,0:),intent(inout)         :: at_data
    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

  !----------------------------------------------------------
  subroutine at_boundariesgrid_dd_1d(t_data,values)

    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

  !----------------------------------------------------------
  subroutine at_boundariesgrid_dn_2d(at_data,values)

    real(DP), dimension(:,0:),intent(inout)         :: at_data
    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

  !----------------------------------------------------------------------
  subroutine at_boundariesgrid_dn_1d(t_data,values)

    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

  !----------------------------------------------------------------------
  !----------------------------------------------------------------------
  subroutine at_boundariesgrid_nd_2d(at_data,values)

    real(DP), dimension(:,0:),intent(inout)         :: at_data
    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

  !----------------------------------------------------------------------
  !----------------------------------------------------------------------
  subroutine at_boundariesgrid_nd_1d(t_data,values)

    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

  !----------------------------------------------------------------------
  !----------------------------------------------------------------------
  subroutine at_boundariesgrid_nn_2d(at_data,values)

    real(DP), dimension(:,0:),intent(inout)         :: at_data
    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

  !----------------------------------------------------------------------
  !----------------------------------------------------------------------
  subroutine at_boundariesgrid_nn_1d(t_data,values)

    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