eq_module

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

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

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

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

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