!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  SPPACK/ISPACK-F90
!
!履歴  2002/02/02  竹広真一 
!
!      制限
!         ・変換する格子点データ, スペクトルデータの配列の大きさは決めうち
!
module sppack_3D
  use snpack_2D, only : im, jm, nm
  use snpack_3D, only : km, sninitial_3D, gg_ep_3D, ep_gg_3D
  use sppack_2D, only : rn, irm, ep_jacobian_ep_ep_2D

  implicit none

  private
 
  public spinitial_3D                         ! 初期化
  public ep_lapla_ep_3D, ep_laplainv_ep_3D    ! ラプラシアンと逆演算
  public ep_dlon_ep_3D                        ! 経度微分
  public gg_gradlon_ep_3D, gg_gradlat_ep_3D   ! 勾配型微分
  public ep_divlon_gg_3D, ep_divlat_gg_3D     ! 発散型微分
  public ep_div_gg_gg_3D                      ! 発散型微分
  public ep_jacobian_ep_ep_3D                 ! ヤコビアン

  contains

  !--------------- 初期化 -----------------
    subroutine spinitial_3D(k_in)

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

      call sninitial_3D(k_in)

    end subroutine spinitial_3D

  !--------------- 微分計算 -----------------
    function ep_lapla_ep_3D(ep)       ! スペクトルに作用する Laplacian
      real(8), intent(in)  :: ep(:,:)
      real(8)              :: 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 演算子
      real(8), intent(in)  :: ep(:,:)
      real(8)              :: 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)        ! スペクトルに作用する経度微分 ∂/∂λ
      real(8), intent(in)  :: ep(:,:)
      real(8)              :: 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φ・∂/∂λ
      real(8), intent(in)  :: ep(:,:)
      real(8)              :: gg_gradlon_ep_3D(im,jm,size(ep,2))

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

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

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

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

      real(8), intent(in)  :: gg(:,:,:)
      real(8)              :: ep_divlon_gg_3D((nm+1)*(nm+1),size(gg,3))

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

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

      real(8), intent(in)  :: gg(:,:,:)
      real(8)              :: ep_divlat_gg_3D((nm+1)*(nm+1),size(gg,3))

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

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

      ep_div_gg_gg_3D = ep_divlon_gg_3D(gg_u) + ep_divlat_gg_3D(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/∂φ

      real(8), intent(in) :: ep_a(:,:)
      real(8), intent(in) :: ep_b(:,:)
      real(8)             :: 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_2D(ep_a(:,k),ep_b(:,k))
      end do
    end function ep_jacobian_ep_ep_3D

  end module sppack_3D
