function AvrYX_yx(yx)    ! 全領域平均
      real(8), dimension(0:jm-1,0:im-1)   :: yx          ! 2 次元格子点
      real(8)                             :: AvrYX_yx    ! 平均値

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

    function IntYX_yx(yx)   ! 全領域積分
      real(8), dimension(0:jm-1,0:im-1)   :: yx          ! 2 次元格子点
      real(8)                             :: IntYX_yx    ! 積分値
      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 ee_JacobianZ_ee(ee_zeta) ! スペクトルに作用するヤコビアン J(psi,zeta)
      real(8), dimension(-lm:lm,-km:km)              :: ee_JacobianZ_ee
      real(8), dimension(-lm:lm,-km:km), intent(in)  :: ee_Zeta

      integer k,l

      call n2ajbs(lm,km,jm,im,ee_Zeta,ee_tmp,ee_work,yx_work3,itj,tj,iti,ti)
     
      do k=-km,km
         do l=-lm,lm
            ee_JacobianZ_ee(l,k)                  = (2*pi/xl)*(2*pi/yl)/((2*pi/xl)**2+(2*pi/yl)**2)                    * ee_tmp(l,k)
         enddo
      enddo

!!$      ee_tmp=ee_lapla_inv(zeta)
!!$      gv=yx_ee(ee_dx_ee(ee_tmp)) 
 gu=yx_ee(ee_dy_ee(ee_tmp))
!!$      suv=ee_yx(gu*gv) 
 suuvv=ee_yx(gv**2-gu**2)
!!$      ee_jacz_ee = ee_dx_ee(ee_dx_ee(suv)) - ee_dy_ee(ee_dy_ee(suv)) 
!!$                - ee_dx_ee(ee_dy_ee(suuvv))
    end function ee_JacobianZ_ee

    function ee_Lapla_ee(ee)   ! スペクトルに作用する Laplacian 演算子
      real(8), dimension(-lm:lm,-km:km)              :: ee_Lapla_ee
      real(8), dimension(-lm:lm,-km:km), intent(in)  :: ee
      integer k,l

      do k=-km,km
         do l=-lm,lm
            ee_Lapla_ee(l,k) = -((2*pi*k/xl)**2+(2*pi*l/yl)**2)*ee(l,k)
         enddo
      enddo
    end function ee_Lapla_ee

    subroutine ee_initial(i,j,k,l,xmin,xmax,ymin,ymax)

      integer,intent(in) :: i, j           ! 格子点の設定(X,Y)
      integer,intent(in) :: k, l           ! 切断波数の設定(X,Y)

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

      integer :: ii, jj

      im = i         
 jm = j
      km = k         
 lm = l
      xl = xmax-xmin 
 yl = ymax-ymin

      allocate(tj(jm*2),ti(im*2))
      allocate(yx_work(jm*im),ee_work((2*km+1)*(2*lm+1)),yx_work3(jm*im*3))
      allocate(yx_tmp(0:jm-1,0:im-1),yx_tmp1(0:jm-1,0:im-1))
      allocate(ee_tmp(-lm:lm,-km:km))

      call n2init(jm,im,itj,tj,iti,ti)

      allocate(x_X(0:im-1), x_X_Weight(0:im-1))
      allocate(y_Y(0:jm-1), y_Y_Weight(0:jm-1))
      allocate(yx_X(0:jm-1,0:im-1), yx_Y(0:jm-1,0:im-1))

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

      do jj=0,jm-1
         y_Y(jj) = ymin + yl/jm*jj
      enddo
      y_Y_Weight = yl/jm

      yx_X = spread(x_X,1,jm)
      yx_Y = spread(y_Y,2,im)

    end subroutine ee_initial

    function yx_ee(ee) ! スペクトル -> 格子
      real(8), dimension(0:jm-1,0:im-1)             :: yx_ee
      real(8), dimension(-lm:lm,-km:km), intent(in) :: ee

      call n2s2ga(lm,km,jm,im,ee,yx_ee,yx_work,itj,tj,iti,ti)
    end function yx_ee