!--
!----------------------------------------------------------------------
!     Copyright 2009-2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ee_module_fftj
!
!      spml/ee_module モジュールは周期境界条件の下での 2 次元矩形領域の
!      流体運動をスペクトル法により数値計算するための Fortran90 関数を
!      提供する. 
!
!      内部で ISPACK/PZPACK(fftj) の Fortran77 サブルーチンを呼んでいる. 
!      スペクトルデータおよび格子点データの格納方法については
!      ISPACK/PZPACK のマニュアルを参照されたい. 
!
!履歴  2009/09/09  竹広真一  ee_module.f90 を pzpack(fftj) 用に改造
!      2016/06/23  竹広真一  ee_JacobianZ_ee の符号を修正
!
!++
module ee_module_fftj
  !
  != ee_module_fftj
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/ee_module モジュールは周期境界条件の下での 2 次元矩形領域の
  ! 流体運動をスペクトル法により数値計算するための Fortran90 関数を
  ! 提供する. 
  !
  ! 内部で ISPACK/PZPACK(fftj) の Fortran77 サブルーチンを呼んでいる. 
  ! スペクトルデータおよび格子点データの格納方法については
  ! ISPACK/PZPACK のマニュアルを参照されたい. 
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !  
  ! * 関数名の先頭 (ee_, yx_, x_, y_) は, 返す値の形を示している.
  !   ee_ :: スペクトルデータ(第 1,2 次元がそれぞれ Y,X 方向波数)
  !   yx_ :: 2 次元格子点データ(第 1,2 次元がそれぞれ Y,X 方向の格子点)
  !   x_  :: X 方向 1 次元格子点データ, y_ : Y 方向 1 次元格子点データ
  !
  ! * 関数名の間の文字列(Dx, Dy, Lapla, LaplaInv, Jacobian)は, 
  !   その関数の作用を表している.
  !
  ! * 関数名の最後 (_ee_ee, _ee, _yx, _x, _y) は, 入力変数の形が
  !   スペクトルデータおよび格子点データであることを示している.
  !   _ee    :: スペクトルデータ
  !   _ee_ee :: 2 つのスペクトルデータ
  !   _yx    :: 2 次元格子点データ
  !   _x     :: X 方向 1 次元格子点データ
  !   _y     :: Y 方向 1 次元格子点データ.
  !  
  !=== 各データの種類の説明
  !  
  ! * yx : 2 次元格子点データ.
  !   * 変数の種類と次元は real(8), dimension(0:jm-1,0:im-1). 
  !   * im, jm はそれぞれ X, Y 座標の格子点数であり, サブルーチン 
  !     ee_initial にてあらかじめ設定しておく.
  !   * 第 1 次元が Y 座標の格子点位置番号, 第 2 次元が X 座標の
  !     格子点位置番号である (X, Y の順ではない)ことに注意.
  !
  ! * ee : スペクトルデータ.
  !   * 変数の種類と次元は real(8), dimension(-lm:lm,-km:km). 
  !   * km, lm はそれぞれ X, Y 方向の最大波数であり, サブルーチン 
  !     ee_initial にてあらかじめ設定しておく.
  !     (X, Y 方向波数の順ではない)ことに注意. 
  !   * スペクトルデータの格納のされ方については...
  !
  ! * x, y : X, Y 方向 1 次元格子点データ.
  !   * 変数の種類と次元はそれぞれ
  !     real(8), dimension(0:im-1) および real(8), dimension(0:jm-1).
  !
  ! * ee_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * yx_ で始まる関数が返す値は 2 次元格子点データに同じ.
  !
  ! * x_, y_ で始まる関数が返す値は 1 次元格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化 
  !
  ! ee_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  ! ee_ChangeResolutionDomain :: 解像度, 領域設定の変更
  !
  !==== 座標変数
  !
  ! x_X, y_Y     ::  格子点座標(X,Y座標)を格納した 1 次元配列
  ! x_X_Weight, y_Y_Weight ::  重み座標を格納した 1 次元配列
  ! yx_X, yx_Y   :: 格子点データの XY 座標(X,Y)(格子点データ型 2 次元配列)
  !
  !==== 基本変換
  !
  ! yx_ee :: スペクトルデータから格子データへの変換
  ! ee_yx :: 格子データからスペクトルデータへの変換
  !
  !==== 微分
  !
  ! ee_Lapla_ee       :: スペクトルデータにラプラシアンを作用させる
  ! ee_LaplaInv_ee    :: スペクトルデータにラプラシアンの逆変換を作用させる
  ! ee_Dx_ee          :: スペクトルデータに X 微分を作用させる
  ! ee_Dy_ee          :: スペクトルデータに Y 微分を作用させる
  ! ee_Jacobian_ee_ee :: 2 つのスペクトルデータからヤコビアンを計算する
  !
  !==== 積分・平均
  !
  ! 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 方向積分および平均
  !
  !==== スペクトル解析
  !
  ! ee_EnergyFromStreamfunc_ee    :: 
  ! 流線関数からエネルギースペクトルを計算する
  !
  ! ee_EnstrophyFromStreamfunc_ee :: 
  ! 流線関数からエンストロフィースペクトルを計算する
  !
  !==== 補間計算
  !
  ! Interpolate_ee       :: 任意の点の値をスペクトルデータから計算する
  !
  use dc_message, only : MessageNotify
  implicit none

  private
  public ee_Initial                                       ! 初期化ルーチン
  public ee_ChangeResolutionDomain                        ! 解像度領域設定変更
  public yx_ee, ee_yx                                     ! 基本変換
  public ee_Lapla_ee, ee_LaplaInv_ee, ee_Dx_ee, ee_Dy_ee  ! 微分
  public ee_JacobianZ_ee, ee_Jacobian_ee_ee               ! 非線形計算
  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 ee_EnergyFromStreamfunc_ee                       ! エネルギー
  public ee_EnstrophyFromStreamfunc_ee                    ! エンストロフィー
  public Interpolate_ee                                   ! 補間
  public x_X, y_Y, x_X_Weight, y_Y_Weight, yx_X, yx_Y     ! 座標変数

  integer, parameter :: ind2imax=11              ! Max(im) = 2048 = 2^{11}
  integer, parameter :: ind2jmax=11              ! Max(jm) = 2048 = 2^{11}

  integer   :: im=32, jm=32                      ! 格子点の設定(X,Y)
  integer   :: km=10, lm=10                      ! 切断波数の設定(X,Y)
  real(8)   :: xl=1.0, yl=1.0                    ! 領域の大きさ

  integer, dimension(:),   pointer :: itj => null()
  real(8), dimension(:),   pointer :: tj => null()
  integer, dimension(:),   pointer :: iti => null()
  real(8), dimension(:),   pointer :: ti => null()

  real(8), dimension(:),   pointer :: x_X => null()
                                         ! 格子点座標(X)
  real(8), dimension(:),   pointer :: y_Y => null()
                                         ! 格子点座標(Y)

  real(8), dimension(:),   pointer :: x_X_Weight => null()
                                         ! 格子点重み(X)
                                         ! X 方向の格子点の間隔が格納してある.
  real(8), dimension(:),   pointer :: y_Y_Weight => null()
                                         ! 格子点重み(Y)
                                         ! Y 方向の格子点の間隔が格納してある.

  real(8), dimension(:,:), pointer :: yx_X => null()
                          ! 格子点(X)座標(2 次元)
                          ! 各格子点(i,j)の位置の X 座標を格納した格子データ
  real(8), dimension(:,:), pointer :: yx_Y => null()
                          ! 格子点(Y)座標(2 次元)
                          ! 各格子点(i,j)の位置の Y 座標を格納した格子データ


  integer, parameter :: nparams_max = 10 ! ee_Initial を呼べる最大回数
  type ee_param                          ! 解像度領域情報構造体
     integer   :: im, jm
     integer   :: km, lm
     real(8)   :: xl, yl
     integer, dimension(:),   pointer :: itj
     real(8), dimension(:),   pointer :: tj
     integer, dimension(:),   pointer :: iti
     real(8), dimension(:),   pointer :: ti
     real(8), dimension(:),   pointer :: x_X
     real(8), dimension(:),   pointer :: y_Y
     real(8), dimension(:),   pointer :: x_X_Weight
     real(8), dimension(:),   pointer :: y_Y_Weight
     real(8), dimension(:,:), pointer :: yx_X
     real(8), dimension(:,:), pointer :: yx_Y 
  end type ee_param
  type(ee_param) :: params(nparams_max)  ! 解像度領域情報
  integer :: nparams                     ! 解像度領域情報の個数

  real(8), parameter  :: pi=3.1415926535897932385D0

  save im, jm, km, lm, itj, tj, iti, ti, xl, yl
  save x_X, y_Y, x_X_Weight, y_Y_Weight, yx_X, yx_Y
  save params, nparams

  contains
  !--------------- 初期化 -----------------
    subroutine ee_Initial(i,j,k,l,xmin,xmax,ymin,ymax,id)
      !
      ! スペクトル変換の格子点数, 波数, 領域の大きさを設定する.
      !
      ! 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
      ! 初期設定をしなければならない.
      !
      ! オプショナル引数 id を用いて異なる解像度, 領域を同時に
      ! 扱うことができる. ee_Initial を解像度, 領域ごとに呼んで
      ! 引数 id をキープし, ee_ChangeResolutionDomain で切替える. 
      !
      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 座標範囲

      integer, intent(out), optional :: id  ! 解像度領域情報番号

      character(len=3) cid
      integer :: ii, jj
      integer :: ind2i, ind2j

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

      if ( nparams .ge. nparams_max ) then
         call MessageNotify('W','ee_initial',&
              'too many call of ee_Initial, nothing was done.')
         if ( present(id) ) id = -1
         return
      end if

      ind2i = 1
      do while ( im /= 2**ind2i .AND. ind2i <= ind2imax ) 
         ind2i=ind2i+1
      enddo
      if ( ind2i > ind2imax ) then
         call MessageNotify('E','ee_Initial', &
              'Number of X grid points IM should be 2^N and <= 2048')
      else if ( km <= 0 ) then
         call MessageNotify('E','ee_Initial', &
              'Number of waves kM should be positive')
      elseif ( km >= im/2 ) then
         call MessageNotify('E','ee_Initial', &
              'KM shoud be less than IM/2')
      endif

      ind2j = 1
      do while ( jm /= 2**ind2j .AND. ind2j <= ind2jmax ) 
         ind2j=ind2j+1
      enddo
      if ( ind2j > ind2jmax ) then
         call MessageNotify('E','ee_Initial', &
              'Number of Y grid points JM should be 2^N and <= 2048')
      else if ( lm <= 0 ) then
         call MessageNotify('E','ee_Initial', &
              'Number of waves LM should be positive')
      elseif ( lm >= jm/2 ) then
         call MessageNotify('E','ee_Initial', &
              'LM shoud be less than JM/2')
      endif


      nparams = nparams + 1

      params(nparams)%im = im
      params(nparams)%jm = jm
      params(nparams)%km = km
      params(nparams)%lm = lm
      params(nparams)%xl = xl
      params(nparams)%yl = yl

      allocate(params(nparams)%itj(4))
      allocate(params(nparams)%iti(4))
      allocate(params(nparams)%tj(jm*6))
      allocate(params(nparams)%ti(im*8))

      allocate(params(nparams)%x_X(0:im-1))
      allocate(params(nparams)%x_X_Weight(0:im-1))
      allocate(params(nparams)%y_Y(0:jm-1))
      allocate(params(nparams)%y_Y_Weight(0:jm-1))
      allocate(params(nparams)%yx_X(0:jm-1,0:im-1))
      allocate(params(nparams)%yx_Y(0:jm-1,0:im-1))

      call ee_ChangeResolutionDomain(nparams)

      call pzinit(jm,im,itj,tj,iti,ti)

      do ii=0,im-1
         x_X(ii) = xmin + xl/im*ii
      enddo
      x_X_Weight = xl/im

      do jj=0,jm-1
         y_Y(jj) = ymin + yl/jm*jj
      enddo
      y_Y_Weight = yl/jm

      yx_X = spread(x_X,1,jm)
      yx_Y = spread(y_Y,2,im)

      if ( present(id) ) id = nparams

      write(cid,'(I3)') nparams
      call MessageNotify('M','ee_initial','ee_module_fftj (2016/06/23) is initialized')
      call MessageNotify('M','ee_initial',&
           'ResolutionDomain ID is '//trim(adjustl(cid)))
    end subroutine ee_initial

  !--------------- id 変更 -----------------
    subroutine ee_ChangeResolutionDomain(id)
      !
      ! 解像度, 領域設定の変更. ee_Initial で設定する際に
      ! かえってくるオプショナル引数 id の値を用いる. 
      ! 
      integer, intent(in) :: id

      if (id .gt. nparams .or. id .lt. 1) then
         write(*,*)"id is invalid"
      end if

      im = params(id)%im
      jm = params(id)%jm
      km = params(id)%km
      lm = params(id)%lm
      xl = params(id)%xl
      yl = params(id)%yl
      itj => params(id)%itj
      tj => params(id)%tj
      iti => params(id)%iti
      ti => params(id)%ti
      x_X => params(id)%x_X
      y_Y => params(id)%y_Y
      x_X_Weight => params(id)%x_X_Weight
      y_Y_Weight => params(id)%y_Y_Weight
      yx_X => params(id)%yx_X
      yx_Y => params(id)%yx_Y

    end subroutine ee_ChangeResolutionDomain

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

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

      real(8), dimension(jm*im)                     :: w
      ! 作業用配列

      call pzs2ga(lm,km,jm,im,ee,yx_ee,w,itj,tj,iti,ti)
    end function yx_ee

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

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

      real(8), dimension(jm*im)                     :: w
      real(8), dimension(0:jm-1,0:im-1)             :: yx_tmp
      ! 作業用配列

      yx_tmp = yx
      call pzg2sa(lm,km,jm,im,yx_tmp,ee_yx,w,itj,tj,iti,ti)

    end function ee_yx

  !--------------- 微分計算 -----------------
    function ee_Lapla_ee(ee)
      !
      ! 入力スペクトルデータにラプラシアン(∂xx+∂yy)を作用する.
      !
      ! スペクトルデータのラプラシアンとは, 対応する格子点データに
      ! ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      ! 実際にはスペクトルデータに全波数 (k**2 + l**2) をかける
      ! 計算を行っている. 
      !
      real(8), dimension(-lm:lm,-km:km)              :: ee_Lapla_ee
      !(out) スペクトルデータのラプラシアン

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

      integer k,l
      ! 作業変数

      do k=-km,km
         do l=-lm,lm
            ee_Lapla_ee(l,k) = -((2*pi*k/xl)**2+(2*pi*l/yl)**2)*ee(l,k)
         enddo
      enddo
    end function ee_Lapla_ee

    function ee_LaplaInv_ee(ee)
      !
      ! 入力スペクトルデータに逆ラプラシアン(∂xx+∂yy)**(-1)を作用する.
      !
      ! スペクトルデータの逆ラプラシアンとは, 対応する格子点データに
      ! 逆ラプラシアンを作用させたデータのスペクトル変換のことである.
      !
      ! 実際にはスペクトルデータに全波数 (k**2 + l**2) で割る
      ! 計算を行っている. k=l=0 成分には 0 を代入している. 
      !
      real(8), dimension(-lm:lm,-km:km)             :: ee_LaplaInv_ee
      !(out) スペクトルデータの逆ラプラシアン

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

      integer k,l

      do k=-km,km
         do l=-lm,lm
            if ( k.ne.0 .or. l.ne.0 ) then
               ee_LaplaInv_ee(l,k) = -ee(l,k)/((2*pi*k/xl)**2+(2*pi*l/yl)**2)
            else
               ee_LaplaInv_ee(l,k) = 0.0
            endif
         enddo
      enddo
    end function ee_LaplaInv_ee

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

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

      integer k,l
      ! 作業変数

      do k=-km,km
         do l=-lm,lm
            ee_Dx_ee(l,k) = -(2*pi*k/xl)*ee(-l,-k)
         enddo
      enddo
    end function ee_Dx_ee

    function ee_Dy_ee(ee)
      !
      ! 入力スペクトルデータに Y 微分(∂y)を作用する.
      !
      ! スペクトルデータの X 微分とは, 対応する格子点データに Y 微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      ! 実際にはスペクトルデータに X 方向波数 l をかけて
      ! sin(ky) <-> cos(ky) 成分に入れ換える計算を行っている.
      !
      real(8), dimension(-lm:lm,-km:km)              :: ee_Dy_ee
      !(out) スペクトルデータの Y 微分

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

      integer k,l
      ! 作業変数

      do k=-km,km
         do l=-lm,lm
            ee_Dy_ee(l,k) = -(2*pi*l/yl)*ee(-l,-k)
         enddo
      enddo

    end function ee_Dy_ee

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

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

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

      real(8), dimension((2*km+1)*(2*lm+1))          :: ws
      real(8), dimension(jm*im*3)                    :: wg
      ! 作業領域

      call pzajcb(lm,km,jm,im,ee_a,ee_b,ee_Jacobian_ee_ee,ws,wg,itj,tj,iti,ti)

      ee_Jacobian_ee_ee  = (2*pi/xl)*(2*pi/yl) * ee_Jacobian_ee_ee

    end function ee_Jacobian_ee_ee

    function ee_JacobianZ_ee(ee_zeta)
      !
      ! 渦度スペクトルデータ ζ から流線関数と渦度のヤコビアン
      !
      !     J(ψ,ζ)=(∂xψ)(∂yζ)-(∂yψ)(∂xζ)
      !
      !  を計算する. ただしψ は (∂xx+∂yy)ψ=ζ を満たす流線関数である.
      !
      real(8), dimension(-lm:lm,-km:km)              :: ee_JacobianZ_ee
      !(out) 流線関数と渦度のヤコビアン

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

      real(8), dimension((2*km+1)*(2*lm+1))          :: ws
      real(8), dimension(jm*im*3)                    :: wg
      ! 作業領域

      call pzajbs(lm,km,jm,im,yl/xl,ee_Zeta,ee_JacobianZ_ee,ws,wg,itj,tj,iti,ti)

      ee_JacobianZ_ee = - (2*pi/xl)*(2*pi/yl)/(2*pi/yl)**2 * ee_JacobianZ_ee

    end function ee_JacobianZ_ee

  !--------------- 積分計算 -----------------
    function IntYX_yx(yx)
      !
      ! 2 次元格子点データの全領域積分および平均.
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:jm-1,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-1
            IntYX_yx = IntYX_yx + yx(j,i) * y_Y_Weight(j) * x_X_Weight(i)
         enddo
      enddo
    end function IntYX_yx

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

      real(8), dimension(0:jm-1)          :: 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 x_IntY_yx(yx)
      !
      ! 2 次元格子点データの Y 方向積分
      !
      ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している. 
      !
      real(8), dimension(0:jm-1,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-1
         x_IntY_yx(:) = x_IntY_yx(:) + yx(j,:) * y_Y_Weight(j)
      enddo
    end function x_IntY_yx

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

      IntY_y = sum(y*y_Y_Weight)
    end function IntY_y

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

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

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

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

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

  !--------------- スペクトル計算 -----------------
    function ee_EnergyFromStreamfunc_ee(ee_StrFunc)
      !
      ! 流線関数からエネルギースペクトルを計算する. 
      !
      !   E_kl = (1/2)(k^2+l^2)|\psi_kl|^2
      !
      ! * E_kl の総和が場の平均運動エネルギーとなる. 
      ! * それに領域の面積をかけると全運動エネルギーとなる. 
      !
      real(8), dimension(-lm:lm,-km:km)    :: ee_EnergyFromStreamfunc_ee
      ! エネルギースペクトル

      real(8), dimension(-lm:lm,-km:km), intent(in) :: ee_StrFunc
      ! 流線関数

      integer k,l
      ! 作業変数

      do k=-km,km
         do l=-lm,lm
            ee_EnergyFromStreamfunc_ee(l,k) &
                 = 0.5 * ( (2*pi*k/xl)**2 + (2*pi*l/yl)**2 ) & 
                          * (ee_StrFunc(l,k)**2 + ee_StrFunc(-l,-k)**2)
         enddo
      enddo

    end function ee_EnergyFromStreamfunc_ee

    function ee_EnstrophyFromStreamfunc_ee(ee_StrFunc)
      !
      ! 流線関数からエンストロフィースペクトルを計算する. 
      !
      !   Q_kl = (1/2)(k^2+l^2)^2|\psi_kl|^2
      !
      ! * Q_kl の総和が場の平均エンストロフィーとなる. 
      ! * それに領域の面積をかけると全エンストロフィーとなる. 
      !

      real(8), dimension(-lm:lm,-km:km)    :: ee_EnstrophyFromStreamfunc_ee
      ! エンストロフィーースペクトル

      real(8), dimension(-lm:lm,-km:km), intent(in) :: ee_StrFunc
      ! 流線関数

      integer k,l
      ! 作業変数

      do k=-km,km
         do l=-lm,lm
            ee_EnstrophyFromStreamfunc_ee(l,k) &
                 = 0.5 * ( (2*pi*k/xl)**2 + (2*pi*l/yl)**2 )**2 & 
                          * (ee_StrFunc(l,k)**2 + ee_StrFunc(-l,-k)**2)
         enddo
      enddo

    end function ee_EnstrophyFromStreamfunc_ee

  !--------------- 補間計算 -----------------
    function Interpolate_ee( ee_Data, x, y )
      real(8), intent(IN)  :: ee_data(-lm:lm,-km:km)  ! スペクトルデータ
      real(8), intent(IN)  :: x                       ! 補間する点の x 座標 
      real(8), intent(IN)  :: y                       ! 補間する点の y 座標 
      real(8)              :: Interpolate_ee          ! 補間した値

      integer :: k, l
      real(8) :: xx, yy

      xx =(2*PI/xl)*(x - x_X(0))
      yy =(2*PI/yl)*(y - y_Y(0))

      Interpolate_ee = ee_Data(0,0)

      ! l=0
      do k=1,km
         Interpolate_ee = Interpolate_ee &
              + 2*( ee_Data(0,k)*cos(k*xx) - ee_Data(0,-k)*sin(k*xx) )
      end do

      ! k=0
      do l=1,lm
         Interpolate_ee = Interpolate_ee &
              + 2*( ee_Data(l,0)*cos(l*yy) - ee_Data(-l,0)*sin(l*yy) )
      end do

      ! k*l > 0
      do l=1,lm
         do k=1,km
            Interpolate_ee = Interpolate_ee &
              + 2*(  ee_Data(l,k)*(   cos(k*xx)*cos(l*yy)   &
                                    - sin(k*xx)*sin(l*yy) ) &
                    -ee_Data(-l,-k)*(   sin(k*xx)*cos(l*yy)   &
                                      + cos(k*xx)*sin(l*yy) ) )
         end do
      end do

      ! k*l < 0
      do l=1,lm
         do k=1,km
            Interpolate_ee = Interpolate_ee &
              + 2*(  ee_Data(-l,k)*(   cos(k*xx)*cos(l*yy)   &
                                     + sin(k*xx)*sin(l*yy) ) &
                    -ee_Data(l,-k)*(   sin(k*xx)*cos(l*yy)   &
                                     - cos(k*xx)*sin(l*yy) ) )
         end do
      end do

    end function Interpolate_ee

end module ee_module_fftj
