ee_mpi_module.f90

Go to the documentation of this file.

!--
!----------------------------------------------------------------------
! Copyright (c) 2016-2020 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ee_mpi_module
!
!      spml/ee_module モジュールは周期境界条件の下での 2 次元矩形領域の
!      流体運動をスペクトル法により数値計算するための Fortran90 関数を
!      提供する. 
!
!      内部で ISPACK3/FXPACK の Fortran77 サブルーチンを呼んでいる.
!      スペクトルデータおよび格子点データの格納方法については...
!
!履歴  2016/08/08  竹広真一  eee_mpi_module より改造
!      2016/09/22  竹広真一  ef_Lapla_ef, ef_Laplainv_ef 初期化
!      2020/04/10  竹広真一  MPI_Send, MPI_Recv 関数を MPI_AlltoAll へ変更
!      2020/05/24  竹広真一  MPI_AlltoAllV へ変更
!
!++
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                  ! エンストロフィー

  public vax_fay, fay_vax

  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(8), dimension(:), allocatable :: itj          ! FFT 用配列(Y)
  real(8),    dimension(:), allocatable :: tj           ! FFT 用配列(Y)
  integer(8), 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 :: iip
    integer :: kp, jp
    integer :: kint, kmod
    integer :: jint, jmod
    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)

    ! x 波数 MPI分割情報
    allocate(ks(0:np-1))
    allocate(ke(0:np-1))
    allocate(kc(0:np-1))
    
    kint = km/np
    kmod = mod(km,np)
    kc   = kint
    do iip=0,kmod
       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

    ! Y 軸 MPI分割情報
    allocate(js(0:np-1))
    allocate(je(0:np-1))
    allocate(jc(0:np-1))

    jint = jm/np
    jmod = mod(jm,np)
    jc   = jint
    do iip=0,jmod-1
       jc(iip) = jc(iip) + 1
    end do

    js(0)=0 ; je(0) = js(0) + jc(0) -1
    do iip=1,np-1
       js(iip) = js(iip-1) + jc(iip-1)
       je(iip) = js(iip) + jc(iip) - 1
    end do
    
    allocate(itj(jm))
    allocate(iti(im/2))
    allocate(tj(jm*2))
    allocate(ti(im*3/2))

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

    call fxrini(int(im,8),iti,ti)
    call fxzini(int(jm,8),itj,tj)


    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 (2021/04/27) 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),2,0:jm-1)          :: fay_work
 
    integer :: k, l

    fay_work = 0.0d0

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

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

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

    call fxztba(int(kc(ip),8),int(jm,8),fay_work,itj,tj)

    vax_work = vax_fay(fay_work)

    call fxrtba(int(jc(ip),8),int(im,8),vax_work,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),2,0:jm-1)          :: fay_work
    integer :: k, l

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

    call fxrtfa(int(jc(ip),8),int(im,8),vax_work,iti,ti)

    fay_work = fay_vax(vax_work)

    if ( ip == 0 ) fay_work(1,2,:) = 0.0d0

    call fxztfa(int(kc(ip),8),int(jm,8),fay_work,itj,tj)

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

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

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

  end function ef_vx

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

    integer :: j, k, iip
    integer :: ierr
    integer :: isndc(0:np-1)
    integer :: isdsp(0:np-1)
    integer :: ircvc(0:np-1)
    integer :: irdsp(0:np-1)
    real(8), allocatable :: aa_work(:,:)
!!$    real(8), allocatable :: xava_work(:,:,:,:)
    integer :: kcm
    integer :: ind1, ind2
    
    vax_fay = 0.0d0

    kcm = maxval(kc)
    allocate(aa_work(kcm*2*jc(ip),0:np-1))
    aa_work = 0.0d0
!!$    kcm = maxval(kc)
!!$    allocate(xava_work(kcm,2,jc(ip),0:np-1))
!!$    xava_work = 0.0D0

    isndc = kc(ip)*2*jc
    ircvc = jc(ip)*2*kc
    do iip=0,np-1
       isdsp(iip) = sum(isndc(0:iip-1))
       irdsp(iip) = kcm*2*jc(ip)*iip
    end do

    call mpi_alltoallv(fay,isndc,isdsp,mpi_real8,&
         aa_work,ircvc,irdsp,mpi_real8,mpi_comm_world,ierr)
!!$    call MPI_AlltoAllV(fay,isndc,isdsp,MPI_REAL8,&
!!$         xava_work,ircvc,irdsp,MPI_REAL8,MPI_COMM_WORLD,IERR)

    do iip=0,np-1
       do k=1,kc(iip)
          do j=1,jc(ip)
             ind1 = k           + kc(iip)*2*(j-1)
             ind2 = k + kc(iip) + kc(iip)*2*(j-1)
             vax_fay(j,1,ks(iip)-1+k) = aa_work(ind1,iip)
             vax_fay(j,2,ks(iip)-1+k) = aa_work(ind2,iip)
!!$             vax_fay(j,:,ks(iip)-1+k) = xava_work(k,:,j,iip)
          end do
       end do
    end do

    deallocate(aa_work)
!!$    deallocate(xava_work)

  end function vax_fay

  function fay_vax(vax)
    real(8), dimension(jc(ip),2,0:im/2-1),intent(IN) :: vax
    real(8), dimension(kc(ip),2,0:jm-1)              :: fay_vax
    
    integer :: j, k, iip
    integer :: ierr
    integer :: isndc(0:np-1)
    integer :: isdsp(0:np-1)
    integer :: ircvc(0:np-1)
    integer :: irdsp(0:np-1)
    real(8), allocatable :: aa_work(:,:)
!!$    real(8), allocatable :: vafa_work(:,:,:,:)
    integer :: jcm
    integer :: ind1, ind2

    fay_vax = 0.0d0

    jcm = maxval(jc)
    allocate(aa_work(jcm*2*kc(ip),0:np-1))
    aa_work = 0.0d0
!!$    allocate(vafa_work(jcm,2,kc(ip),0:np-1))
!!$    vafa_work = 0.0D0

    isndc = jc(ip)*2*kc
    ircvc = kc(ip)*2*jc
    do iip=0,np-1
       isdsp(iip) = sum(isndc(0:iip-1))
       irdsp(iip) = jcm*2*kc(ip)*iip
    end do

    call mpi_alltoallv(vax,isndc,isdsp,mpi_real8,&
         aa_work,ircvc,irdsp,mpi_real8,mpi_comm_world,ierr)
!!$    call MPI_AlltoAllV(vax,isndc,isdsp,MPI_REAL8,&
!!$         vafa_work,ircvc,irdsp,MPI_REAL8,MPI_COMM_WORLD,IERR)

    do iip=0,np-1
       do j=1,jc(iip)
          do k=1,kc(ip)
             ind1 = j           + jc(iip)*2*(k-1)
             ind2 = j + jc(iip) + jc(iip)*2*(k-1)
             fay_vax(k,1,js(iip)+j-1) = aa_work(ind1,iip)
             fay_vax(k,2,js(iip)+j-1) = aa_work(ind2,iip)
!!$             fay_vax(k,:,js(iip)+j-1) = vafa_work(j,:,k,iip)
          end do
       end do
    end do

    deallocate(aa_work)
!!$    deallocate(vafa_work)

  end function fay_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