!= 力学過程 : 水平スペクトル解析
!
!= Dynamical process (horizontal spectral analysis)
!
! Authors::   Shin-ichi Takehiro
! Version::   $Id: dynamics_diagspectrum.f90,v 1.1 2026/03/21 20:00:00 takepiro Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2026. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module dynamics_diag_spectrum
  !
  != 力学過程 : 水平スペクトル解析
  !
  != Dynamical process : horizontal spectral analysis
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 力学過程での水平スペクトル解析を演算するモジュールです. 
  !
  ! This is a horizontal spectral analysis in dynamical process
  !
  !== Procedures List
  !
  ! DynamicsDiagSpAnalysis         :: 水平スペクトル解析
  ! DynamicsDiagSpAnalysisInit     :: 初期化
  ! --------------------------     :: ------------
  ! DynamicsDiagSpAnalysis         :: Horizontal spectral analysis
  ! DynamicsDiagSpAnalysisInit     :: Initialization
  !
  !== NAMELIST
  !
  ! NAMELIST#dynamics_diag_spectrum_nml
  !
  !== References
  !

  ! モジュール引用 ; USE statements
  !


  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &  ! 倍精度実数型. Double precision. 
    &                 TOKEN  ! キーワード.   Keywords. 

  ! メッセージ出力
  ! Message output
  !
  use dc_message, only: MessageNotify

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: DynamicsDiagSpAnalysis
  public :: DynamicsDiagSpAnalysisInit

  ! 非公開変数
  ! Private variables
  !
  logical, save :: dynamics_diag_spectrum_inited = .false.
                              ! 初期設定フラグ. 
                              ! Initialization flag

  character(*), parameter:: module_name = 'dynamics_diag_spectrum'
                              ! モジュールの名称. 
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: dynamics_diag_spectrum.F90,v 1.1 2026/03/21 20:00:00 takepiro Exp $'
                              ! モジュールのバージョン
                              ! Module version
contains

  !--------------------------------------------------------------------------------------

  subroutine DynamicsDiagSpAnalysis( wz_Psi, wz_Chi )
    !
    ! 水平スペクトル解析を実行します.
    ! 運動エネルギースペクトル(渦度, 発散, 全体)と
    ! エンストロフィースペクトル(渦度, 発散, 全体)を計算し
    ! ヒストリー出力します. 
    !
    ! Perform a horizontal spectral analysis.
    ! Calculate the kinetic energy spectra (vorticity, divergence, total) 
    ! and the enstrophy spectra (vorticity, divergence, total).
    ! The results are output to history files.
    !

    ! モジュール引用 ; USE statements
    !

    use constants, only: &
      & RPlanet
                              ! $ a $ [m]. 
                              ! 惑星半径. 
                              ! Radius of planet
    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: nmax, & ! 最大全波数.
                               ! Maximum truncated wavenumber
      &                lmax, & ! スペクトルデータの配列サイズ
                               ! Size of array for spectral data
      &                imax, & ! 経度格子点数. 
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数. 
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数. 
                               ! Number of vertical level
    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $. 
      & TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut
    
    real(DP), intent(IN) :: wz_Psi (lmax, 1:kmax)
                               ! $ \psi $ . 流線関数. Streamline function 
    real(DP), intent(IN) :: wz_Chi (lmax, 1:kmax)
                               ! $ \chi $ . ポテンシャル. Potential

    ! 作業変数
    ! Work variables
    !
    
    real(DP) :: wz_KinEngySpectrum(lmax, 1:kmax)
                              ! 運動エネルギースペクトル
    real(DP) :: wz_KinEngyVorSpectrum(lmax, 1:kmax)
                              ! 運動エネルギースペクトル(渦度)
    real(DP) :: wz_KinEngyDivSpectrum(lmax, 1:kmax)
                              ! 運動エネルギースペクトル(発散)
    real(DP) :: wz_EnstroSpectrum(lmax, 1:kmax)
                              ! エンストロフィースペクトル
    real(DP) :: wz_EnstroVorSpectrum(lmax, 1:kmax)
                              ! エンストロフィー(渦度)スペクトル
    real(DP) :: wz_EnstroDivSpectrum(lmax, 1:kmax)
                              ! エンストロフィー(発散)スペクトル

    ! 実行文 ; Executable statement
    !

    ! 初期化確認
    ! Initialization check
    !
    if ( .not. dynamics_diag_spectrum_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if
    
    ! 計算時間計測開始
    ! Start measurement of computation time
    !
    call TimesetClockStart( module_name )

    ! 運動エネルギー(渦度)スペクトル
    !
    wz_KinEngyVorSpectrum = wa_EnergyFromStreamfunc_wa(wz_Psi/Rplanet)

    ! 運動エネルギー(発散)スペクトル
    wz_KinEngyDivSpectrum = wa_EnergyFromStreamfunc_wa(wz_Chi/Rplanet)

    ! 運動エネルギースペクトル
    wz_KinEngySpectrum = wz_KinEngyVorSpectrum + wz_KinEngyDivSpectrum

    ! エンストロフィー(渦度)スペクトル
    wz_EnstroVorSpectrum = wa_EnstrophyFromStreamfunc_wa(wz_Psi/Rplanet)

    ! エンストロフィー(発散)スペクトル
    wz_EnstroDivSpectrum = wa_EnstrophyFromStreamfunc_wa(wz_Chi/Rplanet)

    ! エンストロフィースペクトル
    wz_EnstroSpectrum = wz_EnstroVorSpectrum + wz_EnstroDivSpectrum 

    ! ヒストリー出力
    ! History output
    !
    call HistoryAutoPut( TimeN, 'KinEngySpVor',wz_KinEngyVorSpectrum)
    call HistoryAutoPut( TimeN, 'KinEngySpDiv',wz_KinEngyDivSpectrum)
    call HistoryAutoPut( TimeN, 'KinEngySp',wz_KinEngySpectrum)

    call HistoryAutoPut( TimeN, 'EnstroSpVor',wz_EnstroVorSpectrum)
    call HistoryAutoPut( TimeN, 'EnstroSpDiv',wz_EnstroDivSpectrum)
    call HistoryAutoPut( TimeN, 'EnstroSp',wz_EnstroSpectrum)

    ! 計算時間計測一時停止
    ! Pause measurement of computation time
    !
    call TimesetClockStop( module_name )
    
  end subroutine DynamicsDiagSpAnalysis
  
  !--------------------------------------------------------------------------------------

  subroutine DynamicsDiagSpAnalysisInit
    !
    ! dynamics_diag_spectrum モジュールの初期化を行います. 
    !
    ! "dynamics_diag_spectrum" module is initialized. 
    !

#ifdef AXISYMMETRY
    use wa_zonal_module, only: nm_l
#elif LIB_MPI
    use ua_mpi_module, only : nm_l => nm_lu
#else
    use wa_module, only: nm_l
#endif

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    !
    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: & 
      &                lmax    ! スペクトルデータの配列サイズ
                               ! Size of array for spectral data
    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
         & AxNameWN, AxNameZ, AxNameT

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN                  ! ステップ $ t $ の時刻. Time of step $ t $. 

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut

    
    integer:: n(lmax)      ! 水平全波数
                           ! Total horizontal wavenumber
    integer:: m(lmax)      ! 経度波数
                           ! longitudinal wavenumber
    integer:: nmary(2)     ! 作業変数(全波数, 経度波数)
                           ! Working variable (total/longitudinal wavenumber)
    integer:: l            ! 水平波数方向に回る DO ループ用作業変数
                           ! Work variables for DO loop in horizontal wavenumber

    ! 実行文 ; Executable statement
    !
    if ( dynamics_diag_spectrum_inited ) return

    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'NSph', &
      & (/ AxNameWN/), &
      & 'total horizontal wavenumber', '1' )
    call HistoryAutoAddVariable( 'MSph', &
      & (/ AxNameWN/), &
      & 'longitudinal wavenumber', '1' )
    do l=1,lmax
       nmary = nm_l(l)
       n(l) = nmary(1) ; m(l) = nmary(2)
    end do
    call HistoryAutoPut( TimeN, 'NSph', n ) 
    call HistoryAutoPut( TimeN, 'MSph', m ) 
    
    call HistoryAutoAddVariable( 'KinEngySpVor', &
      & (/ AxNameWN, AxNameZ, AxNameT /), &
      & 'kinetic energy spectrum (vor)', 'm^2 s-2' )
    call HistoryAutoAddVariable( 'KinEngySpDiv', &
      & (/ AxNameWN, AxNameZ, AxNameT /), &
      & 'kinetic energy spectrum (div)', 'm^2 s-2' )
    call HistoryAutoAddVariable( 'KinEngySp', &
      & (/ AxNameWN, AxNameZ, AxNameT /), &
      & 'kinetic energy spectrum', 'm^2 s-2' )

    call HistoryAutoAddVariable( 'EnstroSpVor', &
      & (/ AxNameWN, AxNameZ, AxNameT /), &
      & 'enstrophy spectrum (vor)', 's-2' )
    call HistoryAutoAddVariable( 'EnstroSpDiv', &
      & (/ AxNameWN, AxNameZ, AxNameT /), &
      & 'enstrophy spectrum (div)', 's-2' )
    call HistoryAutoAddVariable( 'EnstroSp', &
      & (/ AxNameWN, AxNameZ, AxNameT /), &
      & 'enstrophy spectrum', 's-2' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    dynamics_diag_spectrum_inited = .true.
    
  end subroutine DynamicsDiagSpAnalysisInit

  !--------------------------------------------------------------------------------------

  function wa_EnergyFromStreamfunc_wa(wa_Strfunc) result(wa_ESp)
    !
    ! 運動エネルギースペクトル計算
    ! Kinetic energy spectrum calculation

    !
    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: & 
      &                lmax    ! スペクトルデータの配列サイズ
                               ! Size of array for spectral data

#ifdef AXISYMMETRY
    use wa_zonal_module, only: nm_l
#elif LIB_MPI
    use ua_mpi_module, only : nm_l => nm_lu
#else
    use wa_module, only: nm_l
#endif
    
    real(dp), intent(IN) :: wa_Strfunc(:,:)
    real(dp)             :: wa_ESp(size(wa_Strfunc,1),size(wa_Strfunc,2))

    integer :: nmary(2)
    integer :: n, l

    do l=1,lmax
       nmary = nm_l(l)
       n = nmary(1)
       wa_ESp(l,:) = 0.5D0 * n*(n+1) * wa_Strfunc(l,:)**2
    end do
    
  end function wa_EnergyFromStreamfunc_wa

  function wa_EnstrophyFromStreamfunc_wa(wa_Strfunc) result(wa_EnsSp)
    !
    ! エンストフィースペクトル計算
    ! Ensptrophy spectrum calculation

    !
    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: & 
      &                lmax    ! スペクトルデータの配列サイズ
                               ! Size of array for spectral data

#ifdef AXISYMMETRY
    use wa_zonal_module, only: nm_l
#elif LIB_MPI
    use ua_mpi_module, only : nm_l => nm_lu
#else
    use wa_module, only: nm_l
#endif
    
    real(dp), intent(IN) :: wa_Strfunc(:,:)
    real(dp)             :: wa_EnsSp(size(wa_Strfunc,1),size(wa_Strfunc,2))

    integer :: nmary(2)
    integer :: n, l

    do l=1,lmax
       nmary = nm_l(l)
       n = nmary(1)
       wa_EnsSp(l,:) = 0.5D0 * n**2 * (n+1)**2 * wa_Strfunc(l,:)**2
    end do
    
  end function wa_EnstrophyFromStreamfunc_wa

end module dynamics_diag_spectrum
