function AvrYX_yx(yx)    ! 全領域平均
      real(8), dimension(0:jm,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,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
            IntYX_yx = IntYX_yx + yx(j,i) * y_Y_Weight(j) * x_X_Weight(i)
         enddo
      enddo
    end function IntYX_yx

    function ec_Lapla_ec(ec)   ! スペクトル COSY に作用する \lapla 演算子
      real(8), dimension(-km:km,0:lm)                :: ec_Lapla_ec
      real(8), dimension(-km:km,0:lm), intent(in)    :: ec
      integer k,l

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

    function es_Jacobian_es_es(es_a,es_b) !スペクトル SINY に作用するヤコビアン
      real(8), dimension(-km:km,lm)                :: es_Jacobian_es_es
      real(8), dimension(-km:km,lm), intent(in)    :: es_A,es_B

      integer k,l

      call c2ajcb(lm,km,jm,im,es_A,es_B,es_work,ws,wgj,itj,tj,iti,ti)

      do l=1,lm
         do k=-km,km
            es_Jacobian_es_es(k,l) = (2*pi/xl)*(pi/yl)*es_work(k,l)
         enddo
      enddo
    end function es_Jacobian_es_es

    function es_Lapla_es(es)   ! スペクトル SINY に作用する \lapla 演算子
      real(8), dimension(-km:km,lm)                :: es_Lapla_es
      real(8), dimension(-km:km,lm), intent(in)    :: es
      integer k,l

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

    subroutine esc_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*6),ti(im*2))
      allocate(wg((jm+1)*im))
      allocate(ws((2*km+1)*(lm+1)),wgj((jm+1)*im*3))
      allocate(yx_work(0:jm,0:im-1))
      allocate(es_work(-km:km,lm),ec_work(-km:km,0:lm))

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

      allocate(x_X(0:im-1), x_X_Weight(0:im-1))
      allocate(y_Y(0:jm), y_Y_Weight(0:jm))
      allocate(yx_X(0:jm,0:im-1), yx_Y(0:jm,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
         y_Y(jj) = ymin + yl/jm*jj
      enddo
      y_Y_Weight(0) = yl/(2*jm)
      y_Y_Weight(1:jm-1) = yl/jm 
      y_Y_Weight(jm) = yl/(2*jm)

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

    end subroutine esc_initial

    function yx_es(es) ! スペクトル SIN(Y) -> 台形格子
      real(8), dimension(0:jm,0:im-1)              :: yx_es
      real(8), dimension(-km:km,lm), intent(in)    :: es

      call c2s2ga(lm,km,jm,im,es,yx_es,wg,itj,tj,iti,ti,1)
    end function yx_es