et_module

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

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

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

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

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

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

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

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

      IntX_x = sum(x*x_X_Weight)
    end function IntX_x

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

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

      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

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

      IntY_y = sum(y*y_Y_Weight)
    end function IntY_y

    function Interpolate_et(et_data,xval,yval)

      real(8), dimension(-km:km,0:lm), intent(IN) :: et_data
      !(in) 入力フーリエデータ

      real(8), intent(in)                         :: xval, yval
      ! 補間する点の X, Y 座標

      real(8)                                     :: Interpolate_et
      ! 補間した結果の値

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

    end function Interpolate_et

    subroutine et_Boundaries(et,values,cond)
      !
      ! ディリクレ, ノイマン条件の適用. チェビシェフ空間での計算
      !
      ! 実際には中で呼ばれている at_module のサブルーチン at_Boundaries_DD,
      ! at_Boundaries_DN, at_Boundaries_ND, at_Boundaries_NN を用いている.
      ! これらを直接呼ぶことも出来る.
      !
      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
              ! 境界条件. 省略時は 'RR'
              !   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 微分(∂x)を作用する.
      !
      ! スペクトルデータの X 微分とは, 対応する格子点データに X 微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      ! 実際にはスペクトルデータに X 方向波数 k をかけて
      ! sin(kx) <-> cos(kx) 成分に入れ換える計算を行っている.
      !
      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 微分(∂y)を作用する.
      !
      ! スペクトルデータの X 微分とは, 対応する格子点データに Y 微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(-km:km,0:lm)               :: et_Dy_et
      !(out) スペクトルデータの Y 微分

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

      et_Dy_et = at_Dy_at(et)

    end function et_Dy_et

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

      call MessageNotify('M','et_initial','et_module (2013/08/20) is initialized')
    end subroutine et_initial

    function et_Jacobian_et_et(et_a,et_b)
      !
      !  2 つのスペクトルデータからヤコビアン
      !
      !     J(A,B)=(∂xA)(∂yB)-(∂yA)(∂xB)
      !
      !  を計算する.
      !
      !  2 つのスペクトルデータのヤコビアンとは, 対応する 2 つの
      !  格子点データのヤコビアンのスペクトル変換のことである.
      !
      real(8), dimension(-km:km,0:lm)                :: et_Jacobian_et_et
      !(out) 2 つのスペクトルデータのヤコビアン

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

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

      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)
      !
      ! 境界で一様な値を与える条件(ディリクレ条件)下で,
      ! 入力スペクトルデータに逆ラプラシアン(∂xx+∂yy)**(-1)を作用する.
      !
      ! Chebyshev-tau 法による計算
      !
      ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに
      ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      real(8), dimension(-km:km,0:lm),intent(in)  :: et
      !(in) スペクトルデータ

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

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

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

      real(8), dimension(-km:km,0:lm)         :: et_work
      real(8), dimension(-km:km,0:jm)         :: ey_work
      real(8), 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

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

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

      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,cond,new)    ! 渦度から流線を求める.
      !
      ! 渦度から流線を求める.
      !
      ! Chebyshev-tau 法 + Crank Nicolson 法による計算
      ! 渦度 \zeta を与えて流線 \psi を求める.
      !    \nabla^2\psi = \zeta, 
      !    \psi = const. at boundaries.
      ! 粘着条件
      !    \DP{\psi}{y} = 0 at boundaries
      ! 応力なし条件
      !    \DP[2]{\psi}{y} = 0 at boundaries
      !
      real(8), dimension(-km:km,0:lm),intent(in)  :: et
              !(in) 入力渦度分布

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

      real(8), dimension(2), intent(in), optional :: values
              !(in) 流線境界値. 境界で一定なので波数 0 成分のみ
              !     省略時は 0.

      character(len=2), intent(in), optional  :: cond
              ! (in) 境界条件スイッチ. 省略時は 'RR'
              !      RR : 両端粘着条件
              !      RF : 上端粘着下端応力なし条件
              !      FR : 上端応力なし下端粘着条件
              !      FF : 両端応力なし条件

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

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

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

      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
      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: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) = value1*2     ! 運動学的条件. 波数 0 は重み 1/2

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

      et_Vor2Strm_et = lusolve(alu,kp,et_work)

    end function et_Vor2Strm_et

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

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

      et_ey = at_ay(ey)

    end function et_ey

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

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

      real(8), dimension(-km:km,0:jm)              :: ey
      real(8), 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

    function ey_Vor2Strm_ey(ey,values,cond,new)    ! 渦度から流線を求める.
      !
      ! 渦度から流線を求める.
      ! Y 方向格子点空間での Chebyshev-Collocation 法による計算
      !
      ! 渦度 \zeta を与えて流線 \psi を求める.
      !    \nabla^2 \psi = \zeta,
      !    \psi = const. at boundaries.
      ! 粘着条件
      !    \DP{\psi}{y} = 0 at boundaries
      ! 応力なし条件
      !    \DP[2]{\psi}{y} = 0 at boundaries
      !

      real(8), dimension(-km:km,0:jm),intent(in)  :: ey
              !(in) 入力渦度分布

      real(8), dimension(-km:km,0:jm)             :: ey_Vor2Strm_ey
              !(out) 出力流線分布

      real(8), dimension(2), intent(in), optional :: values
              !(in) 流線境界値. 境界で一定なので波数 0 成分のみ
              !     省略時は 0.

      character(len=2), intent(in), optional  :: cond
              ! (in) 境界条件スイッチ. 省略時は 'RR'
              !      RR : 両端粘着条件
              !      RF : 上端粘着下端応力なし条件
              !      FR : 上端応力なし下端粘着条件
              !      FF : 両端応力なし条件

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

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

      real(8), dimension(-km:km,0:jm)         :: ey_work
      real(8), dimension(-km:km,0:jm)         :: ey_I
      real(8)                                 :: 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 
             ey_I(:,j) = 1.0D0
            alu(:,:,j) = ey_et(et_Lapla_et(et_ey(ey_I)))
         enddo

         do j=0,jm
            ey_I = 0
            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

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

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

      ey_et = ay_at(et)

    end function ey_et

    function ey_yx(yx)
      !
      ! X 座標格子データからスペクトルデータへ変換する.
      !
      real(8), dimension(-km:km,0:jm)              :: ey_yx
      !(out) スペクトルデータ

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

      real(8), 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

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

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

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

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

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

      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

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

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

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

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

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

      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

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

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

      real(8), dimension(-km:km,0:jm)              :: ey
      real(8), 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 yx_ey(ey)
      !
      ! X 座標スペクトルデータから格子データへ変換する.
      !
      real(8), dimension(0:jm,0:im-1)              :: yx_ey
      !(out) 格子点データ

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

      real(8), 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