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

  implicit none

  private
 
  public nma_EnergyFromStreamfunc_wa          ! エネルギースペクトル           
                                            ! (水平全波数 n, 帯状波数 m 空間)
  public na_EnergyFromStreamfunc_wa           ! エネルギースペクトル
                                            ! (水平全波数 n 空間) 
  public nma_EnstrophyFromStreamfunc_wa       ! エンストロフィースペクトル     
                                            ! (水平全波数 n, 帯状波数 m 空間)
  public na_EnstrophyFromStreamfunc_wa        ! エンストロフィースペクトル  
                                            !  (水平全波数 n 空間)
  public wa_spectrum_VMiss                   ! 欠損値

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

  contains

  !--------------- エネルギースペクトル計算 -----------------
    function nma_EnergyFromStreamfunc_wa(wa_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 の配列には欠損値が格納される.
      !    欠損値の値はモジュール変数 wa_spectrum_VMiss によって設定できる
      !    (初期値は -999.0)
      !
      real(8), intent(in)   :: wa_Strfunc(:,:)
      !(in) 流線関数(スペクトルデータ)

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

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

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

      integer :: ierr

      nma_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
               nma_Work(n,0,:) = nma_Work(n,0,:) &
                    +  0.5 * n*(n+1) * wa_Strfunc(l,:)**2
            else
               nma_Work(n,m,:) = nma_Work(n,m,:) &
                    +  0.5 * n*(n+1) * wa_Strfunc(l,:)**2
               nma_Work(n,-m,:) = nma_Work(n,-m,:) &
                    +  0.5 * n*(n+1) * wa_Strfunc(l,:)**2
            endif
         endif
      enddo

      CALL MPI_ALLREDUCE(nma_work,nma_EnergyFromStreamfunc_wa,&
                         (nm+1)*(2*nm+1)*size(wa_Strfunc,2), &
                         MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      do m=1,nm
         do n=0,m-1
            nma_EnergyFromStreamfunc_wa(n,m,:)  = wa_spectrum_Vmiss
            nma_EnergyFromStreamfunc_wa(n,-m,:) = wa_spectrum_Vmiss
         enddo
      enddo

    end function nma_EnergyFromStreamfunc_wa

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

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

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

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

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

      integer :: ierr

      na_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
            na_Work(n,:) = na_Work(n,:) &
                 + factor * 0.5 * n*(n+1) * wa_Strfunc(l,:)**2
         endif
      enddo

      CALL MPI_ALLREDUCE(na_Work,na_EnergyFromStreamfunc_wa,&
                         (nm+1)*size(wa_Strfunc,2), &
                         MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

    end function na_EnergyFromStreamfunc_wa

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

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

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

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

      integer :: ierr

      nma_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
               nma_Work(n,0,:) = nma_Work(n,0,:) &
                    +  0.5 * n**2 * (n+1)**2 * wa_Strfunc(l,:)**2
            else
               nma_Work(n,m,:) = nma_Work(n,m,:) &
                    +  0.5 * n**2 * (n+1)**2 * wa_Strfunc(l,:)**2
               nma_Work(n,-m,:) = nma_Work(n,-m,:) &
                    +  0.5 * n**2 * (n+1)**2 * wa_Strfunc(l,:)**2
            endif
         endif
      enddo

      CALL MPI_ALLREDUCE(nma_Work,nma_EnstrophyFromStreamfunc_wa,&
                         (nm+1)*(2*nm+1)*size(wa_Strfunc,2), &
                         MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

      do m=1,nm
         do n=0,m-1
            nma_EnstrophyFromStreamfunc_wa(n,m,:)  = wa_spectrum_Vmiss
            nma_EnstrophyFromStreamfunc_wa(n,-m,:) = wa_spectrum_Vmiss
         enddo
      enddo

    end function nma_EnstrophyFromStreamfunc_wa

    function na_EnstrophyFromStreamfunc_wa(wa_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)      :: wa_Strfunc(:,:)
      !(in) 流線関数(スペクトルデータ)

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

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

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

      integer :: ierr

      na_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
            na_Work(n,:) = na_Work(n,:) &
                 + factor * 0.5 * n**2 * (n+1)**2 * wa_Strfunc(l,:)**2
         endif
      enddo

      CALL MPI_ALLREDUCE(na_Work,na_EnstrophyFromStreamfunc_wa,&
                         (nm+1)*size(wa_Strfunc,2), &
                         MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,IERR)

    end function na_EnstrophyFromStreamfunc_wa

end module wa_spectrum_mpi_module_supack
