!--
!----------------------------------------------------------------------
! Copyright(C) 2020 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  ua_mpi_module_spectrum_mint
!
!  spml/ua_mpi_module_spectrum モジュールは球面上での 2 次元流体運動を
!  球面調和函数を用いたスペクトル法によって数値計算するための
!  モジュール ua_module_mpi の下部モジュールであり,
!  スペクトル解析計算のための Fortran90 関数を提供する.
!
!  スペクトル計算の方法については doc/ua_module_sjpack.tex を参照のこと.
!  内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる.
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK3/SYPACK のマニュアルを参照されたい.
!
!
!履歴  2020/08/04  竹広真一
!
!++
module ua_mpi_module_spectrum_mint
  !
  !=  ua_mpi_module_spectrum_mint
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id: ua_mpi_module_spectrum.f90 590 2013-08-19 08:48:21Z uwabami $
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/ua_spectrum_module_mint モジュールは球面上での 2 次元流体運動を
  ! 球面調和函数を用いたスペクトル法によって数値計算するための
  ! モジュール ua_module の下部モジュールであり,
  ! スペクトル解析計算のための Fortran90 関数を提供する.
  !
  ! スペクトル計算の方法については doc/ua_module.tex を参照のこと.
  ! 内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる.
  ! スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  ! ついては ISPACK3/SYPACK のマニュアルを参照されたい.
  !
  use mpi
  use w_mpi_module_base_mint, only : nm=>nn, nm_l
  use ua_mpi_module_base_mint, only : nc, nmc, km, kc, nm_lu, nms, MPI_COMM_W, MPI_COMM_H

  implicit none

  private

  public nma_EnergyFromStreamfunc_ua          ! エネルギースペクトル
                                              ! (水平全波数 n, 帯状波数 m 空間)
  public na_EnergyFromStreamfunc_ua           ! エネルギースペクトル
                                              ! (水平全波数 n 空間)
  public nma_EnstrophyFromStreamfunc_ua       ! エンストロフィースペクトル
                                              ! (水平全波数 n, 帯状波数 m 空間)
  public na_EnstrophyFromStreamfunc_ua        ! エンストロフィースペクトル
                                              !  (水平全波数 n 空間)
  public nmb_EnergyFromStreamfunc_wb          ! エネルギースペクトル
                                              ! (水平全波数 n, 帯状波数 m 空間)
  public nb_EnergyFromStreamfunc_wb           ! エネルギースペクトル
                                              ! (水平全波数 n 空間)
  public nmb_EnstrophyFromStreamfunc_wb       ! エンストロフィースペクトル
                                              ! (水平全波数 n, 帯状波数 m 空間)
  public nb_EnstrophyFromStreamfunc_wb        ! エンストロフィースペクトル
                                              !  (水平全波数 n 空間)
  public ua_spectrum_VMiss                    ! 欠損値

  real(8) :: ua_spectrum_VMiss = -999.000     ! 欠損値初期値

contains

  !--------------- エネルギースペクトル計算 (ua) -----------------
  function nma_EnergyFromStreamfunc_ua(ua_Strfunc)
    !
    ! 流線関数のスペクトルデータからエネルギーの球面調和函数成分
    ! (スペクトル)を計算する(1 層用).
    !
    !  * 全波数 n, 帯状波数 m の流線関数のスペクトル成分ψ(n,m) から
    !    エネルギースペクトルは (1/2)n(n+1)|ψ(n,m)|^2 と計算される.
    !
    !  * 全てのエネルギースペクトル成分の和に4πをかけたものが球面上での
    !    全エネルギーに等しい.
    !
    !  * スペクトルの (n,m) 成分は全波数 n, 東西波数 m の実部成分,
    !    (n,-m) 成分は全波数 n, 東西波数 m の虚部成分が格納されている.
    !
    !  * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
    !    欠損値の値はモジュール変数 ua_spectrum_VMiss によって設定できる
    !    (初期値は -999.0)
    !
    real(8), intent(in)   :: ua_Strfunc(:,:)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,-nm:nm,size(ua_Strfunc,2)) &
         :: nma_EnergyFromStreamfunc_ua
    !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間)

    real(8), dimension(0:nm,-nm:nm,size(ua_Strfunc,2)) :: nma_work ! 作業領域

    integer :: l, n, m, nm_ary(2)

    integer :: ierr

    nma_work = 0.0D0
    do l=1,nmc
       ! nm_ary = nm_l(l+nms-1)
       nm_ary = nm_lu(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             nma_Work(n,0,:) = nma_Work(n,0,:) &
                  +  0.5 * n*(n+1) * ua_Strfunc(l,:)**2
          else
             nma_Work(n,m,:) = nma_Work(n,m,:) &
                  +  0.5 * n*(n+1) * ua_Strfunc(l,:)**2
             nma_Work(n,-m,:) = nma_Work(n,-m,:) &
                  +  0.5 * n*(n+1) * ua_Strfunc(l,:)**2
          endif
       endif
    enddo

    CALL MPI_ALLREDUCE(nma_work,nma_EnergyFromStreamfunc_ua,&
         (int(nm)+1)*(2*int(nm)+1)*size(ua_Strfunc,2), &
         MPI_REAL8,MPI_SUM,MPI_COMM_W,IERR)

    do m=1,int(nm)
       do n=0,m-1
          nma_EnergyFromStreamfunc_ua(n,m,:)  = ua_spectrum_Vmiss
          nma_EnergyFromStreamfunc_ua(n,-m,:) = ua_spectrum_Vmiss
       enddo
    enddo

  end function nma_EnergyFromStreamfunc_ua

  function na_EnergyFromStreamfunc_ua(ua_Strfunc)
    !
    ! 流線関数のスペクトルデータから各全波数のエネルギー成分(スペクトル)を
    ! 計算する(1 層用).
    !
    !  * 全波数 n の流線関数のスペクトル成分ψ(n,m) から
    !    エネルギースペクトルはΣ[m=-nm]^nm(1/2)n(n+1)|ψ(n,m)|^2
    !    と計算される.
    !
    !  * 全てのエネルギースペクトル成分の和に 4πをかけたものが
    !    球面上での全エネルギーに等しい.
    !

    real(8), intent(in)      :: ua_Strfunc(:,:)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,size(ua_Strfunc,2)) :: na_EnergyFromStreamfunc_ua
    !(out) エネルギースペクトル (水平全波数 n 空間)

    real(8), dimension(0:nm,size(ua_Strfunc,2)) :: na_Work ! 作業領域

    integer :: l, n, m, nm_ary(2)
    real(8) :: factor

    integer :: ierr

    na_work = 0.0D0
    do l=1,nmc
       ! nm_ary = nm_l(l+nms-1)
       nm_ary = nm_lu(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             factor = 1.0D0
          else
             factor = 2.0D0
          endif
          na_Work(n,:) = na_Work(n,:) &
               + factor * 0.5 * n*(n+1) * ua_Strfunc(l,:)**2
       endif
    enddo

    CALL MPI_ALLREDUCE(na_Work,na_EnergyFromStreamfunc_ua,&
         (int(nm)+1)*size(ua_Strfunc,2), &
         MPI_REAL8,MPI_SUM,MPI_COMM_W,IERR)

  end function na_EnergyFromStreamfunc_ua

  !--------------- エンストロフィースペクトル計算 (ua) -----------------
  function nma_EnstrophyFromStreamfunc_ua(ua_Strfunc)
    !
    ! 流線関数のスペクトルデータからエンストロフィーの球面調和函数成分
    ! (スペクトル)を計算する(1 層用).
    !
    ! * 全波数 n, 帯状波数 m の流線関数のスペクトル成分ψ(n,m) から
    !    エンストロフィースペクトルは (1/2)n^2(n+1)^2|ψ(n,m)|^2 と計算される.
    !
    ! * 全てのエンストロフィースペクトル成分の和に4π/R^2をかけたものが
    !   球面上での全エンストロフィーに等しい. ここで R は球面の半径である.
    !
    ! * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
    !   欠損値の値はモジュール変数 ua_spectrum_VMiss によって設定できる
    !   (初期値は -999.0)
    !
    real(8), intent(in)   :: ua_Strfunc(:,:)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,-nm:nm,size(ua_Strfunc,2)) &
         :: nma_EnstrophyFromStreamfunc_ua
    ! エンストロフィースペクトル (水平全波数 n, 帯状波数 m 空間)

    real(8), dimension(0:nm,-nm:nm,size(ua_Strfunc,2)) :: nma_Work ! 作業領域

    integer :: l, n, m, nm_ary(2)

    integer :: ierr

    nma_work = 0.0D0
    do l=1,nmc
       ! nm_ary = nm_l(l+nms-1)
       nm_ary = nm_lu(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             nma_Work(n,0,:) = nma_Work(n,0,:) &
                  +  0.5 * n**2 * (n+1)**2 * ua_Strfunc(l,:)**2
          else
             nma_Work(n,m,:) = nma_Work(n,m,:) &
                  +  0.5 * n**2 * (n+1)**2 * ua_Strfunc(l,:)**2
             nma_Work(n,-m,:) = nma_Work(n,-m,:) &
                  +  0.5 * n**2 * (n+1)**2 * ua_Strfunc(l,:)**2
          endif
       endif
    enddo

    CALL MPI_ALLREDUCE(nma_Work,nma_EnstrophyFromStreamfunc_ua,&
         (int(nm)+1)*(2*int(nm)+1)*size(ua_Strfunc,2), &
         MPI_REAL8,MPI_SUM,MPI_COMM_W,IERR)

    do m=1,int(nm)
       do n=0,m-1
          nma_EnstrophyFromStreamfunc_ua(n,m,:)  = ua_spectrum_Vmiss
          nma_EnstrophyFromStreamfunc_ua(n,-m,:) = ua_spectrum_Vmiss
       enddo
    enddo

  end function nma_EnstrophyFromStreamfunc_ua

  function na_EnstrophyFromStreamfunc_ua(ua_Strfunc)
    !
    ! 流線関数のスペクトルデータから各全波数のエネルギー成分(スペクトル)を
    ! 計算する(1 層用)
    !
    ! * 全波数 n の流線関数のスペクトル成分ψ(n,m) からエンストロフィー
    !   スペクトルはΣ[m=-nm]^nm(1/2)n^2(n+1)^2|ψ(n,m)|^2 と計算される.
    !
    ! * 全てのエネルギースペクトル成分の和に 4π/R^2 をかけたものが
    !   球面上での全エンストフィーに等しい.
    !
    real(8), intent(in)      :: ua_Strfunc(:,:)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,size(ua_Strfunc,2)) :: na_EnstrophyFromStreamfunc_ua
    !(out) エンストロフィースペクトル(水平全波数 n 空間)

    real(8), dimension(0:nm,size(ua_Strfunc,2))     :: na_Work ! 作業領域

    integer :: l, n, m, nm_ary(2)
    real(8) :: factor

    integer :: ierr

    na_work = 0.0D0
    do l=1,nmc
       ! nm_ary = nm_l(l+nms-1)
       nm_ary = nm_lu(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             factor = 1.0D0
          else
             factor = 2.0D0
          endif
          na_Work(n,:) = na_Work(n,:) &
               + factor * 0.5 * n**2 * (n+1)**2 * ua_Strfunc(l,:)**2
       endif
    enddo

    CALL MPI_ALLREDUCE(na_Work,na_EnstrophyFromStreamfunc_ua,&
         (int(nm)+1)*size(ua_Strfunc,2), &
         MPI_REAL8,MPI_SUM,MPI_COMM_W,IERR)

  end function na_EnstrophyFromStreamfunc_ua

  !--------------- エネルギースペクトル計算 (wb) -----------------
  function nmb_EnergyFromStreamfunc_wb(wb_Strfunc)
    !
    ! 流線関数のスペクトルデータからエネルギーの球面調和函数成分
    ! (スペクトル)を計算する(1 層用).
    !
    !  * 全波数 n, 帯状波数 m の流線関数のスペクトル成分ψ(n,m) から
    !    エネルギースペクトルは (1/2)n(n+1)|ψ(n,m)|^2 と計算される.
    !
    !  * 全てのエネルギースペクトル成分の和に4πをかけたものが球面上での
    !    全エネルギーに等しい.
    !
    !  * スペクトルの (n,m) 成分は全波数 n, 東西波数 m の実部成分,
    !    (n,-m) 成分は全波数 n, 東西波数 m の虚部成分が格納されている.
    !
    !  * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
    !    欠損値の値はモジュール変数 wb_spectrum_VMiss によって設定できる
    !    (初期値は -999.0)
    !
    real(8), intent(in)   :: wb_Strfunc(:,:)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,-nm:nm,size(wb_Strfunc,2)) &
         :: nmb_EnergyFromStreamfunc_wb
    !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間)

    real(8), dimension(0:nm,-nm:nm,size(wb_Strfunc,2)) :: nmb_work ! 作業領域

    integer :: l, n, m, nm_ary(2)

    integer :: ierr

    nmb_work = 0.0D0
    do l=1,nc
       nm_ary = nm_l(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             nmb_Work(n,0,:) = nmb_Work(n,0,:) &
                  +  0.5 * n*(n+1) * wb_Strfunc(l,:)**2
          else
             nmb_Work(n,m,:) = nmb_Work(n,m,:) &
                  +  0.5 * n*(n+1) * wb_Strfunc(l,:)**2
             nmb_Work(n,-m,:) = nmb_Work(n,-m,:) &
                  +  0.5 * n*(n+1) * wb_Strfunc(l,:)**2
          endif
       endif
    enddo

    CALL MPI_ALLREDUCE(nmb_work,nmb_EnergyFromStreamfunc_wb,&
         (int(nm)+1)*(2*int(nm)+1)*size(wb_Strfunc,2), &
         MPI_REAL8,MPI_SUM,MPI_COMM_H,IERR)

    do m=1,nm
       do n=0,m-1
          nmb_EnergyFromStreamfunc_wb(n,m,:)  = ua_spectrum_Vmiss
          nmb_EnergyFromStreamfunc_wb(n,-m,:) = ua_spectrum_Vmiss
       enddo
    enddo

  end function nmb_EnergyFromStreamfunc_wb

  function nb_EnergyFromStreamfunc_wb(wb_Strfunc)
    !
    ! 流線関数のスペクトルデータから各全波数のエネルギー成分(スペクトル)を
    ! 計算する(1 層用).
    !
    !  * 全波数 n の流線関数のスペクトル成分ψ(n,m) から
    !    エネルギースペクトルはΣ[m=-nm]^nm(1/2)n(n+1)|ψ(n,m)|^2
    !    と計算される.
    !
    !  * 全てのエネルギースペクトル成分の和に 4πをかけたものが
    !    球面上での全エネルギーに等しい.
    !

    real(8), intent(in)      :: wb_Strfunc(:,:)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,size(wb_Strfunc,2)) :: nb_EnergyFromStreamfunc_wb
    !(out) エネルギースペクトル (水平全波数 n 空間)

    real(8), dimension(0:nm,size(wb_Strfunc,2)) :: nb_Work ! 作業領域

    integer :: l, n, m, nm_ary(2)
    real(8) :: factor

    integer :: ierr

    nb_work = 0.0D0
    do l=1,nc
       nm_ary = nm_l(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             factor = 1.0D0
          else
             factor = 2.0D0
          endif
          nb_Work(n,:) = nb_Work(n,:) &
               + factor * 0.5 * n*(n+1) * wb_Strfunc(l,:)**2
       endif
    enddo

    CALL MPI_ALLREDUCE(nb_Work,nb_EnergyFromStreamfunc_wb,&
         (int(nm)+1)*size(wb_Strfunc,2), &
         MPI_REAL8,MPI_SUM,MPI_COMM_H,IERR)

  end function nb_EnergyFromStreamfunc_wb

  !--------------- エンストロフィースペクトル計算 (wb) -----------------
  function nmb_EnstrophyFromStreamfunc_wb(wb_Strfunc)
    !
    ! 流線関数のスペクトルデータからエンストロフィーの球面調和函数成分
    ! (スペクトル)を計算する(1 層用).
    !
    ! * 全波数 n, 帯状波数 m の流線関数のスペクトル成分ψ(n,m) から
    !    エンストロフィースペクトルは (1/2)n^2(n+1)^2|ψ(n,m)|^2 と計算される.
    !
    ! * 全てのエンストロフィースペクトル成分の和に4π/R^2をかけたものが
    !   球面上での全エンストロフィーに等しい. ここで R は球面の半径である.
    !
    ! * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
    !   欠損値の値はモジュール変数 wb_spectrum_VMiss によって設定できる
    !   (初期値は -999.0)
    !
    real(8), intent(in)   :: wb_Strfunc(:,:)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,-nm:nm,size(wb_Strfunc,2)) &
         :: nmb_EnstrophyFromStreamfunc_wb
    ! エンストロフィースペクトル (水平全波数 n, 帯状波数 m 空間)

    real(8), dimension(0:nm,-nm:nm,size(wb_Strfunc,2)) :: nmb_Work ! 作業領域

    integer :: l, n, m, nm_ary(2)

    integer :: ierr

    nmb_work = 0.0D0
    do l=1,nc
       nm_ary = nm_l(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             nmb_Work(n,0,:) = nmb_Work(n,0,:) &
                  +  0.5 * n**2 * (n+1)**2 * wb_Strfunc(l,:)**2
          else
             nmb_Work(n,m,:) = nmb_Work(n,m,:) &
                  +  0.5 * n**2 * (n+1)**2 * wb_Strfunc(l,:)**2
             nmb_Work(n,-m,:) = nmb_Work(n,-m,:) &
                  +  0.5 * n**2 * (n+1)**2 * wb_Strfunc(l,:)**2
          endif
       endif
    enddo

    CALL MPI_ALLREDUCE(nmb_Work,nmb_EnstrophyFromStreamfunc_wb,&
         (int(nm)+1)*(2*int(nm)+1)*size(wb_Strfunc,2), &
         MPI_REAL8,MPI_SUM,MPI_COMM_H,IERR)

    do m=1,int(nm)
       do n=0,m-1
          nmb_EnstrophyFromStreamfunc_wb(n,m,:)  = ua_spectrum_Vmiss
          nmb_EnstrophyFromStreamfunc_wb(n,-m,:) = ua_spectrum_Vmiss
       enddo
    enddo

  end function nmb_EnstrophyFromStreamfunc_wb

  function nb_EnstrophyFromStreamfunc_wb(wb_Strfunc)
    !
    ! 流線関数のスペクトルデータから各全波数のエネルギー成分(スペクトル)を
    ! 計算する(1 層用)
    !
    ! * 全波数 n の流線関数のスペクトル成分ψ(n,m) からエンストロフィー
    !   スペクトルはΣ[m=-nm]^nm(1/2)n^2(n+1)^2|ψ(n,m)|^2 と計算される.
    !
    ! * 全てのエネルギースペクトル成分の和に 4π/R^2 をかけたものが
    !   球面上での全エンストフィーに等しい.
    !
    real(8), intent(in)      :: wb_Strfunc(:,:)
    !(in) 流線関数(スペクトルデータ)

    real(8), dimension(0:nm,size(wb_Strfunc,2)) :: nb_EnstrophyFromStreamfunc_wb
    !(out) エンストロフィースペクトル(水平全波数 n 空間)

    real(8), dimension(0:nm,size(wb_Strfunc,2))     :: nb_Work ! 作業領域

    integer :: l, n, m, nm_ary(2)
    real(8) :: factor

    integer :: ierr

    nb_work = 0.0D0
    do l=1,nc
       nm_ary = nm_l(l)
       n = nm_ary(1) ; m = abs(nm_ary(2))
       if ( n .ge. 0 ) then
          if ( m == 0 ) then
             factor = 1.0D0
          else
             factor = 2.0D0
          endif
          nb_Work(n,:) = nb_Work(n,:) &
               + factor * 0.5 * n**2 * (n+1)**2 * wb_Strfunc(l,:)**2
       endif
    enddo

    CALL MPI_ALLREDUCE(nb_Work,nb_EnstrophyFromStreamfunc_wb,&
         (int(nm)+1)*size(wb_Strfunc,2), &
         MPI_REAL8,MPI_SUM,MPI_COMM_H,IERR)

  end function nb_EnstrophyFromStreamfunc_wb
  
end module ua_mpi_module_spectrum_mint