function a_Int_ag(ag)
      real(8), dimension(:,0:), intent(in)     :: ag
      real(8), dimension(size(ag,1))           :: a_Int_ag
      integer :: 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 ag_at(at_data)  ! スペクトル -> 台形格子(2 次元)

    double precision, dimension(:,:), intent(in)      :: at_data
    double precision, dimension(size(at_data,1),0:im)  :: ag_at

    real(8), dimension(size(at_data,1)*im) :: y
    integer :: m

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

    ag_at = 0.0D0
    ag_at(:,0:km)=at_data
    call fttctb(m,im,ag_at,y,it,t)

  end function ag_at

  function at_Dx_at(at_data)
    real(8), dimension(:,0:), intent(in)                    :: at_data
    real(8), dimension(size(at_data,1),0:size(at_data,2)-1) :: at_Dx_at

    integer :: m, k
    integer :: nm, kmax

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

    if ( kmax == im ) then
       do m=1,nm
          at_Dx_at(m,kmax)   = 0. 
          at_Dx_at(m,kmax-1) = 2 * km * at_data(m,kmax) /2
       enddo
    else
       do m=1,nm
          at_Dx_at(m,kmax)   = 0. 
          at_Dx_at(m,kmax-1) = 2 * km * at_data(m,kmax) 
          ! スタートはグリッド対応最大波数未満. Factor 1/2 不要
       enddo
    endif

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

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

  end function at_Dx_at

  subroutine at_Initial(i,k,xmin,xmax)
    integer, intent(in) :: i, k
    real(8), intent(in) :: xmin,xmax

    integer :: 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) /= 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(t(3*im))
    call fttcti(im,it,t)

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

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

  end subroutine at_Initial