!--
!----------------------------------------------------------------------
! Copyright (c) 2008-2009 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!
!表題  eee_mpi_module
!
!      spml/eee_module モジュールは周期境界条件の下での 3 次元矩形領域の
!      流体運動をスペクトル法により数値計算するための Fortran90 関数を
!      提供する. 
!
!      内部で ISPACK/P3PACK, P3PACK-MPI の Fortran77 サブルーチンを呼んでいる.
!      スペクトルデータおよび格子点データの格納方法については
!      ISPACK/P3PACK のマニュアルを参照されたい. 
!
!履歴  2008/05/21  竹広真一  eee_module より改造
!      2008/06/03  竹広真一  ee2f_ZetaFromVor_eef_eef_eef bug fix
!      2008/06/03  竹広真一  ESpectralFromZeta bug fix
!      2009/02/23  佐々木洋平  RDoc 用のドキュメントを修正
!      2010/01/07  佐々木洋平  include 'mpif.h' -> use mpi
!
!++
module eee_mpi_module
  !
  != eee_mpi_module
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/eee_module モジュールは周期境界条件の下での 3 次元矩形領域の
  ! 流体運動をスペクトル法により数値計算するための Fortran90 関数を
  ! 提供する.
  !
  ! 内部で ISPACK/P3PACK の Fortran77 サブルーチンを呼んでいる.
  ! スペクトルデータおよび格子点データの格納方法については
  ! ISPACK/P3PACK のマニュアルを参照されたい.
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (eef_, zxv_, x_, v_, z_) は, 返す値の形を示している.
  !   eef_ :: スペクトルデータ
  !           (第 1,2 3 次元がそれぞれ Z,Y,X 方向波数(X について分割配置))
  !   ee2f :: 2 つのスペクトルデータの並んだもの
  !   zxv_ :: 3 次元格子点データ
  !           (第 1,2 3 次元がそれぞれ Z,X,Y 方向の格子点(Y について分割配置)
  !   xv_  :: XY 方向 2 次元格子点データ, zv_ : YZ 方向 2 次元格子点データ
  !   zx_  :: XZ 方向 2 次元格子点データ
  !   x_   :: X 方向 1 次元格子点データ, v_ : Y 方向 1 次元格子点データ
  !   z_   :: Z 方向 1 次元格子点データ
  !
  ! * 関数名の間の文字列(Dx, Dy, Dz, Lapla, LaplaInv)は,
  !   その関数の作用を表している.
  !
  ! * 関数名の最後 (_eee_eee,_eee,_zyx, _x, _v) は, 入力変数の形が
  !   スペクトルデータおよび格子点データであることを示している.
  !    _eef     :: スペクトルデータ
  !    _eef_eef :: 2 つのスペクトルデータ
  !    _ee2f    :: 2 つのスペクトルデータの並んだもの
  !    _zxv     :: 3 次元格子点データ
  !    _x       :: X 方向 1 次元格子点データ
  !    _v       :: Y 方向 1 次元格子点データ.
  !
  !=== 各データの種類の説明
  !
  ! * zxv : 2 次元格子点データ.
  !   * 変数の種類と次元は real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)). 
  !   * im, km はそれぞれ X, Z 座標の格子点数であり, サブルーチン 
  !     eee_mpi_initial にてあらかじめ設定しておく.
  !   * eee_mpi_Initial によってプロセス毎に Y 方向の
  !     分割配置が設定され, その情報が public 変数 js,je,jc に設定される. 
  !   * js が分割配置の最小格子点, je が最大格子点, jc が格子点数である. 
  !   * 第 1 次元が Z 座標の格子点位置番号, 
  !     第 2 次元が X 座標の格子点位置番号, 
  !     第 3 次元が Y 座標であることに注意(X, Y, Zの順ではない).
  !
  ! * eee : スペクトルデータ.
  !   * 変数の種類と次元は real(8), dimension(-nm:nm,-mm:mm,2*lc). 
  !   * mm, nm はそれぞれ Y, Z方向の最大波数であり, サブルーチン 
  !     eee_initial にてあらかじめ設定しておく
  !     (X, Y, Z 方向波数の順ではない)ことに注意. 
  !   * eee_mpi_Initial によってプロセス毎に X 方向波数の
  !     分割配置が設定され, その情報が public 変数 ls,le,lc に設定される. 
  !   * ls が分割配置の最小波数, le が最大波数, jc が波数の数である. 
  !     実際の格納のされ方は ls, ls+1,...,le,-le,-le+1,...,-ls の順である.
  !     ISPACK/P3PACK-MPI マニュアル参照のこと.
  !
  ! * ee2f : 2 つのスペクトルデータのならび. 
  !   * 変数の種類と次元は real(8), dimension(-nm:nm,-mm:mm,2,2*lc). 
  !
  ! * x, v, z : X, Y, Z 方向 1 次元格子点データ.
  !   * 変数の種類と次元はそれぞれ real(8), dimension(0:im-1),
  !     real(8), dimension(js(ip):je(ip)) および real(8), dimension(0:km-1).
  !
  ! * eee_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * zxv_ で始まる関数が返す値は 3 次元格子点データに同じ.
  !
  ! * x_, v_ z_ で始まる関数が返す値は 1 次元格子点データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! eee_Initial :: スペクトル変換の格子点数, 切断波数の設定
  ! eee_ChangeResolution :: 解像度設定の変更
  ! lf_l        :: X 方向波数の位置情報
  !
  !==== 座標変数
  !
  ! x_X, v_Y, z_Z    ::  格子点座標(X,Y座標)を格納した 1 次元配列
  ! x_X_Weight, v_Y_Weight,  z_Z_Weight ::  重み座標を格納した 1 次元配列
  ! zxv_X, zxv_Y, zxv_Z  :: 格子点データの XYZ 座標(X,Y,Z)
  !                         (格子点データ型 3 次元配列)
  !
  !==== 基本変換
  !
  ! zxv_eef   :: スペクトルデータから格子データへの変換
  ! eef_zxv   :: 格子データからスペクトルデータへの変換
  ! xyz_zxv   :: 分散格子点データの集積
  ! eee2_ee2f :: 分散スペクトルデータの集積
  ! ee2f_eee2 :: 集積スペクトルデータの分散
  !
  !==== 微分
  !
  ! eef_Lapla_eef       :: スペクトルデータにラプラシアンを作用させる
  ! eef_LaplaInv_eef    :: スペクトルデータにラプラシアンの逆変換を作用させる
  ! eef_Dx_eef          :: スペクトルデータに X 微分を作用させる
  ! eef_Dy_eef          :: スペクトルデータに Y 微分を作用させる
  ! eef_Dz_eef          :: スペクトルデータに Z 微分を作用させる
  !
  !==== 非線形計算
  !
  ! ee2f_RotVelxVor_ee2f   :: Euler 方程式の非線形項を計算する
  ! eef_VorFromZeta_ee2f   :: 渦度 2 成分(ζ_1, ζ_2)ら渦度の 1 成分を計算する
  ! eef_VelFromZeta_ee2f   :: 渦度 2 成分(ζ_1, ζ_2)ら速度の 1 成分を計算する
  ! ee2f_ZetaFromVor_eef_eef_eef :: 渦度から渦度 2 成分(ζ_1, ζ_2)を計算する
  !
  !==== 積分・平均
  !
  ! IntZXV_zxv, AvrZXV_zxv   :: 3 次元格子点データの全領域積分および平均
  ! z_IntXV_zxv, z_AvrXV_zxv :: 3 次元格子点データの X,Y 方向積分および平均
  ! x_IntZV_zxv, x_AvrZV_zxv :: 3 次元格子点データの Y,Z 方向積分および平均
  ! y_IntZX_zxv, y_AvrZX_zxv :: 3 次元格子点データの Z,X 方向積分および平均
  ! zv_IntX_zxv, zv_AvrX_zxv :: 3 次元格子点データの X 方向積分および平均
  ! zx_IntY_zxv, zx_AvrY_zxv :: 3 次元格子点データの Y 方向積分および平均
  ! vx_IntZ_zxv, vx_AvrZ_zxv :: 3 次元格子点データの Z 方向積分および平均
  !
  ! IntXV_xv,  AvrXV_xv  :: 2 次元(XY)格子点データの X,Y 方向積分および平均
  ! v_IntX_xv, v_AvrX_xv :: 2 次元(XY)格子点データの X 方向積分および平均
  ! x_IntV_xv, x_AvrV_xv :: 2 次元(XY)格子点データの Y 方向積分および平均
  ! IntZX_zx, AvrZX_zx   :: 2 次元(ZX)格子点データの Z,X 方向積分および平均
  ! z_IntX_zx, z_AvrX_zx :: 2 次元(ZX)格子点データの X 方向積分および平均
  ! x_IntZ_zx, x_AvrZ_zx :: 2 次元(ZX)格子点データの Z 方向積分および平均
  ! IntZV_zv, AvrZV_zv   :: 2 次元(YZ)格子点データの Y,Z 方向積分および平均
  ! v_IntZ_zv, v_AvrZ_zv :: 2 次元(YZ)格子点データの Z 方向積分および平均
  ! z_IntV_zv, z_AvrV_zv :: 2 次元(YZ)格子点データの Y 方向積分および平均
  !
  ! IntX_x, AvrX_x       :: 1 次元(X)格子点データの X 方向積分および平均
  ! IntV_v, AvrV_v       :: 1 次元(Y)格子点データの Y 方向積分および平均
  ! IntZ_z, AvrZ_z       :: 1 次元(Z)格子点データの Z 方向積分および平均
  !
  !==== スペクトル解析
  !
  ! EnergyHelicityFromZeta_ee2f  :: 全エネルギーと全ヘリシティーを計算する. 
  ! ESpectralFromZeta            :: エネルギースペクトルを計算する. 
  !
  use dc_message, only : MessageNotify
  use mpi
  implicit none

  private
  public ls, le, lc
  public js, je, jc
  public lf_l
  public eee_mpi_Initial                                  ! 初期化ルーチン
  public eee_ChangeResolution                             ! 解像度設定の変更
  public zxv_eef, eef_zxv                                 ! 基本変換
  public xyz_zxv                                          ! 基本変換
  public eee2_ee2f, ee2f_eee2                             ! 基本変換
  public eef_Dx_eef, eef_Dy_eef, eef_Dz_eef               ! 微分
  public eef_Lapla_eef, eef_LaplaInv_eef                  ! 微分
  public ee2f_RotVelxVor_ee2f                             ! 非線形計算
  public eef_VorFromZeta_ee2f, eef_VelFromZeta_ee2f       ! 渦度速度計算
  public ee2f_ZetaFromVor_eef_eef_eef                     ! 渦度変換

  public IntZXV_zxv, AvrZXV_zxv                           ! 積分・平均
  public z_IntXV_zxv, z_AvrXV_zxv                         ! 積分・平均
  public x_IntZV_zxv, x_AvrZV_zxv                         ! 積分・平均
  public v_IntZX_zxv, v_AvrZX_zxv                         ! 積分・平均
  public zv_IntX_zxv, zv_AvrX_zxv                         ! 積分・平均
  public zx_IntV_zxv, zx_AvrV_zxv                         ! 積分・平均
  public xv_IntZ_zxv, xv_AvrZ_zxv                         ! 積分・平均

  public IntXV_xv,  AvrXV_xv                              ! 積分・平均
  public v_IntX_xv, v_AvrX_xv                             ! 積分・平均
  public x_IntV_xv, x_AvrV_xv                             ! 積分・平均
  public IntZX_zx, AvrZX_zx                               ! 積分・平均
  public z_IntX_zx, z_AvrX_zx                             ! 積分・平均
  public x_IntZ_zx, x_AvrZ_zx                             ! 積分・平均
  public IntZV_zv, AvrZV_zv                               ! 積分・平均
  public v_IntZ_zv, v_AvrZ_zv                             ! 積分・平均
  public z_IntV_zv, z_AvrV_zv                             ! 積分・平均

  public IntX_x, AvrX_x                                   ! 積分・平均
  public IntV_v, AvrV_v                                   ! 積分・平均
  public IntZ_z, AvrZ_z                                   ! 積分・平均

  public EnergyHelicityFromZeta_ee2f                      ! エネルギーヘリシティー
  public ESpectralFromZeta                                ! エネルギースペクトル

  public x_X, v_Y, z_Z                                    ! 座標変数
  public x_X_Weight, v_Y_Weight, z_Z_Weight               ! 座標変数
  public zxv_X, zxv_Y, zxv_Z                              ! 座標変数

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

  integer, dimension(:),   pointer :: le => null()      ! X 波数範囲, 数
  integer, dimension(:),   pointer :: ls => null()      ! X 波数範囲, 数
  integer, dimension(:),   pointer :: lc => null()      ! X 波数範囲, 数
  integer, dimension(:),   pointer :: je => null()      ! Y 座標格子点範囲, 数
  integer, dimension(:),   pointer :: js => null()      ! Y 座標格子点範囲, 数
  integer, dimension(:),   pointer :: jc => null()      ! Y 座標格子点範囲, 数
  integer   :: np, ip                             ! MPI プロセス数, プロセス ID

  integer, dimension(:),   pointer :: itk => null()
  real(8), dimension(:),   pointer :: tk => null()
  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 :: v_Y => null()   ! 格子点座標(Y)
  real(8), dimension(:),   pointer :: z_Z => null()   ! 格子点座標(Y)

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

  real(8), dimension(:,:,:), pointer :: zxv_X => null()
                          ! 格子点(X)座標(3 次元)
                          ! 各格子点(i,j,k)の位置の X 座標を格納した格子データ
  real(8), dimension(:,:,:), pointer :: zxv_Y => null()
                          ! 格子点(Y)座標(3 次元)
                          ! 各格子点(i,j,k)の位置の Y 座標を格納した格子データ
  real(8), dimension(:,:,:), pointer :: zxv_Z => null()
                          ! 格子点(Z)座標(3 次元)
                          ! 各格子点(i,j,k)の位置の Z 座標を格納した格子データ

  real(8), dimension(:), pointer :: w => null(), ws => null(), wg => null()
  real(8), dimension(:), pointer :: sgwork => null()

  integer, parameter :: nparams_max = 10  ! eee_Initial を呼べる最大回数
  type eee_param                          ! 解像度領域情報構造体
     integer   :: im, jm, km
     integer   :: lm, mm, nm
     integer, dimension(:),     pointer :: ls, le, lc
     integer, dimension(:),     pointer :: js, je, jc
     integer, dimension(:),     pointer :: itk
     real(8), dimension(:),     pointer :: tk
     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 :: v_Y
     real(8), dimension(:),     pointer :: z_Z
     real(8), dimension(:),     pointer :: x_X_Weight
     real(8), dimension(:),     pointer :: v_Y_Weight
     real(8), dimension(:),     pointer :: z_Z_Weight
     real(8), dimension(:,:,:), pointer :: zxv_X
     real(8), dimension(:,:,:), pointer :: zxv_Y 
     real(8), dimension(:,:,:), pointer :: zxv_Z
     real(8), dimension(:),     pointer :: w, ws, wg
     real(8), dimension(:),     pointer :: sgwork
  end type eee_param
  type(eee_param) :: params(nparams_max)  ! 解像度領域情報
  integer :: nparams                      ! 解像度領域情報の個数

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

  save im, jm, km, lm, mm, nm, itk, tk, itj, tj, iti, ti
  save ls, le, lc, js, je, jc, np, ip
  save x_X, v_Y, z_Z, x_X_Weight, v_Y_Weight, z_Z_Weight, zxv_X, zxv_Y, zxv_Z
  save params, nparams

  contains
  !--------------- 初期化 -----------------
    subroutine eee_mpi_Initial(i,j,k,l,m,n,id)
      !
      ! スペクトル変換の格子点数, 波数を設定する.
      ! 領域の範囲は [0,2pi]x[0,2pi]x[0,2pi] に決め打ち
      !
      ! 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
      ! 初期設定をしなければならない.
      !
      ! オプショナル引数 id を用いて異なる解像度, 領域を同時に
      ! 扱うことができる. eee_Initial を解像度, 領域ごとに呼んで
      ! 引数 id をキープし, eee_ChangeResolution で切替える.
      !
      integer,intent(in) :: i           ! 格子点数(X)
      integer,intent(in) :: j           ! 格子点数(Y)
      integer,intent(in) :: k           ! 格子点数(Z)
      integer,intent(in) :: l           ! 切断波数(X)
      integer,intent(in) :: m           ! 切断波数(Y)
      integer,intent(in) :: n           ! 切断波数(X)

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

      character(len=3) cid
      integer :: ii, jj, kk
      integer :: lls, lle, llc, jjs, jje, jjc
      integer :: iip
      integer :: lp, jp
      integer :: ierr

      im = i         ; jm = j         ; km = k
      lm = l         ; mm = m         ; nm = n

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

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

      lp=lm/np+1
      lls=lp*ip
      lle=min(lp*(ip+1)-1,lm)
      if(lle.ge.lls) then
         llc=lle-lls+1
      else
         llc=0 ;  lls=0  ; lle=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

      nparams = nparams + 1

      params(nparams)%im = im
      params(nparams)%jm = jm
      params(nparams)%km = km
      params(nparams)%lm = lm
      params(nparams)%mm = mm
      params(nparams)%nm = nm
      allocate(params(nparams)%ls(0:np-1))
      allocate(params(nparams)%le(0:np-1))
      allocate(params(nparams)%lc(0:np-1))
      allocate(params(nparams)%js(0:np-1))
      allocate(params(nparams)%je(0:np-1))
      allocate(params(nparams)%jc(0:np-1))

      allocate(params(nparams)%itk(5))
      allocate(params(nparams)%itj(5))
      allocate(params(nparams)%iti(5))
      allocate(params(nparams)%tk(km*2))
      allocate(params(nparams)%tj(jm*2))
      allocate(params(nparams)%ti(im*2))
      allocate(params(nparams)%w(km*max(im*((jm-1)/np+1),jm*2*(lm/np+1))))
      allocate(params(nparams)%ws((2*mm+1)*(2*nm+1)*2*2*(lm+1)))
      allocate(params(nparams)%wg(4*km*max(im*((jm-1)/np+1),jm*2*(lm/np+1))))
      allocate(params(nparams)%sgwork(km*max(im*((jm-1)/np+1),jm*2*(lm/np+1))))

      allocate(params(nparams)%x_X(0:im-1))
      allocate(params(nparams)%x_X_Weight(0:im-1))
      allocate(params(nparams)%v_Y(jjs:jje))
      allocate(params(nparams)%v_Y_Weight(jjs:jje))
      allocate(params(nparams)%z_Z(0:km-1))
      allocate(params(nparams)%z_Z_Weight(0:km-1))
      allocate(params(nparams)%zxv_X(0:km-1,0:im-1,jjs:jje))
      allocate(params(nparams)%zxv_Y(0:km-1,0:im-1,jjs:jje))
      allocate(params(nparams)%zxv_Z(0:km-1,0:im-1,jjs:jje))

      call eee_ChangeResolution(nparams)

      call p3init(km,jm,im,itk,tk,itj,tj,iti,ti)

      do iip=0,np-1
         lp=lm/np+1
         ls(iip)=lp*iip
         le(iip)=min(lp*(iip+1)-1,lm)
         if(le(iip).ge.ls(iip)) then
            lc(iip)=le(iip)-ls(iip)+1
         else
            lc(iip)=0 ;  ls(iip)=0  ; le(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

      do ii=0,im-1
         x_X(ii) = 2*pi/im*ii
      enddo
      x_X_Weight = 2*pi/im

      do jj=js(ip),je(ip)
         v_Y(jj) = 2*pi/jm*jj
      enddo
      v_Y_Weight = 2*pi/jm

      do kk=0,km-1
         z_Z(kk) = 2*pi/km*kk
      enddo
      z_Z_Weight = 2*pi/km

      zxv_X = spread(spread(x_X,1,km),3,jc(ip))
      zxv_Y = spread(spread(v_Y,1,im),1,km)
      zxv_Z = spread(spread(z_Z,2,im),3,jc(ip))

      if ( present(id) ) id = nparams

      write(cid,'(I3)') nparams
      call MessageNotify('M','eee_mpi_initial','eee_mpi_module is initialized')
      call MessageNotify('M','eee_mpi_initial',&
           'Resolution ID is '//trim(adjustl(cid)))
    end subroutine eee_mpi_Initial

  !--------------- id 変更 -----------------
    subroutine eee_ChangeResolution(id)
      !
      ! 解像度設定の変更. eee_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
      mm = params(id)%mm
      nm = params(id)%nm
      ls => params(id)%ls
      le => params(id)%le
      lc => params(id)%lc
      js => params(id)%js
      je => params(id)%je
      jc => params(id)%jc
      itk => params(id)%itk
      tk  => params(id)%tk
      itj => params(id)%itj
      tj  => params(id)%tj
      iti => params(id)%iti
      ti  => params(id)%ti
      x_X => params(id)%x_X
      v_Y => params(id)%v_Y
      z_Z => params(id)%z_Z
      x_X_Weight => params(id)%x_X_Weight
      v_Y_Weight => params(id)%v_Y_Weight
      z_Z_Weight => params(id)%z_Z_Weight
      zxv_X => params(id)%zxv_X
      zxv_Y => params(id)%zxv_Y
      zxv_Z => params(id)%zxv_Z
      w => params(id)%w
      ws => params(id)%ws
      wg => params(id)%wg
      sgwork => params(id)%sgwork

    end subroutine eee_ChangeResolution

  !--------------- 波数位置問い合わせ -----------------
    function lf_l(l)
      !
      ! スペクトルデータ eef(-nm:nm,-mm,mm,2*lc) での
      ! X 方向波数 L の位置情報(第 3 添え字)を調べる.
      !
      ! X 方向波数については ls,ls+1,... le, -le,-le+1,...,-ls
      ! の順に格納されている(ISPACK/P3PACK マニュアル参照のこと)
      !
      integer, intent(IN)  ::  l        ! X 方向波数
      integer              ::  lf_l     ! スペクトルデータ第 3 添え字

      if ( abs(l) < ls(ip) .OR. le(ip) < abs(l) ) then
         call MessageNotify('E','l_l','Input wavenumber out of range')
      endif
      
      if ( l >= 0 ) then
         lf_l = l-ls(ip)+1
      else
         lf_l = ls(ip)+2*lc(ip)-abs(l)
      endif

    end function lf_l

  !--------------- 基本変換 -----------------
    function zxv_eef(eef)
      !
      ! スペクトルデータから格子データへ変換する.
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip))        :: zxv_eef
                                                        !(out) 格子点データ
      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip)), intent(in) :: eef
                                                        !(in)  スペクトルデータ
      sgwork(1:(2*nm+1)*(2*mm+1)*2*lc(ip)) &
           = reshape(eef,(/(2*nm+1)*(2*mm+1)*2*lc(ip)/))

      call p3smgb(nm,mm,lm,km,jm,im,sgwork,w,itk,tk,itj,tj,iti,ti)

      zxv_eef = reshape(sgwork(1:im*km*jc(ip)),(/km,im,jc(ip)/))
    end function zxv_eef

    function eef_zxv(zxv)
      !
      ! 格子データからスペクトルデータへ変換する.
      !
      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip))                  :: eef_zxv
                                                      !(out)  スペクトルデータ
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(in) :: zxv
                                                      !(in) 格子点データ

      sgwork(1:im*jc(ip)*km) = reshape(zxv,(/im*jc(ip)*km/))

      call p3gmsb(nm,mm,lm,km,jm,im,sgwork,w,itk,tk,itj,tj,iti,ti)

      eef_zxv = reshape(sgwork(1:(2*nm+1)*(2*mm+1)*2*lc(ip)),&
                               (/2*nm+1,2*mm+1,2*lc(ip)/))

    end function eef_zxv

    function eee2_ee2f(ee2f)
      !
      ! 分散スペクトルデータを集積する
      !
      real(8), dimension(-nm:nm,-mm:mm,-lm:lm,2)             :: eee2_ee2f
      !(out)  スペクトルデータ

      real(8), dimension(-nm:nm,-mm:mm,2,2*lc(ip)),intent(IN) :: ee2f
      !(in) 分散スペクトルデータ
      
      real(8), dimension(-nm:nm,-mm:mm,2,0:2*sum(lc)-1)       :: ee2e
      ! 作業データ
      integer, dimension(0:np-1) :: nm2ls
      integer :: ierr, iip
      integer :: l

      nm2ls(0)   = 0
      do iip=1,np-1
         nm2ls(iip) = nm2ls(iip-1) + (2*nm+1)*(2*mm+1)*2*2*lc(iip-1)
      enddo

      call MPI_ALLGATHERV(ee2f,(2*nm+1)*(2*mm+1)*2*2*lc(ip),MPI_REAL8,&
                          ee2e,(2*nm+1)*(2*mm+1)*2*2*lc,nm2ls,&
                          MPI_REAL8,MPI_COMM_WORLD,IERR)

      do l=ls(0),le(0)
         eee2_ee2f(:,:,l,:) = ee2e(:,:,:,l-ls(0))
         if ( l /= 0 ) then
            eee2_ee2f(:,:,-l,:) = ee2e(:,:,:,2*lc(0)-(l-ls(0)+1))
         endif
      enddo

      do iip=1,np-1
         do l=ls(iip),le(iip)
            eee2_ee2f(:,:,l,:) = ee2e(:,:,:,2*sum(lc(0:iip-1))+(l-ls(iip)))
            eee2_ee2f(:,:,-l,:) = ee2e(:,:,:,2*sum(lc(0:iip))-(l-ls(iip)+1))
         enddo
      enddo

    end function eee2_ee2f

    function ee2f_eee2(eee2)
      !
      ! スペクトルデータを分散する
      !
      real(8), dimension(-nm:nm,-mm:mm,2,2*lc(ip))       :: ee2f_eee2
      !(out)  スペクトルデータ

      real(8), dimension(-nm:nm,-mm:mm,-lm:lm,2), intent(IN) :: eee2
      !(in) 分散スペクトルデータ
      
      ! 作業データ
      integer :: l 

      do l=ls(ip),le(ip)
         ee2f_eee2(:,:,:,l-ls(ip)+1) = eee2(:,:,l,:)
         ee2f_eee2(:,:,:,2*lc(ip)-(l-ls(ip))) = eee2(:,:,-l,:)
      enddo
      
    end function ee2f_eee2

    function xyz_zxv(zxv)
      !
      ! 分散格子データの集積
      !
      real(8), dimension(0:im-1,0:jm-1,0:km-1)             :: xyz_zxv
                                                        !(out) 格子点データ

      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(in)  :: zxv
                                                        !(in) 格子点データ
      real(8), dimension(0:km-1,0:im-1,0:jm-1)  :: zxy
      ! 作業変数

      integer :: i,j,k
      integer :: ierr

      call MPI_ALLGATHERV(zxv,km*im*jc(ip),MPI_REAL8,&
                          zxy,km*im*jc,km*im*js,MPI_REAL8,&
                          MPI_COMM_WORLD,IERR)

      do k=0,km-1
         do j=0,jm-1
            do i=0,im-1
               xyz_zxv(i,j,k) = zxy(k,i,j)
            enddo
         enddo
      enddo

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

      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip)), intent(in)  :: eef
      !(in) 入力スペクトルデータ

      integer l,m,n
      ! 作業変数

      eef_Lapla_eef = 0.0

      do l=ls(ip),le(ip)
         do m=-mm,mm
            do n=-nm,nm
               eef_Lapla_eef(n,m,lf_l(l)) = -(l**2+m**2+n**2)*eef(n,m,lf_l(l))
               eef_Lapla_eef(n,m,lf_l(-l)) = -(l**2+m**2+n**2)*eef(n,m,lf_l(-l))
            enddo
         enddo
      enddo

    end function eef_Lapla_eef

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

      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip)), intent(in) :: eef
      !(in) スペクトルデータ

      integer l,m,n

      eef_LaplaInv_eef = 0.0
      do l=ls(ip),le(ip)
         do m=-mm,mm
            do n=-nm,nm
               if ( l.ne.0 .or. m.ne.0 .or. n.ne.0 ) then
                  eef_LaplaInv_eef(n,m,lf_l(l)) = -eef(n,m,lf_l(l))/(l**2+m**2+n**2)
                  eef_LaplaInv_eef(n,m,lf_l(-l)) = -eef(n,m,lf_l(-l))/(l**2+m**2+n**2)
               endif
            enddo
         enddo
      enddo
    end function eef_LaplaInv_eef

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

      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip)), intent(in)  :: eef
      !(in) 入力スペクトルデータ

      integer l,m,n
      ! 作業変数

      eef_Dx_eef = 0.0
      do l=ls(ip),le(ip)
         do m=-mm,mm
            do n=-nm,nm
               eef_Dx_eef(n,m,lf_l(l)) = -l*eef(-n,-m,lf_l(-l))
            enddo
         enddo
      enddo

      do l=-le(ip),-ls(ip)
         do m=-mm,mm
            do n=-nm,nm
               eef_Dx_eef(n,m,lf_l(l)) = -l*eef(-n,-m,lf_l(-l))
            enddo
         enddo
      enddo
    end function eef_Dx_eef

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

      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip)), intent(in)  :: eef
      !(in) 入力スペクトルデータ

      integer l,m,n
      ! 作業変数

      eef_Dy_eef = 0.0
      do l=ls(ip),le(ip)
         do m=-mm,mm
            do n=-nm,nm
               eef_Dy_eef(n,m,lf_l(l)) = -m*eef(-n,-m,lf_l(-l))
            enddo
         enddo
      enddo

      do l=-le(ip),-ls(ip)
         do m=-mm,mm
            do n=-nm,nm
               eef_Dy_eef(n,m,lf_l(l)) = -m*eef(-n,-m,lf_l(-l))
            enddo
         enddo
      enddo

    end function eef_Dy_eef

    function eef_Dz_eef(eef)
      !
      ! 入力スペクトルデータに Z 微分(∂z)を作用する.
      !
      ! スペクトルデータの Z 微分とは, 対応する格子点データに Z 微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      ! 実際にはスペクトルデータに Z 方向波数 n をかけて
      ! 実部 <-> 虚部成分に入れ換える計算を行っている.
      !
      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip))              :: eef_Dz_eef
      !(out) スペクトルデータの X 微分

      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip)), intent(in)  :: eef
      !(in) 入力スペクトルデータ

      integer l,m,n
      ! 作業変数

      eef_Dz_eef = 0.0
      do l=ls(ip),le(ip)
         do m=-mm,mm
            do n=-nm,nm
               eef_Dz_eef(n,m,lf_l(l)) = -n*eef(-n,-m,lf_l(-l))
            enddo
         enddo
      enddo

      do l=-le(ip),-ls(ip)
         do m=-mm,mm
            do n=-nm,nm
               eef_Dz_eef(n,m,lf_l(l)) = -n*eef(-n,-m,lf_l(-l))
            enddo
         enddo
      enddo

    end function eef_Dz_eef

  !--------------- 非線形項計算 -----------------

    function ee2f_RotVelxVor_ee2f(ee2f)
      !
      !  渦度 2 成分(ζ_1, ζ_2)のスペクトルデータから Euler 方程式の非線形項
      !
      !     ▽x(u x ω) 
      !
      !  の 2 成分を計算する.
      !
      !  (ζ_1, ζ_2) と渦度 ω との関係は ISPACK/P3PACK のマニュアルを参照
      !
      real(8), dimension(-nm:nm,-mm:mm,2,2*lc(ip))     :: ee2f_RotVelxVor_ee2f
      !(out) 非線形項のスペクトルデータの 2 つの成分

      real(8), dimension(-nm:nm,-mm:mm,2,2*lc(ip)), intent(in)  :: ee2f
      !(in) 入力スペクトルデータ. 渦度の 2 成分(ζ_1, ζ_2)

      ee2f_RotVelxVor_ee2f = 0.0D0
      call p3emnl(nm,mm,lm,km,jm,im,ee2f,ee2f_RotVelxVor_ee2f, &
                  wg,itk,tk,itj,tj,iti,ti)

    end function ee2f_RotVelxVor_ee2f

    function eef_VorFromZeta_ee2f(ee2f,isw)
      !
      !  渦度 2 成分(ζ_1, ζ_2)のスペクトルデータから渦度のスペクトルの 1 成分
      !  を計算する.
      !
      !  (ζ_1, ζ_2) と渦度 ω との関係は ISPACK/P3PACK のマニュアルを参照
      !
      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip))     :: eef_VorFromZeta_ee2f
      !(out) 非線形項のスペクトルデータの 2 つの成分

      real(8), dimension(-nm:nm,-mm:mm,2,2*lc(ip)), intent(in)  :: ee2f
      !(in) 入力スペクトルデータ. 渦度の 2 成分(ζ_1, ζ_2)

      integer, intent(IN) :: isw
      !(in) 出力する渦度の成分のインデックス(1,2,3)

      eef_VorFromZeta_ee2f = 0.0D0
      call p3gmto(nm,mm,lm,ee2f,eef_VorFromZeta_ee2f,isw)

    end function eef_VorFromZeta_ee2f

    function eef_VelFromZeta_ee2f(ee2f,isw)
      !
      !  渦度 2 成分(ζ_1, ζ_2)のスペクトルデータから速度スペクトルの 1 成分
      !  を計算する.
      !
      !  (ζ_1, ζ_2) と渦度 ω との関係は ISPACK/P3PACK のマニュアルを参照
      !
      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip))     :: eef_VelFromZeta_ee2f
      !(out) 非線形項のスペクトルデータの 2 つの成分

      real(8), dimension(-nm:nm,-mm:mm,2,2*lc(ip)), intent(in)  :: ee2f
      !(in) 入力スペクトルデータ. 渦度の 2 成分(ζ_1, ζ_2)

      integer, intent(IN) :: isw
      !(in) 出力する渦度の成分のインデックス(1,2,3)

      eef_VelFromZeta_ee2f = 0.0D0
      call p3gmtu(nm,mm,lm,ee2f,eef_VelFromZeta_ee2f,isw)

    end function eef_VelFromZeta_ee2f

    function ee2f_ZetaFromVor_eef_eef_eef(eef_1,eef_2,eef_3)
      !
      !  渦度のスペクトルから渦度 2 成分(ζ_1, ζ_2)を計算する.
      !
      !  (ζ_1, ζ_2) と渦度 ω との関係は ISPACK/P3PACK のマニュアルを参照
      !
      real(8), dimension(-nm:nm,-mm:mm,2,2*lc(ip))     :: ee2f_ZetaFromVor_eef_eef_eef
      !(out) 渦度の 2 成分(ζ_1, ζ_2)

      real(8), dimension(-nm:nm,-mm:mm,2*lc(ip)), intent(in)  :: eef_1, eef_2, eef_3
      !(in) 渦度スペクトルデータの各成分

      integer :: l,m,n

      ee2f_ZetaFromVor_eef_eef_eef = 0.0D0
      do l=ls(ip),le(ip)
         do m=-mm,mm
            do n=-nm,nm
               if ( l /= 0 ) then
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,1,lf_l(l)) = eef_2(n,m,lf_l(l))
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,2,lf_l(l)) = eef_3(n,m,lf_l(l))
               elseif( m /= 0 ) then
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,1,lf_l(l)) = eef_3(n,m,lf_l(l))
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,2,lf_l(l)) = eef_1(n,m,lf_l(l))
               else
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,1,lf_l(l)) = eef_1(n,m,lf_l(l))
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,2,lf_l(l)) = eef_2(n,m,lf_l(l))
               endif
            enddo
         enddo
      enddo

      do l=-le(ip),-ls(ip)
         do m=-mm,mm
            do n=-nm,nm
               if ( l /= 0 ) then
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,1,lf_l(l)) = eef_2(n,m,lf_l(l))
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,2,lf_l(l)) = eef_3(n,m,lf_l(l))
               elseif( m /= 0 ) then
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,1,lf_l(l)) = eef_3(n,m,lf_l(l))
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,2,lf_l(l)) = eef_1(n,m,lf_l(l))
               else
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,1,lf_l(l)) = eef_1(n,m,lf_l(l))
                  ee2f_ZetaFromVor_eef_eef_eef(n,m,2,lf_l(l)) = eef_2(n,m,lf_l(l))
               endif
            enddo
         enddo
      enddo

    end function ee2f_ZetaFromVor_eef_eef_eef

  !--------------- 積分計算 -----------------
    function IntZXV_zxv(zxv)
      !
      ! 3 次元格子点データの全領域積分
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight, z_Z_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(IN)   :: zxv
      !(in)  3 次元格子点データ

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

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

      IntZXV_zxv = 0.0d0
      do j=js(ip),je(ip)
         do i=0,im-1
            do k=0,km-1
               IntZXV_zxv = IntZXV_zxv &
                    + zxv(k,i,j) * z_Z_Weight(k) * v_Y_Weight(j) * x_X_Weight(i)
            enddo
         enddo
      end do

      IntZXVTMP=IntZXV_zxv
      CALL MPI_ALLREDUCE(IntZXVTMP,IntZXV_zxv,1,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)

    end function IntZXV_zxv

    function z_IntXV_zxv(zxv)
      !
      ! 3 次元格子点データの X,Y 方向積分
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(IN)   :: zxv
      !(in) 2 次元格子点データ

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

      real(8), dimension(0:km-1)          :: z_IntXVTmp
      integer :: i, j
      integer :: ierr
      ! 作業変数

      z_IntXV_zxv = 0.0d0
      do i=0,im-1
         do j=js(ip),je(ip)
            z_IntXV_zxv(:) = &
                 z_IntXV_zxv(:) + zxv(:,i,j) * x_X_Weight(i)* v_Y_Weight(j)
         enddo
      enddo

      z_IntXVTMP=z_IntXV_zxv
      CALL MPI_ALLREDUCE(z_IntXVTMP,z_IntXV_zxv,km,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)

    end function z_IntXV_zxv

    function v_IntZX_zxv(zxv)
      !
      ! 3 次元格子点データの X,Z 方向積分
      !
      ! 実際には格子点データ各点毎に x_X_Weight, z_Z_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(IN)   :: zxv
      !(in) 2 次元格子点データ

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

      integer :: i, k
      ! 作業変数

      v_IntZX_zxv = 0.0d0
      do i=0,im-1
         do k=0,km-1
            v_IntZX_zxv(:) = &
                 v_IntZX_zxv(:) + zxv(k,i,:) * x_X_Weight(i)* z_Z_Weight(k)
         enddo
      enddo
    end function v_IntZX_zxv

    function x_IntZV_zxv(zxv)
      !
      ! 3 次元格子点データの Y,Z 方向積分
      !
      ! 実際には格子点データ各点毎に y_Y_Weight, z_Z_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(IN)   :: zxv
      !(in) 2 次元格子点データ

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

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

      x_IntZV_zxv = 0.0d0
      do j=js(ip),je(ip)
         do k=0,km-1
            x_IntZV_zxv(:) = &
                 x_IntZV_zxv(:) + zxv(k,:,j) * v_Y_Weight(j)* z_Z_Weight(k)
         enddo
      enddo

      x_IntZVTMP=x_IntZV_zxv
      CALL MPI_ALLREDUCE(x_IntZVTMP,x_IntZV_zxv,im,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)

    end function x_IntZV_zxv

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

      real(8), dimension(0:km-1,js(ip):je(ip))          :: zv_IntX_zxv
      !(out) 積分された 2 次元(ZY)格子点データ

      integer :: i
      ! 作業変数

      zv_IntX_zxv = 0.0d0
      do i=0,im-1
         zv_IntX_zxv(:,:) = zv_IntX_zxv(:,:) + zxv(:,i,:) * x_X_Weight(i)
      enddo

    end function zv_IntX_zxv

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

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

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

      zx_IntV_zxv = 0.0d0
      do j=js(ip),je(ip)
         zx_IntV_zxv(:,:) = zx_IntV_zxv(:,:) + zxv(:,:,j) * v_Y_Weight(j)
      enddo

      zx_IntVTMP=zx_IntV_zxv
      CALL MPI_ALLREDUCE(zx_IntVTMP,zx_IntV_zxv,im*km,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)
    end function zx_IntV_zxv

    function xv_IntZ_zxv(zxv)
      !
      ! 3 次元格子点データの Z 方向積分
      !
      ! 実際には格子点データ各点毎に z_Z_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(IN)   :: zxv
      !(in) 2 次元格子点データ

      real(8), dimension(0:im-1,js(ip):je(ip))          :: xv_IntZ_zxv
      !(out) 積分された 1 次元(YX)格子点データ

      integer :: k 
      ! 作業変数

      xv_IntZ_zxv = 0.0d0
      do k=0,km-1
         xv_IntZ_zxv(:,:) = xv_IntZ_zxv(:,:) + zxv(k,:,:) * z_Z_Weight(k)
      enddo

    end function xv_IntZ_zxv

    function IntXV_xv(xv)
      !
      ! 2 次元(XY)格子点データの全領域積分
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:im-1,js(ip):je(ip)), intent(IN)   :: xv          
      !(in)  2 次元(YX)格子点データ

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

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

      IntXV_xv = 0.0d0
      do j=js(ip),je(ip)
         do i=0,im-1
            IntXV_xv = IntXV_xv + xv(i,j) * v_Y_Weight(j) * x_X_Weight(i)
         enddo
      enddo

      IntXVTMP=IntXV_xv
      CALL MPI_ALLREDUCE(IntXVTMP,IntXV_xv,1,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)

    end function IntXV_xv

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

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

      integer :: i
      ! 作業変数

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

    function x_IntV_xv(xv)
      !
      ! 2 次元(YX)格子点データの Y 方向積分
      !
      ! 実際には格子点データ各点毎に v_Y_Weight をかけた総和を計算している. 
      !
      real(8), dimension(0:im-1,js(ip):je(ip)), intent(IN)   :: xv      
      !(in)  2 次元(XV)格子点データ

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

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

      x_IntV_xv = 0.0d0
      do j=js(ip),je(ip)
         x_IntV_xv(:) = x_IntV_xv(:) + xv(:,j) * v_Y_Weight(j)
      enddo

      x_IntVTMP=x_IntV_xv
      CALL MPI_ALLREDUCE(x_IntVTMP,x_IntV_xv,im,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)
    end function x_IntV_xv

    function IntZV_zv(zv)
      !
      ! 2 次元(ZY)格子点データの全領域積分および平均.
      !
      ! 実際には格子点データ各点毎に z_Z_Weight, y_Y_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:km-1,js(ip):je(ip)), intent(IN)   :: zv          
      !(in)  2 次元(ZY)格子点データ

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

      real(8) :: IntZVTMP
      integer :: j, k
      integer :: ierr
      ! 作業変数

      IntZV_zv = 0.0d0
      do j=js(ip),je(ip)
         do k=0,km-1
            IntZV_zv = IntZV_zv + zv(k,j) * v_Y_Weight(j) * z_Z_Weight(k)
         enddo
      enddo

      IntZVTMP=IntZV_zv
      CALL MPI_ALLREDUCE(IntZVTMP,IntZV_zv,1,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)

    end function IntZV_zv

    function v_IntZ_zv(zv)
      !
      ! 2 次元(ZY)格子点データの Z 方向積分
      !
      ! 実際には格子点データ各点毎に z_Z_Weight をかけた総和を計算している. 
      !
      real(8), dimension(0:km-1,js(ip):je(ip)), intent(IN)   :: zv          
      !(in)  2 次元(ZY)格子点データ

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

      integer :: k
      ! 作業変数

      v_IntZ_zv = 0.0d0
      do k=0,km-1
         v_IntZ_zv(:) = v_IntZ_zv(:) + zv(k,:) * z_Z_Weight(k)
      enddo
    end function v_IntZ_zv

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

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

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

      z_IntV_zv = 0.0d0
      do j=js(ip),je(ip)
         z_IntV_zv(:) = z_IntV_zv(:) + zv(:,j) * v_Y_Weight(j)
      enddo

      z_IntVTMP=z_IntV_zv
      CALL MPI_ALLREDUCE(z_IntVTMP,z_IntV_zv,km,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)

    end function z_IntV_zv

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

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

      integer :: i, k
      ! 作業変数

      IntZX_zx = 0.0d0
      do i=0,im-1
         do k=0,km-1
            IntZX_zx = IntZX_zx + zx(k,i) * x_X_Weight(i) * z_Z_Weight(k)
         enddo
      enddo
    end function IntZX_zx

    function x_IntZ_zx(zx)
      !
      ! 2 次元(ZX)格子点データの Z 方向積分
      !
      ! 実際には格子点データ各点毎に z_Z_Weight をかけた総和を計算している. 
      !
      real(8), dimension(0:km-1,0:im-1), intent(IN)   :: zx          
      !(in)  2 次元(ZX)格子点データ

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

      integer :: k
      ! 作業変数

      x_IntZ_zx = 0.0d0
      do k=0,km-1
         x_IntZ_zx(:) = x_IntZ_zx(:) + zx(k,:) * z_Z_Weight(k)
      enddo
    end function x_IntZ_zx

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

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

      integer :: i
      ! 作業変数

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

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

      IntV_v = sum(v*v_Y_Weight)

      IntVTMP=IntV_v
      CALL MPI_ALLREDUCE(IntVTMP,IntV_v,1,MPI_REAL8, &
                         MPI_SUM,MPI_COMM_WORLD,IERR)
    end function IntV_v

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

      IntZ_z = sum(z*z_Z_Weight)
    end function IntZ_z

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

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

      real(8) :: Vol, VolTMP
      integer :: ierr
      ! 作業変数

      VolTMP=sum(x_X_weight)*sum(v_Y_weight)*sum(z_Z_weight)
      
      CALL MPI_ALLREDUCE(VolTMP,Vol,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      AvrZXV_zxv = IntZXV_zxv(zxv)/Vol

    end function AvrZXV_zxv

    function z_AvrXV_zxv(zxv)
      !
      ! 3 次元格子点データの X,Y 方向平均
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
      ! 総和を計算し, x_X_Weight*y_Y_Weight の総和で割ることで平均している. 
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(IN)   :: zxv
      !(in) 2 次元格子点データ

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

      real(8) :: Area, AreaTMP
      integer :: ierr
      ! 作業変数

      AreaTMP = sum(x_X_weight)*sum(v_Y_weight)
      CALL MPI_ALLREDUCE(AreaTMP,Area,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      z_AvrXV_zxv = z_IntXV_zxv(zxv)/Area

    end function z_AvrXV_zxv

    function v_AvrZX_zxv(zxv)
      !
      ! 3 次元格子点データの X,Z 方向平均
      !
      ! 実際には格子点データ各点毎に x_X_Weight, z_Z_Weight をかけた
      ! 総和を計算し, x_X_Weight*z_Z_Weight の総和で割ることで平均している. 
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(IN)   :: zxv
      !(in) 2 次元格子点データ

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

      v_AvrZX_zxv = v_IntZX_zxv(zxv)/(sum(x_X_weight)*sum(z_Z_weight))

    end function v_AvrZX_zxv

    function x_AvrZV_zxv(zxv)
      !
      ! 3 次元格子点データの Y,Z 方向平均
      !
      ! 実際には格子点データ各点毎に y_Y_Weight, z_Z_Weight をかけた
      ! 総和を計算し, y_Y_Weight*z_Z_Weight の総和で割ることで平均している. 
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(IN)   :: zxv
      !(in) 2 次元格子点データ

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

      real(8) :: Area, AreaTMP
      integer :: ierr
      ! 作業変数

      AreaTMP = sum(z_Z_weight)*sum(v_Y_weight)
      CALL MPI_ALLREDUCE(AreaTMP,Area,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      x_AvrZV_zxv = x_IntZV_zxv(zxv)/Area

    end function x_AvrZV_zxv

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

      real(8), dimension(0:km-1,js(ip):je(ip))          :: zv_AvrX_zxv
      !(out) 積分された 2 次元(ZY)格子点データ

      zv_AvrX_zxv = zv_IntX_zxv(zxv)/sum(x_X_weight)

    end function zv_AvrX_zxv

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

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

      real(8) :: Length, LengthTMP
      integer :: ierr
      ! 作業変数

      LengthTMP = sum(v_Y_weight)
      CALL MPI_ALLREDUCE(LengthTMP,Length,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      zx_AvrV_zxv = zx_IntV_zxv(zxv)/Length

    end function zx_AvrV_zxv

    function xv_AvrZ_zxv(zxv)
      !
      ! 3 次元格子点データの Z 方向平均
      !
      ! 実際には格子点データ各点毎に z_Z_Weight をかけた
      ! 総和を計算し, z_Z_Weight の総和で割ることで平均している. 
      !
      real(8), dimension(0:km-1,0:im-1,js(ip):je(ip)), intent(IN)   :: zxv
      !(in) 2 次元格子点データ

      real(8), dimension(0:im-1,js(ip):je(ip))          :: xv_AvrZ_zxv
      !(out) 積分された 1 次元(YX)格子点データ

      xv_AvrZ_zxv = xv_IntZ_zxv(zxv)/sum(z_Z_weight)

    end function xv_AvrZ_zxv

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

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

      real(8) :: Area, AreaTMP
      integer :: ierr
      ! 作業変数

      AreaTMP = sum(x_X_weight)*sum(v_Y_weight)
      CALL MPI_ALLREDUCE(AreaTMP,Area,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      AvrXV_xv = IntXV_xv(xv)/Area

    end function AvrXV_xv

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

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

      v_AvrX_xv = v_IntX_xv(xv)/sum(x_X_weight)
    end function v_AvrX_xv

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

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

      real(8) :: Length, LengthTMP
      integer :: ierr
      ! 作業変数

      LengthTMP = sum(v_Y_weight)
      CALL MPI_ALLREDUCE(LengthTMP,Length,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      x_AvrV_xv = x_IntV_xv(xv)/Length
    end function x_AvrV_xv

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

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

      AvrZX_zx = IntZX_zx(zx)/(sum(x_X_weight)*sum(z_Z_weight))
    end function AvrZX_zx

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

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

      z_AvrX_zx = z_IntX_zx(zx)/sum(x_X_weight)
    end function z_AvrX_zx

    function x_AvrZ_zx(zx)
      !
      ! 2 次元(ZX)格子点データの Z 方向平均
      !
      ! 実際には格子点データ各点毎に z_Z_Weight をかけた総和を計算し, 
      ! z_Z_Weight の総和で割ることで平均している. 
      !
      real(8), dimension(0:km-1,0:im-1), intent(IN)   :: zx
      !(in) 2 次元(ZX)格子点データ

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

      x_AvrZ_zx = x_IntZ_zx(zx)/sum(z_Z_weight)
    end function x_AvrZ_zx

    function AvrZV_zv(zv)
      !
      ! 2 次元(ZV)格子点データの全領域平均
      !
      ! 実際には格子点データ各点毎に z_Z_Weight, y_Y_Weight をかけた
      ! 総和を計算し, z_Z_Weight*v_Y_Weight の総和で割ることで平均している. 
      !
      real(8), dimension(0:km-1,js(ip):je(ip)), intent(IN)   :: zv
      !(in)  2 次元(ZY)格子点データ

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

      real(8) :: Area, AreaTMP
      integer :: ierr
      ! 作業変数

      AreaTMP = sum(z_Z_weight)*sum(v_Y_weight)
      CALL MPI_ALLREDUCE(AreaTMP,Area,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      AvrZV_zv = IntZV_zv(zv)/Area
    end function AvrZV_zv

    function z_AvrV_zv(zv)
      !
      ! 2 次元(ZY)格子点データの Y 方向平均
      !
      ! 実際には格子点データ各点毎に v_Y_Weight をかけた総和を計算し, 
      ! v_Y_Weight の総和で割ることで平均している. 
      !
      real(8), dimension(0:km-1,js(ip):je(ip)), intent(IN)   :: zv
      !(in) 2 次元(ZV)格子点データ

      real(8), dimension(0:km-1)          :: z_AvrV_zv
      !(out) 平均された 1 次元(Z)格子点データ
      real(8) :: Length, LengthTMP
      integer :: ierr
      ! 作業変数

      LengthTMP = sum(v_Y_weight)
      CALL MPI_ALLREDUCE(LengthTMP,Length,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      z_AvrV_zv = z_IntV_zv(zv)/Length
    end function z_AvrV_zv

    function v_AvrZ_zv(zv)
      !
      ! 2 次元(ZY)格子点データの Z 方向平均
      !
      ! 実際には格子点データ各点毎に z_Z_Weight をかけた総和を計算し, 
      ! z_Z_Weight の総和で割ることで平均している. 
      !
      real(8), dimension(0:km-1,js(ip):je(ip)), intent(IN)   :: zv
      !(in) 2 次元(ZY)格子点データ

      real(8), dimension(js(ip):je(ip))          :: v_AvrZ_zv
      !(out) 平均された 1 次元(X)格子点データ

      v_AvrZ_zv = v_IntZ_zv(zv)/sum(z_Z_weight)
    end function v_AvrZ_zv

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

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

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

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

      real(8) :: Length, LengthTMP
      integer :: ierr
      ! 作業変数

      LengthTMP = sum(v_Y_weight)
      CALL MPI_ALLREDUCE(LengthTMP,Length,1,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      AvrV_v = IntV_v(v)/Length
    end function AvrV_v

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

      AvrZ_z = IntZ_z(z)/sum(z_Z_weight)
    end function AvrZ_z

  !--------------- スペクトル計算 -----------------
    function EnergyHelicityFromZeta_ee2f(ee2f)
      !
      ! 渦度成分(ζ_1, ζ_2)から全領域平均エネルギーとヘリシティーを計算する.
      !
      real(8), dimension(2)                 :: EnergyHelicityFromZeta_ee2f
      ! エネルギーヘリシティー

      real(8), dimension(-nm:nm,-mm:mm,2,2*lc(ip)), intent(in) :: ee2f
      ! 渦度成分(ζ_1, ζ_2)

      real(8) :: E, H           ! エネルギー, ヘリシティー
      
      call p3cmsv(nm,mm,lm,ee2f,E,H)

      EnergyHelicityFromZeta_ee2f(1) = E ; EnergyHelicityFromZeta_ee2f(2) = H

    end function EnergyHelicityFromZeta_ee2f

    subroutine ESpectralFromZeta(esp,ee2f)
      !
      ! 渦度成分(ζ_1, ζ_2)からエネルギースペクトルを計算する. 
      !
      !   * esp のサイズで求められる波数範囲が定められる
      !   * esp の総和が EFFromZeta で求められるエネルギー
      !     (全領域平均値)に等しい
      !
      real(8), dimension(:), intent(OUT)  :: esp
      ! エネルギースペクトル

      real(8), dimension(-nm:nm,-mm:mm,2,2*lc(ip)), intent(in) :: ee2f
      ! 渦度成分(ζ_1, ζ_2)

      integer kmax
      ! 作業変数
      real(8), allocatable :: wesp(:)

      kmax=size(esp)
      allocate(wesp(kmax))
      call p3empt(nm,mm,lm,kmax,ee2f,esp,wesp)
      deallocate(wesp)

    end subroutine ESpectralFromZeta

  end module eee_mpi_module
