!--
!----------------------------------------------------------------------
! Copyright (c) 2016 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ee_mpi_module
!
!      spml/ee_module モジュールは周期境界条件の下での 2 次元矩形領域の
!      流体運動をスペクトル法により数値計算するための Fortran90 関数を
!      提供する. 
!
!      内部で ISPACK/FTPACKI の Fortran77 サブルーチンを呼んでいる.
!      スペクトルデータおよび格子点データの格納方法については...
!
!履歴  2016/08/08  竹広真一  eee_mpi_module より改造
!      2016/09/22  竹広真一  ef_Lapla_ef, ef_Laplainv_ef 初期化
!
!++
module ee_mpi_module
  !
  != ee_mpi_module
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/ee_module モジュールは周期境界条件の下での 2 次元矩形領域の
  ! 流体運動をスペクトル法により数値計算するための Fortran90 関数を
  ! 提供する. 
  !
  ! 内部で ISPACK/FTPACKI の Fortran77 サブルーチンを呼んでいる.
  ! スペクトルデータおよび格子点データの格納方法については...
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (ef_, vx_, x_, v_) は, 返す値の形を示している.
  !   ef_ :: スペクトルデータ
  !           (第 1,2 次元がそれぞれ Y,X 方向波数(X について分割配置))
  !   vx_  :: XY 方向 2 次元格子点データ(Y について分散配置
  !   x_   :: X 方向 1 次元格子点データ
  !   v_   :: Y 方向 1 次元格子点データ(分散配置)
  !
  ! * 関数名の間の文字列(Dx, Dy, Lapla, LaplaInv)は,
  !   その関数の作用を表している.
  !
  ! * 関数名の最後 (_ef_ef,_ef,_vx, _x, _v) は, 入力変数の形が
  !   スペクトルデータおよび格子点データであることを示している.
  !    _ef     :: スペクトルデータ
  !    _ef_ef  :: 2 つのスペクトルデータ
  !    _vx     :: 2 次元格子点データ
  !    _x      :: X 方向 1 次元格子点データ
  !    _v      :: Y 方向 1 次元格子点データ
  !
  !=== 各データの種類の説明
  !
  ! * vx : 2 次元格子点データ.
  !   * 変数の種類と次元は real(8), dimension(,js(ip):je(ip),0:im-1).
  !   * im は X 座標の格子点数であり, サブルーチン ee_mpi_initial にて
  !     Y 方向の全格子点数 jm とともにあらかじめ設定しておく.
  !   * ee_mpi_Initial によってプロセス毎に Y 方向の
  !     分割配置が設定され, その情報が public 変数 js,je,jc に設定される. 
  !   * js が分割配置の最小格子点, je が最大格子点, jc が格子点数である. 
  !
  ! * ef : スペクトルデータ.
  !   * 変数の種類と次元は real(8), dimension(-lm:lm,2*kc). 
  !   * mm は Y 方向の最大波数であり, サブルーチン ee_mpi_initial にて
  !     X 方向の最大波数 km とともにあらかじめ設定しておく. 
  !     (X, Y 方向波数の順ではない)ことに注意. 
  !   * ee_mpi_Initial によってプロセス毎に X 方向波数の
  !     分割配置が設定され, その情報が public 変数 ls,le,lc に設定される. 
  !   * ks が分割配置の最小波数, ke が最大波数, kc が波数の数である. 
  !     実際の格納のされ方は ks, ks+1,...,ke,-ke,-ke+1,...,-ks の順である.
  !   * 分割配置される前のスペクトルデータの格納のされかたは IPSACK の
  !     P2PACK に基づく. 
  !       ee_S( L, K) : s_{kl} の実部
  !       ee_S(-L,-K) : s_{kl} の虚部
  !       ee_S(-L, K) : s_{k(-l)} の実部
  !       ee_S( L,-K) : s_{k(-l)} の虚部
  !       ee_S( L, 0) : s_{0l} の実部
  !       ee_S(-L, 0) : s_{0l} の虚部
  !       ee_S( 0, K) : s_{k0} の実部
  !       ee_S( 0,-K) : s_{k0} の虚部
  !       ee_S( 0, 0) : s_{00} (実数)
  !   * 元の格子点データが実数であるため, s_{-k-l} = s^*_{kl} なる関係がある.
  !     ただし ^* は複素共役を表す. 
  !
  ! * x, v : X, Y 方向 1 次元格子点データ.
  !   * 変数の種類と次元はそれぞれ real(8), dimension(0:im-1),
  !     real(8), dimension(js(ip):je(ip))
  !
  ! * ef_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * vx_ で始まる関数が返す値は 3 次元格子点データに同じ.
  !
  ! * x_, v_ で始まる関数が返す値は 1 次元格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! ee_mpi_Initial :: スペクトル変換の格子点数, 切断波数の設定
  ! kf_k           :: X 方向波数の位置情報
  !
  !==== 座標変数
  !
  ! x_X, v_Y               :: 格子点座標(X,Y座標)を格納した 1 次元配列
  ! x_X_Weight, v_Y_Weight :: 重み座標を格納した 1 次元配列
  ! vx_X, vx_Y_Z           :: 格子点データの XY 座標(X,Y)
  !                           (格子点データ型 3 次元配列)
  !
  !==== 基本変換
  !
  ! vx_ef :: スペクトルデータから格子データへの変換
  ! ef_vx :: 格子データからスペクトルデータへの変換
  ! xy_vx :: 分散格子点データの集積
  ! ee_ef :: 分散スペクトルデータの集積
  ! ef_ee :: 集積スペクトルデータの分散
  !
  !==== 微分
  !
  ! ef_Lapla_ef       :: スペクトルデータにラプラシアンを作用させる
  ! ef_LaplaInv_ef    :: スペクトルデータにラプラシアンの逆変換を作用させる
  ! ef_Dx_ef          :: スペクトルデータに X 微分を作用させる
  ! ef_Dy_ef          :: スペクトルデータに Y 微分を作用させる
  !
  !==== 非線形計算
  !
  ! ef_Joacobian_ef_ef :: 2 つのスペクトルデータからヤコビアンを計算する.
  ! ef_JoacobianZ_ef   :: 渦度スペクトルデータから非発散方程式の非線形項を計算する
  !
  !==== 積分・平均
  !
  ! IntVX_vx,  AvrVX_vx  :: 2 次元(XY)格子点データの X,Y 方向積分および平均
  ! v_IntX_vx, v_AvrX_vx :: 2 次元(XY)格子点データの X 方向積分および平均
  ! x_IntV_vx, x_AvrV_vx :: 2 次元(XY)格子点データの Y 方向積分および平均
  !
  ! IntX_x, AvrX_x       :: 1 次元(X)格子点データの X 方向積分および平均
  ! IntV_v, AvrV_v       :: 1 次元(Y)格子点データの Y 方向積分および平均
  !
  !==== スペクトル解析
  !
  ! ef_EnergyFromStreamfunc_ef    :: 流れ関数からエネルギースペクトルを計算
  ! ef_EnstrophyFromStreamfunc_ef :: 流れ函数からエンストロフィースペクトルを計算
  !
  !==== 補間計算
  !
  ! Interpolate_ef       :: 任意の点の値をスペクトルデータから計算する
  !
  use dc_message, only : MessageNotify
  use mpi
  implicit none

  private
  public ks, ke, kc
  public js, je, jc
  public kf_k
  public ee_mpi_Initial                                 ! 初期化ルーチン
  public vx_ef, ef_vx                                   ! 基本変換
!!$  public xy_vx                                          ! 基本変換
!!$  public ee_ef, ef_ee                                   ! 基本変換
  public ef_Dx_ef, ef_Dy_ef                             ! 微分
  public ef_Lapla_ef, ef_LaplaInv_ef                    ! 微分
  public ef_Jacobian_ef_ef                              ! ヤコビアン
  public ef_JacobianZ_ef                                ! ヤコビアン

  public IntYX_vx,  AvrYX_vx                            ! 積分・平均
  public v_IntX_vx, v_AvrX_vx                           ! 積分・平均
  public x_IntY_vx, x_AvrY_vx                           ! 積分・平均

  public IntX_x, AvrX_x                                 ! 積分・平均
  public IntY_v, AvrY_v                                 ! 積分・平均

  public Interpolate_ef                                 ! 補間

  public x_X, v_Y                                       ! 座標変数
  public x_X_Weight, v_Y_Weight                         ! 座標変数
  public vx_X, vx_Y                                     ! 座標変数

  public ef_EnergyFromStreamfunc_ef                     ! エネルギー
  public ef_EnstrophyFromStreamfunc_ef                  ! エンストロフィー

  integer :: im=32, jm=32                               ! 格子点の設定(X,Y)
  integer :: km=10, lm=10                               ! 切断波数の設定(X,Y)

  integer, dimension(:), allocatable :: ke              ! X 最大波数
  integer, dimension(:), allocatable :: ks              ! X 最小波数
  integer, dimension(:), allocatable :: kc              ! X 波数の数
  integer, dimension(:), allocatable :: je              ! Y 座標最大格子点
  integer, dimension(:), allocatable :: js              ! Y 座標最小格子点
  integer, dimension(:), allocatable :: jc              ! Y 座標格子点数
  integer :: np                                         ! MPI プロセス数
  integer :: ip                                         ! MPI プロセス ID  

  integer, dimension(:), allocatable :: itj             ! FFT 用配列(Y)
  real(8), dimension(:), allocatable :: tj              ! FFT 用配列(Y)
  integer, dimension(:), allocatable :: iti             ! FFT 用配列(X)
  real(8), dimension(:), allocatable :: ti              ! FFT 用配列(X)

  real(8), dimension(:), allocatable :: x_X             ! 格子点座標(X)
  real(8), dimension(:), allocatable :: v_Y             ! 格子点座標(Y)

  real(8), dimension(:), allocatable :: x_X_Weight      ! 格子点重み(X)
  real(8), dimension(:), allocatable :: v_Y_Weight      ! 格子点重み(Y)

  real(8), dimension(:,:), allocatable :: vx_X          ! 格子点(X)座標(2 次元)
  real(8), dimension(:,:), allocatable :: vx_Y          ! 格子点(Y)座標(2 次元)

  real(8) :: xmin, xmax                                 ! X 座標範囲
  real(8) :: ymin, ymax                                 ! Y 座標範囲
  real(8) :: xl, yl                                     ! 領域の大きさ(X,Y)
 
  real(8), parameter  :: pi=3.1415926535897932385D0

  save im, jm, km, lm, itj, tj, iti, ti
  save ks, ke, kc, js, je, jc, np, ip
  save x_X, v_Y, x_X_Weight, v_Y_Weight, vx_X, vx_Y
  save xmin, xmax, ymin, ymax, xl, yl

contains
  !--------------- 初期化 -----------------
  subroutine ee_mpi_Initial(i_in,j_in,k_in,l_in,xmin_in,xmax_in,ymin_in,ymax_in)
    !
    ! スペクトル変換の格子点数, 波数を設定する.
    ! 領域の範囲は [0,2pi]x[0,2pi] に決め打ち
    !
    ! 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
    ! 初期設定をしなければならない.
    !
    integer,intent(in) :: i_in           ! 格子点数(X)
    integer,intent(in) :: j_in           ! 格子点数(Y)
    integer,intent(in) :: k_in           ! 切断波数(X)
    integer,intent(in) :: l_in           ! 切断波数(Y)

    real(8),intent(in) :: xmin_in, xmax_in     ! X 座標範囲
    real(8),intent(in) :: ymin_in, ymax_in     ! Y 座標範囲
 
    integer :: ii, jj
!!$    integer :: kks, kke, kkc, jjs, jje, jjc
    integer :: iip
    integer :: kp, jp
    integer :: kint, kmod
    integer :: ierr

    im = i_in         ; jm = j_in
    km = k_in         ; lm = l_in

    xmin = xmin_in    ; xmax = xmax_in
    ymin = ymin_in    ; ymax = ymax_in
    xl = xmax-xmin    ; yl = ymax-ymin

    CALL MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)
    CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)

!!$    kp=km/np+1
!!$    kks=kp*ip
!!$    kke=min(kp*(ip+1)-1,km)
!!$    if(kke.ge.kks) then
!!$       kkc=kke-kks+1
!!$    else
!!$       kkc=0 ;  kks=0  ; kke=0
!!$    end if
!!$
!!$    jp=(jm-1)/np+1
!!$    jjs=jp*ip
!!$    jje=min(jp*(ip+1)-1,jm-1)
!!$    if(jje.ge.jjs) then
!!$       jjc=jje-jjs+1
!!$    else
!!$       jjc=0 ; jjs=0 ; jje=0
!!$    end if

    allocate(ks(0:np-1))
    allocate(ke(0:np-1))
    allocate(kc(0:np-1))
    allocate(js(0:np-1))
    allocate(je(0:np-1))
    allocate(jc(0:np-1))

    allocate(itj(5))
    allocate(iti(5))
    allocate(tj(jm*2))
    allocate(ti(im*2))

    call fttrui(im,iti,ti)
    call fttzui(jm,itj,tj)

    kint = lm/np
    kmod = mod(lm,np)
    kc   = kint
    do iip=0,kmod-1
       kc(iip) = kc(iip) + 1
    end do

    ks(0)=0 ; ke(0) = ks(0) + kc(0) -1
    do iip=1,np-1
       ks(iip) = ks(iip-1) + kc(iip-1)
       ke(iip) = ks(iip) + kc(iip) - 1
    end do
       
    do iip=0,np-1
!!$       kp=lm/np+1
!!$       ks(iip)=kp*iip
!!$       ke(iip)=min(kp*(iip+1)-1,km)
!!$
!!$
!!$       if(ke(iip).ge.ks(iip)) then
!!$          kc(iip)=ke(iip)-ks(iip)+1
!!$       else
!!$          kc(iip)=0 ;  ks(iip)=0  ; ke(iip)=0
!!$       end if

       jp=(jm-1)/np+1
       js(iip)=jp*iip
       je(iip)=min(jp*(iip+1)-1,jm-1)
       if(je(iip).ge.js(iip)) then
          jc(iip)=je(iip)-js(iip)+1
       else
          jc(iip)=0 ; js(iip)=0 ; je(iip)=0
       end if
    enddo

    allocate(x_X(0:im-1))
    allocate(x_X_Weight(0:im-1))
    allocate(v_Y(jc(ip)))
    allocate(v_Y_Weight(jc(ip)))
    allocate(vx_X(jc(ip),0:im-1))
    allocate(vx_Y(jc(ip),0:im-1))

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

    do jj=1,jc(ip)
       v_Y(jj) = ymin + yl/jm*(js(ip)+jj-1)
    enddo
    v_Y_Weight = yl/jm

    vx_X = spread(x_X,1,jc(ip))
    vx_Y = spread(v_Y,2,im)

    call MessageNotify('M','ee_mpi_initial', &
                       'ee_mpi_module (2020/02/13)is initialized')

  end subroutine ee_mpi_Initial

  !--------------- 波数位置問い合わせ -----------------
  function kf_k(k)
    !
    ! スペクトルデータ eef(-lm,lm,2*kc) での
    ! X 方向波数 L の位置情報(第 2 添え字)を調べる.
    !
    ! X 方向波数については ls,ls+1,... le, -le,-le+1,...,-ls
    ! の順に格納されている.
    !
    integer, intent(IN)  ::  k        ! X 方向波数
    integer              ::  kf_k     ! スペクトルデータ第 2 添え字

    if ( abs(k) < ks(ip) .OR. ke(ip) < abs(k) ) then
       kf_k = -1
    else if ( k >= 0 ) then
       kf_k = k-ks(ip)+1
    else
       kf_k = ks(ip)+2*kc(ip)-abs(k)
    endif

  end function kf_k

  !--------------- 基本変換 -----------------
  function vx_ef(ef)
    !
    ! スペクトルデータから格子データへ変換する.
    !
    real(8), dimension(js(ip):je(ip),0:im-1)        :: vx_ef
                                                        !(out) 格子点データ
    real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef
                                                        !(in)  スペクトルデータ
    real(8), dimension(js(ip):je(ip),2,0:im/2-1) :: vax_work
    real(8), dimension(kc(ip),0:jm-1,2)          :: fya_work
 
    real(8), dimension(jc(ip)*im)       :: y2    
    real(8), dimension(kc(ip)*jm*2)     :: y3

    integer :: k, l

    fya_work = 0.0D0

    do k=1,kc(ip)
       do l=1,lm
          fya_work(k,l,1)    = ef(l,k)
          fya_work(k,l,2)    = ef(-l,2*kc(ip)-k+1)
          fya_work(k,jm-l,1) = ef(-l,k)
          fya_work(k,jm-l,2) = ef(l,2*kc(ip)-k+1)
       enddo
    enddo

    do k=1,kc(ip)
       fya_work(k,0,1) = ef(0,k)
       fya_work(k,0,2) = ef(0,2*kc(ip)-k+1)
    enddo

    if ( ip == 0 ) then
       do l=1,lm
          fya_work(1,l,1) = ef(l,1)
          fya_work(1,l,2) = ef(-l,1)
          fya_work(1,jm-l,1) = ef(l,1)
          fya_work(1,jm-l,2) = -ef(-l,1)
       end do
       fya_work(1,0,1) = ef(0,1)
    endif

    call fttzub(kc(ip),jm,fya_work,y3,itj,tj)

    vax_work = vax_fya(fya_work)

    call fttrub(jc(ip),im,vax_work,y2,iti,ti)

    vx_ef = reshape(vax_work,(/jc(ip),im/))
    
  end function vx_ef

  function ef_vx(vx)
    !
    ! 格子データからスペクトルデータへ変換する.
    !
    real(8), dimension(-lm:lm,2*kc(ip))                  :: ef_vx
                                                      !(out)  スペクトルデータ
    real(8), dimension(js(ip):je(ip),0:im-1), intent(in) :: vx
                                                      !(in) 格子点データ
    real(8), dimension(js(ip):je(ip),2,0:im/2-1) :: vax_work
    real(8), dimension(kc(ip),0:jm-1,2)          :: fya_work
 
    real(8), dimension(jc(ip)*im)       :: y2    
    real(8), dimension(kc(ip)*jm*2)     :: y3

    integer :: k, l

    vax_work = reshape(vx,(/jc(ip),2,im/2/))

    call fttruf(jc(ip),im,vax_work,y2,iti,ti)

    fya_work = fya_vax(vax_work)

    if ( ip == 0 ) fya_work(1,:,2) = 0.0D0
    
    call fttzuf(kc(ip),jm,fya_work,y3,itj,tj)

    do k=1,kc(ip)
       do l=1,lm
          ef_vx(l,k)             = fya_work(k,l,1)
          ef_vx(-l,2*kc(ip)-k+1) = fya_work(k,l,2)
          ef_vx(-l,k)            = fya_work(k,jm-l,1)
          ef_vx(l,2*kc(ip)-k+1)  = fya_work(k,jm-l,2)
       enddo
    enddo

    do k=1,kc(ip)
       ef_vx(0,k)            = fya_work(k,0,1)
       ef_vx(0,2*kc(ip)-k+1) = fya_work(k,0,2)
    enddo

    if ( ip == 0 ) then
       do l=1,lm
          ef_vx(l,1)  = fya_work(1,l,1)
          ef_vx(-l,1) = fya_work(1,l,2)
          ef_vx(l,2*kc(0))  = 0.0D0
          ef_vx(-l,2*kc(0)) = 0.0D0
       end do
       ef_vx(0,1) = fya_work(1,0,1)
       ef_vx(0,2*kc(0)) = 0.0D0
    endif

  end function ef_vx

  function vax_fya(fya)
    real(8), dimension(kc(ip),0:jm-1,2),intent(IN) :: fya
    real(8), dimension(jc(ip),2,0:im/2-1)          :: vax_fya

    real(8), dimension(0:jm-1,2,kc(ip))    :: yaf_work
    real(8), dimension(:,:,:), allocatable :: vaf_work

    integer :: j, k, ips, ipr, ierr, istatus(MPI_STATUS_SIZE)

    vax_fya = 0.0D0

    do k=1,kc(ip)
       do j=0,jm-1
          yaf_work(j,:,k) = fya(k,j,:)
       enddo
    enddo

    do ips=0,np-1
       do ipr=0,np-1
          if ( ipr /= ips ) then
             if ( ip == ips ) then
                allocate(vaf_work(jc(ipr),2,kc(ip)))
                vaf_work = yaf_work(js(ipr):je(ipr),:,:)  
                call MPI_Send(vaf_work,jc(ipr)*2*kc(ips),MPI_REAL8,&
                              ipr,ips,MPI_COMM_WORLD,IERR)
                deallocate(vaf_work)
             endif
             if ( ip == ipr ) then
                call MPI_Recv(vax_fya(1,1,ks(ips)),jc(ipr)*2*kc(ips),MPI_REAL8,&
                              ips,ips,MPI_COMM_WORLD,ISTATUS,IERR)
             endif
          else if ( ip == ipr ) then
             vax_fya(:,:,ks(ip):ke(ip)) = yaf_work(js(ip):je(ip),:,1:kc(ip))
          endif
       enddo
    enddo
    
  end function vax_fya

  function fya_vax(vax)
    real(8), dimension(jc(ip),2,0:im/2-1),intent(IN) :: vax
    real(8), dimension(kc(ip),0:jm-1,2)              :: fya_vax

    real(8), dimension(0:jm-1,2,kc(ip))    :: yaf_work
    real(8), dimension(:,:,:), allocatable :: vaf_work

    integer :: j, k, ips, ipr, ierr, istatus(MPI_STATUS_SIZE)
    
    do ips=0,np-1
       do ipr=0,np-1
          if ( ipr /= ips ) then
             if ( ip == ips ) &
                  call MPI_Send(vax(1,1,ks(ipr)),jc(ips)*2*kc(ipr),MPI_REAL8,&
                                ipr,ips,MPI_COMM_WORLD,IERR)
             if ( ip == ipr ) then
                allocate(vaf_work(jc(ips),2,kc(ip)))
                call MPI_Recv(vaf_work(1,1,1),jc(ips)*2*kc(ipr),MPI_REAL8,&
                                ips,ips,MPI_COMM_WORLD,ISTATUS,IERR)
                yaf_work(js(ips):je(ips),:,:) = vaf_work
                deallocate(vaf_work)
             endif
          else if ( ip == ipr ) then
             yaf_work(js(ip):je(ip),:,1:kc(ip)) = vax(:,:,ks(ip):ke(ip))
          endif
       enddo
    enddo
    
    do k=1,kc(ip)
       do j=0,jm-1
          fya_vax(k,j,:) = yaf_work(j,:,k)
       enddo
    enddo

  end function fya_vax


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

    real(8), dimension(-lm:lm,2*kc(ip)), intent(in)  :: ef
    !(in) 入力スペクトルデータ

    integer kf, k,l
    ! 作業変数

    ef_Lapla_ef = 0.0D0
    do kf=1,kc(ip)
       k = ks(ip) + kf - 1
       do l=-lm,lm
          ef_Lapla_ef(l,kf) = -((2*pi*k/xl)**2+(2*pi*l/yl)**2)*ef(l,kf)
          ef_Lapla_ef(l,2*Kc(ip)-kf+1) &
               = -((2*pi*k/xl)**2+(2*pi*l/yl)**2)*ef(l,2*kc(ip)-kf+1)
       enddo
    enddo
  end function ef_Lapla_ef

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

    real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef
    !(in) スペクトルデータ

    integer kf, k, l

    ef_LaplaInv_ef = 0.0D0
    do kf=1,kc(ip)
       k = ks(ip) + kf - 1
       do l=-lm,lm
          if ( k == 0 .AND. l == 0 ) then
             ef_LaplaInv_ef(l,kf) = 0.0
          else
             ef_LaplaInv_ef(l,kf) = -ef(l,kf)/((2*pi*k/xl)**2+(2*pi*l/yl)**2)
             ef_LaplaInv_ef(l,2*kc(ip)-kf+1) &
                  = -ef(l,2*kc(ip)-kf+1)/((2*pi*k/xl)**2+(2*pi*l/yl)**2)
          endif
       enddo
    enddo
  end function ef_LaplaInv_ef

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

    real(8), dimension(-lm:lm,2*kc(ip)), intent(in)  :: ef
    !(in) 入力スペクトルデータ

    integer kf, k, l
    ! 作業変数

    do kf=1,kc(ip)
       k = ks(ip) + kf - 1
       do l=-lm,lm
          ef_Dx_ef(l,kf) = -(2*pi*k/xl)*ef(-l,2*kc(ip)-kf+1)
          ef_Dx_ef(l,2*kc(ip)-kf+1) = -(2*pi*(-k)/xl)*ef(-l,kf)
       enddo
    enddo
  end function ef_Dx_ef

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

    real(8), dimension(-lm:lm,2*kc(ip)), intent(in)  :: ef
    !(in) 入力スペクトルデータ

    integer kf,l
    ! 作業変数

    do kf=1,kc(ip)
       if ( ip == 0 .AND. kf == 1 ) then
          do l=-lm,lm
             ef_Dy_ef(l,kf) = -(2*pi*l/yl)*ef(-l,kf)
             ef_Dy_ef(l,2*kc(ip)-kf+1) = 0.0D0
          enddo
       else
          do l=-lm,lm
             ef_Dy_ef(l,kf) = -(2*pi*l/yl)*ef(-l,2*kc(ip)-kf+1)
             ef_Dy_ef(l,2*kc(ip)-kf+1) = -(2*pi*l/yl)*ef(-l,kf)   
          enddo
       endif
    enddo
    
  end function ef_Dy_ef

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

    real(8), dimension(-lm:lm,2*kc(ip)), intent(in)  :: ef_a
    !(in) 1つ目の入力スペクトルデータ

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

    ef_Jacobian_ef_ef = ef_vx(   vx_ef(ef_Dx_ef(ef_a))*vx_ef(ef_Dy_ef(ef_b)) &
                               - vx_ef(ef_Dy_ef(ef_a))*vx_ef(ef_Dx_ef(ef_b)) )

  end function ef_Jacobian_ef_ef

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

    real(8), dimension(-lm:lm,2*kc(ip)), intent(in)  :: ef_Zeta
    !(in) 渦度スペクトルデータ

    real(8), dimension(jc(ip),0:im-1)   :: vx_U    ! 速度 X 成分
    real(8), dimension(jc(ip),0:im-1)   :: vx_V    ! 速度 Y 成分
    real(8), dimension(-lm:lm,2*kc(ip)) :: ef_UV   ! 

    vx_U = - vx_ef(ef_Dy_ef(ef_LaplaInv_ef(ef_Zeta)))
    vx_V =   vx_ef(ef_Dx_ef(ef_LaplaInv_ef(ef_Zeta)))
    ef_UV = ef_vx(vx_U*vx_V)
    
    ef_JacobianZ_ef = - ( ef_Dx_ef(ef_Dx_ef(ef_UV))-ef_Dy_ef(ef_Dy_ef(ef_UV)) )&
                      - ef_Dx_ef(ef_Dy_ef(ef_vx(vx_V**2 - vx_U**2)))

  end function ef_JacobianZ_ef

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

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

    real(8) :: IntYXTmp
    integer :: i, j, ierr
    ! 作業変数

    IntYXTmp = 0.0d0
    do i=0,im-1
       do j=1,jc(ip)
          IntYXTmp = IntYXTmp + vx(j,i) * v_Y_Weight(j) * x_X_Weight(i)
       enddo
    enddo

    CALL MPI_ALLREDUCE(IntYXTmp,IntYX_vx,1,MPI_REAL8, &
                       MPI_SUM,MPI_COMM_WORLD,IERR)
    
  end function IntYX_vx

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

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

    integer :: i
    ! 作業変数

    v_IntX_vx = 0.0d0
    do i=0,im-1
       v_IntX_vx(:) = v_IntX_vx(:) + vx(:,i) * x_X_Weight(i)
    enddo
  end function v_IntX_vx

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

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

    real(8), dimension(0:im-1)        :: x_IntYTmp
    integer :: j, ierr
    ! 作業変数

    x_IntYTmp = 0.0d0
    do j=1,jc(ip)
       x_IntYTmp(:) = x_IntYTmp(:) + vx(j,:) * v_Y_Weight(j)
    enddo

    CALL MPI_ALLREDUCE(x_IntYTmp,x_IntY_vx,im,MPI_REAL8, &
                       MPI_SUM,MPI_COMM_WORLD,IERR)

  end function x_IntY_vx

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

    real(8)                      :: IntYTmp     !(out) 積分値
    integer :: ierr
    
    IntYTmp = sum(v*v_Y_Weight)

    CALL MPI_ALLREDUCE(IntYTmp,IntY_v,1,MPI_REAL8, &
                       MPI_SUM,MPI_COMM_WORLD,IERR)
    
  end function IntY_v

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

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

    AvrYX_vx = IntYX_vx(vx)/(xl*yl)

  end function AvrYX_vx

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

    real(8), dimension(jc(ip))          :: v_AvrX_vx
    !(out) 平均された 1 次元(Y)格子点データ

    v_AvrX_vx = v_IntX_vx(vx)/xl

  end function v_AvrX_vx

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

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

    x_AvrY_vx = x_IntY_vx(vx)/yl

  end function x_AvrY_vx

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

  end function AvrX_x

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

    AvrY_v = IntY_v(v)/yl

  end function AvrY_v
  

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

    real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef_StrFunc
    ! 流線関数

    integer kf, k,l
    ! 作業変数

    do kf=1,kc(ip)
       k = ks(ip) + kf -1
       do l=-lm,lm
          ef_EnergyFromStreamfunc_ef(l,kf) &
               = 0.5 * ( (2*pi*k/xl)**2 + (2*pi*l/yl)**2 ) & 
               * (ef_StrFunc(l,kf)**2 + ef_StrFunc(-l,2*kc(ip)-kf+1)**2)
          ef_EnergyFromStreamfunc_ef(l,2*kc(ip)-kf+1) &
               = ef_EnergyFromStreamfunc_ef(l,kf)
       enddo
    enddo

  end function ef_EnergyFromStreamfunc_ef

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

    real(8), dimension(-lm:lm,2*kc(ip))    :: ef_EnstrophyFromStreamfunc_ef
    ! エンストロフィーースペクトル

    real(8), dimension(-lm:lm,2*kc(ip)), intent(in) :: ef_StrFunc
    ! 流線関数

    integer kf, k,l
    ! 作業変数

    do kf=1,kc(ip)
       k = ks(ip) + kf -1
       do l=-lm,lm
          ef_EnstrophyFromStreamfunc_ef(l,kf) &
               = 0.5 * ( (2*pi*k/xl)**2 + (2*pi*l/yl)**2 )**2 & 
               * (ef_StrFunc(l,kf)**2 + ef_StrFunc(-l,2*kc(ip)-kf+1)**2)
          ef_EnstrophyFromStreamfunc_ef(l,2*kc(ip)-kf+1) &
               = ef_EnstrophyFromStreamfunc_ef(l,kf)
       enddo
    enddo

  end function ef_EnstrophyFromStreamfunc_ef


  !--------------- 補間計算 -----------------
  function Interpolate_ef( ef_Data, x, y )
    real(8), intent(IN)  :: ef_data(-lm:lm,2*kc(ip))  ! スペクトルデータ
    real(8), intent(IN)  :: x                         ! 補間する点の x 座標 
    real(8), intent(IN)  :: y                         ! 補間する点の y 座標 
    real(8)              :: Interpolate_ef            ! 補間した値
    real(8)              :: InterpolateTmp            ! 作業変数

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

    xx =(2*PI/xl)*(x - xmin)
    yy =(2*PI/yl)*(y - ymin)

    if ( ip == 0 ) then
       InterpolateTmp = ef_Data(0,1)
    else
       InterpolateTmp = 0.0D0
    endif
    
    ! l=0
    if ( ip == 0 ) then
       do kf=2,kc(ip)
          k = ks(ip) + kf -1
          InterpolateTmp = InterpolateTmp &
               + 2*( ef_Data(0,kf)*cos(k*xx) - ef_Data(0,kf_k(-k))*sin(k*xx) )
       end do
    else
       do kf=1,kc(ip)
          k = ks(ip) + kf -1
          InterpolateTmp = InterpolateTmp &
               + 2*( ef_Data(0,kf)*cos(k*xx) - ef_Data(0,kf_k(-k))*sin(k*xx) )
       end do
    endif
       
    ! k=0
    if ( ip == 0 ) then
       do l=1,lm
          InterpolateTmp = InterpolateTmp &
            + 2*( ef_Data(l,kf_k(0))*cos(l*yy) - ef_Data(-l,kf_k(0))*sin(l*yy) )
       end do
    end if

    ! k*l > 0
    do l=1,lm
       if ( ip == 0 ) then
          do kf=2,kc(ip)
             k = ks(ip) + kf -1
             InterpolateTmp = InterpolateTmp &
                  + 2*(  ef_Data(l,kf)*(   cos(k*xx)*cos(l*yy)   &
                                         - sin(k*xx)*sin(l*yy) ) &
                        -ef_Data(-l,kf_k(-k))*(   sin(k*xx)*cos(l*yy)   &
                                                + cos(k*xx)*sin(l*yy) ) )
          end do
       else
          do kf=1,kc(ip)
             k = ks(ip) + kf - 1
             InterpolateTmp = InterpolateTmp &
                  + 2*(  ef_Data(l,kf)*(   cos(k*xx)*cos(l*yy)   &
                                         - sin(k*xx)*sin(l*yy) ) &
                        -ef_Data(-l,kf_k(-k))*(   sin(k*xx)*cos(l*yy)   &
                                                + cos(k*xx)*sin(l*yy) ) )
          end do
       endif
    end do

    ! k*l < 0
    do l=1,lm
       if ( ip == 0 ) then
          do kf=2,kc(ip)
             k = ks(ip) + kf - 1
             InterpolateTmp = InterpolateTmp &
                  + 2*(  ef_Data(-l,kf)*(   cos(k*xx)*cos(l*yy)   &
                                          + sin(k*xx)*sin(l*yy) ) &
                        -ef_Data(l,kf_k(-k))*(   sin(k*xx)*cos(l*yy)   &
                                               - cos(k*xx)*sin(l*yy) ) )
          end do
       else
          do kf=1,kc(ip)
             k = ks(ip) + kf - 1
             InterpolateTmp = InterpolateTmp &
                  + 2*(  ef_Data(-l,kf)*(   cos(k*xx)*cos(l*yy)   &
                                          + sin(k*xx)*sin(l*yy) ) &
                        -ef_Data(l,kf_k(-k))*(   sin(k*xx)*cos(l*yy)   &
                                               - cos(k*xx)*sin(l*yy) ) )
          end do
       endif
    end do

    CALL MPI_ALLREDUCE(InterpolateTmp,Interpolate_ef,1,MPI_REAL8, &
                       MPI_SUM,MPI_COMM_WORLD,IERR)

  end function Interpolate_ef

end module ee_mpi_module
