!--
!----------------------------------------------------------------------
!     Copyright (c) 2024 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  wa_zonal_module_interpolate
!
!  spml/wa_zonal_module_interpolate モジュールは球面上での流体運動を
!  球面調和函数を用いたスペクトル法によって東西 1 波数数値計算するための
!  モジュール wa_zonal_module の下部モジュールであり, スペクトル法による
!  補間計算のための Fortran90 関数を提供する. 
!
!  球面上の 1 層モデル用 w_zonal_module_interpolate モジュールを
!  多層モデル用に拡張したものであり, 同時に複数個のスペクトルデータを
!  用いた補間計算が行える.
!
!  補間計算の方法については doc/w_module.tex を参照のこと. 
!  このサブルーチンは内部で ISPACK3 の LXACK のサブルーチンを
!  呼んでいない. 
!
!
!履歴  2024/02/17  竹広真一 wa_wave_module_interpolate より改造
!
!++
module wa_zonal_module_interpolate
  !
  != wa_zonal_module_interpolate
  !
  !== 概要
  !
  !  spml/wa_zonal_module_interpolate モジュールは球面上での流体運動を
  !  球面調和函数を用いたスペクトル法によって東西 1 波数数値計算するための
  !  モジュール wa_zonal_module の下部モジュールであり, スペクトル法による
  !  補間計算のための Fortran90 関数を提供する. 
  !
  !  球面上の 1 層モデル用 w_zonal_module_interpolate モジュールを
  !  多層モデル用に拡張したものであり, 同時に複数個のスペクトルデータを
  !  用いた補間計算が行える.
  !
  !  補間計算の方法については doc/w_module.tex を参照のこと. 
  !  このサブルーチンは内部で ISPACK3 の LXACK のサブルーチンを
  !  呼んでいない. 
  !
  use dc_message, only : MessageNotify
  use spml_utils, only : gammaln
  use w_zonal_module_base, only : nm=>nn, l_nm
  implicit none
  private

  public a_Interpolate_wa                        ! 補間関数

  interface a_Interpolate_wa
     !
     ! 緯度 alon, 経度 alat における関数値を
     ! 球面調和変換係数 wa_data から補間計算する
     !
     ! 入力する緯度経度座標は１点の 1 種類
     !
     module procedure a_Interpolate_array00_wa
  end interface

contains

  !--------------- 補間計算 -----------------    
  function a_Interpolate_array00_wa(wa_data,alon,alat)
    !
    ! 緯度 alat, 経度 alon における関数値を
    ! その球面調和変換係数 w_data から補間計算する
    !
    real(8), intent(IN) :: wa_data(:,:)             ! スペクトルデータ
    real(8), intent(IN) :: alon                     ! 補間する位置(経度)
    real(8), intent(IN) :: alat                     ! 補間する位置(緯度)
    real(8)             :: a_Interpolate_array00_wa(size(wa_data,2)) ! 補間した値

    real(8) :: mu
    real(8), dimension(size(wa_data,2)) :: y0, y1, y2
    integer :: k

    mu = sin(alat)
    a_Interpolate_array00_wa = 0.0D0

    !---- Σa_n^0 L_n^0 の計算
    y2 = 0 ; y1 = 0
    do k=nm,1,-1
       y0 = alpha(k,mu) * y1 + beta(k+1)*y2 + wa_data(l_nm(k,0),:)
       y2 = y1 ; y1 = y0
    enddo
    a_Interpolate_array00_wa = beta(1) * y2 + mu*sqrt(3.0D0) * y1 + wa_data(l_nm(0,0),:) 

  end function a_Interpolate_array00_wa

  !--------------- 下部ルーチン -----------------
  function alpha(n,x)
    !
    !  漸化式の P_n の係数
    !
    integer, intent(IN) :: n
    real(8), intent(IN) :: x
    real(8)             :: alpha

    alpha = sqrt( (2.0D0*n+3)*(2.0D0*n+1)/((n+1)*(n+1)) ) * x
  end function alpha

  function beta(n)
    !
    !  漸化式の P_{n-1} の係数
    !
    integer, intent(IN) :: n
    real(8)             :: beta

    beta = - sqrt( (2.0D0*n+3)*n*n/((2*n-1)*(n+1)*(n+1)) )
  end function beta

end module wa_zonal_module_interpolate