eea_module

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

      integer :: n
      ! 作業変数

      do n=1,size(xa,2)
         a_AvrX_xa(n) = AvrX_x(xa(:,n))
      end do

    end function a_AvrX_xa

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

      real(8), dimension(size(yxa,3)) :: a_AvrYX_yxa
      !(out) 平均値

      integer :: n
      ! 作業変数

      do n=1,size(yxa,3)
         a_AvrYX_yxa(n) = AvrYX_yx(yxa(:,:,n))
      end do
    end function a_AvrYX_yxa

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

      integer :: n
      ! 作業変数

      do n=1,size(ya,2)
         a_AvrY_ya(n) = AvrY_y(ya(:,n))
      end do

    end function a_AvrY_ya

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

      real(8), dimension(size(xa,2)) :: a_IntX_xa     !(out) 積分値
      !(out) 積分された 1 次元格子点データ

      integer :: n
      ! 作業変数

      do n=1,size(xa,2)
         a_IntX_xa(n) = IntX_x(xa(:,n))
      enddo

    end function a_IntX_xa

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

      real(8), dimension(size(yxa,3)) :: a_IntYX_yxa
      !(out) 積分値

      integer :: n
      ! 作業変数

      do n=1,size(yxa,3)
         a_IntYX_yxa(n) = IntYX_yx(yxa(:,:,n))
      enddo
    end function a_IntYX_yxa

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

      real(8), dimension(size(ya,2)) :: a_IntY_ya   !(out) 積分値
      !(out) 積分された 1 次元格子点データ

      integer :: n
      ! 作業変数

      do n=1,size(ya,2)
         a_IntY_ya(n) = IntY_y(ya(:,n))
      enddo

    end function a_IntY_ya

    function a_Interpolate_eea( eea_Data, xa, ya )
      real(8), intent(IN)  :: eea_data(-lm:,-km:,:)  ! スペクトルデータ
      real(8), intent(IN)  :: xa(:)                   ! 補間する点の x 座標 
      real(8), intent(IN)  :: ya(:)                   ! 補間する点の y 座標 
      real(8)              :: a_Interpolate_eea(size(xa)) ! 補間した値

      integer :: n

      do n=1,size(xa)
         a_Interpolate_eea(n) = Interpolate_ee(eea_Data(:,:,n),xa(n),ya(n))
      end do

    end function a_Interpolate_eea

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

      real(8), dimension(-lm:lm,-km:km,size(eea,3)) :: eea_Dx_eea
      !(out) スペクトルデータの X 微分

      integer n
      ! 作業変数

      do n=1,size(eea,3)
         eea_Dx_eea(:,:,n) = ee_Dx_ee(eea(:,:,n))
      enddo
    end function eea_Dx_eea

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

      real(8), dimension(-lm:lm,-km:km,size(eea,3)) :: eea_Dy_eea
      !(out) スペクトルデータの Y 微分

      integer n
      ! 作業変数

      do n=1,size(eea,3)
         eea_Dy_eea(:,:,n) = ee_Dy_ee(eea(:,:,n))
      enddo

    end function eea_Dy_eea

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

      real(8), dimension(-lm:lm,-km:km,size(eea_StrFunc,3))  :: eea_EnergyFromStreamfunc_eea
      ! エネルギースペクトル

      integer n
      ! 作業変数

      do n=1,size(eea_StrFunc,3)
         eea_EnergyFromStreamfunc_eea(:,:,n) = ee_EnergyFromStreamfunc_ee(eea_StrFunc(:,:,n))
      enddo

    end function eea_EnergyFromStreamfunc_eea

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

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

      integer n
      ! 作業変数

      do n=1,size(eea_StrFunc,3)
         eea_EnstrophyFromStreamfunc_eea(:,:,n) = ee_EnstrophyFromStreamfunc_ee(eea_StrFunc(:,:,n))
      enddo

    end function eea_EnstrophyFromStreamfunc_eea

    subroutine eea_Initial(i,j,k,l,xmin,xmax,ymin,ymax)
      !
      ! スペクトル変換の格子点数, 波数, 領域の大きさを設定する.
      !
      ! 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
      ! 初期設定をしなければならない.
      !
      integer,intent(in) :: i           ! 格子点数(X)
      integer,intent(in) :: j           ! 格子点数(Y)
      integer,intent(in) :: K           ! 切断波数(X)
      integer,intent(in) :: l           ! 切断波数(Y)

      real(8),intent(in) :: xmin, xmax     ! X 座標範囲
      real(8),intent(in) :: ymin, ymax     ! Y 座標範囲

      im = i 
       jm = j 
       km = k 
       lm = l

      call ee_Initial(i,j,k,l,xmin,xmax,ymin,ymax)

      allocate(xy_X(0:im-1,0:jm-1),xy_Y(0:im-1,0:jm-1))
      xy_X = transpose(yx_X) 
        xy_Y = transpose(yx_Y)

      call MessageNotify('M','eea_initial','eae_module (2011/12/10) is initialized')
    end subroutine eea_Initial

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

      real(8), dimension(-lm:lm,-km:km,size(eea_Zeta,3)) :: eea_JacobianZ_eea
      !(out) 流線関数と渦度のヤコビアン

      integer n
      ! 作業変数

      do n=1,size(eea_Zeta,3)
         eea_JacobianZ_eea(:,:,n) = ee_JacobianZ_ee(eea_Zeta(:,:,n))
      end do

    end function eea_JacobianZ_eea

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

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

      real(8), dimension(-lm:lm,-km:km,size(eea_a,3)) :: eea_Jacobian_eea_eea
      !(out) 2 つのスペクトルデータのヤコビアン

      integer n
      ! 作業変数

      do n=1,size(eea_a,3)
         eea_Jacobian_eea_eea(:,:,n) = ee_Jacobian_ee_ee(eea_a(:,:,n),eea_b(:,:,n))
      end do

    end function eea_Jacobian_eea_eea

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

      real(8), dimension(-lm:lm,-km:km,size(eea,3)) :: eea_LaplaInv_eea
      !(out) スペクトルデータの逆ラプラシアン

      integer n

      do n=1,size(eea,3)
         eea_LaplaInv_eea(:,:,n) = ee_LaplaInv_ee(eea(:,:,n))
      enddo
    end function eea_LaplaInv_eea

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

      real(8), dimension(-lm:lm,-km:km,size(eea,3))  :: eea_Lapla_eea
      !(out) スペクトルデータのラプラシアン

      integer n
      ! 作業変数

      do n=1,size(eea,3)
         eea_Lapla_eea(:,:,n) = ee_Lapla_ee(eea(:,:,n))
      enddo
    end function eea_Lapla_eea

    function eea_yxa(yxa)
      !
      ! 格子データからスペクトルデータへ変換する.
      !
      real(8), dimension(0:,0:,:), intent(in)  :: yxa
      !(in) 格子点データ

      real(8), dimension(-lm:lm,-km:km,size(yxa,3))              :: eea_yxa
      !(out)  スペクトルデータ

      integer :: n
      ! DO LOOP 変数

      do n = 1, size(yxa,3)
         eea_yxa(:,:,n) = ee_yx(yxa(:,:,n))
      end do

    end function eea_yxa

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

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

      integer :: n
      ! 作業変数

      do n=1,size(yxa,3)
         xa_AvrY_yxa(:,n) = x_AvrY_yx(yxa(:,:,n))
      end do

    end function xa_AvrY_yxa

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

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

      integer :: n
      ! 作業変数

      do n=1,size(yxa,3)
         xa_IntY_yxa(:,n) = x_IntY_yx(yxa(:,:,n))
      enddo
    end function xa_IntY_yxa

    function xya_yxa(yxa)
      !
      ! 格子データの次元入換
      !
      real(8), dimension(0:,0:,:), intent(in)          :: yxa
      !(in) 格子点データ

      real(8), dimension(0:im-1,0:jm-1,size(yxa,3))    :: xya_yxa
      !(out)  スペクトルデータ

      integer :: n
      ! DO LOOP 変数

      do n = 1, size(yxa,3)
         xya_yxa(:,:,n) = transpose(yxa(:,:,n))
      end do

    end function xya_yxa

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

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

      integer :: n
      ! 作業変数

      do n=1,size(yxa,3)
         ya_AvrX_yxa(:,n) = y_AvrX_yx(yxa(:,:,n))
      end do
    end function ya_AvrX_yxa

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

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

      integer :: n
      ! 作業変数

      do n=1,size(yxa,3)
         ya_IntX_yxa(:,n) = y_IntX_yx(yxa(:,:,n))
      enddo
    end function ya_IntX_yxa

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

      real(8), dimension(0:jm-1,0:im-1,size(eea,3))             :: yxa_eea
      !(out) 格子点データ

      integer :: n
      ! DO LOOP 変数

      do n = 1, size(eea,3)
         yxa_eea(:,:,n) = yx_ee(eea(:,:,n))
      end do

    end function yxa_eea

    function yxa_xya(xya)
      !
      ! 格子データの次元入換
      !
      real(8), dimension(0:,0:,:), intent(in)        :: xya
      !(in) 格子点データ

      real(8), dimension(0:jm-1,0:im-1,size(xya,3))  :: yxa_xya
      !(out)  スペクトルデータ

      integer :: n
      ! DO LOOP 変数

      do n = 1, size(xya,3)
         yxa_xya(:,:,n) = transpose(xya(:,:,n))
      end do

    end function yxa_xya