!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  SPPACK/ISPACK-F90
!
!履歴  2002/02/02  竹広真一 
!
!      制限
!         ・変換する格子点データ, スペクトルデータの配列の大きさは決めうち
!
module sppack_multi
  use snpack_multi
  use sppack, only : spinitial, ep_lapla_ep, ep_laplainv_ep, ep_dlon_ep, &
                     gg_gradlon_ep, gg_gradlat_ep, &
                     ep_divlon_gg, ep_divlat_gg, ep_div_gg_gg, &
                     ep_jacobian_ep_ep

  implicit none

  real(8), allocatable  :: rn(:,:)            ! ラプラシアン演算用配列
  integer, allocatable  :: irm(:,:)           ! 経度微分演算用配列

  interface ep_lapla_ep
     module procedure ep_lapla_ep, ep_lapla_ep_3D
  end interface

  interface ep_laplainv_ep
     module procedure ep_laplainv_ep, ep_laplainv_ep_3D
  end interface

  interface ep_dlon_ep
     module procedure ep_dlon_ep, ep_dlon_ep_3D
  end interface

  interface gg_gradlon_ep
     module procedure gg_gradlon_ep, gg_gradlon_ep_3D
  end interface

  interface gg_gradlat_ep
     module procedure gg_gradlat_ep, gg_gradlat_ep_3D
  end interface

  interface ep_divlon_gg
     module procedure ep_divlon_gg, ep_divlon_gg_3D
  end interface

  interface ep_divlat_gg
     module procedure ep_divlat_gg, ep_divlat_gg_3D
  end interface

  interface ep_div_gg_gg
     module procedure ep_div_gg_gg, ep_div_gg_gg_3D
  end interface

  interface ep_jacobian_ep_ep
     module procedure ep_jacobian_ep_ep, ep_jacobian_ep_ep_3D
  end interface

  private

  public longitude, latitude                  ! 格子座標
  public lon_weight, lat_weight               ! 格子座標重み
  public gg_lon, gg_lat                       ! 格子座標(im,jm)
  public gg_ep, ep_gg, l_nm, nm_l             ! 変換関数
  
  public spminitial                           ! 初期化
  public ep_lapla_ep, ep_laplainv_ep          ! ラプラシアンと逆演算
  public ep_dlon_ep                           ! 経度微分
  public gg_gradlon_ep, gg_gradlat_ep         ! 勾配型微分
  public ep_divlon_gg, ep_divlat_gg           ! 発散型微分
  public ep_div_gg_gg                         ! 発散型微分
  public ep_jacobian_ep_ep                    ! ヤコビアン

  save rn, irm

  contains

  !--------------- 初期化 -----------------
    subroutine spminitial(n_in,i_in,j_in,k_in)

      integer,intent(in) :: n_in              ! 切断波数の設定
      integer,intent(in) :: i_in, j_in        ! 格子点の設定(東西, 南北)
      integer,intent(in) :: k_in              ! 最大データ数(層数)を記憶

      call spinitial(n_in,i_in,j_in)
      call snminitial(k_in)

      allocate(rn((nm+1)*(nm+1),2))           ! ラプラシアン演算用配列
      allocate(irm((nm+1)*(nm+1),2))          ! 経度微分演算用配列
      call spnini(nm,rn)
      call spmini(nm,irm)

    end subroutine spminitial

  !--------------- 微分計算 -----------------
    function ep_lapla_ep_3D(ep)       ! スペクトルに作用する Laplacian
      double precision, intent(in)  :: ep(:,:)
      double precision              :: ep_lapla_ep_3D((nm+1)*(nm+1),size(ep,2))

      integer :: l,k

      do k=1,size(ep,2)
         do l=1,(nm+1)*(nm+1)
            ep_lapla_ep_3D(l,k) = rn(l,1)*ep(l,k)
         enddo
      enddo
    end function ep_lapla_ep_3D

    function ep_laplainv_ep_3D(ep)    ! スペクトルに作用する逆 Laplacian 演算子
      double precision, intent(in)  :: ep(:,:)
      double precision              :: ep_laplainv_ep_3D((nm+1)*(nm+1),size(ep,2))

      integer :: l,k

      do k=1,size(ep,2)
         do l=1,(nm+1)*(nm+1)
            ep_laplainv_ep_3D(l,k) = rn(l,2)*ep(l,k)
         enddo
      enddo
    end function ep_laplainv_ep_3D

    function ep_dlon_ep_3D(ep)        ! スペクトルに作用する経度微分 ∂/∂λ
      double precision, intent(in)  :: ep(:,:)
      double precision              :: ep_dlon_ep_3D((nm+1)*(nm+1),size(ep,2))

      integer :: l,k

      do k=1,size(ep,2)
         do l=1,(nm+1)*(nm+1)
            ep_dlon_ep_3D(irm(l,1),k) = irm(l,2)*ep(l,k)
         enddo
      enddo
    end function ep_dlon_ep_3D

    function gg_gradlon_ep_3D(ep) ! スペクトルに作用する勾配型経度微分
                                  ! 1/cosφ・∂/∂λ
      double precision, intent(in)  :: ep(:,:)
      double precision              :: gg_gradlon_ep_3D(im,jm,size(ep,2))

      gg_gradlon_ep_3D = gg_ep(ep,ipow=1,iflag=-1)
    end function gg_gradlon_ep_3D

    function gg_gradlat_ep_3D(ep) ! スペクトルに作用する勾配型緯度微分, ∂/∂φ
      double precision, intent(in)  :: ep(:,:)
      double precision              :: gg_gradlat_ep_3D(im,jm,size(ep,2))

      gg_gradlat_ep_3D = gg_ep(ep,ipow=1,iflag=1)
    end function gg_gradlat_ep_3D

    function ep_divlon_gg_3D(gg)   ! 格子に作用する発散型経度微分 
                                   ! 1/cosφ・∂/∂λ

      double precision, intent(in)  :: gg(:,:,:)
      double precision              :: ep_divlon_gg_3D((nm+1)*(nm+1),size(gg,3))

      ep_divlon_gg_3D = ep_gg(gg,ipow=1,iflag=-1)
    end function ep_divlon_gg_3D

    function ep_divlat_gg_3D(gg)   ! 格子に作用する発散型緯度微分
                                   ! 1/cosφ・∂(f cosφ)/∂φ

      double precision, intent(in)  :: gg(:,:,:)
      double precision              :: ep_divlat_gg_3D((nm+1)*(nm+1),size(gg,3))

      ep_divlat_gg_3D = ep_gg(gg,ipow=1,iflag=1)
    end function ep_divlat_gg_3D

    function ep_div_gg_gg_3D(gg_u,gg_v)   ! 格子に作用する発散
      double precision, intent(in)  :: gg_u(:,:,:)   ! ベクトル経度成分
      double precision, intent(in)  :: gg_v(:,:,:)   ! ベクトル緯度成分
      double precision              :: ep_div_gg_gg_3D((nm+1)*(nm+1),size(gg_u,3))

      ep_div_gg_gg_3D = ep_divlon_gg(gg_u) + ep_divlat_gg(gg_v)
    end function ep_div_gg_gg_3D

    function ep_jacobian_ep_ep_3D(ep_a,ep_b)  ! スペクトルに作用するヤコビアン
                             ! J(f,g) = ∂f/∂λ・∂g/∂μ - ∂g/∂λ・∂f/∂μ
                             !        = ∂f/∂λ・1/cosφ・∂g/∂φ
                             !           - ∂g/∂λ・1/cosφ・∂f/∂φ

      double precision, intent(in) :: ep_a(:,:)
      double precision, intent(in) :: ep_b(:,:)
      double precision             :: ep_jacobian_ep_ep_3D((nm+1)*(nm+1),size(ep_a,2))

      integer :: k

      do k=1,size(ep_a,2)
         ep_jacobian_ep_ep_3D(:,k) = ep_jacobian_ep_ep(ep_a(:,k),ep_b(:,k))
      end do
    end function ep_jacobian_ep_ep_3D

  end module sppack_multi
