!--
!----------------------------------------------------------------------
!     Copyright (c) 2016 SPMODEL Development Group.
!----------------------------------------------------------------------
!表題  wa_interpolate_mpi_module_supack
!
!  spml/wa_interpolate_mpi_module_supack モジュールは球面上での 2 次元流体運動を
!  球面調和函数を用いたスペクトル法によって数値計算するための 
!  モジュール wa_mpi_module_supack の下部モジュールであり, スペクトル法による
!  補間計算のための Fortran90 関数を提供する. 
!
!  補間計算の方法については doc/wa_module.tex を参照のこと. 
!  内部で ISPACK2 の SUPACK の Fortran77 サブルーチンを呼んでいる. 
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK2/SUPACK のマニュアルを参照されたい.
!
!履歴  2016/03/16  竹広真一 w_interpolate_mpi_module_supack より改造
!
!++
module wa_interpolate_mpi_module_supack
  !
  !=  wa_interpolate_mpi_module_supack
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id: wa_interpolate_mpi_module_supack.f90 590 2013-08-19 08:48:21Z uwabami $
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/wa_interplata_mpi_module_supack モジュールは球面上での 2 次元流体運動を
  ! 球面調和函数を用いたスペクトル法によって数値計算するための 
  ! モジュール wa_module の下部モジュールであり, スペクトル法による
  ! 補間計算のための Fortran90 関数を提供する. 
  !
  ! 内部で ISPACK2 の SUPACK の Fortran77 サブルーチンを呼んでいる. 
  ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  ! ついては ISPACK2/SUPACK のマニュアルを参照されたい.
  !
  use dc_message, only : MessageNotify
  use w_base_mpi_module_supack, only : nm=>nn, nc, l_nm
  use mpi
  implicit none
  private

  public a_Interpolate_wa                        ! 補間関数

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

  interface alpha
     module procedure alpha_array0
  end interface

  interface Pmm
     module procedure Pmm_array0
  end interface

  contains

  !--------------- 補間計算 -----------------
    function a_Interpolate_array00_wa(wa_data,alon,alat)
      !
      ! 緯度 alat, 経度 alon における関数値を
      ! その球面調和変換係数 wa_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, AnmPnm
      real(8),  dimension(size(wa_data,2)) :: gdata
      logical :: LFlagAnd, LFlagOr
      integer :: k,m
      integer :: ierr
      logical,parameter :: wdata_check = .false.

      mu = sin(alat)
      gdata = 0.0D0

      if ( wdata_check ) then
         !---- 波数チェック
         LFlagAnd = .true. ; LFlagOr  = .false.
         do k=nm,0,-1
            LFlagAnd = (l_nm(k,0) /= 0) .and. LFlagAnd
            LFlagOr = (l_nm(k,0) /= 0) .or. LFlagOr
         enddo
         if ( (.not. LFlagAnd) .and. LFlagOr ) then
            call MessageNotify &
                 ('E','wa_Interpolate_array00','Invalid spectral data for m=0')
         endif

         do m=1,nm
            LFlagAnd = .true. ; LFlagOr  = .false.
            do k=nm,m,-1
               LFlagAnd = (l_nm(k,m) /= 0) .and. LFlagAnd
               LFlagOr = (l_nm(k,m) /= 0) .or. LFlagOr
            enddo
            if ( (.not. LFlagAnd) .and. LFlagOr ) then
               call MessageNotify &
                    ('E','wa_Interpolate_array00','Invalid spectral data for m nonzedo')
            endif
         end do

         do m=1,nm
            LFlagAnd = .true. ; LFlagOr  = .false.
            do k=nm,m,-1
               LFlagAnd = (l_nm(k,-m) /= 0) .and. LFlagAnd
               LFlagOr = (l_nm(k,-m) /= 0) .or. LFlagOr
            enddo
            if ( (.not. LFlagAnd) .and. LFlagOr ) then
               call MessageNotify &
                    ('E','wa_Interpolate_array00','Invalid spectral data for m nonzedo')
            endif
         end do
      endif

      !---- Σa_n^0 P_n^0 の計算
      if ( l_nm(0,0) /= 0 ) then
         y2 = 0 ; y1 = 0
         do k=nm,1,-1
            y0 = alpha(k,0,mu) * y1 + beta(k+1,0)*y2 + wa_data(l_nm(k,0),:)
            y2 = y1 ; y1 = y0
         enddo
         gdata = gdata + (  beta(1,0) * y2 + mu*sqrt(3.0D0) * y1 &
                          + wa_data(l_nm(0,0),:)  ) * Pmm(0,mu)
      endif
         
      !----  Σ Re[s_n^m] P_n^m exp(imλ)の計算
      do m=1,nm
         if ( l_nm(m,m) /= 0 ) then
            y2 = 0 ; y1 = 0
            do k=nm,m+1,-1
               y0 = alpha(k,m,mu) * y1 + beta(k+1,m) * y2 + wa_data(l_nm(k,m),:)
               y2 = y1 ; y1 = y0
            enddo

            AnmPnm =(wa_data(l_nm(m,m),:) + beta(m+1,m)*y2 &
                       + mu*sqrt(2.0D0*m+3)*y1 ) * Pmm(m,mu)

            gdata = gdata + AnmPnm * 2 * cos(m*alon)
         endif
      end do

      !----  Σ Im[s_n^m] P_n^m exp(imλ)の計算
      do m=1,nm
         if ( l_nm(m,-m) /= 0 ) then
            y2 = 0 ; y1 = 0
            do k=nm,m+1,-1
               y0 = alpha(k,m,mu) * y1 + beta(k+1,m)*y2 + wa_data(l_nm(k,-m),:)
               y2 = y1 ; y1 = y0
            enddo

            AnmPnm =(wa_data(l_nm(m,-m),:) + beta(m+1,m)*y2 &
                      + mu*sqrt(2.0D0*m+3)*y1 ) * Pmm(m,mu)


            gdata = gdata - AnmPnm * 2 * sin(m*alon)
         endif
      end do

      CALL MPI_ALLREDUCE(gdata,a_Interpolate_array00_wa,size(wa_data,2),&
                         MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)
      
    end function a_Interpolate_array00_wa

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

      alpha_array0 = sqrt( (2.0D0*n+3)*(2.0D0*n+1)/((n-m+1)*(n+m+1)) ) * x
    end function alpha_array0

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

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

    function Pmm_array0(m,x)
      !
      ! ルジャンドル陪函数 P_m^m(x) の計算
      !         sqrt( factrl(2*m+1) )/(2.0D0**m * factrl(m)) &
      !          * (1-x**2)**(m/2.0D0)
      !
      ! 係数が発散しないように対数で計算する
      !
      integer, intent(IN) :: m           ! ルジャンドル陪函数の次数
      real(8), intent(IN) :: x           ! 引き数
      real(8)             :: Pmm_array0  ! ルジャンドル陪函数の値
      real(8)             :: gammaln     ! Γ関数の対数

      if ( m < 0 ) call MessageNotify('E','Pmm in wa_Intepolate_module',&
                                'order m should be larger equal to zero')

      Pmm_array0 = 1.0
      if ( m > 0 )then
            Pmm_array0 = &
                 exp(gammaln(2.0D0*m+2.0)/2.0D0 - m*log(2.0D0) - gammaln(m+1.0D0)) &
                 * (1-x**2)**(m/2.0D0)
      endif
    end function Pmm_array0

end module wa_interpolate_mpi_module_supack
