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

  integer               :: im=64, jm=32     ! 格子点の設定(東西, 南北)
  integer               :: nm=21            ! 切断波数の設定

  integer               :: it(6)            ! 変換用配列
  real(8), allocatable  :: t(:)             ! 変換用配列
  integer, allocatable  :: ip(:)            ! 変換用配列
  real(8), allocatable  :: p(:), r(:)       ! 変換用配列
  integer, allocatable  :: ia(:)            ! 変換用配列
  real(8), allocatable  :: a(:)             ! 変換用配列
  real(8), allocatable  :: y(:,:)           ! 変換用配列

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

  real(8), allocatable  :: longitude(:), latitude(:)     ! 緯度経度
  real(8), allocatable  :: lon_weight(:), lat_weight(:)  ! 座標重み
  real(8), allocatable  :: gg_lon(:,:), gg_lat(:,:)

  real(8), parameter    :: pi=3.1415926535897932385D0

  private

  public im, jm, nm                           ! 格子点数, 切断波数, 半径
  public it, t, y, ip, p, r, ia, a            ! 変換用作業配列

  public sninitial_2D                         ! 初期化サブルーチン
  public longitude, latitude                  ! 格子座標
  public lon_weight, lat_weight               ! 格子座標重み
  public gg_lon, gg_lat                       ! 格子座標(im,jm)
  public l_nm, nm_l                           ! 波数格納位置
  public gg_ep_2D, ep_gg_2D                   ! 変換関数

  save im, jm, nm                             ! 格子点数, 切断波数, 半径を記憶
  save it, t, y, ip, p, r, ia, a              ! 変換用配列を記憶

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

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

      integer :: iw, i, j

      im = i_in  ; jm = j_in  ; nm = n_in

      allocate(t(im*2))                       ! 変換用配列
      allocate(ip(((nm+1)/2+nm+1)*2))         ! 変換用配列
      allocate(p(((nm+1)/2+nm+1)*jm))         ! 変換用配列
      allocate(r(((nm+1)/2*2+3)*(nm/2+1)))    ! 変換用配列
      allocate(ia((nm+1)*(nm+1)*4))           ! 変換用配列
      allocate(a((nm+1)*(nm+1)*6))            ! 変換用配列
      allocate(y(jm/2,4))                     ! 変換用配列

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

      allocate(longitude(im) ,latitude(jm))   ! 格子点座標格納配列
      allocate(lon_weight(im),lat_weight(jm)) ! 格子点座標格納配列
      allocate(gg_lon(im,jm),gg_lat(im,jm))   ! 格子点座標格納配列

      call sninit(nm,im,jm,it,t,y,ip,p,r,ia,a)

      do i=1,im
         longitude(i)  = 2*pi/im*(i-1)        ! 経度座標
         lon_weight(i) = 2*pi/im              ! 経度座標重み
      enddo

      do j=1,jm/2
         latitude(jm/2+j)   =  asin(y(j,1))   ! 緯度座標
         latitude(jm/2-j+1) = -asin(y(j,1))   ! 緯度座標
         lat_weight(jm/2+j)   = y(j,2)        ! 緯度重み(Gauss grid)
         lat_weight(jm/2-j+1) = y(j,2)        ! 緯度重み(Gauss grid)
      enddo

      do j=1,jm
         gg_lon(:,j) = longitude
      enddo

      do i=1,im
         gg_lat(i,:) = latitude
      enddo

    end subroutine sninitial_2D

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

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

      call snnm2l(n,m,l_nm)
    end function l_nm

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

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

    function gg_ep_2D(ep,ipow,iflag)    ! 球面調和関数スペクトル -> 格子点
      real(8)               :: gg_ep_2D(im,jm)    ! 格子点
      real(8), intent(in)   :: ep((nm+1)*(nm+1))  ! スペクトル
      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

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

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

      call snts2g(nm,im,im,jm,jm,1,ep,gg_ep_2D,&
           it,t,y,ip,p,r,ia,a,q,ws,ww,ipval,ifval)

    end function gg_ep_2D

    function ep_gg_2D(gg,ipow,iflag) ! 格子点 -> 球面調和関数スペクトル
      real(8)               :: ep_gg_2D((nm+1)*(nm+1)) ! スペクトル
      real(8), intent(in)   :: gg(im,jm)               ! 格子点
      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

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

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

      call sntg2s(nm,im,im,jm,jm,1,gg,ep_gg_2D,&
           it,t,y,ip,p,r,ia,a,q,ws,ww,ipval,ifval)
    end function ep_gg_2D

  end module snpack_2D

