! -*- mode: f90; coding: utf-8 -*-
!----------------------------------------------------------------------
! Copyright(c) 2023-2025 SPMDODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @copyright Copyright (C) SPMODEL Development Group, 2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @brief 2 次元水路領域: Fourier 展開 + Chebyshev 展開法
!>
!> @details
!>
!> spml/et_mpi_module モジュールは 2 次元水路領域での流体運動を
!> スペクトル法により数値計算を実行するための Fortran90 関数を提供する.
!> 周期的な境界条件を扱うための X 方向へのフーリエ変換と境界壁を扱うための
!> Y 方向のチェビシェフ変換を用いる場合のスペクトル計算のためのさまざまな
!> 関数を提供する.
!>
!> 内部で ae_module, at_module を用いている.
!> 最下部ではフーリエ変換およびチェビシェフ変換のエンジンとして
!> ISPACK3/FXPACKを用いている.
!>
!> ## 関数・変数の名前と型について
!>
!> ### 命名法
!>
!> * 関数名の先頭 (et_, yx_, x_, y_) は, 返す値の形を示している.
!>   et_ :: 2次元スペクトルデータ
!>   yx_ :: 2 次元格子点データ
!>   x_  :: X 方向 1 次元格子点データ
!>   y_  :: Y 方向 1 次元格子点データ
!>
!> * 関数名の間の文字列(Dx, Dy, Lapla, LaplaInv, Jacobian)は,
!>   その関数の作用を表している.
!>
!> * 関数名の最後 (_et_et, _et, _yx, _x, _y) は, 入力変数のスペクトルデータ
!>   および格子点データであることを示している.
!>   _et    :: 2次元スペクトルデータ
!>   _et_et :: 2 つの2次元スペクトルデータ
!>   _yx    :: 2 次元格子点データ
!>   _x     :: X 方向 1 次元格子点データ
!>   _y     :: Y 方向 1 次元格子点データ
!>
!> ### 各データの種類の説明
!>
!> * yx : 2 次元格子点データ.
!>   * 変数の種類と次元は real(DP), dimension(0:jm,0:im-1).
!>   * im, jm はそれぞれ X, Y 座標の格子点数であり, サブルーチン
!>     et_initial にてあらかじめ設定しておく.
!>   * 第 1 次元が Y 座標の格子点位置番号, 第 2 次元が X 座標の
!>     格子点位置番号である (X, Y の順ではない)ことに注意.
!>
!> * et : 2 次元スペクトルデータ.
!>   * 変数の種類と次元は real(DP), dimension(-km:km,0:lm).
!>   * km, lm はそれぞれ X, Y 方向の最大波数であり, サブルーチン
!>     et_initial にてあらかじめ設定しておく.
!>   * スペクトルデータの格納のされ方については...
!>
!> * x, y : X, Y 方向 1 次元格子点データ.
!>   * 変数の種類と次元はそれぞれ real(DP), dimension(0:im-1)
!>     および real(DP), dimension(0:jm).
!>
!> * e, t : 1 次元スペクトルデータ.
!>   * 変数の種類と次元は real(DP), dimension(-km:km)
!>     および real(DP), dimension(-lm:lm).
!>
!> * ax, ay : 1 次元格子点データの並んだ 2 次元配列.
!>   * 変数の種類と次元は real(DP), dimension(:,0:im-1)
!>     および real(DP), dimension(:,0:jm).
!>
!> * ae, at : 1 次元スペクトルデータの並んだ 2 次元配列.
!>   * 変数の種類と次元は real(DP), dimension(:,-km:km)
!>     および real(DP), dimension(:,0:lm).
!>
!> * et_ で始まる関数が返す値はスペクトルデータに同じ.
!>
!> * yx_ で始まる関数が返す値は 2 次元格子点データに同じ.
!>
!> * x_, y_ で始まる関数が返す値は 1 次元格子点データに同じ.
!>
!> * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
!>   微分などを作用させたデータをスペクトル変換したものことである.
!----------------------------------------------------------------------
!履歴  2023/03/04  竹広真一
!      2025/11/20  竹広真一  f_e, e_f, y_v を追加
!----------------------------------------------------------------------
module et_mpi_module
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! et_mpi_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  ! kf_k           :: X 方向波数の位置情報
  !
  !==== 座標変数
  !
  ! x_X, y_Y     ::  格子点座標(X,Y座標)を格納した 1 次元配列
  ! x_X_Weight, y_Y_Weight ::  重み座標を格納した 1 次元配列
  ! yx_X, yx_Y   :: 格子点データの XY 座標(X,Y)(格子点データ型 2 次元配列)
  !
  ! v_Y          ::  格子点座標(Y座標, 分割格子)を格納した 1 次元配列
  ! v_Y_Weight   ::  重み座標を格納した 1 次元配列(分割格子)
  ! vx_X, vx_Y   :: 格子点データの XY 座標(X,Y)(格子点データ型 2 次元配列, 分割格子)
  !
  !==== 基本変換
  !
  ! vx_ft :: スペクトルデータから格子データへの変換
  ! ft_vx :: 格子データからスペクトルデータへの変換
  ! vx_ft, x_e :: X 方向のスペクトルデータから格子データへの変換
  ! fy_ft, y_t :: Y 方向のスペクトルデータから格子データへの変換
  ! fy_yx, e_x :: X 方向の格子点データからスペクトルデータへの変換
  ! ft_fy, t_y :: Y 方向の格子点データからスペクトルデータへの変換
  !
  !==== 微分
  !
  ! ft_Lapla_ft  :: スペクトルデータにラプラシアンを作用させる
  ! ft_Dx_ft, ae_Dx_ae, e_Dx_e :: スペクトルデータに X 微分を作用させる
  ! ft_Dy_ft, at_Dy_at, t_Dy_t :: スペクトルデータに Y 微分を作用させる
  ! ft_Jacobian_ft_ft :: 2 つのスペクトルデータからヤコビアンを計算する
  !
  !==== 境界値問題
  !
  ! ft_Boundaries  :: ディリクレ, ノイマン境界条件の適用
  ! ft_LaplaInv_ft :: スペクトルデータにラプラシアンの逆変換を作用させる
  ! ft_Vor2Strm_ft :: 渦度から流線を計算する
  ! fy_Vor2Strm_fy :: 渦度から流線を計算する
  !
  !==== 積分・平均
  !
  ! IntYX_vx, AvrYX_vx   :: 2 次元格子点データの全領域積分および平均
  ! v_IntX_vx, v_AvrX_vx :: 2 次元格子点データの X 方向積分および平均
  ! IntX_x, AvrX_x       :: 1 次元(X)格子点データの X 方向積分および平均
  ! x_IntY_vx, x_AvrY_vx :: 2 次元格子点データの Y 方向積分および平均
  ! IntY_y, AvrY_v       :: 1 次元(Y)格子点データの Y 方向積分および平均
  !
  !==== 補間
  !
  ! Interpolate_ft       :: スペクトルデータからの特定の座標の値の補間
  !
  use dc_types, only: DP
  use dc_message, only: MessageNotify
  use mpi
  use lumatrix
  use spml_utils, only: pi
  use ae_module, only : x_X => g_X, x_X_weight => g_X_Weight, &
    ae_initial, &
    e_x => e_g, ae_ax => ae_ag, &
    x_e => g_e, ax_ae => ag_ae, &
    Interpolate_e, a_Interpolate_ae
  use at_module, only : y_Y => g_X, y_Y_Weight => g_X_Weight, &
    at_initial, &
    at_ay => at_ag, t_y => t_g, &
    ay_at => ag_at, y_t => g_t, &
    t_Dy_t => t_Dx_t, at_Dy_at => at_Dx_at, &
    at_Boundaries_DD, at_Boundaries_DN, at_Boundaries_ND, at_Boundaries_NN, &
    ft_BoundariesGrid => at_BoundariesGrid, &
    Interpolate_t, a_Interpolate_at
  implicit none

  private

  public ks, ke, kc
  public js, je, jc
  public kf_k

  public et_mpi_Initial                                   ! 初期化

  public x_X, y_Y, x_X_Weight, y_Y_Weight, yx_X, yx_Y     ! 座標変数
  public v_Y, v_Y_Weight, vx_X, vx_Y                      ! 座標変数

  public vx_ft, ft_vx                                     ! 基本変換
  public ve_fy, fy_ve                                     ! 基本変換
  public ft_fy, fy_ft                                     ! 基本変換
  public vx_fy, fy_vx                                     ! 基本変換
  public ae_af, af_ae                                     ! 基本変換
  public ay_av, av_ay                                     ! 基本変換
  public e_f, f_e                                         ! 基本変換
  public y_v                                              ! 基本変換
!!$  public vx_ey, ey_yx                                     ! 基本変換
!!$  public ey_et, et_ey                                     ! 基本変換
  public e_x, x_e, ae_ax, ax_ae                           ! 基本変換
  public t_y, y_t, at_ay, ay_at                           ! 基本変換

  public ft_Dx_ft                       ! 微分
  public ft_Dy_ft                       ! 微分
  public ft_Lapla_ft                                      ! 微分
!!$  public ft_Dx_ft, f_Dx_e, ae_Dx_ae                       ! 微分
!!$  public ft_Dy_ft, t_Dy_t, at_Dy_at                       ! 微分

  public ft_Jacobian_ft_ft                                ! 非線形計算

  public ft_Boundaries                                    ! 境界値問題
!!$  public at_Boundaries_DD, at_Boundaries_DN               ! 境界値問題
!!$  public at_Boundaries_ND, at_Boundaries_NN               ! 境界値問題
  public ft_BoundariesGrid                                ! 境界値問題
  public ft_LaplaInv_ft
  public fy_Vor2Strm_fy                                   ! 境界値問題
  public ft_Vor2Strm_ft                                   ! 境界値問題
  public ft_Vor2Strm2_ft                                  ! 境界値問題
  public ft_Vor2Strm3_ft                                  ! 境界値問題

  public IntYX_vx, v_IntX_vx, x_IntY_vx, IntX_x, IntY_v   ! 積分
  public AvrYX_vx, v_AvrX_vx, x_AvrY_vx, AvrX_x, AvrY_v   ! 平均

  public Interpolate_ft                                   ! 補間

  !> 格子点データの Y 座標: 格子点データ型分割(Y)
  real(DP), dimension(:), allocatable :: v_Y
  real(DP), dimension(:), allocatable :: v_Y_Weight

  !> 格子点データの X 座標: 格子点データ型 2 次元配列(X,Y)
  real(DP), dimension(:,:), allocatable :: yx_X
  real(DP), dimension(:,:), allocatable :: vx_X
  !> 格子点データの Y 座標: 格子点データ型 2 次元配列(X,Y)
  real(DP), dimension(:,:), allocatable :: yx_Y
  real(DP), dimension(:,:), allocatable :: vx_Y

  !> X方向格子点
  integer            :: im
  !> Y方向格子点
  integer            :: jm
  !> X方向切断波数
  integer            :: km
  !> Y方向切断波数
  integer            :: lm
  !> X 方向計算領域
  real(DP)           :: xl
  !> Y 方向計算領域
  real(DP)           :: yl

  real(DP) :: xmin, xmax                                 ! X 座標範囲
  real(DP) :: ymin, ymax                                 ! Y 座標範囲

  integer, dimension(:), allocatable :: ke              ! X 最大波数
  integer, dimension(:), allocatable :: ks              ! X 最小波数
  integer, dimension(:), allocatable :: kc              ! X 波数の数
  integer, dimension(:), allocatable :: je              ! Y 座標最大格子点
  integer, dimension(:), allocatable :: js              ! Y 座標最小格子点
  integer, dimension(:), allocatable :: jc              ! Y 座標格子点数
  integer :: np                                         ! MPI プロセス数
  integer :: ip                                         ! MPI プロセス ID

  save im, jm, km, lm, xl, yl
  save ks, ke, kc, js, je, jc, np, ip
  save xmin, xmax, ymin, ymax

contains
  !----------------------------------------------------------
  !> スペクトル変換の格子点数, 波数, 領域の大きさを設定する.
  !>
  !> 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
  !> 初期設定をしなければならない.
  !----------------------------------------------------------
  subroutine et_mpi_Initial(i,j,k,l,xmin,xmax,ymin,ymax)
    !> 格子点数(X)
    integer,intent(in) :: i
    !> 格子点数(Y)
    integer,intent(in) :: j
    !> 切断波数(X)
    integer,intent(in) :: k
    !> 切断波数(Y)
    integer,intent(in) :: l
    !> X 座標範囲
    real(DP),intent(in) :: xmin, xmax
    !> Y 座標範囲
    real(DP),intent(in) :: ymin, ymax

    integer :: kint, kmod
    integer :: jint, jmod
    integer :: ierr, iip

    im = i
    jm = j
    km = k
    lm = l
    xl = xmax-xmin
    yl = ymax-ymin

    CALL MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)
    CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)

    ! x 波数 MPI分割情報
    allocate(ks(0:np-1))
    allocate(ke(0:np-1))
    allocate(kc(0:np-1))

    kint = km/np
    kmod = mod(km,np)
    kc   = kint
    do iip=0,kmod
       kc(iip) = kc(iip) + 1
    end do

    ks(0)=0 ; ke(0) = ks(0) + kc(0) -1
    do iip=1,np-1
       ks(iip) = ks(iip-1) + kc(iip-1)
       ke(iip) = ks(iip) + kc(iip) - 1
    end do

    ! Y 軸 MPI分割情報
    allocate(js(0:np-1))
    allocate(je(0:np-1))
    allocate(jc(0:np-1))

    jint = (jm+1)/np
    jmod = mod(jm+1,np)
    jc   = jint
    do iip=0,jmod-1
       jc(iip) = jc(iip) + 1
    end do

    js(0)=0 ; je(0) = js(0) + jc(0) -1
    do iip=1,np-1
       js(iip) = js(iip-1) + jc(iip-1)
       je(iip) = js(iip) + jc(iip) - 1
    end do

    call ae_initial(im,km,xmin,xmax)
    call at_initial(jm,lm,ymin,ymax)

    allocate(v_Y(jc(ip)))
    allocate(v_Y_Weight(jc(ip)))

    allocate(yx_X(0:jm,0:im-1),yx_Y(0:jm,0:im-1))
    allocate(vx_X(jc(ip),0:im-1),vx_Y(jc(ip),0:im-1))

    v_Y        = y_Y(js(ip):je(ip))
    v_Y_Weight = y_Y_Weight(js(ip):je(ip))

    yx_X = spread(x_X,1,jm+1)
    yx_Y = spread(y_Y,2,im)

    vx_X = spread(x_X,1,jc(ip))
    vx_Y = spread(v_Y,2,im)

    call MessageNotify('M','et_mpi_initial', &
      & 'et_mpi_module (2025/11/20) is initialized')
  end subroutine et_mpi_Initial

  !--------------- 波数位置問い合わせ -----------------
  function kf_k(k)
    !
    ! スペクトルデータ ft(2*kc,0:lm) での
    ! X 方向波数 k の位置情報(第 1 添え字)を調べる.
    !
    ! X 方向波数については ks,ks+1,... ke, -ke,-ke+1,...,-ks
    ! の順に格納されている.
    !
    integer, intent(IN)  ::  k        ! X 方向波数
    integer              ::  kf_k     ! スペクトルデータ第 2 添え字

    if ( abs(k) < ks(ip) .OR. ke(ip) < abs(k) ) then
       kf_k = -1
    else if ( k >= 0 ) then
       kf_k = k-ks(ip)+1
    else
       kf_k = ks(ip)+2*kc(ip)-abs(k)
    endif

  end function kf_k

  !----------------------------------------------------------
  !> スペクトルデータから格子データへ変換する.
  !----------------------------------------------------------
  function vx_ft(ft)
    !> スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm), intent(in)  :: ft
    !> 格子点データ
    real(DP), allocatable, dimension(:,:)           :: vx_ft

!!     real(DP), dimension(2*kc(ip),0:jm)              :: fy
!!     real(DP), dimension(js(ip):je(ip),-km:km)       :: ve
!!
!!     fy = ay_at(ft)
!!     ve = ve_fy(fy)
!!     vx_ft = ax_ae(ve)
    allocate(vx_ft(js(ip):je(ip),0:im-1))
    vx_ft = ax_ae(ve_fy(ay_at(ft)))

  end function vx_ft

  !----------------------------------------------------------
  ! 格子データからスペクトルデータへ変換する.
  !----------------------------------------------------------
  function ft_vx(vx)
    !> 格子点データ
    real(DP), dimension(js(ip):je(ip),0:im-1), intent(in) :: vx
    !> スペクトルデータ
    real(DP), allocatable, dimension(:,:)                 :: ft_vx

!!     real(DP), dimension(2*kc(ip),0:jm)              :: fy
!!     real(DP), dimension(js(ip):je(ip),-km:km)       :: ve
!!
!!     ve = ae_ax(vx)
!!     fy = fy_ve(ve)
!!     ft_vx = at_ay(fy)
    allocate(ft_vx(2*kc(ip),0:lm))
    ft_vx = at_ay(fy_ve(ae_ax(vx)))
  end function ft_vx

  !----------------------------------------------------------
  !> X スペクトル・Y 格子の転置
  !----------------------------------------------------------
  function ve_fy(fy)
    !> スペクトルデータ
    real(DP), dimension(2*kc(ip),0:jm), intent(in)  :: fy

    !> 格子点データ
    real(DP), dimension(js(ip):je(ip),-km:km)       :: ve_fy

    integer :: j, k, iip
    integer :: ierr
    integer :: isndc(0:np-1)
    integer :: isdsp(0:np-1)
    integer :: ircvc(0:np-1)
    integer :: irdsp(0:np-1)
    real(DP), allocatable :: aa_work(:,:)
    integer :: kcm
    integer :: ind1, ind2

    kcm = maxval(kc)
    allocate(aa_work(kcm*2*jc(ip),0:np-1))
    aa_work = 0.0D0

    isndc = kc(ip)*2*jc
    ircvc = jc(ip)*2*kc
    do iip=0,np-1
       isdsp(iip) = sum(isndc(0:iip-1))
       irdsp(iip) = kcm*2*jc(ip)*iip
    end do

    call MPI_AlltoAllV(fy,isndc,isdsp,MPI_REAL8,&
         aa_work,ircvc,irdsp,MPI_REAL8,MPI_COMM_WORLD,IERR)

    do iip=0,np-1
       do k=1,kc(iip)
          do j=1,jc(ip)
             ind1 = k           + kc(iip)*2*(j-1)
             ind2 = k + kc(iip) + kc(iip)*2*(j-1)
             ve_fy(js(ip)+j-1,-ke(iip)-1+k) = aa_work(ind2,iip)
             ve_fy(js(ip)+j-1, ks(iip)-1+k) = aa_work(ind1,iip)
          end do
       end do
    end do

    deallocate(aa_work)

  end function ve_fy

  function fy_ve(ve)
    real(DP), dimension(jc(ip),-km:km),intent(IN) :: ve
    real(DP), dimension(2*kc(ip),0:jm)            :: fy_ve

    integer :: j, k, iip
    integer :: ierr
    integer :: isndc(0:np-1)
    integer :: isdsp1(0:np-1)
    integer :: isdsp2(0:np-1)
    integer :: ircvc(0:np-1)
    integer :: irdsp(0:np-1)
    real(DP), allocatable :: aa_work1(:,:)
    real(DP), allocatable :: aa_work2(:,:)
    integer :: jcm
    integer :: ind

    jcm = maxval(jc)
    allocate(aa_work1(jcm*kc(ip),0:np-1))
    allocate(aa_work2(jcm*kc(ip),0:np-1))
    aa_work1 = 0.0D0 ; aa_work2 = 0.0D0

    isndc = jc(ip)*kc
    ircvc = kc(ip)*jc
    do iip=0,np-1
       isdsp1(iip) = (km + ks(iip))*jc(ip)
       isdsp2(iip) = (km - ke(iip))*jc(ip)
       irdsp(iip) = jcm*kc(ip)*iip
    end do

    call MPI_AlltoAllV(ve,isndc,isdsp1,MPI_REAL8,&
         aa_work1,ircvc,irdsp,MPI_REAL8,MPI_COMM_WORLD,IERR)
    call MPI_AlltoAllV(ve,isndc,isdsp2,MPI_REAL8,&
         aa_work2,ircvc,irdsp,MPI_REAL8,MPI_COMM_WORLD,IERR)

    do iip=0,np-1
       do j=1,jc(iip)
          do k=1,kc(ip)
             ind = j           + jc(iip)*(k-1)
             fy_ve(k,js(iip)+j-1)        = aa_work1(ind,iip)
             fy_ve(kc(ip)+k,js(iip)+j-1) = aa_work2(ind,iip)
          end do
       end do
    end do

    deallocate(aa_work1)
    deallocate(aa_work2)

  end function fy_ve

  !----------------------------------------------------------
  !> X スペクトルの分散・統合
  !----------------------------------------------------------
  function ae_af(af)
    real(DP), dimension(:,:),intent(IN)       :: af
    real(DP), dimension(size(af,1),-km:km)    :: ae_af

    integer :: icm
    integer :: ircvc(0:np-1)
    integer :: irdsp1(0:np-1)
    integer :: irdsp2(0:np-1)
    integer :: ierr, ipp

    icm = size(af,1)
    do ipp=0,np-1
       ircvc(ipp) = icm*kc(ipp)
    end do
    irdsp1(0) = 0
    do ipp=1,np-1
       irdsp1(ipp) = sum(ircvc(0:ipp-1))
    end do
    do ipp=0,np-1
       irdsp2(ipp) = icm*(km+1) - sum(ircvc(0:ipp))
    end do

    call MPI_AllGatherV(af(:,kc(ip)+1:2*kc(ip)),icm*kc(ip),MPI_REAL8,&
         ae_af(:,-km:0),ircvc,irdsp2,MPI_REAL8,MPI_COMM_WORLD,IERR)
    call MPI_AllGatherV(af(:,1:kc(ip)),icm*kc(ip),MPI_REAL8,&
         ae_af(:,0:km),ircvc,irdsp1,MPI_REAL8,MPI_COMM_WORLD,IERR)

  end function ae_af

  function af_ae(ae)
    real(DP), dimension(:,-km:),intent(IN)      :: ae
    real(DP), dimension(size(ae,1),2*kc(ip))    :: af_ae

    integer :: k, kk

    do kk=1,kc(ip)
       k = ks(ip)+kk-1
       af_ae(:,kk) = ae(:,k)
       af_ae(:,2*kc(ip)-kk+1) = ae(:,-k)
    end do

  end function af_ae

  function e_f(f)
    real(DP), dimension(2*kc(ip)),intent(IN) :: f
    real(DP), dimension(-km:km)              :: e_f

    real(DP), dimension(1,2*kc(ip)) :: af
    real(DP), dimension(1,-km:km)   :: ae

    af(1,:) = f
    ae = ae_af(af)
    e_f = ae(1,:)
  end function e_f

  function f_e(e)
    real(DP), dimension(-km:km),intent(IN) :: e
    real(DP), dimension(2*kc(ip))          :: f_e

    real(DP), dimension(1,2*kc(ip)) :: af
    real(DP), dimension(1,-km:km)   :: ae

    ae(1,:) = e
    af = af_ae(ae)
    f_e = af(1,:)
  end function f_e

  !----------------------------------------------------------
  !> Y 格子の分散・統合
  !----------------------------------------------------------
  function ay_av(av)
    real(DP), dimension(:,:),intent(IN)       :: av
    real(DP), dimension(size(av,1),0:jm)      :: ay_av

    integer :: icm
    integer :: ircvc(0:np-1)
    integer :: irdsp(0:np-1)
    integer :: ierr, ipp

    icm = size(av,1)
    do ipp=0,np-1
       ircvc(ipp) = icm*jc(ipp)
    end do
    irdsp(0) = 0
    do ipp=1,np-1
       irdsp(ipp) = sum(ircvc(0:ipp-1))
    end do

    call MPI_AllGatherV(av,icm*jc(ip),MPI_REAL8,&
         ay_av,ircvc,irdsp,MPI_REAL8,MPI_COMM_WORLD,IERR)

  end function ay_av

  function av_ay(ay)
    real(DP), dimension(:,0:),intent(IN)      :: ay
    real(DP), dimension(size(ay,1),jc(ip))    :: av_ay

    av_ay = ay(:,js(ip):je(ip))

  end function av_ay

  function y_v(v)
    real(DP), dimension(jc(ip)),intent(IN) :: v
    real(DP), dimension(0:jm)              :: y_v

    real(DP), dimension(1,0:jm)    :: ay
    real(DP), dimension(1,jc(ip))  :: av

    av(1,:) = v
    ay = ay_av(av)
    y_v = ay(1,:)
  end function y_v

  !----------------------------------------------------------
  !> Y 座標の格子データから格子データへ変換する.
  !----------------------------------------------------------
  function vx_fy(fy)
    !> 格子点データ
    real(DP), dimension(2*kc(ip),0:jm), intent(in)  :: fy
    !> スペクトルデータ
    real(DP), dimension(jc(ip),0:im-1)              :: vx_fy
    vx_fy = ax_ae(ve_fy(fy))

  end function vx_fy

  !----------------------------------------------------------
  !> Y 座標の格子データから格子データへ変換する.
  !----------------------------------------------------------
  function fy_vx(vx)
    !> 格子点データ
    real(DP), dimension(jc(ip),0:im-1), intent(in)  :: vx
    !> スペクトルデータ
    real(DP), dimension(2*kc(ip),0:jm)              :: fy_vx
    fy_vx = fy_ve(ae_ax(vx))

  end function fy_vx

  !----------------------------------------------------------
  !> Y 座標の格子データからスペクトルデータへ変換する.
  !----------------------------------------------------------
  function ft_fy(fy)
    !> 格子点データ
    real(DP), dimension(2*kc(ip),0:jm), intent(in)  :: fy
    !> スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm)              :: ft_fy
    ft_fy = at_ay(fy)

  end function ft_fy

  !----------------------------------------------------------
  !> Y 座標スペクトルデータから格子データへ変換する.
  !----------------------------------------------------------
  function fy_ft(ft)
    !> スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm), intent(in)  :: ft
    !> 格子点データ
    real(DP), dimension(2*kc(ip),0:jm)              :: fy_ft

    fy_ft = ay_at(ft)

  end function fy_ft

  !----------------------------------------------------------
  !>  入力スペクトルデータに X 微分 \f$\partial_x\f$ を作用する.
  !----------------------------------------------------------
  function ft_Dx_ft(ft)
    !> 入力スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm), intent(in)    :: ft
    real(DP), dimension(2*kc(ip),0:lm)                :: ft_Dx_ft

    integer k, kp, km
    do kp=1,kc(ip)
       km = 2*kc(ip)-kp+1
       k = ks(ip)+kp-1
       ft_Dx_ft(kp,:) = (-2*pi*k/xl)*ft(km,:)
       ft_Dx_ft(km,:) = ( 2*pi*k/xl)*ft(kp,:)
    enddo

  end function ft_Dx_ft

  !----------------------------------------------------------
  !>  入力スペクトルデータに Y 微分 \f$\partial_y\f$ を作用する.
  !----------------------------------------------------------
  function ft_Dy_ft(ft)
    !> 入力スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm), intent(in)   :: ft
    real(DP), dimension(2*kc(ip),0:lm)               :: ft_Dy_ft

    ft_Dy_ft = at_Dy_at(ft)

  end function ft_Dy_ft

  !----------------------------------------------------------
  !> 入力スペクトルデータにラプラシアン
  !> \f$(\partial_{xx}+\partial_{yy})\f$ を作用する.
  !----------------------------------------------------------
  function ft_Lapla_ft(ft)
    !> 入力スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm), intent(in)    :: ft
    real(DP), dimension(2*kc(ip),0:lm)                :: ft_Lapla_ft

    integer k, kp, km
    do kp=1,kc(ip)
       km = 2*kc(ip)-kp+1
       k = ks(ip)+kp-1
       ft_Lapla_ft(kp,:) = -(2*pi*k/xl)**2*ft(kp,:)
       ft_Lapla_ft(km,:) = -(2*pi*k/xl)**2*ft(km,:)
    enddo

    ft_Lapla_ft = ft_Lapla_ft + ft_Dy_ft(ft_Dy_ft(ft))

  end function ft_Lapla_ft

  !----------------------------------------------------------
  !> 2 つのスペクトルデータからヤコビアン
  !> \f$ J(A,B) = (\partial_{x}A)(\partial_{y}B)
  !>            -(\partial_{y}A)(\partial_{x}B)\f$ を計算する.
  !----------------------------------------------------------
  function ft_Jacobian_ft_ft(ft_a,ft_b)
    !> 1つ目の入力スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm), intent(in)    :: ft_a
    !> 2つ目の入力スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm), intent(in)    :: ft_b
    !> 2 つのスペクトルデータのヤコビアン
    real(DP), dimension(2*kc(ip),0:lm)                :: ft_Jacobian_ft_ft
    ft_Jacobian_ft_ft = ft_vx(                           &
      &  vx_ft(ft_Dx_ft(ft_a)) * vx_ft(ft_Dy_ft(ft_b))   &
      & - vx_ft(ft_Dy_ft(ft_a)) * vx_ft(ft_Dx_ft(ft_b)) )

  end function ft_Jacobian_ft_ft


  !----------------------------------------------------------
  !> ディリクレ, ノイマン条件の適用. チェビシェフ空間での計算
  !>
  !> 実際には中で呼ばれている at_module のサブルーチン at_Boundaries_DD,
  !> at_Boundaries_DN, at_Boundaries_ND, at_Boundaries_NN を用いている.
  !> これらを直接呼ぶことも出来る.
  !----------------------------------------------------------
  subroutine ft_Boundaries(ft,values,cond)
    !> 境界での 値/勾配 分布を水平スペクトル変換したものを与える.
    !> 省略時は値/勾配 0 となる.
    real(DP), dimension(2*kc(ip),2), intent(in), optional :: values

    !> 境界条件を適用するデータ. 修正された値を返す.
    real(DP), dimension(2*kc(ip),0:lm),intent(inout)      :: ft

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

    if (.not. present(cond)) then
       if (present(values)) then
          call at_Boundaries_DD(ft,values)
       else
          call at_Boundaries_DD(ft)
       endif
       return
    endif

    select case(cond)
    case ('NN')
       if (present(values)) then
          call at_Boundaries_NN(ft,values)
       else
          call at_Boundaries_NN(ft)
       endif
    case ('DN')
       if (present(values)) then
          call at_Boundaries_DN(ft,values)
       else
          call at_Boundaries_DN(ft)
       endif
    case ('ND')
       if (present(values)) then
          call at_Boundaries_ND(ft,values)
       else
          call at_Boundaries_ND(ft)
       endif
    case ('DD')
       if (present(values)) then
          call at_Boundaries_DD(ft,values)
       else
          call at_Boundaries_DD(ft)
       endif
    case default
       call MessageNotify('E','ft_Boundaries','B.C. not supported')
    end select

  end subroutine ft_Boundaries

  !----------------------------------------------------------
  !> 境界で一様な値を与える条件(ディリクレ条件)下で
  !> 入力スペクトルデータに逆ラプラシアン
  !> \f$(\partial_{xx}+\partial_{yy})^{-1}\f$ を作用する
  !> (Chebyshev-tau法)
  !----------------------------------------------------------
  function ft_LaplaInv_ft(ft,values)
    !> スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm),intent(in)  :: ft
    !> 境界値. 省略時は 0 が設定される.
    real(DP), dimension(2*kc(ip),2), intent(in), optional :: values
    !> スペクトルデータの逆ラプラシアン
    real(DP), dimension(2*kc(ip),0:lm)             :: ft_LaplaInv_ft

    real(DP), dimension(:,:,:), allocatable  :: alu
    integer,  dimension(:,:), allocatable    :: kp
    real(DP), dimension(2*kc(ip),0:lm)       :: ft_work
    real(DP), dimension(2*kc(ip),0:jm)       :: fy_work
    real(DP), dimension(2*kc(ip))            :: value1, value2   ! 境界値

    logical :: first = .true.
    integer :: l
    save    :: alu, kp, first

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

    if ( first ) then
      first = .false.

      allocate(alu(2*kc(ip),0:lm,0:lm),kp(2*kc(ip),0:lm))

      do l=0,lm
        ft_work=0.0D0
        ft_work(:,l) = 1.0D0
        alu(:,:,l) = ft_Lapla_ft(ft_work)
        fy_work = ay_at(ft_work)
        alu(:,lm-1,l) = fy_work(:,0)
        alu(:,lm,l)   = fy_work(:,jm)
      enddo

      call ludecomp(alu,kp)
    endif

    ft_work = ft
    ft_work(:,lm-1) = value1
    ft_work(:,lm)   = value2
    ft_LaplaInv_ft = lusolve(alu,kp,ft_work)

  end function ft_LaplaInv_ft

  !----------------------------------------------------------
  !> @brief
  !> Y 方向格子点空間での Chebyshev-Collocation 法により，
  !> 渦度から流線を求める.
  !>
  !> @details
  !> 渦度 \f$\zeta\f$ を与えて流線 \f$\psi\f$ を求める.
  !> \f{align*}{
  !>    \nabla^2 \psi & = \zeta, \\
  !>    \psi & = \mbox{const. at boundaries.}
  !> \f}
  !> 粘着条件
  !> \f[ \partial_{y} \psi  = 0 \, \mbox{at boundaries} \f]
  !> 応力なし条件
  !> \f[ (\partial_{y})^{2} \psi = 0\, \mbox{at boundaries} \f]
  !----------------------------------------------------------
  function fy_Vor2Strm_fy(fy,values,cond,new)
    !> 入力渦度分布
    real(DP), dimension(2*kc(ip),0:jm),intent(in)  :: fy
    !> 流線境界値. 境界で一定なので波数 0 成分のみ. 省略時は 0.
    real(DP), dimension(2), intent(in), optional :: values
    !> 境界条件スイッチ. 省略時は 'RR'
    !> * RR : 両端粘着条件
    !> * RF : 上端粘着下端応力なし条件
    !> * FR : 上端応力なし下端粘着条件
    !> * FF : 両端応力なし条件
    character(len=2), intent(in), optional  :: cond
    !> true だと境界条件計算用行列を強制的に新たに作る. default は false.
    logical, intent(IN), optional :: new
    !> 出力流線分布
    real(DP), dimension(2*kc(ip),0:jm)             :: fy_Vor2Strm_fy

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

    real(DP), dimension(2*kc(ip),0:jm)         :: fy_work
    real(DP), dimension(2*kc(ip),0:jm)         :: fy_I
    real(DP)                                 :: value1, value2   ! 境界値
    logical                                  :: rigid1 = .true.
    logical                                  :: rigid2 = .true.
    logical :: first = .true.
    logical :: new_matrix = .false.
    integer :: j
    save    :: alu, kp, first

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

    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','fy_Vor2Strm_fy','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*kc(ip),0:jm,0:jm),kp(2*kc(ip),0:jm))

      do j=0,jm
        fy_I = 0.0D0
        fy_I(:,j) = 1.0D0
        alu(:,:,j) = fy_ft(ft_Lapla_ft(ft_fy(fy_I)))
      enddo

      do j=0,jm
        fy_I = 0.0D0
        fy_I(:,j) = 1.0D0

        ! 運動学的条件. 流線は境界で一定
        alu(:,0,j)   = fy_I(:,0)
        alu(:,jm,j)  = fy_I(:,jm)

        if ( rigid1 ) then                               ! 粘着条件(top)
          fy_work=fy_ft(ft_Dy_ft(ft_fy(fy_I)))
        else
          fy_work=fy_ft(ft_Dy_ft(ft_Dy_ft(ft_fy(fy_I)))) ! すべり条件(top)
        endif

        alu(:,1,j) = fy_work(:,0)

        if ( rigid2 ) then                               ! 粘着条件(bottom)
          fy_work=fy_ft(ft_Dy_ft(ft_fy(fy_I)))
        else
          fy_work=fy_ft(ft_Dy_ft(ft_Dy_ft(ft_fy(fy_I)))) ! すべり条件(bottom)
        endif
        alu(:,jm-1,j) = fy_work(:,jm)
      enddo

      call ludecomp(alu,kp)

      call MessageNotify('M','fy_Vor2Strm_fy',&
        'Matrix for stream func. calc. produced')
    endif

    fy_work = fy
    fy_work(:,1)    = 0.0D0               ! 力学的条件
    fy_work(:,jm-1) = 0.0D0               ! 力学的条件

    fy_work(:,0) = 0.0D0                  ! 運動学的条件. 波数 0 以外は 0
    if ( kf_k(0) > 0 ) then
       fy_work(kf_k(0),0) = value1               ! 運動学的条件.
       fy_work(2*kc(ip)-kf_k(0)+1,0) = value1    ! 運動学的条件.
    endif

    fy_work(:,jm)   = 0.0D0               ! 運動学的条件. 波数 0 以外は 0
    if ( kf_k(0) > 0 ) then
       fy_work(kf_k(0),jm)   = value2                 ! 運動学的条件.
       fy_work(2*kc(ip)-kf_k(0)+1,jm)   = value2      ! 運動学的条件.
    endif

    fy_Vor2Strm_fy = lusolve(alu,kp,fy_work)

  end function fy_Vor2Strm_fy

  !----------------------------------------------------------
  !> @brief Chebyshev-tau 法により，渦度から流線を求める(境界値は定数)
  !> @details
  !> 渦度 \f$\zeta\f$ を与えて流線 \f$\psi\f$ を求める.
  !> \f{align*}{
  !>    \nabla^2\psi & = \zeta, \\
  !>    \psi & = \mbox{const. at boundaries.}
  !> \f}
  !> 粘着条件
  !> \f[   \partial_{y} \psi = 0 \, \mbox{at boundaries} \f]
  !> 応力なし条件
  !> \f[  (\partial_{y})^{2} \psi = 0 \, \mbox{at boundaries} \f]
  !----------------------------------------------------------
  function ft_Vor2Strm_ft(ft,values,cond,new)
    !> 入力渦度分布
    real(DP), dimension(2*kc(ip),0:lm),intent(in)  :: ft
    !> 流線境界値. 境界で一定なので波数 0 成分のみ. 省略時は 0.
    real(DP), dimension(2), intent(in), optional :: values
    !> 境界条件スイッチ. 省略時は 'RR'
    !> * RR : 両端粘着条件
    !> * RF : 上端粘着下端応力なし条件
    !> * FR : 上端応力なし下端粘着条件
    !> * FF : 両端応力なし条件
    character(len=2), intent(in), optional  :: cond
    !> true だと境界条件計算用行列を強制的に新たに作る.
    !> default は false.
    logical, intent(IN), optional :: new
    !> 出力流線分布
    real(DP), dimension(2*kc(ip),0:lm)       :: ft_Vor2Strm_ft

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

    real(DP), dimension(2*kc(ip),0:lm)       :: ft_work
    real(DP), dimension(2*kc(ip),0:jm)       :: fy_work
    real(DP)                                 :: value1
    real(DP)                                 :: value2
    logical                                  :: rigid1 = .true.
    logical                                  :: rigid2 = .true.


    logical :: first = .true.
    logical :: new_matrix = .false.
    integer :: l
    save    :: alu, kp, first

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

    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','ft_Vor2Strm_ft','B.C. not supported')
      end select
    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*kc(ip),0:lm,0:lm),kp(2*kc(ip),0:lm))

      alu = 0.0D0
      do l=0,lm-2
        ft_work(:,:)= 0.0D0
        ft_work(:,l)= 1.0D0

        alu(:,:,l) = ft_Lapla_ft(ft_work)
      enddo

      do l=0,lm
        ft_work(:,:)= 0.0D0
        ft_work(:,l)= 1.0D0

        ! 運動学的条件. 流線は境界で一定
        fy_work = ay_at(ft_work)
        alu(:,lm,l) = fy_work(:,0)
        alu(:,lm-1,l) = fy_work(:,jm)

        if ( rigid1 ) then                   ! 力学的条件粘着条件(Top)
          fy_work=ay_at(ft_Dy_ft(ft_work))
          alu(:,lm-2,l) = fy_work(:,0)
        else                                 ! 力学的条件すべり条件(Top)
          fy_work=ay_at(ft_Dy_ft(ft_Dy_ft(ft_work)))
          alu(:,lm-2,l) = fy_work(:,0)
        endif
        if ( rigid2 ) then                   ! 力学的条件粘着条件(bottom)
          fy_work=ay_at(ft_Dy_ft(ft_work))
          alu(:,lm-3,l) = fy_work(:,jm)
        else                                 ! 力学的条件すべり条件(bottom)
          fy_work=ay_at(ft_Dy_ft(ft_Dy_ft(ft_work)))
          alu(:,lm-3,l) = fy_work(:,jm)
        endif
      enddo

      call ludecomp(alu,kp)

      call MessageNotify('M','ft_Vor2Strm_ft',&
        'Matrix for stream func. calc. produced')
    endif

    ! 内部領域計算
    ft_work = ft

    ft_work(:,lm-2) = 0.0D0           ! 力学的条件
    ft_work(:,lm-3) = 0.0D0           ! 力学的条件

    ft_work(:,lm-1) = 0.0D0        ! 運動学的条件. 波数 0 以外は 0
    !!$      ft_work(0,lm-1) = value2*2     ! 運動学的条件. 波数 0 は重み 1/2
    if ( kf_k(0) > 0 ) then
       ft_work(kf_k(0),lm-1) = value2       ! 運動学的条件. 波数 0 は重み 1/2
       ft_work(2*kc(ip)-kf_k(0)+1,lm-1) = value2       ! 運動学的条件. 波数 0 は重み 1/2
    endif

    ft_work(:,lm)   = 0.0D0        ! 運動学的条件. 波数 0 以外は 0
!!$      ft_work(0,lm)   = value1*2     ! 運動学的条件. 波数 0 は重み 1/2
    if ( kf_k(0) > 0 ) then
       ft_work(kf_k(0),lm)   = value1     ! 運動学的条件. 波数 0 は重み 1/2
       ft_work(2*kc(ip)-kf_k(0)+1,lm)   = value1     ! 運動学的条件. 波数 0 は重み 1/2
    endif
    ft_Vor2Strm_ft = lusolve(alu,kp,ft_work)

  end function ft_Vor2Strm_ft

  !----------------------------------------------------------
  !> @brief Chebyshev-tau 法により，渦度から流線を求める．
  !>        境界値として X方向の スペクトルを与えることが可能
  !> @details
  !> 渦度 \f$\zeta\f$ を与えて流線 \f$\psi\f$ を求める.
  !> \f{align*}{
  !>    \nabla^2\psi & = \zeta, \\
  !>    \psi & = \mbox{const. at boundaries.}
  !> \f}
  !> 粘着条件
  !> \f[   \partial_{y} \psi = 0 \, \mbox{at boundaries} \f]
  !> 応力なし条件
  !> \f[  (\partial_{y})^{2} \psi = 0 \, \mbox{at boundaries} \f]
  !----------------------------------------------------------
  function ft_Vor2Strm2_ft(ft,fa,cond,new)
    !> 入力渦度分布
    real(DP), dimension(2*kc(ip),0:lm),intent(in)  :: ft
    !> 流線境界値. 省略時は 0. 1:Top,2:Btm
    real(DP), dimension(2*kc(ip),2), intent(in), optional :: fa
    !> 境界条件スイッチ. 省略時は 'RR'
    !> * RR : 両端粘着条件
    !> * RF : 上端粘着下端応力なし条件
    !> * FR : 上端応力なし下端粘着条件
    !> * FF : 両端応力なし条件
    character(len=2), intent(in), optional  :: cond
    !> true だと境界条件計算用行列を強制的に新たに作る.
    !> default は false.
    logical, intent(IN), optional :: new
    !> 出力流線分布
    real(DP), dimension(2*kc(ip),0:lm)         :: ft_Vor2Strm2_ft

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

    real(DP), dimension(2*kc(ip),0:lm)         :: ft_work
    real(DP), dimension(2*kc(ip),0:jm)         :: fy_work
    real(DP), dimension(2*kc(ip))              :: f_Top
    real(DP), dimension(2*kc(ip))              :: f_Btm
    logical                                  :: rigid1 = .true.
    logical                                  :: rigid2 = .true.

    logical :: first = .true.
    logical :: new_matrix = .false.
    integer :: l
    save    :: alu, kp, first

    if (.not. present(fa)) then
      f_Top = 0.0D0
      f_Btm = 0.0D0
    else
      f_Top = fa(:,1)
      f_Btm = fa(:,2)
    endif

    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','ft_Vor2Strm_ft','B.C. not supported')
      end select
    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*kc(ip),0:lm,0:lm),kp(2*kc(ip),0:lm))

      alu = 0.0D0
      do l=0,lm-2
        ft_work(:,:)= 0.0D0
        ft_work(:,l)= 1.0D0

        alu(:,:,l) = ft_Lapla_ft(ft_work)
      enddo

      do l=0,lm
        ft_work(:,:)= 0.0D0
        ft_work(:,l)= 1.0D0

        ! 運動学的条件. 流線は境界で一定
        fy_work = ay_at(ft_work)
        alu(:,lm,l) = fy_work(:,0)
        alu(:,lm-1,l) = fy_work(:,jm)

        if ( rigid1 ) then                   ! 力学的条件粘着条件(Top)
          fy_work=ay_at(ft_Dy_ft(ft_work))
          alu(:,lm-2,l) = fy_work(:,0)
        else                                 ! 力学的条件すべり条件(Top)
          fy_work=ay_at(ft_Dy_ft(ft_Dy_ft(ft_work)))
          alu(:,lm-2,l) = fy_work(:,0)
        endif
        if ( rigid2 ) then                   ! 力学的条件粘着条件(bottom)
          fy_work=ay_at(ft_Dy_ft(ft_work))
          alu(:,lm-3,l) = fy_work(:,jm)
        else                                 ! 力学的条件すべり条件(bottom)
          fy_work=ay_at(ft_Dy_ft(ft_Dy_ft(ft_work)))
          alu(:,lm-3,l) = fy_work(:,jm)
        endif
      enddo

      call ludecomp(alu,kp)

      call MessageNotify('M','ft_Vor2Strm_ft',&
        'Matrix for stream func. calc. produced')
    endif

    ! 内部領域計算
    ft_work = ft

    ft_work(:,lm-2) = 0.0D0           ! 力学的条件
    ft_work(:,lm-3) = 0.0D0           ! 力学的条件

    ft_work(:,lm-1) = f_Btm           ! 運動学的条件. 波数 0 以外は 0
    ft_work(:,lm)   = f_Top           ! 運動学的条件. 波数 0 以外は 0

    ft_Vor2Strm2_ft = lusolve(alu,kp,ft_work)

  end function ft_Vor2Strm2_ft

  !----------------------------------------------------------
  !> @brief Chebyshev-tau 法により，渦度から流線を求める．
  !>        境界での値と勾配のスペクトルを与えることが可能.
  !> @details
  !> 渦度 \f$\zeta\f$ を与えて流線 \f$\psi\f$ を求める.
  !> \f{align*}{
  !>    \nabla^2\psi & = \zeta, \\
  !>    \psi & = \mbox{const. at boundaries.}
  !> \f}
  !> 粘着条件
  !> \f[   \partial_{y} \psi = 0 \, \mbox{at boundaries} \f]
  !> 応力なし条件
  !> \f[  (\partial_{y})^{2} \psi = 0 \, \mbox{at boundaries} \f]
  !----------------------------------------------------------
  function ft_Vor2Strm3_ft(ft,fa,cond,new)
    !> 入力渦度分布
    real(DP), dimension(2*kc(ip),0:lm),intent(in)  :: ft
    !> 流線/微分境界値. 省略時は 0.
    !> 1. Top 値,
    !> 2. Bottom 値
    !> 3. Top 微分(勾配)
    !> 4. Bottom 微分(勾配)
    real(DP), dimension(2*kc(ip),4), intent(in), optional :: fa
    !> 境界条件スイッチ. 省略時は 'RR'
    !> * RR : 両端粘着条件
    !> * RF : 上端粘着下端応力なし条件
    !> * FR : 上端応力なし下端粘着条件
    !> * FF : 両端応力なし条件
    character(len=2), intent(in), optional  :: cond
    !> true だと境界条件計算用行列を強制的に新たに作る.
    !> default は false.
    logical, intent(IN), optional :: new
    !> 出力流線分布
    real(DP), dimension(2*kc(ip),0:lm)             :: ft_Vor2Strm3_ft

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

    real(DP), dimension(2*kc(ip),0:lm)        :: ft_work
    real(DP), dimension(2*kc(ip),0:jm)        :: fy_work
    real(DP), dimension(2*kc(ip))             :: f_TopValue
    real(DP), dimension(2*kc(ip))             :: f_BtmValue
    real(DP), dimension(2*kc(ip))             :: f_TopDiff
    real(DP), dimension(2*kc(ip))             :: f_BtmDiff
    logical                                 :: rigid1 = .true.
    logical                                 :: rigid2 = .true.
    logical :: first = .true.
    logical :: new_matrix = .false.
    integer :: l
    save    :: alu, kp, first

    if (.not. present(fa)) then
      f_TopValue = 0.0D0
      f_BtmValue = 0.0D0
      f_TopDiff  = 0.0D0
      f_BtmDiff  = 0.0D0
    else
      f_TopValue = fa(:,1)
      f_BtmValue = fa(:,2)
      f_TopDiff  = fa(:,3)
      f_BtmDiff  = fa(:,4)
    endif

    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','ft_Vor2Strm_ft','B.C. not supported')
      end select
    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*kc(ip),0:lm,0:lm),kp(2*kc(ip),0:lm))

      alu = 0.0D0
      do l=0,lm-2
        ft_work(:,:)= 0.0D0
        ft_work(:,l)= 1.0D0

        alu(:,:,l) = ft_Lapla_ft(ft_work)
      enddo

      do l=0,lm
        ft_work(:,:)= 0.0D0
        ft_work(:,l)= 1.0D0

        ! 運動学的条件. 流線は境界で一定
        fy_work = ay_at(ft_work)
        alu(:,lm,l) = fy_work(:,0)
        alu(:,lm-1,l) = fy_work(:,jm)

        if ( rigid1 ) then                   ! 力学的条件粘着条件(Top)
          fy_work=ay_at(ft_Dy_ft(ft_work))
          alu(:,lm-2,l) = fy_work(:,0)
        else                                 ! 力学的条件すべり条件(Top)
          fy_work=ay_at(ft_Dy_ft(ft_Dy_ft(ft_work)))
          alu(:,lm-2,l) = fy_work(:,0)
        endif
        if ( rigid2 ) then                   ! 力学的条件粘着条件(bottom)
          fy_work=ay_at(ft_Dy_ft(ft_work))
          alu(:,lm-3,l) = fy_work(:,jm)
        else                                 ! 力学的条件すべり条件(bottom)
          fy_work=ay_at(ft_Dy_ft(ft_Dy_ft(ft_work)))
          alu(:,lm-3,l) = fy_work(:,jm)
        endif
      enddo

      call ludecomp(alu,kp)

      call MessageNotify('M','ft_Vor2Strm_ft',&
        'Matrix for stream func. calc. produced')
    endif

    ! 内部領域計算
    ft_work = ft

    ft_work(:,lm-2) = f_TopDIff         ! 力学的条件
    ft_work(:,lm-3) = f_BtmDiff         ! 力学的条件

    ft_work(:,lm-1) = f_BtmValue        ! 運動学的条件. 波数 0 以外は 0
    ft_work(:,lm)   = f_TopValue        ! 運動学的条件. 波数 0 以外は 0

    ft_Vor2Strm3_ft = lusolve(alu,kp,ft_work)

  end function ft_Vor2Strm3_ft

  !----------------------------------------------------------
  !> 2 次元格子点データの全領域積分および平均.
  !----------------------------------------------------------
  function IntYX_vx(vx)
    !> 2 次元格子点データ
    real(DP), dimension(jc(ip),0:im-1)   :: vx
    !> 積分値
    real(DP)                           :: IntYX_vx

    real(DP) :: IntYXTmp
    integer :: i, j
    integer :: ierr

    IntYXTmp = 0.0d0
    do i=0,im-1
       do j=1,jc(ip)
          IntYXTmp = IntYXTmp + vx(j,i) * v_Y_Weight(j) * x_X_Weight(i)
       enddo
    enddo

    CALL MPI_ALLREDUCE(IntYXTmp,IntYX_vx,1,MPI_REAL8, &
                       MPI_SUM,MPI_COMM_WORLD,IERR)

  end function IntYX_vx

  function v_IntX_vx(vx)
    !
    ! 2 次元格子点データの X 方向積分
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
    !
    real(DP), dimension(jc(ip),0:im-1)   :: vx
    !(in) 2 次元格子点データ

    real(DP), dimension(jc(ip))          :: v_IntX_vx
    !(out) 積分された 1 次元(Y)格子点データ

    integer :: i
    ! 作業変数

    v_IntX_vx = 0.0d0
    do i=0,im-1
       v_IntX_vx(:) = v_IntX_vx(:) + vx(:,i) * x_X_Weight(i)
    enddo
  end function v_IntX_vx

  function x_IntY_vx(vx)
    !
    ! 2 次元格子点データの Y 方向積分
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している.
    !
    real(DP), dimension(jc(ip),0:im-1)   :: vx
    !(in)  2 次元格子点データ

    real(DP), dimension(0:im-1)        :: x_IntY_vx
    !(out) 積分された 1 次元(X)格子点データ

    real(DP), dimension(0:im-1)        :: x_IntYTmp
    integer :: j, ierr
    ! 作業変数

    x_IntYTmp = 0.0d0
    do j=1,jc(ip)
       x_IntYTmp(:) = x_IntYTmp(:) + vx(j,:) * v_Y_Weight(j)
    enddo

    CALL MPI_ALLREDUCE(x_IntYTmp,x_IntY_vx,im,MPI_REAL8, &
                       MPI_SUM,MPI_COMM_WORLD,IERR)

  end function x_IntY_vx

  function IntX_x(x)
    !
    ! 1 次元(X)格子点データの X 方向積分
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
    !
    real(DP), dimension(0:im-1)   :: x          !(in)  1 次元格子点データ
    real(DP)                      :: IntX_x     !(out) 積分値

    IntX_x = sum(x*x_X_Weight)
  end function IntX_x

  function IntY_v(v)      ! Y 方向積分
    !
    ! 1 次元(Y)格子点データの Y 方向積分
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している.
    !
    real(DP), dimension(jc(ip))   :: v          !(in)  1 次元格子点データ
    real(DP)                      :: IntY_v     !(out) 積分値

    real(DP)                      :: IntYTmp     !(out) 積分値
    integer :: ierr

    IntYTmp = sum(v*v_Y_Weight)

    CALL MPI_ALLREDUCE(IntYTmp,IntY_v,1,MPI_REAL8, &
                       MPI_SUM,MPI_COMM_WORLD,IERR)

  end function IntY_v

  !--------------- 平均計算 -----------------
  function AvrYX_vx(vx)
    !
    ! 2 次元格子点データの全領域平均
    !
    ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
    ! 総和を計算し, x_X_Weight*y_Y_Weight の総和で割ることで平均している.
    !
    real(DP), dimension(jc(ip),0:im-1)   :: vx
    !(in)  2 次元格子点データ

    real(DP)                             :: AvrYX_vx
    !(out) 平均値

    AvrYX_vx = IntYX_vx(vx)/(xl*yl)

  end function AvrYX_vx

  function v_AvrX_vx(vx)
    !
    ! 2 次元格子点データの X 方向平均
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し,
    ! x_X_Weight の総和で割ることで平均している.
    !
    real(DP), dimension(jc(ip),0:im-1)   :: vx
    !(in) 2 次元格子点データ

    real(DP), dimension(jc(ip))          :: v_AvrX_vx
    !(out) 平均された 1 次元(Y)格子点データ

    v_AvrX_vx = v_IntX_vx(vx)/xl

  end function v_AvrX_vx

  function x_AvrY_vx(vx)
    !
    ! 2 次元格子点データの Y 方向平均
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算し,
    ! y_Y_Weight の総和で割ることで平均している.
    !
    real(DP), dimension(jc(ip),0:im-1)   :: vx
    !(in) 2 次元格子点データ

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

    x_AvrY_vx = x_IntY_vx(vx)/yl

  end function x_AvrY_vx

  function AvrX_x(x)
    !
    ! 1 次元(X)格子点データの X 方向平均
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し,
    ! x_X_Weight の総和で割ることで平均している.
    !
    real(DP), dimension(0:im-1)   :: x          !(in)  1 次元格子点データ
    real(DP)                      :: AvrX_x     !(out) 平均値

    AvrX_x = IntX_x(x)/xl

  end function AvrX_x

  function AvrY_v(v)
    !
    ! 1 次元(Y)格子点データの Y 方向平均
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算し,
    ! y_Y_Weight の総和で割ることで平均している.
    !
    real(DP), dimension(jc(ip)) :: v          !(in)  1 次元格子点データ
    real(DP)                    :: AvrY_v     !(out) 平均値

    AvrY_v = IntY_v(v)/yl

  end function AvrY_v


  !----------------------------------------------------------
  !> 補間関数
  !----------------------------------------------------------
  function Interpolate_ft(ft_data,xval,yval)
    !> 入力スペクトルデータ
    real(DP), dimension(2*kc(ip),0:lm), intent(IN) :: ft_data
    !> 補間する点の X 座標
    real(DP), intent(in)                         :: xval
    !> 補間する点の Y 座標
    real(DP), intent(in)                         :: yval
    !> 補間結果
    real(DP)                                     :: Interpolate_ft
    real(DP)                                     :: af_work(1,2*kc(ip))
    real(DP)                                     :: a_work(1)

    af_work(1,:) = a_Interpolate_at(ft_data,yval)
    a_work = a_Interpolate_ae(ae_af(af_work),xval)
    Interpolate_ft = a_work(1)

  end function Interpolate_ft

end module et_mpi_module