! -*- mode: f90; coding: utf-8 -*-
!----------------------------------------------------------------------
! Copyright (c) 2019 SPMODEL Development Group. All rights reserved.
!
! 履歴  2016/03/09  竹広真一
!       2019/10/04  佐々木洋平
!       2019/10/09  佐々木洋平
!       2020/11/06  竹広真一   セクター計算オプション導入
!----------------------------------------------------------------------
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @copyright Copyright (C) SPMODEL Development Group, 2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @brief 2次元球面: 球面調和関数展開
!>
!> @details
!> spml/w_module モジュールは球面上での 2 次元流体運動を
!> 球面調和函数を用いたスペクトル法によって数値計算するための
!> Fortran90 関数を提供する.
!> w_module は実際には以下の下部モジュールからなる．
!> * w_module_base: 基本変換
!> * w_module_deriv: 微分計算
!> * w_module_integral: 積分・平均計算
!> * w_module_spectrum: スペクトル解析
!> * w_module_interpolate: 補間計算
!>
!> 内部で ISPACK3 の SXPACKサブルーチンを呼んでいる. スペクトルデータお
!> よび格子点データの格納方法や変換の詳しい計算法については
!> ISPACK3/SXPACK のマニュアルを参照されたい.
!>
!> @warning
!> 後方非互換である点に注意されたい
!> * 旧版と異なり (ISPACK/SXPACKのレベルで) 球面調和関数の正規化が異なる．
!> * 旧版と異なり,
!>   [l_nm](@ref w_module_base#l_nm),
!>   [nm_l](@ref w_module_base#nm_l)
!>   は初期化したのちにしか使うことができない.
!>
!> ## 関数・変数の名前と型について
!>
!> ### 命名法
!>
!>  * 関数名の先頭 (`w_`, `nm_`, `n_`, `xy_`, `x_`, `y_`) は,
!>  返す値の形を示している.
!>    * `w_`  : スペクトルデータ
!>    * `xy_` : 2 次元格子点データ
!>    * `nm_` : スペクトルデータの並んだ 3 次元配列
!>    * `n_`  : スペクトルデータの並んだ 2 次元配列
!>    * `x_`  : 経度方向 1 次元格子点データ
!>    * `y_`  : 緯度方向 1 次元格子点データ
!>
!>  * 関数名の間の文字列(`DLon`, `GradLat`, `GradLat`, `DivLon`, `DivLat`,
!>    `Lapla`, `LaplaInv`, `Jacobian`)は, その関数の作用を表している.
!>
!>  * 関数名の最後 (`_w_w`, `_w`, `_xy`, `_x`, `_y`) は,
!>    入力変数の形スペクトルデータおよび格子点データであることを示している.
!>    * `_w`   :: スペクトルデータ
!>    * `_w_w` :: 2 つのスペクトルデータ
!>    * `_xy`  :: 2 次元格子点データ
!>    * `_x`   :: 経度方向 1 次元格子点データ
!>    * `_y`   :: 緯度方向 1 次元格子点データ
!>
!> ### 各データの種類の説明
!>
!>  * 以下, `DP = kind(1.0D0)`
!>  * `xy` : 2 次元格子点データ.
!>    * 変数の種類と次元は `real(DP), dimension(0:im-1,1:jm)`.
!>    * `im`, `jm` はそれぞれ経度, 緯度座標の格子点数であり, サブルーチン
!>      [w_Initial](@ref w_initial) にてあらかじめ設定しておく.
!>
!>  * `w` : スペクトルデータ.
!>    * 変数の種類と次元は `real(DP), dimension((nm+1)*(nm+1))`.
!>    * nm は球面調和函数の最大全波数であり, サブルーチン w_Initial にて
!>      あらかじめ設定しておく.
!>    * スペクトルデータの格納のされ方は関数
!>      [l_nm](@ref w_module_base#l_nm),
!>      [nm_l](@ref w_module_base#nm_l) によって
!>      調べることができる.
!>
!>  * `nm` : スペクトルデータの並んだ 2 次元配列.
!>    * 変数の種類と次元は `real(DP), dimension(0:nm,-nm:nm)`.
!>      第 1 次元が水平全波数,  第 2 次元が帯状波数を表す.
!>    * `nm` は球面調和函数の最大全波数であり,
!>      サブルーチン [w_Initial](@ref w_initial)
!>      にてあらかじめ設定しておく.
!>
!>  * `n` : スペクトルデータの並んだ 1 次元配列.
!>    * 変数の種類と次元は `real(DP), dimension(0:nm)`.
!>    * 第 1 次元が水平全波数を表す. nm は球面調和函数の最大全波数であり,
!>      サブルーチン [w_Initial](@ref w_initial)
!>      にてあらかじめ設定しておく.
!>
!>  * `x`, `y` : 経度, 緯度方向 1 次元格子点データ.
!>    * 変数の種類と次元はそれぞれ `real(DP), dimension(0:im-1)`
!>      および `real(DP), dimension(1:jm)`.
!>
!>  * `w_` で始まる関数が返す値はスペクトルデータに同じ.
!>
!>  * `xy_` で始まる関数が返す値は 2 次元格子点データに同じ.
!>
!>  * `x_`, `y_` で始まる関数が返す値は 1 次元格子点データに同じ.
!>
!>  * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
!>    微分などを作用させたデータをスペクトル変換したものことである.
!>
!> ## 変数・手続き群の要約
!>
!> ### 初期化, 終了処理
!>
!> w_module_base にて定義
!>
!> * [w_Initial]
!>   (@ref w_initial)
!>   : スペクトル変換の格子点数, 波数, 領域の大きさの設定
!> * [w_Finalize]
!>   (@ref w_finalize)
!>   : モジュールの終了処理(割り付け配列の解放)をおこなう.
!>
!> ### 座標変数
!>
!> * [x_Lon]
!>   (@ref w_module_base#x_lon)
!>   : 経度座標格子点を格納した 1 次元配列
!> * [y_Lat]
!>   (@ref w_module_base#y_lat)
!>   : 緯度座標格子点を格納した 1 次元配列
!> * [x_Lon_Weight]
!>   (@ref w_module_base#x_lon_weight)
!>   : 重み座標を格納した 1 次元配列
!> * [y_Lat_Weight]
!>   (@ref w_module_base#y_lat_weight)
!>   : 重み座標を格納した 1 次元配列
!> * [xy_Lon]
!>   (@ref w_module_base#xy_lon)
!>   : 経度座標(X,Y)の格子点データ(2 次元配列)
!> * [xy_Lat]
!>   (@ref w_module_base#xy_lat)
!>   : 緯度座標(X,Y)の格子点データ(2 次元配列)
!>
!> ### 基本変換
!>
!> * [xy_w]
!>   (@ref w_module_base#xy_w)
!>   : スペクトルデータから格子データへの変換
!> * [w_xy]
!>   (@ref w_module_base#w_xy)
!>   : 格子データからスペクトルデータへの変換
!> * [l_nm]
!>   (@ref w_module_base#l_nm),
!>   [nm_l]
!>   (@ref w_module_base#nm_l)
!>   : スペクトルデータの格納位置と全波数・帯状波数の変換
!> * [w_StreamPotential2Vector]
!>   (@ref w_module_base#w_streampotential2vector)
!>   : 流線ポテンシャルから速度場計算
!> * [w_Vector2VorDiv]
!>   (@ref w_module_base#w_vector2vordiv)
!>   : 速度場から渦度発散を計算
!> * [w_VectorCosLat2VorDiv]
!>   (@ref w_module_base#w_vectorcoslat2vordiv)
!>   : 速度場から渦度発散を計算
!>
!> ### 微分
!>
!> w_module_deriv にて定義
!>
!> * [w_Lapla_w]
!>   (@ref w_module_deriv#w_lapla_w)
!>   : スペクトルデータにラプラシアンを作用させる
!> * [rn]
!>   (@ref w_module_deriv#rn)
!>   : スペクトルデータのラプラシアンを計算するための係数.
!> * [w_LaplaInv_w]
!>   (@ref w_module_deriv#w_laplainv_w)
!>   : スペクトルデータにラプラシアンの逆変換を作用させる
!> * [w_DLon_w]
!>   (@ref w_module_deriv#w_dlon_w)
!>   : スペクトルデータに経度微分\f$∂/∂λ\f$を作用させる
!> * [xy_GradLon_w]
!>   (@ref w_module_deriv#xy_gradlon_w)
!>   : スペクトルデータに勾配型経度微分
!>   \f$1/\cosφ・∂/∂λ\f$を作用させる
!> * [xy_GradLat_w]
!>   (@ref w_module_deriv#xy_gradlat_w)
!>   : スペクトルデータに勾配型緯度微分\f$∂/∂φ\f$を作用させる
!> * [w_DivLon_xy]
!>   (@ref w_module_deriv#w_divlon_xy)
!>   : 格子データに発散型経度微分
!>   \f$1/\cosφ・∂/∂λ\f$を作用させる
!> * [w_DivLat_xy]
!>   (@ref w_module_deriv#w_divlat_xy)
!>   : 格子データに発散型緯度微分
!>   \f$1/\cosφ・∂(g \cosφ)/∂φ\f$を作用させる
!> * [w_Div_xy_xy]
!>   (@ref w_module_deriv#w_div_xy_xy)
!>   : ベクトル成分である 2 つの格子データに発散を作用させる
!> * [w_Jacobian_w_w]
!>   (@ref w_module_deriv#w_jacobian_w_w)
!>   : 2 つのスペクトルデータからヤコビアンを計算する
!>
!>  \f$(λ,μ)\f$ 座標 \f$(μ=\sinφ)\f$ での微分
!>
!> * [xy_GradLambda_w]
!>   (@ref w_module_deriv#xy_gradlambda_w)
!>   : スペクトルデータに勾配型経度微分\f$∂/∂λ\f$を作用させる
!> * [xy_GradMu_w]
!>   (@ref w_module_deriv#xy_gradmu_w)
!>   : スペクトルデータに勾配型緯度微分\f$(1-μ^2)∂/∂μ\f$を作用させる
!> * [w_DivLambda_xy]
!>   (@ref w_module_deriv#w_divlambda_xy)
!>   : 格子データに発散型経度微分\f$1/(1-μ^2)・∂/∂λ\f$を作用させる
!> * [w_DivMu_xy]
!>   (@ref w_module_deriv#w_divmu_xy)
!>   : 格子データに発散型緯度微分\f$∂/∂μ\f$を作用させる
!>
!> ### 補間
!>
!> w_module_interpolate にて定義
!>
!> * [Interpolate_w]
!>   (@ref w_module_interpolate#interpolate_w)
!>   : スペクトルデータから任意の点での値を求める.
!>
!> ### 積分・平均
!>
!> w_module_integral にて定義
!>
!> * [IntLonLat_xy]
!>   (@ref w_module_integral#intlonlat_xy)
!>   : 2 次元格子点データの全領域積分
!> * [AvrLonLat_xy]
!>   (@ref w_module_integral#avrlonlat_xy)
!>   : 2 次元格子点データの全領域平均
!> * [y_IntLon_xy]
!>   (@ref w_module_integral#y_intlon_xy)
!>   : 2 次元格子点データの経度方向積分
!> * [y_AvrLon_xy]
!>   (@ref w_module_integral#y_avrlon_xy)
!>   : 2 次元格子点データの経度方向平均
!> * [IntLon_x]
!>   (@ref w_module_integral#intlon_x)
!>   : 1 次元(X)格子点データの経度方向積分
!> * [AvrLon_x]
!>   (@ref w_module_integral#avrlon_x)
!>   : 1 次元(X)格子点データの経度方向平均
!> * [x_IntLat_xy]
!>   (@ref w_module_integral#x_intlat_xy)
!>   : 2 次元格子点データの緯度方向積分
!> * [x_AvrLat_xy]
!>   (@ref w_module_integral#x_avrlat_xy)
!>   : 2 次元格子点データの緯度方向平均
!> * [IntLat_y]
!>   (@ref w_module_integral#intlat_y)
!>   : 1 次元(Y)格子点データの緯度方向積分
!> * [AvrLat_y]
!>   (@ref w_module_integral#avrlat_y)
!>   : 1 次元(Y)格子点データの緯度方向平均
!>
!> ### スペクトル解析
!>
!> w_module_spectrum にて定義
!>
!> * [nm_EnergyFromStreamfunc_w]
!>   (@ref w_module_spectrum#nm_energyfromstreamfunc_w)
!>   : 流線関数からエネルギースペクトルを計算する(水平全波数 n, 帯状波数 m 空間)
!> * [n_EnergyFromStreamfunc_w]
!>   (@ref w_module_spectrum#n_energyfromstreamfunc_w)
!>   : 流線関数からエネルギースペクトルを計算する(水平全波数 n 空間)
!> * [nm_EnstrophyFromStreamfunc_w]
!>   (@ref w_module_spectrum#nm_enstrophyfromstreamfunc_w)
!>   : 流線関数からエンストロフィースペクトルを計算する (水平全波数 n, 帯状波数 m 空間)
!> * [n_EnstrophyFromStreamfunc_w]
!>   (@ref w_module_spectrum#n_enstrophyfromstreamfunc_w)
!>   : 流線関数からエンストロフィースペクトルを計算する (水平全波数 n 空間)
!> * [w_spectrum_VMiss]
!>   (@ref w_module_spectrum#w_spectrum_vmiss)
!>   :  欠損値
!>
module w_module_mint
  use dc_message, only : MessageNotify
  use w_module_base_mint, only :    &
    & x_Lon, y_Lat,                 &
    & x_Lon_weight, y_Lat_Weight,   &
    & xy_Lon, xy_Lat,               &
    & xy_w, w_xy, l_nm, nm_l,       &
    & w_StreamPotential2Vector,     &
    & w_Vector2VorDiv,              &
    & w_VectorCosLat2VorDiv,        &
    & w_base_initial,               &
    & w_base_finalize
  use w_module_deriv_mint, only:    &
    & w_Lapla_w, w_LaplaInv_w,      &
    & rn,                           &
    & w_DLon_w,                     &
    & xy_GradLon_w, xy_GradLat_w,   &
    & w_DivLon_xy, w_DivLat_xy,     &
    & w_Div_xy_xy,                  &
    & w_Jacobian_w_w,               &
    & xy_GradLambda_w, xy_GradMu_w, &
    & w_DivLambda_xy, w_DivMu_xy,   &
    & w_deriv_initial,              &
    & w_deriv_Finalize
  use w_module_integral_mint, only: &
    & IntLonLat_xy,                 &
    & y_IntLon_xy, IntLon_x,        &
    & x_IntLat_xy, IntLat_y,        &
    & AvrLonLat_xy,                 &
    & y_AvrLon_xy, AvrLon_x,        &
    & x_AvrLat_xy, AvrLat_y
  use w_module_spectrum_mint, only: &
    & nm_EnergyFromStreamfunc_w,    &
    & n_EnergyFromStreamfunc_w,     &
    & nm_EnstrophyFromStreamfunc_w, &
    & n_EnstrophyFromStreamfunc_w,  &
    & w_spectrum_VMiss
  use w_module_interpolate_mint, only:  &
    & Interpolate_w

  private
  public w_Initial
  public w_Finalize
  ! w_base_module
  public x_Lon
  public y_Lat
  public x_Lon_weight
  public y_Lat_Weight
  public xy_Lon
  public xy_Lat
  public xy_w
  public w_xy
  public l_nm
  public nm_l
  public w_StreamPotential2Vector
  public w_Vector2VorDiv
  public w_VectorCosLat2VorDiv
  ! w_derive_module
  public w_Lapla_w
  public w_LaplaInv_w
  public rn
  public w_DLon_w
  public xy_GradLon_w
  public xy_GradLat_w
  public w_DivLon_xy
  public w_DivLat_xy
  public w_Div_xy_xy
  public w_Jacobian_w_w
  public xy_GradLambda_w
  public xy_GradMu_w
  public w_DivLambda_xy
  public w_DivMu_xy
  ! w_integral_module
  public IntLonLat_xy
  public y_IntLon_xy
  public IntLon_x
  public x_IntLat_xy
  public IntLat_y
  public AvrLonLat_xy
  public y_AvrLon_xy
  public AvrLon_x
  public x_AvrLat_xy
  public AvrLat_y
  ! w_spectrum_module
  public nm_EnergyFromStreamfunc_w
  public n_EnergyFromStreamfunc_w
  public nm_EnstrophyFromStreamfunc_w
  public n_EnstrophyFromStreamfunc_w
  public w_spectrum_VMiss
  ! w_interpolate_module
  public Interpolate_w
contains

  !---------------------------------------------------------------------
  !> @brief
  !> スペクトル変換の格子点数, 波数および OPENMP 使用時の
  !> 最大スレッド数を設定する.
  !>
  !> @details
  !> 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
  !> しなければならない.
  !>
  !> np_in に 1 より大きな値を指定すれば ISPACK の球面調和函数変換
  !> OPENMP 並列計算ルーチンが用いられる. 並列計算を実行するには,
  !> 実行時に環境変数 OMP_NUM_THREADS を np_in 以下の数字に設定する等の
  !> システムに応じた準備が必要となる.
  !>
  !> np_in に 1 より大きな値を指定しなければ並列計算ルーチンは呼ばれない.
  !>
  subroutine w_initial(n_in,i_in,j_in,np,mint)
    !> 格子点数(東西)
    integer,intent(in) :: i_in
    !> 格子点数(南北)
    integer,intent(in) :: j_in
    !> 切断波数の設定
    integer,intent(in) :: n_in
    !> OPENMP での最大スレッド数
    integer,intent(in), optional :: np
    !> 経度方向対称性
    integer,intent(in), optional :: mint

    if ( present (np) ) then
      call MessageNotify('M','w_Initial', &
        & 'Optional argnument np is dummy.')
    endif

    if ( present (mint) ) then
       call w_base_initial(n_in,i_in,j_in,mint=mint)
    else
       call w_base_initial(n_in,i_in,j_in)
    endif
    call w_deriv_initial

    call MessageNotify('M','w_initial',&
      'w_module_mint (2020/11/18) is initialized')

  end subroutine w_initial

  !---------------------------------------------------------------------
  !> @brief
  !> モジュールの終了処理(割り付け配列の解放)をおこなう.
  !>
  !> @warning
  !> 解像度を変更する際にはこのサブルーチンを呼んで終了処理を
  !> おこなったのちに再度 w_Initial で初期設定しなければ
  !> ならない.
  !>
  subroutine w_Finalize
    call w_base_Finalize
    call w_deriv_Finalize

    call MessageNotify('M','w_Finalize',&
      'w_module_mint (2020/11/18) is finalized')

  end subroutine w_Finalize

end module w_module_mint
