!----------------------------------------------------------------------
!     Copyright (c) 2001 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  N2PACK/ISPACK-F90
!
!履歴  2001/10/07 竹広真一
!      2001/12/25 竹広真一  関数, 変数命名法変更
!
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  :: gg_work, ee_work, gg_work3
  double precision, dimension(:,:),allocatable  :: gg_tmp, gg_tmp1, ee_tmp
  double precision, parameter  :: pi=3.1415926535897932385D0

  private
  public n2initial
  public gg_ee, ee_gg
  public ee_lapla_ee, ee_laplainv_ee, ee_dx_ee, ee_dy_ee
  public ee_jacobianz_ee, ee_jacobian_ee_ee

  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(gg_work(jm*im),ee_work((2*km+1)*(2*lm+1)),gg_work3(jm*im*3))
      allocate(gg_tmp(0:jm-1,0:im-1),gg_tmp1(0:jm-1,0:im-1))
      allocate(ee_tmp(-lm:lm,-km:km))

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

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

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

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

      gg_tmp = gg
      call n2g2sa(lm,km,jm,im,gg_tmp,ee_gg,gg_work,itj,tj,iti,ti)
    end function ee_gg

  !--------------- 微分計算 -----------------
    function ee_lapla_ee(ee)   ! スペクトルに作用する Laplacian 演算子
      double precision, dimension(-lm:lm,-km:km)              :: ee_lapla_ee
      double precision, 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

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

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

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

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

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

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

    end function ee_dy_ee

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

      !ee_jacobian_ee_ee =   ee_gg(  gg_ee(ee_dx_ee(sa))*gg_ee(ee_dy_ee(sb)) &
      !                 - gg_ee(ee_dy_ee(sa))*gg_ee(ee_dx_ee(sb)))

      gg_tmp1=gg_ee(ee_a)
      ee_jacobian_ee_ee =   ee_dx_ee(ee_gg(gg_tmp1*gg_ee(ee_dy_ee(ee_b)))) &
                 - ee_dy_ee(ee_gg(gg_tmp1*gg_ee(ee_dx_ee(ee_b))))

    end function ee_jacobian_ee_ee

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

      integer k,l

      call n2ajbs(lm,km,jm,im,ee_zeta,ee_tmp,ee_work,gg_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=gg_ee(ee_dx_ee(ee_tmp)) ; gu=gg_ee(ee_dy_ee(ee_tmp))
!!$      suv=ee_gg(gu*gv) ; suuvv=ee_gg(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
   
  end module n2pack
