! -*- mode: f90; coding: utf-8 -*-
!----------------------------------------------------------------------
! Copyright(c) 2002-2019 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_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 次元格子点データに同じ.
!>
!> * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
!>   微分などを作用させたデータをスペクトル変換したものことである.
!----------------------------------------------------------------------
!履歴  2002/01/27  竹広真一
!      2002/03/30  竹広真一  モジュール名変更
!      2002/08/19  竹広真一  格子データ添字を gg -> xy, eg -> ey に変更
!                            積分・平均関数追加
!      2005/01/09  竹広真一  msgdmp -> MessageNotify に変更
!      2005/03/15  竹広真一  xy -> yx に接頭子を変更
!      2006/03/04  竹広真一  余計な変数を削除(et_Jacobian_et_et)
!      2006/03/06  竹広真一  コメントを RDoc 用に修正
!      2007/11/20  竹広真一  渦度-流線関数計算に行列設定スイッチ導入
!      2007/11/21  竹広真一  初期化サブルーチンメッセージ出力
!      2009/01/09  竹広真一  et_Initial メッセージに日付を追加
!      2009/01/29  佐々木洋平 コメントを RDoc 用に修正
!      2009/07/31  竹広真一  境界条件計算用配列を threadprivate 指定(OpenMP)
!      2010/03/10  佐々木洋平  threadprivate 削除(コンパイラ依存)
!      2011/03/06  佐々木洋平 境界条件スイッチを修正, RR をデフォルトに.
!      2011/12/03  竹広真一  補間関数を追加
!      2012/07/08  竹広真一  et_Vor2Strm_et を追加
!      2013/08/20  竹広真一  gnu fortran 対応
!      2015/12/11  竹広真一  et_Vort2Strm_et の境界値オプションを追加
!      2016/09/26  竹広真一  et_BoundariesGrid 追加
!----------------------------------------------------------------------
module et_module
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! et_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  !
  !==== 座標変数
  !
  ! x_X, y_Y     ::  格子点座標(X,Y座標)を格納した 1 次元配列
  ! x_X_Weight, y_Y_Weight ::  重み座標を格納した 1 次元配列
  ! yx_X, yx_Y   :: 格子点データの XY 座標(X,Y)(格子点データ型 2 次元配列)
  !
  !==== 基本変換
  !
  ! yx_et :: スペクトルデータから格子データへの変換
  ! et_yx :: 格子データからスペクトルデータへの変換
  ! ax_ae, x_e :: X 方向のスペクトルデータから格子データへの変換
  ! ay_at, y_t :: Y 方向のスペクトルデータから格子データへの変換
  ! ae_ax, e_x :: X 方向の格子点データからスペクトルデータへの変換
  ! at_ay, t_y :: Y 方向の格子点データからスペクトルデータへの変換
  !
  !==== 微分
  !
  ! et_Lapla_et  :: スペクトルデータにラプラシアンを作用させる
  ! et_Dx_et, ae_Dx_ae, e_Dx_e :: スペクトルデータに X 微分を作用させる
  ! et_Dy_et, at_Dy_at, t_Dy_t :: スペクトルデータに Y 微分を作用させる
  ! et_Jacobian_et_et :: 2 つのスペクトルデータからヤコビアンを計算する
  !
  !==== 境界値問題
  !
  ! et_Boundaries  :: ディリクレ, ノイマン境界条件の適用
  ! et_LaplaInv_et :: スペクトルデータにラプラシアンの逆変換を作用させる
  ! et_Vor2Strm_et :: 渦度から流線を計算する
  ! ey_Vor2Strm_ey :: 渦度から流線を計算する
  !
  !==== 積分・平均
  !
  ! IntYX_yx, AvrYX_yx   :: 2 次元格子点データの全領域積分および平均
  ! y_IntX_yx, y_AvrX_yx :: 2 次元格子点データの X 方向積分および平均
  ! IntX_x, AvrX_x       :: 1 次元(X)格子点データの X 方向積分および平均
  ! x_IntY_yx, x_AvrY_yx :: 2 次元格子点データの Y 方向積分および平均
  ! IntY_y, AvrY_y       :: 1 次元(Y)格子点データの Y 方向積分および平均
  !
  !==== 補間
  !
  ! Interpolate_et       :: スペクトルデータからの特定の座標の値の補間
  !
  use dc_types, only: DP
  use dc_message, only: MessageNotify
  use lumatrix
  use spml_utils, only: pi
  use ae_module, x_X => g_X, x_X_weight => g_X_Weight, &
    e_x => e_g, ae_ax => ae_ag, &
    x_e => g_e, ax_ae => ag_ae
  use at_module, y_Y => g_X, y_Y_Weight => g_X_Weight, &
    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, &
    et_BoundariesGrid => at_BoundariesGrid
  implicit none
  private

  public et_Initial                                       ! 初期化
  public x_X, y_Y, x_X_Weight, y_Y_Weight, yx_X, yx_Y     ! 座標変数

  public yx_et, et_yx                                     ! 基本変換
  public yx_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 et_Dx_et, e_Dx_e, ae_Dx_ae                       ! 微分
  public et_Dy_et, t_Dy_t, at_Dy_at                       ! 微分
  public et_Lapla_et                                      ! 微分

  public et_Jacobian_et_et                                ! 非線形計算

  public et_Boundaries                                    ! 境界値問題
  public at_Boundaries_DD, at_Boundaries_DN               ! 境界値問題
  public at_Boundaries_ND, at_Boundaries_NN               ! 境界値問題
  public et_BoundariesGrid                                ! 境界値問題
  public et_LaplaInv_et, ey_Vor2Strm_ey                   ! 境界値問題
  public et_Vor2Strm_et                                   ! 境界値問題
  public et_Vor2Strm2_et                                  ! 境界値問題
  public et_Vor2Strm3_et                                  ! 境界値問題

  public IntYX_yx, y_IntX_yx, x_IntY_yx, IntX_x, IntY_y   ! 積分
  public AvrYX_yx, y_AvrX_yx, x_AvrY_yx, AvrX_x, AvrY_y   ! 平均

  public Interpolate_et                                   ! 補間

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

  !> X方向格子点
  integer            :: im=32
  !> Y方向格子点
  integer            :: jm=8
  !> X方向切断波数
  integer            :: km=10
  !> Y方向切断波数
  integer            :: lm=5
  !> X 方向計算領域
  real(DP)           :: xl=2.0D0
  !> Y 方向計算領域
  real(DP)           :: yl=1.0D0

  save im, jm, km, lm, xl, yl

contains
  !----------------------------------------------------------
  !> スペクトル変換の格子点数, 波数, 領域の大きさを設定する.
  !>
  !> 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
  !> 初期設定をしなければならない.
  !----------------------------------------------------------
  subroutine et_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

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

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

    allocate(yx_X(0:jm,0:im-1),yx_Y(0:jm,0:im-1))
    yx_X = spread(x_X,1,jm+1)
    yx_Y = spread(y_Y,2,im)

    call MessageNotify('M','et_initial', &
      & 'et_module (2019/09/30) is initialized')
  end subroutine et_initial


  !----------------------------------------------------------
  !> スペクトルデータから格子データへ変換する.
  !----------------------------------------------------------
  function yx_et(et)
    !> スペクトルデータ
    real(DP), dimension(-km:km,0:lm), intent(in)  :: et
    !> 格子点データ
    real(DP), dimension(0:jm,0:im-1)              :: yx_et

    real(DP), dimension(-km:km,0:jm)              :: ey
    real(DP), dimension(0:jm,-km:km)              :: ye

    integer :: j, k

    ! yx_et = ax_ae(transpose(ay_at(et)))
    !
    ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように
    ! 書き直す
    !
    ey = ay_at(et)
    do j=0,jm
      do k=-km,km
        ye(j,k) = ey(k,j)
      enddo
    enddo
    yx_et = ax_ae(ye)

  end function yx_et

  !----------------------------------------------------------
  ! 格子データからスペクトルデータへ変換する.
  !----------------------------------------------------------
  function et_yx(yx)
    !> 格子点データ
    real(DP), dimension(0:jm,0:im-1), intent(in)  :: yx
    !> スペクトルデータ
    real(DP), dimension(-km:km,0:lm)              :: et_yx


    real(DP), dimension(-km:km,0:jm)              :: ey
    real(DP), dimension(0:jm,-km:km)              :: ye

    integer :: j, k

    ! et_yx = at_ay(transpose(ae_ax(yx)))
    !
    ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
    !
    ye = ae_ax(yx)
    do j=0,jm
      do k=-km,km
        ey(k,j) = ye(j,k)
      enddo
    enddo
    et_yx = at_ay(ey)

  end function et_yx

  !----------------------------------------------------------
  !> X 座標スペクトルデータから格子データへ変換する.
  !----------------------------------------------------------
  function yx_ey(ey)
    !> スペクトルデータ
    real(DP), dimension(-km:km,0:jm), intent(in)  :: ey

    !> 格子点データ
    real(DP), dimension(0:jm,0:im-1)              :: yx_ey

    real(DP), dimension(0:jm,-km:km)              :: ye
    integer :: j, k

    ! yx_ey = ax_ae(transpose(ey))
    !
    ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書
    ! き直す
    !
    do j=0,jm
      do k=-km,km
        ye(j,k) = ey(k,j)
      enddo
    enddo
    yx_ey = ax_ae(ye)

  end function yx_ey

  !----------------------------------------------------------
  ! X 座標格子データからスペクトルデータへ変換する.
  !----------------------------------------------------------
  function ey_yx(yx)
    !> 格子点データ
    real(DP), dimension(0:jm,0:im-1), intent(in)  :: yx
    !> スペクトルデータ
    real(DP), dimension(-km:km,0:jm)              :: ey_yx

    real(DP), dimension(0:jm,-km:km)              :: ye
    integer :: j, k

    ! ey_yx = transpose(ae_ax(yx))
    !
    ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書
    ! き直す
    !
    ye = ae_ax(yx)
    do j=0,jm
      do k=-km,km
        ey_yx(k,j) = ye(j,k)
      enddo
    enddo

  end function ey_yx

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

  end function et_ey

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

    ey_et = ay_at(et)

  end function ey_et


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

    integer k
    do k=-km,km
      et_Dx_et(k,:)  =  (-2*pi*k/xl)*et(-k,:)
    enddo
  end function et_Dx_et

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

    et_Dy_et = at_Dy_at(et)

  end function et_Dy_et

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

    integer k
    do k=-km,km
      et_Lapla_et(k,:) = -(2*pi*k/xl)**2*et(k,:)
    enddo

    et_Lapla_et = et_Lapla_et + et_Dy_et(et_Dy_et(et))

  end function et_Lapla_et

  !----------------------------------------------------------
  !> 2 つのスペクトルデータからヤコビアン
  !> \f$ J(A,B) = (\partial_{x}A)(\partial_{y}B)
  !>            -(\partial_{y}A)(\partial_{x}B)\f$ を計算する.
  !----------------------------------------------------------
  function et_Jacobian_et_et(et_a,et_b)
    !> 1つ目の入力スペクトルデータ
    real(DP), dimension(-km:km,0:lm), intent(in)    :: et_a
    !> 2つ目の入力スペクトルデータ
    real(DP), dimension(-km:km,0:lm), intent(in)    :: et_b
    !> 2 つのスペクトルデータのヤコビアン
    real(DP), dimension(-km:km,0:lm)                :: et_Jacobian_et_et
    et_Jacobian_et_et = et_yx(                           &
      &  yx_et(et_Dx_et(et_a)) * yx_et(et_Dy_et(et_b))   &
      & - yx_et(et_Dy_et(et_a)) * yx_et(et_Dx_et(et_b)) )

  end function et_Jacobian_et_et


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

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

    !> 境界条件. 省略時は '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(et,values)
      else
        call at_Boundaries_DD(et)
      endif
      return
    endif

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

  end subroutine et_Boundaries

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

    real(DP), dimension(:,:,:), allocatable  :: alu
    integer,  dimension(:,:), allocatable    :: kp
    real(DP), dimension(-km:km,0:lm)         :: et_work
    real(DP), dimension(-km:km,0:jm)         :: ey_work
    real(DP), dimension(-km:km)              :: 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(-km:km,0:lm,0:lm),kp(-km:km,0:lm))

      do l=0,lm
        et_work=0.0D0
        et_work(:,l) = 1.0D0
        alu(:,:,l) = et_Lapla_et(et_work)
        ey_work = ay_at(et_work)
        alu(:,lm-1,l) = ey_work(:,0)
        alu(:,lm,l)   = ey_work(:,jm)
      enddo

      call ludecomp(alu,kp)
    endif

    et_work = et
    et_work(:,lm-1) = value1
    et_work(:,lm)   = value2
    et_LaplaInv_et = lusolve(alu,kp,et_work)

  end function et_LaplaInv_et

  !----------------------------------------------------------
  !> @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 ey_Vor2Strm_ey(ey,values,cond,new)
    !> 入力渦度分布
    real(DP), dimension(-km:km,0:jm),intent(in)  :: ey
    !> 流線境界値. 境界で一定なので波数 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(-km:km,0:jm)             :: ey_Vor2Strm_ey

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

    real(DP), dimension(-km:km,0:jm)         :: ey_work
    real(DP), dimension(-km:km,0:jm)         :: ey_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','ey_Vor2Strm_ey','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(-km:km,0:jm,0:jm),kp(-km:km,0:jm))

      do j=0,jm
        ey_I = 0.0D0
        ey_I(:,j) = 1.0D0
        alu(:,:,j) = ey_et(et_Lapla_et(et_ey(ey_I)))
      enddo

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

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

        if ( rigid1 ) then                               ! 粘着条件(top)
          ey_work=ey_et(et_Dy_et(et_ey(ey_I)))
        else
          ey_work=ey_et(et_Dy_et(et_Dy_et(et_ey(ey_I)))) ! すべり条件(top)
        endif

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

        if ( rigid2 ) then                               ! 粘着条件(bottom)
          ey_work=ey_et(et_Dy_et(et_ey(ey_I)))
        else
          ey_work=ey_et(et_Dy_et(et_Dy_et(et_ey(ey_I)))) ! すべり条件(bottom)
        endif
        alu(:,jm-1,j) = ey_work(:,jm)
      enddo

      call ludecomp(alu,kp)

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

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

    ey_work(:,0) = 0.0D0                  ! 運動学的条件. 波数 0 以外は 0
    ey_work(0,0) = value1*2.0D0           ! 運動学的条件. 波数 0 は重み 1/2

    ey_work(:,jm)   = 0.0D0               ! 運動学的条件. 波数 0 以外は 0
    ey_work(0,jm)   = value2*2.0D0        ! 運動学的条件. 波数 0 は重み 1/2

    ey_Vor2Strm_ey = lusolve(alu,kp,ey_work)

  end function ey_Vor2Strm_ey

  !----------------------------------------------------------
  !> @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 et_Vor2Strm_et(et,values,cond,new)
    !> 入力渦度分布
    real(DP), dimension(-km:km,0:lm),intent(in)  :: et
    !> 流線境界値. 境界で一定なので波数 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(-km:km,0:lm)             :: et_Vor2Strm_et

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

    real(DP), dimension(-km:km,0:lm)         :: et_work
    real(DP), dimension(-km:km,0:jm)         :: ey_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','et_Vor2Strm_et','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(-km:km,0:lm,0:lm),kp(-km:km,0:lm))

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

        alu(:,:,l) = et_Lapla_et(et_work)
      enddo

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

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

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

      call ludecomp(alu,kp)

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

    ! 内部領域計算
    et_work = et

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

    et_work(:,lm-1) = 0.0D0        ! 運動学的条件. 波数 0 以外は 0
    !!$      et_work(0,lm-1) = value2*2     ! 運動学的条件. 波数 0 は重み 1/2
    et_work(0,lm-1) = value2       ! 運動学的条件. 波数 0 は重み 1/2

    et_work(:,lm)   = 0.0D0        ! 運動学的条件. 波数 0 以外は 0
    !!$      et_work(0,lm)   = value1*2     ! 運動学的条件. 波数 0 は重み 1/2
    et_work(0,lm)   = value1     ! 運動学的条件. 波数 0 は重み 1/2

    et_Vor2Strm_et = lusolve(alu,kp,et_work)

  end function et_Vor2Strm_et

  !----------------------------------------------------------
  !> @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 et_Vor2Strm2_et(et,ea,cond,new)
    !> 入力渦度分布
    real(DP), dimension(-km:km,0:lm),intent(in)  :: et
    !> 流線境界値. 省略時は 0. 1:Top,2:Btm
    real(DP), dimension(-km:km,2), intent(in), optional :: ea
    !> 境界条件スイッチ. 省略時は 'RR'
    !> * RR : 両端粘着条件
    !> * RF : 上端粘着下端応力なし条件
    !> * FR : 上端応力なし下端粘着条件
    !> * FF : 両端応力なし条件
    character(len=2), intent(in), optional  :: cond
    !> true だと境界条件計算用行列を強制的に新たに作る.
    !> default は false.
    logical, intent(IN), optional :: new
    !> 出力流線分布
    real(DP), dimension(-km:km,0:lm)         :: et_Vor2Strm2_et

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

    real(DP), dimension(-km:km,0:lm)         :: et_work
    real(DP), dimension(-km:km,0:jm)         :: ey_work
    real(DP), dimension(-km:km)              :: e_Top
    real(DP), dimension(-km:km)              :: e_Btm
    logical                                  :: rigid1 = .true.
    logical                                  :: rigid2 = .true.

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

    if (.not. present(ea)) then
      e_Top = 0.0D0
      e_Btm = 0.0D0
    else
      e_Top = ea(:,1)
      e_Btm = ea(:,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','et_Vor2Strm_et','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(-km:km,0:lm,0:lm),kp(-km:km,0:lm))

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

        alu(:,:,l) = et_Lapla_et(et_work)
      enddo

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

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

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

      call ludecomp(alu,kp)

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

    ! 内部領域計算
    et_work = et

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

    et_work(:,lm-1) = e_Btm           ! 運動学的条件. 波数 0 以外は 0
    et_work(:,lm)   = e_Top           ! 運動学的条件. 波数 0 以外は 0

    et_Vor2Strm2_et = lusolve(alu,kp,et_work)

  end function et_Vor2Strm2_et

  !----------------------------------------------------------
  !> @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 et_Vor2Strm3_et(et,ea,cond,new)
    !> 入力渦度分布
    real(DP), dimension(-km:km,0:lm),intent(in)  :: et
    !> 流線/微分境界値. 省略時は 0.
    !> 1. Top 値,
    !> 2. Bottom 値
    !> 3. Top 微分(勾配)
    !> 4. Bottom 微分(勾配)
    real(DP), dimension(-km:km,4), intent(in), optional :: ea
    !> 境界条件スイッチ. 省略時は 'RR'
    !> * RR : 両端粘着条件
    !> * RF : 上端粘着下端応力なし条件
    !> * FR : 上端応力なし下端粘着条件
    !> * FF : 両端応力なし条件
    character(len=2), intent(in), optional  :: cond
    !> true だと境界条件計算用行列を強制的に新たに作る.
    !> default は false.
    logical, intent(IN), optional :: new
    !> 出力流線分布
    real(DP), dimension(-km:km,0:lm)             :: et_Vor2Strm3_et

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

    real(DP), dimension(-km:km,0:lm)        :: et_work
    real(DP), dimension(-km:km,0:jm)        :: ey_work
    real(DP), dimension(-km:km)             :: e_TopValue
    real(DP), dimension(-km:km)             :: e_BtmValue
    real(DP), dimension(-km:km)             :: e_TopDiff
    real(DP), dimension(-km:km)             :: e_BtmDiff
    logical                                 :: rigid1 = .true.
    logical                                 :: rigid2 = .true.
    logical :: first = .true.
    logical :: new_matrix = .false.
    integer :: l
    save    :: alu, kp, first

    if (.not. present(ea)) then
      e_TopValue = 0.0D0
      e_BtmValue = 0.0D0
      e_TopDiff  = 0.0D0
      e_BtmDiff  = 0.0D0
    else
      e_TopValue = ea(:,1)
      e_BtmValue = ea(:,2)
      e_TopDiff  = ea(:,3)
      e_BtmDiff  = ea(:,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','et_Vor2Strm_et','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(-km:km,0:lm,0:lm),kp(-km:km,0:lm))

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

        alu(:,:,l) = et_Lapla_et(et_work)
      enddo

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

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

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

      call ludecomp(alu,kp)

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

    ! 内部領域計算
    et_work = et

    et_work(:,lm-2) = e_TopDIff         ! 力学的条件
    et_work(:,lm-3) = e_BtmDiff         ! 力学的条件

    et_work(:,lm-1) = e_BtmValue        ! 運動学的条件. 波数 0 以外は 0
    et_work(:,lm)   = e_TopValue        ! 運動学的条件. 波数 0 以外は 0

    et_Vor2Strm3_et = lusolve(alu,kp,et_work)

  end function et_Vor2Strm3_et

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

    integer :: i, j

    IntYX_yx = 0.0d0
    do i=0,im-1
      do j=0,jm
        IntYX_yx = IntYX_yx + yx(j,i) * y_Y_Weight(j) * x_X_Weight(i)
      enddo
    enddo
  end function IntYX_yx

  !----------------------------------------------------------
  !> 2 次元格子点データの X 方向積分
  !----------------------------------------------------------
  function y_IntX_yx(yx)  ! X 方向積分
    !> 2 次元格子点データ
    real(DP), dimension(0:jm,0:im-1)   :: yx
    !> 積分された 1 次元(Y)格子点データ
    real(DP), dimension(0:jm)          :: y_IntX_yx

    integer :: i

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

  !----------------------------------------------------------
  !> 2 次元格子点データの Y 方向積分
  !----------------------------------------------------------
  function x_IntY_yx(yx)  ! Y 方向積分
    !> 2 次元格子点データ
    real(DP), dimension(0:jm,0:im-1)   :: yx
    !> 積分された 1 次元(X)格子点データ
    real(DP), dimension(0:im-1)        :: x_IntY_yx

    integer :: j
    ! 作業変数

    x_IntY_yx = 0.0d0
    do j=0,jm
      x_IntY_yx(:) = x_IntY_yx(:) + yx(j,:) * y_Y_Weight(j)
    enddo
  end function x_IntY_yx

  !----------------------------------------------------------
  !> 1 次元(X)格子点データの X 方向積分
  !----------------------------------------------------------
  function IntX_x(x)      ! X 方向積分
    !> 1 次元格子点データ
    real(DP), dimension(0:im-1)   :: x
    !> 積分値
    real(DP)                      :: IntX_x

    IntX_x = sum(x*x_X_Weight)
  end function IntX_x

  !----------------------------------------------------------
  !> 1 次元(Y)格子点データの Y 方向積分
  !----------------------------------------------------------
  function IntY_y(y)      ! Y 方向積分
    !> 1 次元格子点データ
    real(DP), dimension(0:jm)   :: y
    !> 積分値
    real(DP)                    :: IntY_y

    IntY_y = sum(y*y_Y_Weight)
  end function IntY_y

  !----------------------------------------------------------
  !> 2 次元格子点データの全領域平均
  !----------------------------------------------------------
  function AvrYX_yx(yx)
    !> 2 次元格子点データ
    real(DP), dimension(0:jm,0:im-1)   :: yx
    !> 平均値
    real(DP)                           :: AvrYX_yx

    AvrYX_yx = IntYX_yx(yx)/(sum(x_X_weight)*sum(y_Y_weight))
  end function AvrYX_yx

  !----------------------------------------------------------
  !>  2 次元格子点データの X 方向平均
  !----------------------------------------------------------
  function y_AvrX_yx(yx)
    !> 2 次元格子点データ
    real(DP), dimension(0:jm,0:im-1)   :: yx
    !> 平均された 1 次元(Y)格子点
    real(DP), dimension(0:jm)          :: y_AvrX_yx

    y_AvrX_yx = y_IntX_yx(yx)/sum(x_X_weight)
  end function y_AvrX_yx

  !----------------------------------------------------------
  ! 2 次元格子点データの Y 方向平均
  !----------------------------------------------------------
  function x_AvrY_yx(yx)
    !> 2 次元格子点データ
    real(DP), dimension(0:jm,0:im-1)   :: yx
    !> 平均された 1 次元(X)格子点
    real(DP), dimension(0:im-1)        :: x_AvrY_yx

    x_AvrY_yx = x_IntY_yx(yx)/sum(y_Y_weight)
  end function x_AvrY_yx

  !----------------------------------------------------------
  !> 1 次元(X)格子点データの X 方向平均
  !----------------------------------------------------------
  function AvrX_x(x)
    !> 1 次元格子点データ
    real(DP), dimension(0:im-1)   :: x
    !> 平均値
    real(DP)                      :: AvrX_x

    AvrX_x = IntX_x(x)/sum(x_X_weight)
  end function AvrX_x

  !----------------------------------------------------------
  !> 1 次元(Y)格子点データの Y 方向平均
  !----------------------------------------------------------
  function AvrY_y(y)
    !> 1 次元格子点データ
    real(DP), dimension(0:jm)   :: y
    !> 平均値
    real(DP)                    :: AvrY_y

    AvrY_y = IntY_y(y)/sum(y_Y_weight)
  end function AvrY_y

  !----------------------------------------------------------
  !> 補間関数
  !----------------------------------------------------------
  function Interpolate_et(et_data,xval,yval)
    !> 入力スペクトルデータ
    real(DP), dimension(-km:km,0:lm), intent(IN) :: et_data
    !> 補間する点の X 座標
    real(DP), intent(in)                         :: xval
    !> 補間する点の Y 座標
    real(DP), intent(in)                         :: yval
    !> 補間結果
    real(DP)                                     :: Interpolate_et

    Interpolate_et = Interpolate_e(a_Interpolate_at(et_data,yval),xval)
  end function Interpolate_et

end module et_module
