!表題  STPACK/ISPACK-F90
!
!履歴  2001/10/07  竹広真一
!
module stpack
  implicit none

  integer            :: im=64, jm=32     ! 格子点の設定(東西, 南北)
  integer            :: nm=21            ! 切断波数の設定
  double precision   :: rad              ! 球の大きさ

  integer,dimension(5)                       :: it
  double precision,dimension(:),allocatable  :: t
  double precision,dimension(:),allocatable  :: q, r
  double precision,dimension(:),allocatable  :: x, y

  double precision, dimension(:),   allocatable  :: wp
  double precision, dimension(:,:), allocatable  :: wgg,wss,wsc
  double precision, parameter  :: pi=3.141592653589793

  private
  public stinitial
  public g_s, s_g, l_nm, nm_l, xgrid, ygrid
  public s_lapla_s, s_laplainv_s, s_jac_ss
  public v_grad_s, v_ktimesnabla_s, s_div_v, s_rot_v

  save im, jm, nm, rad, it, t, q, r, x, y

  contains
  !--------------- 初期化 -----------------
    subroutine stinitial(n,i,j,radius)

      integer,intent(in) :: i, j              ! 格子点の設定(東西, 南北)
      integer,intent(in) :: n                 ! 切断波数の設定

      double precision,intent(in) :: radius   ! 領域の大きさ

      im = i  ; jm = j  ; nm = n
      rad = radius

      allocate(t(im*2),q(jm*(nm+1)),r((nm+1)*(nm+1)))   ! 変換用配列
      allocate(x(im),y(jm))                             ! 格子点座標格納配列

      allocate(wp(jm*im*4))                             ! 作業用配列
      !allocate(ws((2*km+1)*(lm+1)),wgj((jm+1)*im*3))
      !allocate(wss(-km:km,lm),wsc(-km:km,0:lm))

      call stinit(nm,jm,im,q,r,it,t)

      call stogrd(jm,im,y,x,q)
    end subroutine stinitial

  !--------------- 基本変換 -----------------
    function xgrid()    ! 格子点データの座標(X)
      double precision, dimension(im)     :: xgrid
      xgrid = x
    end function xgrid

    function ygrid()    ! 格子点データの座標(Y)
      double precision, dimension(jm)     :: ygrid
      ygrid = y
    end function ygrid

    function l_nm(n,m) ! スペクトルデータの格納位置
      integer               :: l_nm     ! スペクトルデータ位置
      integer, intent(in)   :: n, m     ! 全波数, 帯状波数

      call stnm2l(nm,n,m,l_nm)
    end function l_nm

    function nm_l(l)   ! スペクトルデータの格納位置から波数を求める
      integer, dimension(2) :: nm_l     ! 全波数, 帯状波数
      integer, intent(in)   :: l        ! スペクトルデータの格納位置

      call stl2nm(nm,nm_l(1),nm_l(2),l)
    end function nm_l

    function g_s(s) ! スペクトル -> 格子点
      double precision, dimension(jm,im)                     :: g_s  ! 格子点
      double precision, dimension((nm+1)*(nm+1)), intent(in) :: s    ! スペクトル

      call sts2ga(nm,jm,im,s,g_s,wp,q,r,it,t)
    end function g_s

    function s_g(g) ! 格子点 -> スペクトル
      double precision, dimension((nm+1)*(nm+1))      :: s_g    ! スペクトル
      double precision, dimension(jm,im), intent(in)  :: g      ! 格子点

      call stg2sa(nm,jm,im,g,s_g,wp,q,r,it,t)
    end function s_g

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

      call stclfa(nm,s,s_lapla_s)
    end function s_lapla_s

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

      call stclba(nm,s,s_laplainv_s)
    end function s_laplainv_s

    function s_jac_ss(sa,sb) ! スペクトルに作用するヤコビアン J(sa,sb)
      double precision, dimension((nm+1)*(nm+1))              :: s_jac_ss
      double precision, dimension((nm+1)*(nm+1)), intent(in)  :: sa, sb

      call stajba(nm,jm,im,sa,sb,s_jac_ss,wp,q,r,it,t)
    end function s_jac_ss

  !--------------- ベクトル解析計算 -----------------
    function v_grad_s(s)    ! スペクトルデータから勾配ベクトル格子点データを計算
      double precision, dimension(jm,im,2)                    :: v_grad_s ! ベクトルデータ
      double precision, dimension((nm+1)*(nm+1)), intent(in)  :: s        ! 入力スペクトル

      call sts2va(nm,jm,im,s,v_grad_s(1,1,1),v_grad_s(1,1,2),wp,q,r,it,t)
    end function v_grad_s

    function v_ktimesnabla_s(s)  ! スペクトルデータから回転的勾配ベクトルを計算
      double precision, dimension(jm,im,2)                    :: v_ktimesnabla_s ! ベクトルデータ
      double precision, dimension((nm+1)*(nm+1)), intent(in)  :: s        ! 入力スペクトル

      call stsrva(nm,jm,im,s,v_ktimesnabla_s(1,1,1),v_ktimesnabla_s(1,1,2),wp,q,r,it,t)
    end function v_ktimesnabla_s

    function s_div_v(v)    ! ベクトル格子点データの発散
      double precision, dimension((nm+1)*(nm+1))        :: s_div_v ! スペクトル
      double precision, dimension(jm,im,2), intent(in)  :: v       ! ベクトルデータ

      call stv2sa(nm,jm,im,v(1,1,1),v(1,1,2),s_div_v,wp,q,r,it,t)
    end function s_div_v

    function s_rot_v(v)    ! ベクトル格子点データの回転
      double precision, dimension((nm+1)*(nm+1))        :: s_rot_v ! スペクトル
      double precision, dimension(jm,im,2), intent(in)  :: v       ! ベクトルデータ

      call stvrsa(nm,jm,im,v(1,1,1),v(1,1,2),s_rot_v,wp,q,r,it,t)
    end function s_rot_v

  end module stpack
