! -*- mode: f90; coding: utf-8 -*-
!-------------------------------------------------------------------------
! Copyright (c) 2019 SPMODEL Development Group. All rights reserved.
!
! 履歴  2019/03/16  竹広真一
!       2019/10/04  佐々木洋平
!
!----------------------------------------------------------------------
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @copyright Copyright (C) SPMODEL Development Group, 2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @brief
!> @ref w_module の下部モジュール: スペクトル解析
!> @note
!> スペクトル計算の方法については doc/w_module.tex を参照のこと.
!> @warning
!> 本モジュールを直接呼ぶ事は想定されていないことに注意されたい．
!> ユーザは w_module を呼ぶこと．
!>
module w_module_spectrum_mint
  use dc_types, only: DP
  use w_module_base_mint, only : mm, nn, mi, l_nm

  implicit none

  private

  public nm_EnergyFromStreamfunc_w
  public n_EnergyFromStreamfunc_w
  public nm_EnstrophyFromStreamfunc_w
  public n_EnstrophyFromStreamfunc_w
  public w_spectrum_VMiss

  !> 欠損値. 初期値は -999.000
  real(DP) :: w_spectrum_VMiss = -999.000

contains

  !----------------------------------------------------------
  !> 流線関数のスペクトルデータからエネルギーの球面調和函数成分
  !> (スペクトル)を計算する(1 層用).
  !>
  !>  * 全波数 n, 帯状波数 m の流線関数のスペクトル成分\f$ψ(n,m)\f$
  !>    からエネルギースペクトルは
  !>    \f$(1/2)n(n+1)|ψ(n,m)|^2\f$ と計算される.
  !>
  !>  * 全てのエネルギースペクトル成分の和に\f$4π\f$をかけたものが
  !>    球面上での全エネルギーに等しい.
  !>
  !>  * スペクトルの (n,m) 成分は全波数 n, 東西波数 m の実部成分,
  !>    (n,-m) 成分は全波数 n, 東西波数 m の虚部成分が格納されている.
  !>
  !>  * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納
  !>    される. 欠損値の値はモジュール変数
  !>    [w_spectrum_VMiss](@ref w_spectrum_vmiss) によって
  !>    設定できる(初期値は -999.0)
  !>
  function nm_EnergyFromStreamfunc_w(w_Strfunc)
    !> 流線関数(スペクトルデータ)
    real(DP), intent(in)   :: w_Strfunc(:)
    !> エネルギースペクトル(水平全波数 n, 帯状波数 m 空間)
    real(DP), dimension(0:nn,-mm:mm) :: nm_EnergyFromStreamfunc_w

    integer(8) :: n
    integer(8) :: m

    nm_EnergyFromStreamfunc_w = w_spectrum_VMiss

    do n=0,nn
       nm_EnergyFromStreamfunc_w(n,0) &
           & = 0.5D0 * n*(n+1) * w_Strfunc(l_nm(n,0_8))**2
    end do

    do m=1,mm
       if ( mod(m,mi).eq.0) then
          do n=m,nn
             nm_EnergyFromStreamfunc_w(n,m)        &
                  = 0.5D0 * n*(n+1)                &
                  * (   w_Strfunc(l_nm(n,m))**2    &
                  &   + w_Strfunc(l_nm(n,-m))**2 )
             nm_EnergyFromStreamfunc_w(n,-m)       &
                  = nm_EnergyFromStreamfunc_w(n,m)
          end do
       else
          do n=m,nn
             nm_EnergyFromStreamfunc_w(n,m) = 0.0D0
             nm_EnergyFromStreamfunc_w(n,-m) = 0.0D0
          enddo
       endif
    end do
  end function nm_EnergyFromStreamfunc_w

  !----------------------------------------------------------
  !> 流線関数のスペクトルデータから各全波数のエネルギー成分(スペクトル)を
  !> 計算する(1 層用).
  !>
  !>  * 全波数 n の流線関数のスペクトル成分\f$ψ(n,m)\f$から
  !>    エネルギースペクトルは
  !>    \f[\sum_{m=-nm}^{nm}(1/2)n(n+1)|ψ(n,m)|^2\f]
  !>    と計算される.
  !>
  !>  * 全てのエネルギースペクトル成分の和に \f$4π\f$ をかけたものが
  !>    球面上での全エネルギーに等しい.
  !>
  function n_EnergyFromStreamfunc_w(w_Strfunc)
    !> 流線関数(スペクトルデータ)
    real(DP), intent(in)      :: w_Strfunc(:)
    !> エネルギースペクトル (水平全波数 n 空間)
    real(DP), dimension(0:nn) :: n_EnergyFromStreamfunc_w

    integer(8) :: n
    integer(8) :: m

    n_EnergyFromStreamfunc_w = 0.0D0
    do n=0,nn
       n_EnergyFromStreamfunc_w(n) = &
            0.5D0 * n*(n+1_8) * w_StrFunc(l_nm(n,0_8))**2_8
    enddo

    do m=mi,mm,mi
       do n=m,nn
          n_EnergyFromStreamfunc_w(n) =    &
               & n_EnergyFromStreamfunc_w(n)  &
               !!           & + 2 * 0.5 * n*(n+1) &
               & + n*(n+1)  &
               & * (w_StrFunc(l_nm(n,m))**2_8 + w_StrFunc(l_nm(n,-m))**2_8)
       end do
    end do
  end function n_EnergyFromStreamfunc_w

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

      integer(8) :: n
      integer(8) :: m

      nm_EnstrophyFromStreamfunc_w = w_spectrum_VMiss

      do n=0,nn
         nm_EnstrophyFromStreamfunc_w(n,0) &
              &  = 0.5D0 * n**2_8 * (n+1_8)**2_8 * w_Strfunc(l_nm(n,0_8))**2_8
      enddo
      
      do m=1,mm
         if ( mod(m,mi).eq.0) then
            do n=m,nn
               nm_EnstrophyFromStreamfunc_w(n,0) &
                    &  = 0.5D0 * n**2_8 * (n+1_8)**2_8 * w_Strfunc(l_nm(n,0_8))**2_8
               nm_EnstrophyFromStreamfunc_w(n,m)    &
                    &  = 0.5D0 * n**2_8 * (n+1_8)**2_8 &
                 &    * ( w_Strfunc(l_nm(n,m))**2_8 +w_Strfunc(l_nm(n,-m))**2_8 )
               nm_EnstrophyFromStreamfunc_w(n,-m)   &
                    &  = nm_EnstrophyFromStreamfunc_w(n,m)
            enddo
         else
            do n=m,nn
               nm_EnstrophyFromStreamfunc_w(n,m) = 0.0D0
               nm_EnstrophyFromStreamfunc_w(n,-m) = 0.0D0
            enddo
         endif
      enddo
    end function nm_EnstrophyFromStreamfunc_w

    !--------------------------------------------------------------------
    !> 流線関数のスペクトルデータから各全波数のエネルギー成分(スペクトル)を
    !> 計算する(1 層用)
    !>
    !> * 全波数 n の流線関数のスペクトル成分\f$ψ(n,m)\f$から
    !>   エンストロフィースペクトルは
    !>   \f[\sum{m=-nm}^{nm} (1/2)n^2(n+1)^2|ψ(n,m)|^2\f]
    !>   と計算される.
    !>
    !> * 全てのエネルギースペクトル成分の和に \f$4π/R^2\f$
    !>   をかけたものが球面上での全エンストフィーに等しい.
    !>   ここで\f$R\f$ は球面の半径である.
    !>
    function n_EnstrophyFromStreamfunc_w(w_Strfunc)
      !> 流線関数(スペクトルデータ)
      real(DP), intent(in)      :: w_Strfunc(:)
      !> エンストロフィースペクトル(水平全波数 n 空間)
      real(DP), dimension(0:nn) :: n_EnstrophyFromStreamfunc_w

      integer(8) :: n
      integer(8) :: m

      do n=0,nn
         n_EnstrophyFromStreamfunc_w(n)  &
              &  = 0.5D0 * n**2_8 * (n+1_8)**2_8 * w_StrFunc(l_nm(n,0_8))**2_8
      enddo
      do m=mi,mm,mi
         do n=m,nn
            n_EnstrophyFromStreamfunc_w(n) &
                 & =   n_EnstrophyFromStreamfunc_w(n) &
                 !!             &   + 2 * 0.5D0 * n**2_8 * (n+1)**2_8 &
                 &   + n**2_8 * (n+1)**2_8 &
                 &      * (w_StrFunc(l_nm(n,m))**2_8 + w_StrFunc(l_nm(n,-m))**2_8)
         end do
      end do
    end function n_EnstrophyFromStreamfunc_w

end module w_module_spectrum_mint
