!-------------------------------------------------
!  SHTlib Module
!-------------------------------------------------
module shtlib
  use dcl_common

  type work
    integer                     :: mm,jm,im
    real, dimension(:), pointer :: array
  end type work

  type(work), dimension(10) :: wk
  private :: work, wk
  contains
!-------------------------------------------------
  subroutine DclInitSHT(mm,jm,im,index)        !初期化ルーチン． 
    integer,   intent(in)           :: mm      !入力．切断波数
    integer,   intent(in)           :: jm      !入力．南北分割数の1/2
    integer,   intent(in)           :: im      !入力．東西分割数の1/2
    integer,   intent(in), optional :: index   !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitSHT', &
                   & 'Working area has been allocated already.')
    wk(idx)%mm = mm
    wk(idx)%jm = jm
    wk(idx)%im = im
    len = (JM+1)*(4*JM+5*MM+14) + (MM+1)*(MM+1) +MM + 2 + 6*IM + 15
    allocate(wk(idx)%array(len))
    call shtint(mm,jm,im,wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocSHT(index)           ! 作業領域解放
    integer, intent(in), optional :: index  ! 作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if
    deallocate(wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclGetSpectrumNumber(mm,n,m,lr,li) !スペクトルデータの格納位置を求める．
    integer,   intent(in) :: mm                 !入力．切断波数
    integer,   intent(in) :: n                  !入力．全波数
    integer,   intent(in) :: m                  !入力．帯状波数(
    integer,   intent(out) :: lr                !出力．
    integer,   intent(out) :: li                !出力．
    call shtnml(mm,n,m,lr,li)
  end subroutine
!-------------------------------------------------
  subroutine DclOperateLaplacian(mm,ind,a,b)    !スペクトルデータに対してラプラシアンを演算する．
    integer,   intent(in) :: mm                 !入力．切断波数(m).
    integer,   intent(in) :: ind                !入力．ラプラシアンの演算形式を指定する
    real,      intent(in), dimension(*) :: a    !入力．
    real,      intent(out), dimension(*) :: b   !出力．
    call shtlap(mm,ind,a,b)
  end subroutine
!-------------------------------------------------
  subroutine DclGridToWave(g,w,index)   !グリッドデータからウエーブデータへの変換
    real,    intent(in),  dimension(*) :: g      !入力．変換の種類の指定
    real,    intent(out), dimension(*) :: w      !出力．ウエーブデータ
    integer, intent(in),  optional     :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    call shtg2w(wk(idx)%mm,wk(idx)%jm,wk(idx)%im,g,w,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclWaveToSpectrum(isw,w,s,index)  !ウエーブデータからスペクトルデータへの変換  
    integer, intent(in)                :: isw  !入力．変換の種類の指定
    real,    intent(in),  dimension(*) :: w    !入力．ウエーブデータ
    real,    intent(out), dimension(*) :: s    !出力．スペクトルデータ
    integer, intent(in),  optional     :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    call shtw2s(wk(idx)%mm,wk(idx)%jm,isw,w,s,wk(idx)%array) 
  end subroutine
!-------------------------------------------------
  subroutine DclGridToSpectrum(isw,g,w,s,index)  !グリッドデータからスペクトルデータへの変換
    integer, intent(in)                :: isw  !入力．変換の種類の指定
    real,    intent(in),  dimension(*) :: g    !入力．グリッドデータ
    real,    intent(out), dimension(*) :: w    !出力．ウエーブデータ
    real,    intent(out), dimension(*) :: s    !出力．スペクトルデータ
    integer, intent(in),  optional     :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    call shtg2s(wk(idx)%mm,wk(idx)%jm,wk(idx)%im,isw,g,w,s,wk(idx)%array)
  end subroutine
!---------------------------------------------------------- 
!下位ルーチン     
  subroutine DclSpectrumToWave(isw,s,w,min,max,index)   !SHTS2Wの下位ルーチン
    integer, intent(in)                :: isw       !入力．変換の種類の指定
    real,    intent(in),  dimension(*) :: s         !入力．スペクトルデータ
    real,    intent(out), dimension(*) :: w         !出力．ウエーブデータ
    integer, intent(in),  optional     :: min, max  !入力．波数区間の最小最大値
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1        ; end if

    if(present(min)) then; m1 = min
                     else; m1 = 0           ; end if

    if(present(max)) then; m2 = max
                     else; m1 = wk(idx)%mm  ; end if
                       
    call shtswa(wk(idx)%mm,wk(idx)%jm,isw,m1,m2,s,w,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclWaveToGrid(w,g,min,max,index)    !SHTW2Gの下位ルーチン   
    real,    intent(in),  dimension(*) :: w         !入力．ウエーブデータ
    real,    intent(out), dimension(*) :: g         !出力．グリッドデータ
    integer, intent(in),  optional     :: min, max  !入力．波数区間の最小最大値
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    if(present(min)) then; m1 = min
                     else; m1 = 0           ; end if

    if(present(max)) then; m2 = max
                     else; m1 = wk(idx)%mm  ; end if
                       
    call shtwga(wk(idx)%mm,wk(idx)%jm,wk(idx)%im,m1,m2,w,g,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclSpectrumToGrid(isw,s,w,g,min,max,index) !SHTS2Gの下位ルーチン．
    integer, intent(in)                :: isw !入力．変換の種類の指定
    real,    intent(in),  dimension(*) :: s   !入力．スペクトルデータ
    real,    intent(out), dimension(*) :: w   !出力．ウエーブデータ
    real,    intent(out), dimension(*) :: g   !出力．グリッドデータ

    integer, intent(in),  optional     :: min, max  !入力．波数区間の最小最大値
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

    if(present(min))   then; m1 = min
                       else; m1 = 0           ; end if

    if(present(max))   then; m2 = max
                       else; m1 = wk(idx)%mm  ; end if
                       
    call shtsga(wk(idx)%mm,wk(idx)%jm,wk(idx)%im,isw,m1,m2,s,w,g,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclSpectrumToWaveForWave(m,isw,s,wr,wi,index)   !SHTS2Wの下位ルーチン 
    integer, intent(in)                :: m     !入力．変換する波数
    integer, intent(in)                :: isw   !入力．変換の種類の指定
    real,    intent(in),  dimension(*) :: s     !入力．スペクトルデータ
    real,    intent(out), dimension(*) :: wr    !出力． w^m（φ）の実数部分
    real,    intent(out), dimension(*) :: wi    !出力． w^m（φ）の虚数部分

    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

    call shtswm(wk(idx)%mm,wk(idx)%jm,m,isw,s,wr,wi,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclWaveToGridForWave(m,wr,wi,g,index)    !SHTW2Gの下位ルーチン 
    integer, intent(in)                :: m    !入力．変換する波数
    real,    intent(in),  dimension(*) :: wr   !入力． w^m（φ）の実数部分
    real,    intent(in),  dimension(*) :: wi   !入力． w^m（φ）の虚数部分
    real,    intent(out), dimension(*) :: g    !出力．グリッドデータ
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

    call shtwgm(wk(idx)%mm,wk(idx)%jm,wk(idx)%im,m,wr,wi,g,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclSpectrumToGridForWave(m,isw,s,wr,wi,g,index) !SHTS2Gの下位ルーチ
    integer, intent(in)                :: m    !入力．変換する波数
    integer, intent(in)                :: isw  !入力．変換の種類の指定
    real,    intent(in),  dimension(*) :: s    !入力．スペクトルデータ
    real,    intent(out), dimension(*) :: wr   !出力． w^m（φ）の実数部分
    real,    intent(out), dimension(*) :: wi   !出力． w^m（φ）の虚数部分
    real,    intent(out), dimension(*) :: g    !出力．グリッドデータ
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

    call shtsgm(wk(idx)%mm,wk(idx)%jm,wk(idx)%im,m,isw,s,wr,wi,g,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclSpectrumToWaveForZonal(isw,s,wz,index)        !SHTS2Wの下位ルーチン
    integer, intent(in)                :: isw  !入力．変換の種類の指定
    real,    intent(in),  dimension(*) :: s    !入力．スペクトルデータ
    real,    intent(out), dimension(*) :: wz   !出力． w^0（φ）が格納される
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if
    call shtswz(wk(idx)%mm,wk(idx)%jm,isw,s,wz,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclWaveToGridForZonal(jm,im,wz,g)  !SHTW2Gの下位ルーチン
    integer,   intent(in) :: jm                 !入力．南北分割数の1/2(j)
    integer,   intent(in) :: im                 !入力．東西分割数の1/2(i)
    real,      intent(in), dimension(*) :: wz   !入力．ウエーブデータ(w^0（φ）)
    real,      intent(out), dimension(*) :: g   !出力．グリッドデータ
    call shtwgz(jm,im,wz,g) 
  end subroutine
!-------------------------------------------------
  subroutine DclSpectrumToGridForZonal(isw,s,wz,g,index)   !SHTS2Gの下位ルーチン s
    integer, intent(in)                :: isw  !入力．変換の種類の指定
    real,    intent(in),  dimension(*) :: s    !入力．スペクトルデータ
    real,    intent(out), dimension(*) :: wz   !出力．ウエーブデータ(w^0（φ）)
    real,    intent(out), dimension(*) :: g    !出力．グリッドデータ
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if
                       
    call shtsgz(wk(idx)%mm,wk(idx)%jm,wk(idx)%im,isw,s,wz,g,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclSpectrumToWaveForLatitude(isw,j,s,wj,min,max,index) !SHTS2Wの下位ルーチン
    integer, intent(in) :: isw                !入力．変換の種類の指定
    integer, intent(in) :: j                  !入力．変換を行う緯度円の指定
    real,    intent(in), dimension(*) :: s    !入力．スペクトルデータ
    real,    intent(out), dimension(*) :: wj  !出力． w^m（φ_j) が格納される． 長さ2*mm+1の配列

    integer, intent(in),  optional     :: min, max  !入力．波数区間の最小最大値
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

    if(present(min))   then; m1 = min
                       else; m1 = 0           ; end if

    if(present(max))   then; m2 = max
                       else; m1 = wk(idx)%mm  ; end if
                       
                       
    call shtswj(wk(idx)%mm,wk(idx)%jm,isw,j,m1,m2,s,wj,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclWaveToGridForLatitude(wj,gj,min,max,index)     !SHTW2Gの下位ルーチン．
    real,    intent(in),  dimension(*) :: wj   !入力． w^m（φ_j)
    real,    intent(out), dimension(*) :: gj   !出力．グリッドデータ．長さ2*im+1の配列

    integer, intent(in),  optional     :: min, max  !入力．波数区間の最小最大値
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

    if(present(min))   then; m1 = min
                       else; m1 = 0           ; end if

    if(present(max))   then; m2 = max
                       else; m1 = wk(idx)%mm  ; end if
                       
    call shtwgj(wk(idx)%mm,wk(idx)%im,m1,m2,wj,gj,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclSpectrumToGridForLatitude(isw,j,s,wj,gj,min,max,index) !SHTS2Gの下位  
    integer, intent(in) :: isw                !入力．変換の種類の指定
    integer, intent(in) :: j                  !入力．変換を行う緯度円の指定
    real,    intent(in), dimension(*) :: s    !入力．スペクトルデータ
    real,    intent(out), dimension(*) :: wj  !出力． w^m（φ_j)
    real,    intent(out), dimension(*) :: gj  !出力．グリッドデータ

    integer, intent(in),  optional     :: min, max  !入力．波数区間の最小最大値
    integer, intent(in),  optional     :: index     !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

    if(present(min))   then; m1 = min
                       else; m1 = 0           ; end if

    if(present(max))   then; m2 = max
                       else; m1 = wk(idx)%mm  ; end if

    call shtsgj(wk(idx)%mm,wk(idx)%jm,wk(idx)%im,isw,j,m1,m2,s,wj,gj,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclGetLegendreFunctions(m,fun,index)           !ルジャンドル陪関数を計算する．
    integer, intent(in)                :: m      !入力．求めるルジャンドル陪関数の帯状波数(m).
    real,    intent(out), dimension(*) :: fun    !p^m_n(φ）が格納される長さ (2*jm+1)*(mm-m+1)の配列
    integer, intent(in),  optional     :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

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

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

    call shtlfw(wk(idx)%mm,wk(idx)%jm,m,isw,wm,sm,wk(idx)%array)
  end subroutine
!-------------------------------------------------
  subroutine DclLegendreTransform_B(m,isw,sm,wm,index)     !ルジャンドル逆変換を行う．
    integer, intent(in)                :: isw  !入力．変換の種類の指定
    integer, intent(in)                :: m    !入力．変換を行う帯状波数(m).
    real,    intent(in),  dimension(*) :: sm   !入力．スペクトルデータ．長さmm-m+1の配列
    real,    intent(out), dimension(*) :: wm   !出力．ウエーブデータ． 長さ(2*jm+1)の配列
    integer, intent(in),  optional     :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1          ; end if

    call shtlbw(wk(idx)%mm,wk(idx)%jm,m,isw,sm,wm,wk(idx)%array) 
  end subroutine
!-------------------------------------------------
end module
