!表題  N2PACK/ISPACK-F90
!
!履歴  2001/10/07 竹広真一
!
module n2pack
  implicit none

  integer            :: im=32, jm=32     ! 格子点の設定(X,Y)
  integer            :: km=10, lm=10     ! 切断波数の設定(X,Y)
  double precision   :: xl=1.0, yl=1.0   ! 領域の大きさ

  integer,dimension(5)                       :: itj
  double precision,dimension(:),allocatable  :: tj
  integer,dimension(5)                       :: iti
  double precision,dimension(:),allocatable  :: ti

  double precision, dimension(:),  allocatable  :: wg, ws, wgj
  double precision, dimension(:,:),allocatable  :: gtmp,gtmp1,stmp
  double precision, parameter  :: pi=3.141592653589793

  private
  public n2initial
  public g_s, s_g
  public s_lapla_s, s_laplainv_s, s_dx_s, s_dy_s
  public s_jacz_s, s_jac_ss

  save im, jm, km, lm, itj, tj, iti, ti, xl, yl

  contains
  !--------------- 初期化 -----------------
    subroutine n2initial(i,j,k,l,xlength,ylength)

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

      double precision,intent(in)   :: xlength, ylength   ! 領域の大きさ

      im = i       ; jm = j
      km = k       ; lm = l
      xl = xlength ; yl = ylength

      allocate(tj(jm*2),ti(im*2))
      allocate(wg(jm*im),ws((2*km+1)*(2*lm+1)),wgj(jm*im*3))
      allocate(gtmp(0:jm-1,0:im-1),gtmp1(0:jm-1,0:im-1))
      allocate(stmp(-lm:lm,-km:km))

      call n2init(jm,im,itj,tj,iti,ti)
    end subroutine n2initial

  !--------------- 基本変換 -----------------
    function g_s(s) ! スペクトル -> 格子
      double precision, dimension(0:jm-1,0:im-1)             :: g_s
      double precision, dimension(-lm:lm,-km:km), intent(in) :: s

      call n2s2ga(lm,km,jm,im,s,g_s,wg,itj,tj,iti,ti)
    end function g_s

    function s_g(g)  ! 格子 -> スペクトル
      double precision, dimension(-lm:lm,-km:km)              :: s_g
      double precision, dimension(0:jm-1,0:im-1), intent(in)  :: g

      gtmp = g
      call n2g2sa(lm,km,jm,im,gtmp,s_g,wg,itj,tj,iti,ti)
    end function s_g

  !--------------- 微分計算 -----------------
    function s_lapla_s(s)   ! スペクトルに作用する Laplacian 演算子
      double precision, dimension(-lm:lm,-km:km)              :: s_lapla_s
      double precision, dimension(-lm:lm,-km:km), intent(in)  :: s
      integer k,l

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

    function s_laplainv_s(s)   ! スペクトル SINY に作用する逆 \nabla 演算子
      double precision, dimension(-lm:lm,-km:km)             :: s_laplainv_s
      double precision, dimension(-lm:lm,-km:km), intent(in) :: s
      integer k,l

      do k=-km,km
         do l=-lm,lm
            if ( k.ne.0 .or. l.ne.0 ) then
               s_laplainv_s(l,k) = -s(l,k)/((2*pi*k/xl)**2+(2*pi*l/yl)**2)
            else
               s_laplainv_s(l,k) = 0.0
            endif
         enddo
      enddo
    end function s_laplainv_s

    function s_dx_s(s)   ! スペクトルに作用する x 微分演算子
      double precision, dimension(-lm:lm,-km:km)              :: s_dx_s
      double precision, dimension(-lm:lm,-km:km), intent(in)  :: s
      integer k,l

      do k=-km,km
         do l=-lm,lm
            s_dx_s(l,k) = -(2*pi*k/xl)*s(-l,-k)
         enddo
      enddo
    end function s_dx_s

    function s_dy_s(s)   ! スペクトルに作用する y 微分演算子
      double precision, dimension(-lm:lm,-km:km)              :: s_dy_s
      double precision, dimension(-lm:lm,-km:km), intent(in)  :: s
      integer k,l

      do k=-km,km
         do l=-lm,lm
            s_dy_s(l,k) = -(2*pi*l/yl)*s(-l,-k)
         enddo
      enddo

    end function s_dy_s

    function s_jac_ss(sa,sb)   ! スペクトルに作用するヤコビアン J(sa,sb)
      double precision, dimension(-lm:lm,-km:km)              :: s_jac_ss
      double precision, dimension(-lm:lm,-km:km), intent(in)  :: sa
      double precision, dimension(-lm:lm,-km:km), intent(in)  :: sb
      !integer k,l

      !s_jac_ss =   s_g(  g_s(s_dx_s(sa))*g_s(s_dy_s(sb)) &
      !                 - g_s(s_dy_s(sa))*g_s(s_dx_s(sb)))

      gtmp1=g_s(sa)
      s_jac_ss =   s_dx_s(s_g(gtmp1*g_s(s_dy_s(sb)))) &
                 - s_dy_s(s_g(gtmp1*g_s(s_dx_s(sb))))

    end function s_jac_ss

    function s_jacz_s(zeta) ! スペクトルに作用するヤコビアン J(psi,zeta)
      double precision, dimension(-lm:lm,-km:km)              :: s_jacz_s
      double precision, dimension(-lm:lm,-km:km), intent(in)  :: zeta

      integer k,l

      call n2ajbs(lm,km,jm,im,zeta,stmp,ws,wgj,itj,tj,iti,ti)
     
      do k=-km,km
         do l=-lm,lm
            s_jacz_s(l,k) = (2*pi/xl)*(2*pi/yl)/((2*pi/xl)**2+(2*pi/yl)**2) &
                            * stmp(l,k)
         enddo
      enddo
!!$
!!$      stmp=s_lapla_inv(zeta)
!!$      gv=g_s(s_dx_s(stmp)) ; gu=g_s(s_dy_s(stmp))
!!$      suv=s_g(gu*gv) ; suuvv=s_g(gv**2-gu**2)
!!$      s_jacz_s = s_dx_s(s_dx_s(suv)) - s_dy_s(s_dy_s(suv)) &
!!$                - s_dx_s(s_dy_s(suuvv))
    end function s_jacz_s

  end module n2pack
