!--
!----------------------------------------------------------------------
! Copyright(c) 2008-2013 SPMDODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!
!表題  eq_module
!      2 次元円盤領域問題, Fourier 展開 + 多項式展開法
!
!      spml/eq_module モジュールは 2 次元円盤領域での流体運動を
!      スペクトル法により数値計算を実行するための Fortran90 関数を提供する.
!      周期的な境界条件を扱うための方位角方向へのフーリエ変換と
!      境界壁を扱うための動径方向の多項式変換を用いる場合の
!      スペクトル計算のためのさまざまな関数を提供する.
!
!      内部で ae_module, aq_module を用いている.
!      最下部ではフーリエ変換およびチェビシェフ変換のエンジンとして
!      ISPACK/FTPACK の Fortran77 サブルーチンを用いている.
!
!      Matsushima and Marcus (1994) の多項式に関する説明は
!      doc/spectral_radial.tex を参照のこと.
!
!
!履歴  2008/04/11  竹広真一  et_module より改造
!      2008/05/02  竹広真一  コメント修正
!      2008/10/29  竹広真一  eq_Vor2Strm_eq 行列生成ルーチン改良
!      2009/01/09  竹広真一  eq_Initial メッセージに日付を追加
!      2009/01/29  佐々木洋平 コメントを RDoc 用に修正
!      2009/07/31  竹広真一  境界条件計算用配列を threadprivate 指定(OpenMP)
!      2010/03/10  佐々木洋平  threadprivate 削除(コンパイラ依存)
!      2011/03/07  佐々木洋平  境界条件用論理変数を初期化, if 文削除
!      2013/08/20  竹広真一  gnu fortran 対応
!
!++
module eq_module
  !
  != eq_module
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/eq_module モジュールは 2 次元円盤領域での流体運動を
  ! スペクトル法により数値計算を実行するための Fortran90 関数を提供する.
  ! 周期的な境界条件を扱うための方位角方向へのフーリエ変換と
  ! 境界壁を扱うための動径方向の多項式変換を用いる場合の
  ! スペクトル計算のためのさまざまな関数を提供する.
  !
  ! 内部で ae_module, aq_module を用いている.
  ! 最下部ではフーリエ変換およびチェビシェフ変換のエンジンとして
  ! ISPACK/FTPACK の Fortran77 サブルーチンを用いている.
  !
  ! Matsushima and Marcus (1994) の多項式に関する説明は
  ! 動径座標のスペクトル法(spectral_radial.pdf[link:./spectral_radial.pdf])を
  ! 参照のこと.
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (eq_, rp_, r_, p_) は, 返す値の形を示している.
  !   eq_ :: 2次元スペクトルデータ
  !   rp_ :: 2 次元格子点データ
  !   r_  :: 動径方向 1 次元格子点データ
  !   p_  :: 方位角方向 1 次元格子点データ
  !
  ! * 関数名の間の文字列(Dr, Dp, Lapla, LaplaInv, Jacobian)は,
  !   その関数の作用を表している.
  !
  ! * 関数名の最後 (_eq_eq,_eq,_rp, _r, _p) は, 入力変数のスペクトルデータ
  !   および格子点データであることを示している.
  !   _eq    :: 2次元スペクトルデータ
  !   _eq_eq :: 2 つの2次元スペクトルデータ
  !   _rp    :: 2 次元格子点データ
  !   _r     :: 動径方向 1 次元格子点データ
  !   _p     :: 方位角方向 1 次元格子点データ
  !
  !=== 各データの種類の説明
  !
  ! * rp : 2 次元格子点データ.
  !   * 変数の種類と次元は real(8), dimension(jm,0:im-1).
  !   * im, jm はそれぞれ方位角, 動径座標の格子点数であり,
  !     サブルーチン eq_initial にてあらかじめ設定しておく.
  !   * 第 1 次元が動径座標の格子点位置番号, 第 2 次元が方位角座標の
  !     格子点位置番号である.
  !
  ! * eq : 2 次元スペクトルデータ.
  !   * 変数の種類と次元は real(8), dimension(-km:km,0:lm).
  !   * km, lm はそれぞれ方位角, 動径方向の最大波数であり,
  !     サブルーチン eq_initial にてあらかじめ設定しておく.
  !   * 動径スペクトルデータの格納のされ方については
  !     aq_module.f90 を参照のこと.
  !
  ! * p, r : X, Y 方向 1 次元格子点データ.
  !   * 変数の種類と次元はそれぞれ real(8), dimension(0:im-1)
  !     および real(8), dimension(jm).
  !
  ! * e, q : 1 次元スペクトルデータ.
  !   * 変数の種類と次元は real(8), dimension(-km:km)
  !     および real(8), dimension(0:lm).
  !
  ! * ap, ar : 1 次元格子点データの並んだ 2 次元配列.
  !   * 変数の種類と次元は real(8), dimension(:,0:im-1)
  !     および real(8), dimension(:,jm).
  !
  ! * ae, aq : 1 次元スペクトルデータの並んだ 2 次元配列.
  !   * 変数の種類と次元は real(8), dimension(:,-km:km)
  !     および real(8), dimension(:,0:lm).
  !
  ! * eq_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * rp_ で始まる関数が返す値は 2 次元格子点データに同じ.
  !
  ! * p_, p_ で始まる関数が返す値は 1 次元格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! eq_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  !
  !==== 座標変数
  !
  ! p_Phi, r_Rad               ::  格子点座標(X,Y座標)を格納した 1 次元配列
  ! p_Phi_Weight, r_Rad_Weight ::  重み座標を格納した 1 次元配列
  ! rp_Phi, rp_Rad             ::  格子点データの XY 座標(X,Y)
  !                                (格子点データ型 2 次元配列)
  !
  !==== 基本変換
  !
  ! rp_eq :: スペクトルデータから格子データへの変換
  ! eq_rp :: 格子データからスペクトルデータへの変換
  ! ap_ae, p_e :: 方位角方向のスペクトルデータから格子データへの変換
  ! ar_aq, r_q :: 動径方向のスペクトルデータから格子データへの変換
  ! ae_ap, e_p :: 方位角方向の格子点データからスペクトルデータへの変換
  ! aq_ar, q_r :: 動径方向の格子点データからスペクトルデータへの変換
  !
  !==== 微分
  !
  ! eq_Lapla_eq  :: スペクトルデータにラプラシアンを作用させる
  ! eq_DPhi_eq, ae_DPhi_ae, e_DPhi_e :: スペクトルデータに
  !                                     方位角微分を作用させる
  ! eq_RadDRad_eq, aq_RadDRad_aq, q_RadDRad_q :: スペクトルデータに
  !                                              動径微分を作用させる
  ! eq_Jacobian_eq_eq :: 2 つのスペクトルデータからヤコビアンを計算する
  !
  !==== 境界値問題
  !
  ! eq_Boundary    :: ディリクレ, ノイマン境界条件の適用
  ! eq_LaplaInv_eq :: スペクトルデータにラプラシアンの逆変換を作用させる
  ! eq_Vor2Strm_eq :: 渦度から流線を計算する
  !
  !==== 積分・平均
  !
  ! IntRadPhi_rp, AvrRadPhi_rp   :: 2 次元格子点データの全領域積分および平均
  ! r_IntPhi_rp, r_AvrPhi_rp :: 2 次元格子点データの方位角方向積分および平均
  ! IntPhi_p, AvrPhi_p       :: 1 次元(X)格子点データの方位角方向積分および平均
  ! p_IntRad_rp, p_AvrRad_rp :: 2 次元格子点データの動径方向積分および平均
  ! IntRad_r, AvrRad_r       :: 1 次元(Y)格子点データの動径方向積分および平均
  !
  use dc_message
  use lumatrix
  use ae_module, p_Phi => g_X, p_Phi_weight => g_X_Weight, &
                 e_p => e_g, ae_ap => ae_ag, &
                 p_e => g_e, ap_ae => ag_ae, &
                 ae_DPhi_ae => ae_Dx_ae, e_DPhi_e => e_Dx_e
  use aq_module, r_Rad => g_R, r_Rad_Weight => g_R_Weight, &
                 aq_ar => aq_ag, q_r => q_g, &
                 ar_aq => ag_aq, r_q => g_q, &
                 q_RadDRad_q => q_rDr_q, aq_RadDRad_aq => aq_rDr_aq

  implicit none
  private

  public eq_Initial                                       ! 初期化

  public p_Phi, p_Phi_Weight, rp_Phi                      ! 座標変数
  public r_Rad, r_Rad_Weight, rp_Rad, er_Rad              ! 座標変数

  public rp_eq, eq_rp                                     ! 基本変換
  public er_eq, eq_er                                     ! 基本変換
  public er_rp, rp_er                                     ! 基本変換
  public e_p, p_e, ae_ap, ap_ae                           ! 基本変換
  public q_r, r_q, aq_ar, ar_aq                           ! 基本変換

  public eq_DPhi_eq, e_DPhi_e, ae_DPhi_ae                 ! 微分
  public eq_RadDRad_eq, q_RadDRad_q, aq_RadDRad_aq        ! 微分
  public er_Lapla_eq, eq_Lapla_eq                         ! 微分

  public eq_Jacobian_eq_eq                                ! 非線形計算

  public eq_Boundary                                      ! 境界値問題
  public aq_Boundary_D, aq_Boundary_N                     ! 境界値問題
  public eq_LaplaInv_eq, eq_Vor2Strm_eq                   ! 境界値問題

  public IntRadPhi_rp, r_IntPhi_rp, p_IntRad_rp, IntPhi_p, IntRad_r   ! 積分
  public AvrRadPhi_rp, r_AvrPhi_rp, p_AvrRad_rp, AvrPhi_p, AvrRad_r   ! 平均

  integer            :: im=32, jm=8      ! 格子点の設定(Phi,Rad)
  integer            :: km=10, lm=5      ! 切断波数の設定(Phi,Rad)
  real(8)            :: ra=1.0           ! 領域の大きさ
  real(8), parameter :: pi=3.1415926535897932385D0

  real(8), parameter :: alpha = 1.0D0        ! 展開多項式パラメター  0 < α <= 1
  real(8), parameter :: beta  = 1.0D0        ! 展開多項式パラメター  0 < β

  real(8), dimension(:,:), allocatable :: rp_Phi, rp_Rad
  real(8), dimension(:,:), allocatable :: er_Rad
  integer, dimension(:), allocatable   :: md

  save im, jm, km, lm, ra, md

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

      real(8),intent(in) :: ra_in       ! 半径

      integer :: kk

      im = i       ; jm = j
      km = k       ; lm = l
      ra = ra_in

      allocate(md(-km:km))

      do kk=-km,km
         md(kk) = abs(kk)
      enddo

      call ae_initial(im,km,0.0D0,2*pi)
      call aq_Initial(jm,lm,ra,alpha,beta,md)

      allocate(rp_Phi(jm,0:im-1),rp_Rad(jm,0:im-1))
      rp_Phi = spread(p_Phi,1,jm)
      rp_Rad = spread(r_Rad,2,im)

      allocate(er_Rad(-km:km,jm))
      er_Rad = spread(r_Rad,1,2*km+1)

      call MessageNotify('M','eq_initial','eq_module (2013/08/20) is initialized')
    end subroutine eq_initial

  !--------------- 基本変換 -----------------

    function rp_eq(eq)
      !
      ! スペクトルデータから格子データへ変換する.
      !
      real(8), dimension(jm,0:im-1)                :: rp_eq
      !(out) 格子点データ

      real(8), dimension(-km:km,0:lm), intent(in)  :: eq
      !(in) スペクトルデータ

      real(8), dimension(-km:km,jm)                :: er
      real(8), dimension(jm,-km:km)                :: re

      integer :: j, k
      
      !rp_eq = ap_ae(transpose(ar_aq(eq)))
      !
      ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
      !
      er = ar_aq(eq)
      do j=1,jm
         do k=-km,km
            re(j,k) = er(k,j)
         enddo
      enddo
      rp_eq = ap_ae(re)

    end function rp_eq

    function eq_rp(rp)
      !
      ! 格子データからスペクトルデータへ変換する.
      !
      real(8), dimension(-km:km,0:lm)              :: eq_rp
      !(out) スペクトルデータ

      real(8), dimension(jm,0:im-1), intent(in)    :: rp
      !(in) 格子点データ

      real(8), dimension(-km:km,jm)                :: er
      real(8), dimension(jm,-km:km)                :: re

      integer :: j, k
      
      !eq_rp = aq_ar(transpose(ae_ap(rp)))
      !
      ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
      !
      re = ae_ap(rp)
      do j=1,jm
         do k=-km,km
            er(k,j) = re(j,k)
         enddo
      enddo
      eq_rp = aq_ar(er)

    end function eq_rp

    function eq_er(er)
      !
      ! 動径格子データからスペクトルデータへ変換する.
      !
      real(8), dimension(-km:km,0:lm)              :: eq_er
      !(out) スペクトルデータ

      real(8), dimension(-km:km,jm), intent(in)    :: er
      !(in) 格子点データ

      eq_er = aq_ar(er)

    end function eq_er

    function er_eq(eq)
      !
      ! 動径格子データからスペクトルデータへ変換する.
      !
      real(8), dimension(-km:km,jm)                :: er_eq
      !(out) 格子点データ

      real(8), dimension(-km:km,0:lm), intent(in)  :: eq
      !(out) スペクトルデータ

      er_eq = ar_aq(eq)

    end function er_eq

    function rp_er(er)
      !
      ! スペクトルデータから格子データへ変換する.
      !
      real(8), dimension(jm,0:im-1)              :: rp_er
      !(out) 格子点データ

      real(8), dimension(-km:km,jm), intent(in)  :: er
      !(in) スペクトルデータ

      real(8), dimension(jm,-km:km)              :: re

      integer :: j, k

      !rp_er = ap_ae(transpose(er))
      !
      ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
      !
      do j=1,jm
         do k=-km,km
            re(j,k) = er(k,j)
         enddo
      enddo
      rp_er = ap_ae(re)

    end function rp_er

    function er_rp(rp)
      !
      ! 格子データからスペクトルデータへ変換する.
      !
      real(8), dimension(-km:km,jm)                :: er_rp
      !(out) スペクトルデータ

      real(8), dimension(jm,0:im-1), intent(in)    :: rp
      !(in) 格子点データ

      real(8), dimension(jm,-km:km)                :: re

      integer :: j, k

      !er_rp = transpose(ae_ap(rp))
      !
      ! 本来上の書き方でいいはずだが, gfortran で通らないので下のように書き直す
      !
      re = ae_ap(rp)
      do j=1,jm
         do k=-km,km
            er_rp(k,j) = re(j,k)
         enddo
      enddo

    end function er_rp

  !--------------- 微分計算 -----------------

    function eq_DPhi_eq(eq)
      !
      ! 入力スペクトルデータに方位角微分(∂φ)を作用する.
      !
      ! スペクトルデータのφ微分とは, 対応する格子点データにφ微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      ! 実際にはスペクトルデータに X 方向波数 k をかけて
      ! sin(kx) <-> cos(kx) 成分に入れ換える計算を行っている.
      !
      real(8), dimension(-km:km,0:lm)                :: eq_DPhi_eq
      real(8), dimension(-km:km,0:lm), intent(in)    :: eq
      integer k

      do k=-km,km
         eq_DPhi_eq(k,:)  =  -k*eq(-k,:)
      enddo
    end function eq_DPhi_eq

    function eq_RadDRad_eq(eq)
      !
      ! 入力スペクトルデータに動径微分(r∂r)を作用する.
      !
      ! スペクトルデータの動径微分とは, 対応する格子点データに動径微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(-km:km,0:lm)               :: eq_RadDRad_eq
      !(out) スペクトルデータの動径微分

      real(8), dimension(-km:km,0:lm), intent(in)   :: eq
      !(in) 入力スペクトルデータ

      eq_RadDRad_eq = aq_RadDRad_aq(eq)

    end function eq_RadDRad_eq

    function er_Lapla_eq(eq)
      !
      ! 入力スペクトルデータにラプラシアン
      !  (1/r)(∂r(r∂r)+ (1/r^2) ∂φφ を作用する.
      !
      real(8), dimension(-km:km,jm)                :: er_Lapla_eq
      !(out) スペクトルデータのラプラシアン

      real(8), dimension(-km:km,0:lm), intent(in)  :: eq
      !(in) 入力スペクトルデータ

      real(8), dimension(-km:km,0:lm)              :: eq_work

      integer k

      do k=-km,km
         eq_work(k,:) = -k**2 * eq(k,:)
      enddo

      er_Lapla_eq = er_eq(eq_work + eq_RadDRad_eq(eq_RadDRad_eq(eq)))/er_Rad**2

    end function er_Lapla_eq

    function eq_Lapla_eq(eq)
      !
      ! 入力スペクトルデータにラプラシアン
      !  (1/r)(∂r(r∂r)+ (1/r^2) ∂φφ を作用する.
      !
      ! スペクトルデータのラプラシアンとは,
      ! 対応する格子点データにラプラシアンを作用させたデータの
      ! スペクトル変換のことである.
      !
      real(8), dimension(-km:km,0:lm)              :: eq_Lapla_eq
      !(out) スペクトルデータのラプラシアン

      real(8), dimension(-km:km,0:lm), intent(in)  :: eq
      !(in) 入力スペクトルデータ

      eq_Lapla_eq = eq_er(er_Lapla_eq(eq))

    end function eq_Lapla_eq

    function eq_Jacobian_eq_eq(eq_a,eq_b)
      !
      !  2 つのスペクトルデータからヤコビアン
      !
      !     J(A,B)=1/r[(∂rA)(∂φB)-(∂φA)(∂rB)]
      !
      !  を計算する. 1/r のファクターがついていることに注意.
      !
      !  2 つのスペクトルデータのヤコビアンとは, 対応する 2 つの
      !  格子点データのヤコビアンのスペクトル変換のことである.
      !
      real(8), dimension(-km:km,0:lm)                :: eq_Jacobian_eq_eq
      !(out) 2 つのスペクトルデータのヤコビアン

      real(8), dimension(-km:km,0:lm), intent(in)    :: eq_a
      !(in) 1つ目の入力スペクトルデータ

      real(8), dimension(-km:km,0:lm), intent(in)    :: eq_b
      !(in) 2つ目の入力スペクトルデータ

      eq_Jacobian_eq_eq = eq_rp(&
           (  rp_eq(eq_RadDRad_eq(eq_a)) * rp_eq(eq_DPhi_eq(eq_b))   &
             -rp_eq(eq_DPhi_eq(eq_a)) * rp_eq(eq_RadDRad_eq(eq_b)) ) &
           /rp_Rad**2)

    end function eq_Jacobian_eq_eq


  !--------------- 境界値問題 -----------------

    subroutine eq_Boundary(eq,value,cond)
      !
      ! ディリクレ, ノイマン条件の適用. チェビシェフ空間での計算
      !
      ! 実際には中で呼ばれている aq_module のサブルーチン
      ! aq_Boundary_D,, aq_Boundary_N を用いている.
      ! これらを直接呼ぶことも出来る.
      !
      real(8), dimension(-km:km,0:lm),intent(inout)      :: eq
              ! 境界条件を適用するデータ. 修正された値を返す.

      real(8), dimension(-km:km), intent(in), optional   :: value
              ! 境界での 値/勾配 分布を水平スペクトル変換したものを与える.
              ! 省略時は値/勾配 0 となる.

      character(len=1), intent(in), optional             :: cond
              ! 境界条件. 省略時は 'D'
              !   D : 両端ディリクレ
              !   N : 両端ノイマン

      if (.not. present(cond)) then
         if (present(value)) then
            call aq_Boundary_D(eq,value)
         else
            call aq_Boundary_D(eq)
         endif
         return
      endif

      select case(cond)
      case ('N')
         if (present(value)) then
            call aq_Boundary_N(eq,value)
         else
            call aq_Boundary_N(eq)
         endif
      case ('D')
         if (present(value)) then
            call aq_Boundary_D(eq,value)
         else
            call aq_Boundary_D(eq)
         endif
      case default
         call MessageNotify('E','eq_Boundaries','B.C. not supported')
      end select

    end subroutine eq_Boundary

    function eq_LaplaInv_eq(eq,value)
      !
      ! 境界で値を与える条件(ディリクレ条件)下で,
      ! 入力スペクトルデータに逆ラプラシアン
      ! [(1/r)(∂r(r∂r)+ (1/r^2) ∂φφ]^{-1} を作用する.
      !
      ! タウ法による計算
      !
      ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに
      ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(-km:km,0:lm),intent(in)  :: eq
      !(in) スペクトルデータ

      real(8), dimension(-km:km,0:lm)             :: eq_LaplaInv_eq
      !(out) スペクトルデータの逆ラプラシアン

      real(8), dimension(-km:km), intent(in), optional :: value
      !(in) 境界値. 省略時は 0 が設定される.

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

      real(8), dimension(-km:km,0:lm)         :: eq_work
      real(8), dimension(-km:km,jm)           :: er_work
      real(8), dimension(-km:km)              :: value1       ! 境界値

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

      if (.not. present(value)) then
         value1=0
      else
         value1 = value
      endif

      if ( first ) then
         first = .false.

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

         do l=0,lm
            eq_work = 0.0 ; eq_work(:,l) = 1.0
            alu(:,:,l) = eq_er(er_Lapla_eq(eq_work))
         enddo

         ! 0 成分のところを 1 で埋める.
         do k=-km,km
            do l=0,md(k)-1
               alu(k,l,l) = 1.0D0
            enddo
            do l=md(k)+1,lm,2
               alu(k,l,l) = 1.0D0
            enddo
         enddo

         ! 境界条件 r=ra で値を与える.
         do k=-km,km
            do l=0,lm
               eq_work=0 ; eq_work(k,l)=1.0
               er_work=er_eq(eq_work)
               if ( mod(md(k),2) .eq. mod(lm,2) ) then
                  alu(k,lm,l) = er_work(k,jm)
               else
                  alu(k,lm-1,l) = er_work(k,jm)
               endif
            enddo
         enddo

         call ludecomp(alu,kp)
      endif

      eq_work = eq
      do k=-km,km
         if ( mod(md(k),2) .eq. mod(lm,2) ) then
            eq_work(k,lm)   = value1(k)
         else
            eq_work(k,lm-1) = value1(k)
         endif
      enddo
      eq_LaplaInv_eq = lusolve(alu,kp,eq_work)

    end function eq_LaplaInv_eq

    function eq_Vor2Strm_eq(eq,value,cond,new)
      !
      ! 渦度から流線を求める.
      !
      ! Chebyshev-tau 法による計算
      ! 渦度 \zeta を与えて流線 \psi を求める.
      !    \nabla^2 \psi = \zeta,
      !    \psi = const. at the boundary
      ! 粘着条件
      !    \DP{\psi}{r} = 0 at the boundary
      ! 応力なし条件
      !    r\DP{}{r}(1/r\DP{\psi}{r})  = 0 at the boundary
      !
      ! l=0,lm 成分の式の代わりに境界条件を与える.
      !
      real(8), dimension(-km:km,0:lm),intent(in)  :: eq
              !(in) 入力渦度分布

      real(8), dimension(-km:km,0:lm)             :: eq_Vor2Strm_eq
              !(out) 出力流線関数分布

      real(8), intent(in), optional               :: value
              ! 流線境界値. 境界で一定なので波数 0 成分のみ

      character(len=1), intent(in), optional  :: cond
              !(in) 境界条件スイッチ. 省略時は 'R'
              !     R    : 上側粘着条件
              !     F    : 上側応力なし条件

      logical, intent(IN), optional :: new
              !(in) true だと境界条件計算用行列を強制的に新たに作る.
              !     default は false.

      real(8), dimension(:,:,:), allocatable  :: alu, alub
      integer, dimension(:,:), allocatable    :: kp, kpb

      real(8), dimension(-km:km,0:lm)         :: eq_work
      real(8), dimension(-km:km,jm)           :: er_work
      real(8)                                 :: value1          ! 境界値
      logical                                 :: rigid = .true.

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

      if (.not. present(value)) then
         value1 = 0
      else
         value1 = value
      endif

      select case (cond)
      case ('R')
        rigid = .TRUE.
      case ('F')
        rigid = .FALSE.
      case default
        call MessageNotify('E','eq_Vor2Strm_eq','B.C. not supported')
      end select

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

         if ( allocated(alub) ) deallocate(alub)
         if ( allocated(kpb) ) deallocate(kpb)
         allocate(alub(-km:km,0:lm,0:lm),kpb(-km:km,0:lm))

         ! 内部条件
         do l=0,lm
            eq_work = 0.0 ; eq_work(:,l) = 1.0
            alu(:,:,l) = eq_er(er_Lapla_eq(eq_work))
         enddo

         ! 0 成分のところを 1 で埋める.
         do k=-km,km
            do l=0,md(k)-1
               alu(k,l,l) = 1.0D0
            enddo
            do l=md(k)+1,lm,2
               alu(k,l,l) = 1.0D0
            enddo
         enddo

         ! alu(:,:,nd(k)) 列は 0 なので 1 をいれておく.
         ! l=md(k) 成分は境界条件で決める.
         do k=-km,km
            if ( mod(md(k),2) .eq. mod(lm,2) ) then
               alu(k,lm,md(k)) = 1.0D0
            else
               alu(k,lm-1,md(k)) = 1.0D0
            endif
         enddo

         call ludecomp(alu,kp)

         !---- 境界条件計算用行列 -----
         alub = 0.0
         do l=0,lm
            alub(:,l,l) = 1.0D0
         enddo

         ! 運動学的条件. 流線は境界で一定
         !     l=nd(n) 成分を境界条件で決める.
         do l=0,lm
            eq_work = 0.0 ; eq_work(:,l)=1.0D0
            er_work = er_eq(eq_work)
            do k=-km,km
               alub(k,md(k),l) = er_work(k,jm)
            enddo
         enddo

         ! 力学的条件粘着条件
         !     l=lm or lm-1 成分を境界条件で決める.
         if ( rigid ) then
            do l=0,lm
               eq_work = 0.0 ;  eq_work(:,l)=1.0D0
               er_work=er_eq(eq_RadDRad_eq(eq_work))/er_Rad
               do k=-km,km
                  if ( mod(md(k),2) .eq. mod(lm,2) ) then
                     alub(k,lm,l) = er_work(k,jm)
                  else
                     alub(k,lm-1,l) = er_work(k,jm)
                  endif
               end do
            enddo
         else
            do l=0,lm
               eq_work = 0.0 ; eq_work(:,l)=1.0D0
               er_work=er_eq(eq_RadDRad_eq(eq_RadDRad_eq(eq_work)) &
                            -2*eq_RadDRad_eq(eq_work))/er_Rad**2
               do k=-km,km
                  if ( mod(md(k),2) .eq. mod(lm,2) ) then
                     alub(k,lm,l) = er_work(k,jm)
                  else
                     alub(k,lm-1,l) = er_work(k,jm)
                  endif
               end do
            enddo
         endif

         call ludecomp(alub,kpb)

         if ( rigid ) then
            call MessageNotify('M','eq_Vor2Strm_eq',&
                              'Matrix to apply rigid b.c. newly produced.')
         else
            call MessageNotify('M','eq_Vor2Strm_eq',&
                              'Matrix to apply stress-free b.c. newly produced.')
         endif
      endif

      ! 内部領域計算
      eq_work = eq
      eq_work = lusolve(alu,kp,eq_work)

      ! 境界条件計算
      do k=-km,km
         eq_work(k,md(k)) = 0.0D0
         if ( mod(md(k),2) .eq. mod(lm,2) ) then
            eq_work(k,lm)   = 0.0D0
         else
            eq_work(k,lm-1) = 0.0D0
         endif
      enddo
      eq_work(0,0)   = value1*2.0D0     ! 運動学的条件. 波数 0 は重み 1/2

      eq_Vor2Strm_eq = lusolve(alub,kpb,eq_work)

    end function eq_Vor2Strm_eq

  !--------------- 積分計算 -----------------
    function IntRadPhi_rp(rp)
      !
      ! 2 次元格子点データの全領域積分および平均.
      !
      ! 実際には格子点データ各点毎に p_Phi_Weight, r_Rad_Weight をかけた
      ! 総和を計算している.
      !
      real(8), dimension(jm,0:im-1)   :: rp
      !(in)  2 次元格子点データ

      real(8)                         :: IntRadPhi_rp
      !(out) 積分値

      integer :: i, j

      IntRadPhi_rp = 0.0d0
      do i=0,im-1
         do j=1,jm
            IntRadPhi_rp = IntRadPhi_rp + rp(j,i) * r_Rad_Weight(j) * p_Phi_Weight(i)
         enddo
      enddo
    end function IntRadPhi_rp

    function r_IntPhi_rp(rp)
      !
      ! 2 次元格子点データの Phi 方向積分
      !
      ! 実際には格子点データ各点毎に p_Phi_Weight をかけた総和を計算している.
      !
      real(8), dimension(jm,0:im-1)   :: rp
      !(in) 2 次元格子点データ

      real(8), dimension(jm)          :: r_IntPhi_rp
      !(out) 積分された 1 次元(Rad)格子点データ

      integer :: i
      ! 作業変数

      r_IntPhi_rp = 0.0d0
      do i=0,im-1
         r_IntPhi_rp(:) = r_IntPhi_rp(:) + rp(:,i) * p_Phi_Weight(i)
      enddo
    end function r_IntPhi_rp

    function p_IntRad_rp(rp)
      !
      ! 2 次元格子点データの Rad 方向積分
      !
      ! 実際には格子点データ各点毎に r_Rad_Weight をかけた総和を計算している.
      !
      real(8), dimension(jm,0:im-1)   :: rp
      !(in)  2 次元格子点データ

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

      integer :: j
      ! 作業変数

      p_IntRad_rp = 0.0d0
      do j=1,jm
         p_IntRad_rp(:) = p_IntRad_rp(:) + rp(j,:) * r_Rad_Weight(j)
      enddo
    end function p_IntRad_rp

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

      IntPhi_p = sum(p*p_Phi_Weight)
    end function IntPhi_p

    function IntRad_r(r)
      !
      ! 1 次元(Rad)格子点データの Rad 方向積分
      !
      ! 実際には格子点データ各点毎に r_Rad_Weight をかけた総和を計算している.
      !
      real(8), dimension(jm)   :: r          !(in)  1 次元格子点データ
      real(8)                  :: IntRad_r     !(out) 積分値

      IntRad_r = sum(r*r_Rad_Weight)
    end function IntRad_r

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

      real(8)                         :: AvrRadPhi_rp
      !(out) 平均値

      AvrRadPhi_rp = IntRadPhi_rp(rp)/(sum(p_Phi_weight)*sum(r_Rad_weight))
    end function AvrRadPhi_rp

    function r_AvrPhi_rp(rp)
      !
      ! 2 次元格子点データの Phi 方向平均
      !
      ! 実際には格子点データ各点毎に p_Phi_Weight をかけた総和を計算し,
      ! p_Phi_Weight の総和で割ることで平均している.
      !
      real(8), dimension(jm,0:im-1)   :: rp
      !(in) 2 次元格子点データ

      real(8), dimension(jm)          :: r_AvrPhi_rp
      !(out) 平均された 1 次元(Rad)格子点

      r_AvrPhi_rp = r_IntPhi_rp(rp)/sum(p_Phi_weight)
    end function r_AvrPhi_rp

    function p_AvrRad_rp(rp)
      !
      ! 2 次元格子点データの Rad 方向平均
      !
      ! 実際には格子点データ各点毎に r_Rad_Weight をかけた総和を計算し,
      ! r_Rad_Weight の総和で割ることで平均している.
      !
      real(8), dimension(jm,0:im-1)   :: rp
      !(in) 2 次元格子点データ

      real(8), dimension(0:im-1)      :: p_AvrRad_rp
      !(out) 平均された 1 次元(Phi)格子点

      p_AvrRad_rp = p_IntRad_rp(rp)/sum(r_Rad_weight)
    end function p_AvrRad_rp

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

      AvrPhi_p = IntPhi_p(p)/sum(p_Phi_weight)
    end function AvrPhi_p

    function AvrRad_r(r)
      !
      ! 1 次元(Rad)格子点データの Rad 方向平均
      !
      ! 実際には格子点データ各点毎に r_Rad_Weight をかけた総和を計算し,
      ! r_Rad_Weight の総和で割ることで平均している.
      !
      real(8), dimension(jm)   :: r            !(in)  1 次元格子点データ
      real(8)                  :: AvrRad_r     !(out) 平均値

      AvrRad_r = IntRad_r(r)/sum(r_Rad_weight)
    end function AvrRad_r

end module eq_module
