!-------------------------------------------------
!  SHTlib Module
!-------------------------------------------------
module shtrlib
!  use dcl_common
!  use sht_interface

  type list !作業領域などを格納するための連結リスト
    integer :: idx             ! 位置のインデックス
    type(list),pointer :: next ! 次のリストへのポインタ
    integer :: mm,jm,im        ! 格納する定数
    real, dimension(:), pointer :: array ! 格納する作業領域
  end type list

  type(list),pointer :: top,now,pre

  public
  private :: list,top,now,pre,seek
  save
  
contains
!------------------------
  function seek(idx) ! 連結リスト検索の関数
    ! idx に対応する位置があれば, 真となり, now がその位置へのポインタとなる.
    ! また, その位置が先頭でなければ, pre が前の位置へのポインタとなっている.
    ! idx に対応する位置がなければ, 偽となる.
    integer,intent(in) :: idx
    logical :: seek

    now => top

    do
      if(.not. associated(now)) then
        seek = .false.
        return
      else if(now%idx == idx) then
        seek = .true.
        return
      end if
      pre => now
      now => now%next
    end do
    
  end function seek
!-------------------------------------------------
  subroutine DclInitSHT(mm,jm,im,idx) ! 初期化ルーチン． 
    integer,intent(in)          :: mm    ! 切断波数
    integer,intent(in)          :: jm    ! 南北分割数の1/2
    integer,intent(in)          :: im    ! 東西分割数の1/2
    integer,intent(in),optional :: idx   ! 作業領域番号
    integer :: idxl, length

    call prcopn('DclInitSHT')
    idxl = 1; if(present(idx)) idxl = idx

    if(seek(idxl)) then
      call msgdmp('E','DclInitSHT','The working area has been allocated already.')
    else
      nullify(now)
      allocate(now)
      now%next => top
      top => now
      now%idx = idxl
      now%mm = mm
      now%jm = jm
      now%im = im
      length = (jm+1)*(4*jm+5*mm+14) + (mm+1)*(mm+1) +mm + 2 + 6*im + 15
      allocate(now%array(length))
      call shtint(mm,jm,im,now%array)
    end if
    call prccls('DclInitSHT')

  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocSHT(idx) ! 作業領域解放
    integer,optional :: idx    ! 作業領域番号
    integer :: idxl
    
    call prcopn('DclDeallocSHT')
    idxl = 1; if(present(idx)) idxl = idx

    if(associated(top)) then
      if(top%idx == idxl) then
        deallocate(top%array)
        now => top%next
        deallocate(top)
        top => now
      else
        if(seek(idxl)) then
          deallocate(now%array)
          pre%next => now%next
          deallocate(now)
        end if
      end if
    end if

    call prccls('DclDeallocSHT')

  end subroutine
!-------------------------------------------------
  function DclGetSpectrumNumber(n,m,idx) !スペクトルデータの格納位置を求める．
    integer :: DclGetSpectrumNumber
    integer,intent(in) :: n            ! 全波数
    integer,intent(in) :: m            ! 帯状波数
    integer,intent(in),optional :: idx ! 作業領域番号
    integer :: idxl, lr, li

    call prcopn('DclGetSpectrumNumber')
    idxl = 1; if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclGetSpectrumNumber','Working area has not been allocated yet.')
    end if

    call shtnml(now%mm,n,abs(m),lr,li)

    if(m >= 0) then              
      DclGetSpectrumNumber = lr   ! m >= 0 のときは実部の位置
    else
      DclGetSpectrumNumber = li   ! m <  0 のときは虚部の位置
    end if
    call prccls('DclGetSpectrumNumber')

  end function
!-------------------------------------------------
  subroutine DclOperateLaplacian(a,b,ind,idx)   !スペクトルデータに対してラプラシアンを演算する．
    real,intent(in),dimension(*) :: a   ! 入力配列
    real,intent(out),dimension(*) :: b  ! 出力配列
    integer,intent(in),optional :: ind  ! ラプラシアンの演算形式を指定
    integer,intent(in),optional :: idx  ! 作業領域番号
    integer :: idxl, indl

    call prcopn('DclOperateLaplacian')
    idxl = 1
    if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclOperateLaplacian','Working area has not been allocated yet.')
    end if

    indl = 1
    if(present(ind)) indl = ind
    
    call shtlap(now%mm,indl,a,b)
    call prccls('DclOperateLaplacian')
  end subroutine
!-------------------------------------------------
  subroutine DclSpectrumToGrid(s,w,g,isw,idx,m1,m2) ! スペクトルデータからグリッドデータへの変換
    real,intent(in),dimension(*),optional :: s   ! スペクトルデータ
    real,dimension(*) :: w                       ! ウエーブデータ
    real,intent(out),dimension(*),optional :: g  ! グリッドデータ
    integer,intent(in),optional :: isw           ! 変換の種類の指定
    integer,intent(in),optional :: idx           ! 作業領域番号
    integer,intent(in),optional :: m1, m2        ! 波数区間の最小最大値
    integer :: idxl, m1l, m2l, iswl

    call prcopn('DclSpectrumToGrid')
    idxl = 1; if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclSpectrumToGrid','Working area has not been allocated yet.')
    end if

    m1l = 0; if(present(m1)) m1l = m1
    m2l = now%mm; if(present(m2)) m2l = m2
    iswl = 0; if(present(isw)) iswl = isw

    if(present(s)) then
      if(present(g)) then
        call shtsga(now%mm,now%jm,now%im,iswl,m1l,m2l,s,w,g,now%array)
      else
        call shtswa(now%mm,now%jm,iswl,m1l,m2l,s,w,now%array)
      end if
    else
      if(present(g)) then
        call shtwga(now%mm,now%jm,now%im,m1l,m2l,w,g,now%array)
      else
        call msgdmp('E','DclSpectrumToGrid','Either S or G must be specified.')
      end if
    end if
    call prccls('DclSpectrumToGrid')
    
  end subroutine
!---------------------------------------------------------------------------
  subroutine DclGridToSpectrum(g,w,s,isw,idx) ! グリッドデータからスペクトルデータへの変換
    real,intent(in),dimension(*),optional :: g  ! グリッドデータ
    real,intent(out),dimension(*) :: w          ! ウエーブデータ
    real,intent(out),dimension(*),optional :: s ! スペクトルデータ
    integer,intent(in),optional :: isw          ! 変換の種類の指定    
    integer,intent(in),optional :: idx          ! 作業領域番号
    integer :: idxl, iswl

    call prcopn('DclGridToSpectrum')
    idxl = 1; if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclGridToSpectrum','Working area has not been allocated yet.')
    end if

    iswl = 0; if(present(isw)) iswl = isw
    
    if(present(g)) then
      if(present(s)) then
        call shtg2s(now%mm,now%jm,now%im,iswl,g,w,s,now%array)
      else
        call shtg2w(now%mm,now%jm,now%im,g,w,now%array)
      end if
    else
      if(present(s)) then
        call shtw2s(now%mm,now%jm,iswl,w,s,now%array)
      else
        call msgdmp('E','DclGridToSpectrum','Either G or S must be specified.')
      end if
    end if
    call prccls('DclGridToSpectrum')

  end subroutine
!---------------------------------------------------------------------------
  subroutine DclSpectrumToGridForWave(m,s,wr,wi,g,isw,idx) ! スペクトルデータからグリッドデータへの変換(波数指定)
    integer,intent(in) :: m                     ! 変換する波数
    real,intent(in),dimension(*),optional :: s  ! スペクトルデータ
    real,dimension(*) :: wr                     ! w^m(φ)の実数部分
    real,dimension(*) :: wi                     ! w^m(φ)の虚数部分
    real,intent(out),dimension(*),optional :: g ! グリッドデータ
    integer,intent(in),optional :: isw          ! 変換の種類の指定
    integer,intent(in),optional :: idx          ! 作業領域番号
    integer :: idxl, iswl

    call prcopn('DclSpectrumToGridForWave')
    idxl = 1; if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclSpectrumToGridForWave','Working area has not been allocated yet.')
    end if

    iswl = 0; if(present(isw)) iswl = isw

    if(present(s)) then
      if(present(g)) then
        call shtsgm(now%mm,now%jm,now%im,m,iswl,s,wr,wi,g,now%array)
      else
        call shtswm(now%mm,now%jm,m,iswl,s,wr,wi,now%array)
      end if
    else
      if(present(g)) then
        call shtwgm(now%mm,now%jm,now%im,m,wr,wi,g,now%array)        
      else
        call msgdmp('E', 'DclSpectrumToGridForWave', &
&         'Either S or G must be specified.')
      end if
    end if
    call prccls('DclSpectrumToGridForWave')

  end subroutine
!--------------------------------------------------------------------------
  subroutine DclSpectrumToGridForZonal(s,wz,g,isw,idx) ! スペクトルデータからグリッドデータへの変換(帯状成分)
    real,intent(in),dimension(*),optional :: s  ! スペクトルデータ
    real,dimension(*) :: wz                     ! w^0(φ)
    real,intent(out),dimension(*),optional :: g ! グリッドデータ
    integer,intent(in),optional :: isw          ! 変換の種類の指定
    integer,intent(in),optional :: idx          ! 作業領域番号
    integer :: idxl, iswl

    call prcopn('DclSpectrumToGridForZonal')
    idxl = 1; if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclSpectrumToGridForZonal','Working area has not been allocated yet.')
    end if

    iswl = 0; if(present(isw)) iswl = isw

    if(present(s)) then
      if(present(g)) then
        call shtsgz(now%mm,now%jm,now%im,iswl,s,wz,g,now%array)
      else
        call shtswz(now%mm,now%jm,iswl,s,wz,now%array)
      end if
    else
      if(present(g)) then
        call shtwgz(now%jm,now%im,wz,g)
      else
        call msgdmp('E','DclSpectrumToGridForZonal','Either S or G must be specified.')
      end if
    end if
    call prccls('DclSpectrumToGridForZonal')
    
  end subroutine
!--------------------------------------------------------------------------
  subroutine DclSpectrumToGridForLatitude(j,s,wj,gj,isw,idx,m1,m2) ! スペクトルデータからグリッドデータへの変換(緯度円指定)
    integer,intent(in) :: j                      ! 変換を行う緯度円の指定
    real,intent(in), dimension(*),optional :: s  ! スペクトルデータ
    real,dimension(*) :: wj                      ! w^m(φ_j)
    real,intent(out),dimension(*),optional :: gj ! グリッドデータ
    integer,intent(in),optional :: isw           ! 変換の種類の指定
    integer,intent(in),optional :: idx           ! 作業領域番号
    integer,intent(in),optional :: m1, m2        ! 波数区間の最小最大値
    integer :: idxl, iswl, m1l, m2l

    call prcopn('DclSpectrumToGridForLatitude')
    idxl = 1; if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclSpectrumToGridForLatitude','Working area has not been allocated yet.')
    end if

    iswl = 0; if(present(isw)) iswl = isw

    m1l = 0; if(present(m1)) m1l = m1
    m2l = now%mm; if(present(m2)) m2l = m2

    if(present(s)) then
      if(present(gj)) then
        call shtsgj(now%mm,now%jm,now%im,iswl,j,m1l,m2l,s,wj,gj,now%array)
      else
        call shtswj(now%mm,now%jm,iswl,j,m1l,m2l,s,wj,now%array)
      end if
    else
      if(present(gj)) then
        call shtwgj(now%mm,now%im,m1l,m2l,wj,gj,now%array)
      else
        call msgdmp('E','DclSpectrumToGridForLatitude','Either S or G must be specified.')
      end if
    end if
    call prccls('DclSpectrumToGridForLatitude')

  end subroutine
!------------------------------------------------------------------------
  subroutine DclGetLegendreFunctions(m,fun,idx) ! ルジャンドル陪関数の計算
    integer,intent(in) :: m              ! 帯状波数
    real,intent(out),dimension(*) :: fun ! ルジャンドル陪関数が格納される配列
    integer,intent(in),optional :: idx   ! 作業領域番号
    integer :: idxl

    call prcopn('DclGetLegendreFunctions')
    idxl = 1; if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclGetLegendreFunctions','Working area has not been allocated yet.')
    end if

    call shtfun(now%mm,now%jm,m,fun,now%array)
    call prccls('DclGetLegendreFunctions')
  end subroutine
!-------------------------------------------------
  subroutine DclLegendreTransform_F(m,wm,sm,isw,idx) ! ルジャンドル正変換
    integer,intent(in) :: m             ! 変換を行う帯状波数).
    real,intent(in),dimension(*) :: wm  ! ウエーブデータ
    real,intent(out),dimension(*) :: sm ! スペクトルデータ
    integer,intent(in),optional :: isw  ! 変換の種類の指定
    integer,intent(in),optional :: idx  ! 作業領域番号
    integer :: idxl, iswl

    call prcopn('DclLegendreTransform_F')
    idxl = 1; if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclLegendreTransform_F','Working area has not been allocated yet.')
    end if

    iswl = 0; if(present(isw)) iswl = isw

    call shtlfw(now%mm,now%jm,m,iswl,wm,sm,now%array)
    call prccls('DclLegendreTransform_F')
  end subroutine
!-------------------------------------------------
  subroutine DclLegendreTransform_B(m,sm,wm,isw,idx) ! ルジャンドル逆変換
    integer,intent(in) :: m             ! 変換を行う帯状波数.
    real,intent(in),dimension(*) :: sm  ! スペクトルデータ
    real,intent(out),dimension(*) :: wm ! ウエーブデータ
    integer,intent(in),optional :: isw  ! 変換の種類の指定
    integer,intent(in),optional :: idx  ! 作業領域番号
    integer :: idxl, iswl

    call prcopn('DclLegendreTransform_B')
    idxl = 1; if(present(idx)) idxl = idx

    if(.not. seek(idxl)) then
      call msgdmp('E','DclLegendreTransform_B','Working area has not been allocated yet.')
    end if

    iswl = 0; if(present(isw)) iswl = isw

    call shtlbw(now%mm,now%jm,m,iswl,sm,wm,now%array) 
    call prccls('DclLegendreTransform_B')
  end subroutine
!-------------------------------------------------
end module
