eee_module

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

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

      AvrYX_yx = IntYX_yx(yx)/(sum(x_X_weight)*sum(y_Y_weight))
    end function AvrYX_yx

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

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

    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 AvrZYX_zyx(zyx)
      !
      ! 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:jm-1,0:im-1), intent(IN)   :: zyx
      !(in)  3 次元格子点データ

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

      AvrZYX_zyx = IntZYX_zyx(zyx)/(sum(x_X_weight)*sum(y_Y_weight)*sum(z_Z_weight))

    end function AvrZYX_zyx

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

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

      AvrZY_zy = IntZY_zy(zy)/(sum(y_Y_weight)*sum(z_Z_weight))
    end function AvrZY_zy

    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

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

      real(8), dimension(-nm:nm,-mm:mm,-lm:lm,2), intent(in) :: eee2
      ! 渦度成分(ζ_1, ζ_2)

      integer kmax
      ! 作業変数

      kmax=size(esp)
      call p3espt(nm,mm,lm,kmax,eee2,esp)

    end subroutine ESpectralFromZeta

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

      real(8), dimension(-nm:nm,-mm:mm,-lm:lm,2), intent(in) :: eee2
      ! 渦度成分(ζ_1, ζ_2)

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

      EnergyHelicityFromZeta_eee2(1) = E 
       EnergyHelicityFromZeta_eee2(2) = H

    end function EnergyHelicityFromZeta_eee2

    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 IntYX_yx(yx)
      !
      ! 2 次元(YX)格子点データの全領域積分
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:jm-1,0:im-1), intent(IN)   :: yx          
      !(in)  2 次元(YX)格子点データ

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

      integer :: i, j
      ! 作業変数

      IntYX_yx = 0.0d0
      do i=0,im-1
         do j=0,jm-1
            IntYX_yx = IntYX_yx + yx(j,i) * y_Y_Weight(j) * x_X_Weight(i)
         enddo
      enddo
    end function IntYX_yx

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

      IntY_y = sum(y*y_Y_Weight)
    end function IntY_y

    function 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 IntZYX_zyx(zyx)
      !
      ! 2 次元格子点データの全領域積分
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight, z_Z_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:km-1,0:jm-1,0:im-1), intent(IN)   :: zyx
      !(in)  3 次元格子点データ

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

      integer :: i, j, k
      ! 作業変数

      IntZYX_zyx = 0.0d0
      do i=0,im-1
         do j=0,jm-1
            do k=0,km-1
               IntZYX_zyx = IntZYX_zyx + zyx(k,j,i) * z_Z_Weight(k) * y_Y_Weight(j) * x_X_Weight(i)
            enddo
         enddo
      end do
    end function IntZYX_zyx

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

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

      integer :: j, k
      ! 作業変数

      IntZY_zy = 0.0d0
      do j=0,jm-1
         do k=0,km-1
            IntZY_zy = IntZY_zy + zy(k,j) * y_Y_Weight(j) * z_Z_Weight(k)
         enddo
      enddo
    end function IntZY_zy

    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 eee2_RotVelxVor_eee2(eee2)
      !
      !  渦度 2 成分(ζ_1, ζ_2)のスペクトルデータから Euler 方程式の非線形項
      !
      !     ▽x(u x ω) 
      !
      !  の 2 成分を計算する.
      !
      !  (ζ_1, ζ_2) と渦度 ω との関係は ISPACK/P3PACK のマニュアルを参照
      !
      real(8), dimension(-nm:nm,-mm:mm,-lm:lm,2)     :: eee2_RotVelxVor_eee2
      !(out) 非線形項のスペクトルデータの 2 つの成分

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

      real(8), dimension((2*lm+1)*(2*mm+1)*(2*nm+1))      :: ws
      real(8), dimension(km*jm*im*4)                      :: wg

      call p3elnl(nm,mm,lm,km,jm,im,eee2,eee2_RotVelxVor_eee2, ws,wg,itk,tk,itj,tj,iti,ti)

    end function eee2_RotVelxVor_eee2

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

      real(8), dimension(-nm:nm,-mm:mm,-lm:lm), intent(in)  :: eee_1, eee_2, eee_3
      !(in) 渦度スペクトルデータの各成分

      integer :: l,m,n

      do l=-lm,lm
         do m=-mm,mm
            do n=-nm,nm
               if ( l /= 0 ) then
                  eee2_ZetaFromVor_eee_eee_eee(n,m,l,1) = eee_2(n,m,l)
                  eee2_ZetaFromVor_eee_eee_eee(n,m,l,2) = eee_3(n,m,l)
               elseif( m /= 0 ) then
                  eee2_ZetaFromVor_eee_eee_eee(n,m,l,1) = eee_3(n,m,l)
                  eee2_ZetaFromVor_eee_eee_eee(n,m,l,2) = eee_1(n,m,l)
               else
                  eee2_ZetaFromVor_eee_eee_eee(n,m,l,1) = eee_1(n,m,l)
                  eee2_ZetaFromVor_eee_eee_eee(n,m,l,2) = eee_2(n,m,l)
               endif
            enddo
         enddo
      enddo

    end function eee2_ZetaFromVor_eee_eee_eee

    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
      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
      y_Y => params(id)%y_Y
      z_Z => params(id)%z_Z
      x_X_Weight => params(id)%x_X_Weight
      y_Y_Weight => params(id)%y_Y_Weight
      z_Z_Weight => params(id)%z_Z_Weight
      zyx_X => params(id)%zyx_X
      zyx_Y => params(id)%zyx_Y
      zyx_Z => params(id)%zyx_Z

    end subroutine eee_ChangeResolution

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

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

      integer l,m,n
      ! 作業変数


      do l=-lm,lm
         do m=-mm,mm
            do n=-nm,nm
               eee_Dx_eee(n,m,l) = -l*eee(-n,-m,-l)
            enddo
         enddo
      enddo
    end function eee_Dx_eee

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

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

      integer l,m,n
      ! 作業変数


      do l=-lm,lm
         do m=-mm,mm
            do n=-nm,nm
               eee_Dy_eee(n,m,l) = -m*eee(-n,-m,-l)
            enddo
         enddo
      enddo
    end function eee_Dy_eee

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

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

      integer l,m,n
      ! 作業変数


      do l=-lm,lm
         do m=-mm,mm
            do n=-nm,nm
               eee_Dz_eee(n,m,l) = -n*eee(-n,-m,-l)
            enddo
         enddo
      enddo
    end function eee_Dz_eee

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

      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

      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)%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)%x_X(0:im-1))
      allocate(params(nparams)%x_X_Weight(0:im-1))
      allocate(params(nparams)%y_Y(0:jm-1))
      allocate(params(nparams)%y_Y_Weight(0:jm-1))
      allocate(params(nparams)%z_Z(0:km-1))
      allocate(params(nparams)%z_Z_Weight(0:km-1))
      allocate(params(nparams)%zyx_X(0:km-1,0:jm-1,0:im-1))
      allocate(params(nparams)%zyx_Y(0:km-1,0:jm-1,0:im-1))
      allocate(params(nparams)%zyx_Z(0:km-1,0:jm-1,0:im-1))

      call eee_ChangeResolution(nparams)

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

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

      do jj=0,jm-1
         y_Y(jj) = 2*pi/jm*jj
      enddo
      y_Y_Weight = 2*pi/jm

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

      zyx_X = spread(spread(x_X,1,jm),1,km)
      zyx_Y = spread(spread(y_Y,2,im),1,km)
      zyx_Z = spread(spread(z_Z,2,jm),3,im)

      if ( present(id) ) id = nparams

      write(cid,'(I3)') nparams
      call MessageNotify('M','eee_initial','eee_module (2009/07/31) is initialized')
      call MessageNotify('M','eee_initial', 'Resolution ID is '//trim(adjustl(cid)))
    end subroutine eee_initial

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

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

      integer l,m,n

      do l=-lm,lm
         do m=-mm,mm
            do n=-nm,nm
               if ( l.ne.0 .or. m.ne.0 .or. n.ne.0 ) then
                  eee_LaplaInv_eee(n,m,l) = -eee(n,m,l)/(l**2+m**2+n**2)
               else
                  eee_LaplaInv_eee(n,m,l) = 0.0
               endif
            enddo
         enddo
      enddo
    end function eee_LaplaInv_eee

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

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

      integer l,m,n
      ! 作業変数

         do l=-lm,lm
            do m=-mm,mm
               do n=-nm,nm
                  eee_Lapla_eee(n,m,l) = -(l**2+m**2+n**2)*eee(n,m,l)
               enddo
            enddo
         enddo
    end function eee_Lapla_eee

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

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

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

      call p3getu(nm,mm,lm,eee2,eee_VelFromZeta_eee2,isw)

    end function eee_VelFromZeta_eee2

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

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

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

      call p3geto(nm,mm,lm,eee2,eee_VorFromZeta_eee2,isw)

    end function eee_VorFromZeta_eee2

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

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

      real(8), dimension(km*jm*im)                         :: w
      real(8), dimension(0:km-1,0:jm-1,0:im-1)             :: zyx_tmp
      ! 作業領域

      zyx_tmp = zyx
      call p3g2sa(nm,mm,lm,km,jm,im,zyx_tmp,eee_zyx,w,itk,tk,itj,tj,iti,ti)
    end function eee_zyx

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

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

      x_AvrY_yx = x_IntY_yx(yx)/sum(y_Y_weight)
    end function x_AvrY_yx

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

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

      x_AvrZY_zyx = x_IntZY_zyx(zyx)/(sum(y_Y_weight)*sum(z_Z_weight))

    end function x_AvrZY_zyx

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

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

      integer :: j
      ! 作業変数

      x_IntY_yx = 0.0d0
      do j=0,jm-1
         x_IntY_yx(:) = x_IntY_yx(:) + yx(j,:) * y_Y_Weight(j)
      enddo
    end function x_IntY_yx

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

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

      integer :: j, k
      ! 作業変数

      x_IntZY_zyx = 0.0d0
      do j=0,jm-1
         do k=0,km-1
            x_IntZY_zyx(:) = x_IntZY_zyx(:) + zyx(k,j,:) * y_Y_Weight(j)* z_Z_Weight(k)
         enddo
      enddo
    end function x_IntZY_zyx

    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 xyz_zyx(zyx)
      !
      ! 格子データの添え字順序変更
      !
      real(8), dimension(0:im-1,0:jm-1,0:km-1)             :: xyz_zyx
                                                        !(out) 格子点データ

      real(8), dimension(0:km-1,0:jm-1,0:im-1), intent(in)  :: zyx
                                                      !(in) 格子点データ
      integer :: i,j,k

      do k=0,km-1
         do j=0,jm-1
            do i=0,im-1
               xyz_zyx(i,j,k) = zyx(k,j,i)
            enddo
         enddo
      enddo

    end function xyz_zyx

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

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

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

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

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

      y_AvrZX_zyx = y_IntZX_zyx(zyx)/(sum(x_X_weight)*sum(z_Z_weight))

    end function y_AvrZX_zyx

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

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

      y_AvrZ_zy = y_IntZ_zy(zy)/sum(z_Z_weight)
    end function y_AvrZ_zy

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

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

      integer :: i
      ! 作業変数

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

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

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

      integer :: i, k
      ! 作業変数

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

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

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

      integer :: k
      ! 作業変数

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

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

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

      yx_AvrZ_zyx = yx_IntZ_zyx(zyx)/sum(z_Z_weight)

    end function yx_AvrZ_zyx

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

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

      integer :: k 
      ! 作業変数

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

    end function yx_IntZ_zyx

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

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

      z_AvrYX_zyx = z_IntYX_zyx(zyx)/(sum(x_X_weight)*sum(y_Y_weight))

    end function z_AvrYX_zyx

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

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

      z_AvrY_zy = z_IntY_zy(zy)/sum(y_Y_weight)
    end function z_AvrY_zy

    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 z_IntYX_zyx(zyx)
      !
      ! 3 次元格子点データの X,Y 方向積分
      !
      ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
      ! 総和を計算している. 
      !
      real(8), dimension(0:km-1,0:jm-1,0:im-1), intent(IN)   :: zyx
      !(in) 2 次元格子点データ

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

      integer :: i, j
      ! 作業変数

      z_IntYX_zyx = 0.0d0
      do i=0,im-1
         do j=0,jm-1
            z_IntYX_zyx(:) = z_IntYX_zyx(:) + zyx(:,j,i) * x_X_Weight(i)* y_Y_Weight(j)
         enddo
      enddo
    end function z_IntYX_zyx

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

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

      integer :: j
      ! 作業変数

      z_IntY_zy = 0.0d0
      do j=0,jm-1
         z_IntY_zy(:) = z_IntY_zy(:) + zy(:,j) * y_Y_Weight(j)
      enddo
    end function z_IntY_zy

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

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

      zx_AvrY_zyx = zx_IntY_zyx(zyx)/sum(y_Y_weight)

    end function zx_AvrY_zyx

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

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

      integer :: j
      ! 作業変数

      zx_IntY_zyx = 0.0d0
      do j=0,jm-1
         zx_IntY_zyx(:,:) = zx_IntY_zyx(:,:) + zyx(:,j,:) * y_Y_Weight(j)
      enddo
    end function zx_IntY_zyx

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

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

      zy_AvrX_zyx = zy_IntX_zyx(zyx)/sum(x_X_weight)

    end function zy_AvrX_zyx

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

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

      integer :: i
      ! 作業変数

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

    end function zy_IntX_zyx

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

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

      real(8), dimension(km*jm*im)                         :: w
      ! 作業領域

      call p3s2ga(nm,mm,lm,km,jm,im,eee,zyx_eee,w,itk,tk,itj,tj,iti,ti)
    end function zyx_eee