function AvrYX_yx(yx)    ! 全領域平均
      real(8), dimension(0:jm,0:im-1)   :: yx          ! 2 次元格子点
      real(8)                           :: AvrYX_yx    ! 平均値

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

    function IntYX_yx(yx)   ! 全領域積分
      real(8), dimension(0:jm,0:im-1)   :: yx          ! 2 次元格子点
      real(8)                           :: 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

    subroutine et_Boundaries(et,values,cond)   ! ディリクレ, ノイマン条件
      ! Chebyshev 空間での境界条件適用

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

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

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

      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

    function et_Dx_et(et)   ! スペクトルに作用する x 微分演算子
      real(8), dimension(-km:km,0:lm)                :: et_Dx_et
      real(8), dimension(-km:km,0:lm), intent(in)    :: et
      integer k

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

    function et_Dy_et(et)   ! スペクトルに作用する y 微分演算子
      real(8), dimension(-km:km,0:lm)               :: et_Dy_et
      real(8), dimension(-km:km,0:lm), intent(in)   :: et

      et_Dy_et = at_Dy_at(et)

    end function et_Dy_et

    function et_Jacobian_et_et(et_a,et_b) ! スペクトルに作用する Jacobian
      real(8), dimension(-km:km,0:lm)                :: et_Jacobian_et_et
      real(8), dimension(-km:km,0:lm), intent(in)    :: et_a
      real(8), dimension(-km:km,0:lm), intent(in)    :: et_b
      integer k

      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

    function et_LaplaInv_et(et,values) ! 逆 Laplacian, Dirichlet 型境界条件
      ! Chebyshev-tau 法による計算

      real(8), dimension(-km:km,0:lm),intent(in)  :: et
      real(8), dimension(-km:km,0:lm)             :: et_LaplaInv_et
      real(8), dimension(-km:km,2), intent(in), optional :: values   ! 境界値

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

      real(8), dimension(-km:km,0:lm)         :: et_work
      real(8), dimension(0:lm,0:lm)           :: tt_work
      real(8), dimension(0:lm,0:jm)           :: ty_work
      real(8), dimension(-km:km)              :: value1, value2   ! 境界値

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

      if (.not. present(values)) then
         value1=0 
 value2=0
      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))

         tt_work=0
         do l=0,lm
            tt_work(l,l)=1
         enddo
         ty_work=ay_at(tt_work)

         do k=-km,km
            alu(k,:,:) = transpose(at_Dy_at(at_Dy_at(tt_work))                                    - (2*pi*k/xl)**2*tt_work)
            alu(k,lm-1,:) = ty_work(:,0)
            alu(k,lm,:)   = ty_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

    function et_Lapla_et(et)   ! スペクトルに作用する Laplacian 演算子
      real(8), dimension(-km:km,0:lm)                :: et_Lapla_et
      real(8), dimension(-km:km,0:lm), intent(in)    :: 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

    function et_Vor2Strm_et(et,values,rigid) ! 渦度から流線を求める. 
                                             ! デフォルトは粘着条件
      ! Chebyshev-tau 法による計算
      ! 渦度 \zeta を与えて流線 \psi を求める.
      !    \nabla^2 \psi = \zeta, 
      !    \psi = const. at boundaries.
      ! 粘着条件
      !    \DP{\psi}{y} = 0 at boundaries
      ! 応力なし条件
      !    \DP[2]{\psi}{y} = 0 at boundaries
      !
      ! l=0,1,lm-1,lm 成分の式の代わりに境界条件を与える. 
      ! 渦度の低次成分を無視することは
      ! \nabla^4 \psi = \zeta^2 を解いていることに相当. 
      ! 4 階の微分方程式にすることで境界条件の数とつじつまがあう. 

      real(8), dimension(-km:km,0:lm),intent(in)  :: et
      real(8), dimension(-km:km,0:lm)             :: et_Vor2Strm_et

      ! 流線境界値. 境界で一定なので波数 0 成分のみ
      real(8), dimension(2), intent(in), optional :: values
      ! 境界条件スイッチ
      logical, dimension(2), intent(in), optional :: rigid

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

      real(8), dimension(-km:km,0:lm)         :: et_work
      real(8), dimension(0:lm,0:lm)           :: tt_work
      real(8), dimension(0:lm,0:jm)           :: ty_work
      real(8)                                 :: value1, value2   ! 境界値
      logical                                 :: rigid1, rigid2   ! 境界条件

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

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

      if (.not. present(rigid)) then
         rigid1=.true. 
 rigid2=.true.
      else
         rigid1 = rigid(1) 
 rigid2 = rigid(2)
      endif

      if ( first ) then
         first = .false.

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

         tt_work=0
         do l=0,lm
            tt_work(l,l)=1
         enddo
         do k=-km,km
            alu(k,:,:) = transpose(at_Dy_at(at_Dy_at(tt_work))                                    - (2*pi*k/xl)**2*tt_work)
         enddo

         ! 運動学的条件. 流線は境界で一定
         ty_work=ay_at(tt_work)
         do k=-km,km
            alu(k,lm-1,:) = ty_work(:,0)
            alu(k,lm,:)   = ty_work(:,jm)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid1 ) then
            ty_work=ay_at(at_Dy_at(tt_work))
         else
            ty_work=ay_at(at_Dy_at(at_Dy_at(tt_work)))
         endif
         do k=-km,km
            alu(k,0,:) = ty_work(:,0)
         enddo

         ! 力学的条件粘着条件 
         if ( rigid2 ) then
            ty_work=ay_at(at_Dy_at(tt_work))
         else
            ty_work=ay_at(at_Dy_at(at_Dy_at(tt_work)))
         endif
         do k=-km,km
            alu(k,1,:) = ty_work(:,jm)
         enddo

         call ludecomp(alu,kp)
      endif

      et_work = et
      et_work(:,0) = 0               ! 力学的条件
      et_work(:,1) = 0               ! 力学的条件

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

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

      et_Vor2Strm_et = lusolve(alu,kp,et_work)

    end function et_Vor2Strm_et

    subroutine et_initial(i,j,k,l,xmin,xmax,ymin,ymax)

      integer,intent(in) :: i, j           ! 格子点の設定(X,Y)
      integer,intent(in) :: k, l           ! 切断波数の設定(X,Y)

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

      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)

    end subroutine et_initial

    function yx_et(et)  ! スペクトル -> 台形格子
      real(8), dimension(0:jm,0:im-1)              :: yx_et
      real(8), dimension(-km:km,0:lm), intent(in)  :: et

      yx_et = ax_ae(transpose(ay_at(et)))

    end function yx_et