! -*- mode: f90; coding: utf-8 -*-
!-----------------------------------------------------------------------
! Copyright (c) 2019-2021 SPMODEL Development Group. All rights reserved.
!
! 履歴  2019/03/21  竹広真一
!       2019/10/09  佐々木洋平
!       2021/02/18  竹広真一
!
!-----------------------------------------------------------------------
!> @copyright Copyright (C) SPMODEL Development Group, 2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @brief
!> @ref w_mpi_module の下部モジュール: 補間
!>
!> @details
!> 補間計算の方法については doc/w_module.tex を参照のこと.
!> 内部で ISPACK3 の SYPACK を呼んでいる.
!> スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!> ついては ISPACK3/SYPACK のマニュアルを参照されたい.
!>
!> @warning
!> 本モジュールを直接呼ぶ事は想定されていないことに注意されたい．
!> ユーザは w_mpi_module を呼ぶこと．
module w_mpi_module_interpolate_mint
  use dc_message, only : MessageNotify
  use dc_types, only: DP
  use spml_utils, only: gammaln
  use w_mpi_module_base_mint, only : &
    & nm=>nn, mi, nc, l_nm, icom
  use mpi
  implicit none
  private

  public Interpolate_w

  !----------------------------------------------------------
  !>  入力緯度 `alon`, 入力経度 `alat` における関数値を
  !>  その球面調和変換係数 w_data から補間計算する
  !>
  !>  入力を受けつける緯度経度座標は,
  !>  * 1点,
  !>  * 経度1点, 緯度複数点,
  !>  * 経度複数点, 緯度1点,
  !>  * 経度複数点, 緯度複数点
  !>  の 4 種類
  !>
  interface Interpolate_w
    module procedure Interpolate_array00_w
    module procedure Interpolate_array10_w
  end interface Interpolate_w

  interface alpha
    module procedure alpha_array0
  end interface alpha

  interface Pmm
    module procedure Pmm_array0
  end interface Pmm

contains

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

    real(DP) :: mu
    real(DP) :: y0, y1, y2, AnmPnm
    real(DP) :: 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=int(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','w_Interpolate_array00','Invalid spectral data for m=0')
      endif

      do m=int(mi),int(nm),int(mi)
        LFlagAnd = .true. ; LFlagOr  = .false.
        do k=int(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','w_Interpolate_array00',          &
            & 'Invalid spectral data for m nonzedo')
        endif
      end do

      do m=int(mi),int(nm),int(mi)
        LFlagAnd = .true. ; LFlagOr  = .false.
        do k=int(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','w_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=int(nm),1,-1
        y0 = alpha(k,0,mu) * y1 + beta(k+1,0)*y2 + w_data(l_nm(k,0))
        y2 = y1 ; y1 = y0
      enddo
      gdata = gdata + (  beta(1,0) * y2 + mu*sqrt(3.0D0) * y1 &
        + w_data(l_nm(0,0))  ) * Pmm(0,mu)
    endif

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

        AnmPnm =(w_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=int(mi),int(nm),int(mi)
      if ( l_nm(m,-m) /= 0 ) then
        y2 = 0 ; y1 = 0
        do k=int(nm),m+1,-1
          y0 = alpha(k,m,mu) * y1 + beta(k+1,m)*y2 + w_data(l_nm(k,-m))
          y2 = y1 ; y1 = y0
        enddo

        AnmPnm =(w_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,Interpolate_array00_w,1,&
      MPI_REAL8,MPI_SUM,int(icom),IERR)

  end function Interpolate_array00_w

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

    real(DP) :: mu
    real(DP) :: y0, y1, y2, AnmPnm
    real(DP) :: gdata(size(alon))
    logical :: LFlagAnd, LFlagOr
    integer :: i,k,m
    integer :: ierr
    logical,parameter :: wdata_check = .false.

    mu = sin(alat)
    !$omp parallel
    !$omp workshare
    gdata = 0.0D0
    !$omp end workshare
    !$omp end parallel

    if ( wdata_check ) then
      !---- 波数チェック
      LFlagAnd = .true. ; LFlagOr  = .false.
      do k=int(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','w_Interpolate_array00','Invalid spectral data for m=0')
      endif

      do m=int(mi),int(nm),int(mi)
        LFlagAnd = .true. ; LFlagOr  = .false.
        do k=int(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','w_Interpolate_array00','Invalid spectral data for m nonzedo')
        endif
      end do

      do m=int(mi),int(nm),int(mi)
        LFlagAnd = .true. ; LFlagOr  = .false.
        do k=int(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','w_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=int(nm),1,-1
        y0 = alpha(k,0,mu) * y1 + beta(k+1,0)*y2 + w_data(l_nm(k,0))
        y2 = y1 ; y1 = y0
      enddo

      AnmPnm = (  beta(1,0) * y2 + mu*sqrt(3.0D0) * y1 &
        + w_data(l_nm(0,0))  ) * Pmm(0,mu)

      !$omp parallel
      !$omp workshare
      gdata = gdata + AnmPnm
      !$omp end workshare
      !$omp end parallel
    endif

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

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

        !$omp parallel do default(shared)
        do i=1,size(alon)
          gdata(i) = gdata(i) + AnmPnm * 2 * cos(m*alon(i))
        enddo
        !$omp end parallel do
      endif
    end do

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

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

        !$omp parallel do default(shared)
        do i=1,size(alon)
          gdata(i) = gdata(i) - AnmPnm * 2 * sin(m*alon(i))
        enddo
        !$omp end parallel do
      endif
    end do

    CALL MPI_ALLREDUCE(gdata,Interpolate_array10_w,size(alon),&
      MPI_REAL8,MPI_SUM,int(icom),IERR)

  end function Interpolate_array10_w

  !---------------------------------------------------------------------
  !>  漸化式の P_n^m の係数
  function alpha_array0(n,m,x)
    integer, intent(IN) :: n,m
    real(DP), intent(IN) :: x
    real(DP)             :: alpha_array0

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

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

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

  !---------------------------------------------------------------------
  !> ルジャンドル陪函数 P_m^m(x) の計算
  !>
  !> 係数が発散しないように対数で計算する
  function Pmm_array0(m,x)
    integer, intent(IN)  :: m           ! ルジャンドル陪函数の次数
    real(DP), intent(IN) :: x           ! 引き数
    real(DP)             :: Pmm_array0  ! ルジャンドル陪函数の値

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

    Pmm_array0 = 1.0D0
    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 w_mpi_module_interpolate_mint