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

  integer            :: im=32, jm=8     ! 格子点の設定(X,Y)
  integer            :: km=10, lm=5     ! 切断波数の設定(X,Y)
  double precision   :: xl=2.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  :: wgg,wss,wsc
  double precision, parameter  :: pi=3.141592653589793

  private
  public c2initial
  public g_s, g_c, s_g, c_g
  public s_lapla_s, s_laplainv_s, s_dx_s, c_dy_s, s_jac_ss
  public c_lapla_c, c_dx_c, s_dy_c, c_jac_sc

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

  contains
  !--------------- 初期化 -----------------
    subroutine c2initial(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*6),ti(im*2))
      allocate(wg((jm+1)*im),wgg(0:jm,0:im-1))
      allocate(ws((2*km+1)*(lm+1)),wgj((jm+1)*im*3))
      allocate(wss(-km:km,lm),wsc(-km:km,0:lm))

      call c2init(jm,im,itj,tj,iti,ti)
    end subroutine c2initial

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

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

    function g_c(c)  ! スペクトル COS(Y) -> 台形格子
      double precision, dimension(0:jm,0:im-1)              :: g_c
      double precision, dimension(-km:km,0:lm), intent(in)  :: c

      call c2s2ga(lm,km,jm,im,c,g_c,wg,itj,tj,iti,ti,2)
    end function g_c

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

      wgg = g
      call c2g2sa(lm,km,jm,im,wgg,s_g,wg,itj,tj,iti,ti,1)
    end function s_g

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

      wgg = g
      call c2g2sa(lm,km,jm,im,wgg,c_g,wg,itj,tj,iti,ti,2)
    end function c_g

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

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

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

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

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

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

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

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

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

      do l=0,lm
         do k=-km,km
            c_dx_c(k,l)  =  (-2*pi*k/xl)*c(-k,l)
         enddo
      enddo
    end function c_dx_c

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

      do k=-km,km
         c_dy_s(k,0)  =  0.0
      enddo
      do l=1,lm
         do k=-km,km
            c_dy_s(k,l)  =  (pi*l/yl)*s(k,l)
         enddo
      enddo
    end function c_dy_s

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

      do l=1,lm
         do k=-km,km
            s_dy_c(k,l)  =  -(pi*l/yl)*c(k,l)
         enddo
      enddo
    end function s_dy_c

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

      integer k,l

      call c2ajcb(lm,km,jm,im,sa,sb,wss,ws,wgj,itj,tj,iti,ti)

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

    function c_jac_sc(s,c)    ! スペクトル COS(Y) に作用するヤコビアン
      double precision, dimension(-km:km,0:lm)              :: c_jac_sc
      double precision, dimension(-km:km,lm), intent(in)    :: s
      double precision, dimension(-km:km,0:lm), intent(in)  :: c
      integer k,l

      call c2ajcc(lm,km,jm,im,s,c,wsc,ws,wgj,itj,tj,iti,ti)

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

  end module c2pack
