!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  SNPACK/ISPACK-F90
!
!履歴  2002/02/02  竹広真一  多層用に改造
!
!         ・変換する格子点データ, スペクトルデータの配列の大きさは決めうち
!
module snpack_3D

  use snpack_2D, only : im, jm, nm, it, t, y, ip, p, r, ia, a
  implicit none

  integer               :: km=16         ! 同時に処理する最大データ数(層の数)

  integer, allocatable  :: ipk(:,:)            ! 変換用配列(多層用)
  real(8), allocatable  :: pk(:,:), rk(:,:)    ! 変換用配列(多層用)

  real(8), allocatable  :: q(:)                ! 作業配列
  real(8), allocatable  :: ww(:), ws(:)        ! 作業配列

  real(8), parameter    :: pi=3.14159265358979

  private

  public km                                    ! 層数
  public sninitial_3D                          ! 初期化サブルーチン
  public gg_ep_3D, ep_gg_3D                    ! 変換関数

  save km                                      ! 最大データ数(層数)を記憶
  save ipk, pk, rk                             ! 変換用配列を記憶

  contains
  !--------------- 初期化 -----------------
    subroutine sninitial_3D(k_in)

      integer,intent(in) :: k_in               ! 最大データ数(層数)を記憶

      integer :: iw

      km = k_in

      allocate(ipk(km,((nm+1)/2+nm+1)*2))      ! 変換用配列(多層用)
      allocate(pk(km,((nm+1)/2+nm+1)*jm))      ! 変換用配列(多層用)
      allocate(rk(km,((nm+1)/2*2+3)*(nm/2+1))) ! 変換用配列(多層用)

      allocate(q(km*((nm+1)/2+nm+1)*jm))       ! 作業配列(多層用)
      iw=km*(im+3*(nm+1))*jm
      allocate(ws(iw),ww(iw))                  ! 作業用配列(多層用)

      call snkini(nm,jm,km,ip,p,r,ipk,pk,rk)

    end subroutine sninitial_3D

  !--------------- 基本変換 -----------------

    function gg_ep_3D(ep,ipow,iflag)    ! 球面調和関数スペクトル -> 格子点
      real(8), intent(in)   :: ep(:,:)                    ! スペクトル
      real(8)               :: gg_ep_3D(im,jm,size(ep,2)) ! 格子点
      integer, intent(in), optional  :: ipow      ! 作用させる 1/cosφ の次数
      integer, intent(in), optional  :: iflag     ! 変換の種類

      integer, parameter  :: ipow_default  = 0
      integer, parameter  :: iflag_default = 0

      integer ipval, ifval
      integer k

      if (present(ipow)) then
         ipval = ipow
      else
         ipval = ipow_default
      endif

      if (present(iflag)) then
         ifval = iflag
      else
         ifval = iflag_default
      endif

      k= size(ep,2)
      if  ( k > km ) then
         call msgdmp('E','gg_ep','Size of 3rd dimension invalid.')
      else
         call snts2g(nm,im,im,jm,jm,size(ep,2),ep,gg_ep_3D,&
              it,t,y,ipk(1:k,:),pk(1:k,:),rk(1:k,:),ia,a,q,ws,ww,ipval,ifval)
      endif

    end function gg_ep_3D

    function ep_gg_3D(gg,ipow,iflag) ! 格子点 -> 球面調和関数スペクトル
      real(8), intent(in)   :: gg(:,:,:)      ! 格子点(im,jm,*)
      real(8)               :: ep_gg_3D((nm+1)*(nm+1),size(gg,3))  ! スペクトル
      integer, intent(in), optional  :: ipow        ! 作用させる 1/cosφ の次数
      integer, intent(in), optional  :: iflag       ! 変換の種類

      integer, parameter  :: ipow_default  = 0      ! スイッチデフォルト値
      integer, parameter  :: iflag_default = 0      ! スイッチデフォルト値

      integer ipval, ifval
      integer k

      if (present(ipow)) then
         ipval = ipow
      else
         ipval = ipow_default
      endif

      if (present(iflag)) then
         ifval = iflag
      else
         ifval = iflag_default
      endif

      k = size(gg,3)
      if ( k > km ) then
         call msgdmp('E','ep_gg','Size of 3rd dimension invalid.')
      else
         call sntg2s(nm,im,im,jm,jm,k,gg,ep_gg_3D,&
              it,t,y,ipk(1:k,:),pk(1:k,:),rk(1:k,:),ia,a,q,ws,ww,ipval,ifval)
      endif
    end function ep_gg_3D

  end module snpack_3D
