!----------------------------------------------------------------------
!     Copyright (c) 2001 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  CTPACK/ISPACK-F90
!      水路領域問題, Fourier 展開 + Chebyshev 展開法
!
!履歴  2002/01/18  竹広真一
!
module ctpack
  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  :: gg_work,et_work
  double precision, parameter  :: pi=3.1415926535897932385D0
  double precision, dimension(:),allocatable  :: g_x, g_y

  private
  public ctinitial
  public gg_x, gg_y
  public gg_et, et_gg
  public et_lapla_et, et_dx_et, et_dy_et

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

  contains
  !--------------- 初期化 -----------------
    subroutine ctinitial(i,j,k,l,xmin,xmax,ymin,ymax)

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

      double precision,intent(in)  :: xmin, xmax         ! Y 座標範囲
      double precision,intent(in)  :: ymin, ymax         ! Y 座標範囲

      integer ix, jy                       ! DO 変数

      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(gg_work(0:jm,0:im-1))
      allocate(et_work(-km:km,0:lm))

      allocate(g_x(0:im-1),g_y(0:jm))

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

      do ix=0,im-1
         g_x(ix) = xmin + (xmax-xmin)/im*ix
      enddo

      do jy=0,jm
         g_y(jy) = (ymin+ymax)/2 + (ymax-ymin)/2 * cos(pi*jy/jm)
      enddo

    end subroutine ctinitial

  !--------------- 基本変換 -----------------
    function gg_x()                    ! 格子点データの座標(経度)
      double precision :: gg_x(0:jm,0:im-1)

      integer :: j

      do j=0,jm
         gg_x(j,:) = g_x
      enddo
    end function gg_x

    function gg_y()                    ! 格子点データの座標(緯度)
      double precision :: gg_y(0:jm,0:im-1)

      integer          :: i

      do i=0,im-1
         gg_y(:,i) = g_y
      enddo
    end function gg_y

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

      ! Y 波数 0 成分の和が Factor 1/2 であることの処理. 
      ! FTTCTB では重み 1/2 の和なのだが, C2S2GA 内部で修正されている. 
      ! それを再び修正する. 
      ! 最大波数成分に対する和の重みも Factor 1/2 であるはずだが, 
      ! 実は LM < JM なる条件からグリッド数に対応する最大波数が
      ! 切断されているので考慮しなくて良い. 

      et_work=et
      et_work(:,0)=et(:,0)/2.0D0

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

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

      gg_work = gg
      call c2g2sa(lm,km,jm,im,gg_work,et_work,wg,itj,tj,iti,ti,2)

      ! Y 波数 0 成分の和が Factor 1/2 であることの処理. 
      ! FTTCTF では重み 1/2 の和なのだが, C2G2SSA 内部で修正されている. 
      ! それを再び修正する. 
      ! 最大波数成分に対する和の重みも Factor 1/2 であるはずだが, 
      ! 実は LM < JM なる条件からグリッド数に対応する最大波数が
      ! 切断されているので考慮しなくて良い. 

      et_gg=et_work
      et_gg(:,0)=et_work(:,0)*2.0D0

    end function et_gg

  !--------------- 微分計算 -----------------

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

      do l=0,lm
         do k=-km,km
            et_lapla_et(k,l) = -(2*pi*k/xl)**2*et(k,l)
         enddo
      enddo
      et_lapla_et = et_lapla_et + et_dy_et(et_dy_et(et))
    end function et_lapla_et

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

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

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

      ! 本来最大波数成分に対する和の重みが Factor 1/2 であることに
      ! 注意せねばならないが, 
      ! 実は LM < JM なる条件からグリッド数に対応する最大波数が
      ! 切断されているので考慮しなくて良い. 

      do k=-km,km
         et_dy_et(k,lm)   = 0. 
         et_dy_et(k,lm-1) = 2 * lm * et(k,lm) 
         ! スタートはグリッド対応最大波数未満. Factor 1/2 不要
      end do

      do l=lm-2,0,-1
         do k=-km,km
            et_dy_et(k,l) = et_dy_et(k,l+2) + 2*(l+1)*et(k,l+1)
         enddo
      enddo

      do l=0,lm                       !" 領域の factor
         do k=-km,km
            et_dy_et(k,l) = 2/yl * et_dy_et(k,l) 
         enddo
      enddo

    end function et_dy_et

  end module ctpack
