!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  SPPACK/ISPACK-F90
!
!履歴  2001/12/08  竹広真一
!      2001/12/26  竹広真一 関数,変数名前変更
!
!      制限
!         ・変換する格子点データ, スペクトルデータの配列の大きさは決めうち
!         ・とりあえず 1 層用
!
module sppack
  use snpack
  implicit none

  real(8), allocatable  :: rn(:,:)            ! ラプラシアン演算用配列
  integer, allocatable  :: irm(:,:)           ! 経度微分演算用配列
  integer, allocatable  :: ip2(:), ip3(:)     ! ヤコビアン計算用配列
  real(8), allocatable  :: p2(:), p3(:)       ! ヤコビアン計算用配列
  real(8), allocatable  :: r2(:), r3(:)       ! ヤコビアン計算用配列

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

  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 spinitial                            ! 初期化
  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, ip2, ip3, p2, p3, r2, r3

  contains

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

      integer,intent(in) :: i_in, j_in        ! 格子点の設定(東西, 南北)
      integer,intent(in) :: n_in              ! 切断波数の設定

      integer iw

      call sninitial(n_in,i_in,j_in)

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

      allocate(ip2(2*((nm+1)/2+nm+1)*2))      ! ヤコビアン計算用配列
      allocate(p2(2*((nm+1)/2+nm+1)*jm))      ! ヤコビアン計算用配列
      allocate(r2(2*((nm+1)/2*2+3)*(nm/2+1))) ! ヤコビアン計算用配列
      allocate(ip3(3*((nm+1)/2+nm+1)*2))      ! ヤコビアン計算用配列
      allocate(p3(3*((nm+1)/2+nm+1)*jm))      ! ヤコビアン計算用配列
      allocate(r3(3*((nm+1)/2*2+3)*(nm/2+1))) ! ヤコビアン計算用配列
      call snkini(nm,jm,2,ip,p,r,ip2,p2,r2)
      call snkini(nm,jm,3,ip,p,r,ip3,p3,r3)

      allocate(q(3*((nm+1)/2+nm+1)*jm))       ! 作業用配列
      iw=3*max( ((nm+1)/2*2+3)*(nm/2+2)*2, &
                jm*((nm+1)/2+nm+1)*2, jm*jm )
      allocate(ws(iw),ww(iw))                 ! 作業用配列
    end subroutine spinitial

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

      call spclap(nm,ep,ep_lapla_ep,rn(1,1))
    end function ep_lapla_ep

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

      call spclap(nm,ep,ep_laplainv_ep,rn(1,2))
    end function ep_laplainv_ep

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

      call spclam(nm,ep,ep_dlon_ep,irm)
    end function ep_dlon_ep

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

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

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

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

    function ep_divlon_gg(gg)   ! 格子に作用する発散型経度微分 1/cosφ・∂/∂λ
      double precision              :: ep_divlon_gg((nm+1)*(nm+1))
      double precision, intent(in)  :: gg(im,jm)

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

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

      double precision              :: ep_divlat_gg((nm+1)*(nm+1))
      double precision, intent(in)  :: gg(im,jm)

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

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

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

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

      double precision             :: ep_jacobian_ep_ep((nm+1)*(nm+1))
      double precision, intent(in) :: ep_a((nm+1)*(nm+1))
      double precision, intent(in) :: ep_b((nm+1)*(nm+1))

      call spnjcb(nm,im,im,jm,jm,ep_a,ep_b,ep_jacobian_ep_ep,&
           it,t,y,ip2,p2,r2,ip3,p3,r3,ia,a,q,ws,ww)
    end function ep_jacobian_ep_ep

  end module sppack
