! -*- mode: f90; coding: utf-8 -*-
!-----------------------------------------------------------------------
! Copyright (c) 2019-2020 SPMODEL Development Group. All rights reserved.
!-----------------------------------------------------------------------
!> @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次元球面: 球面調和関数展開，MPI並列
!> @details
!> spml/w_mpi_module モジュールは球面上での 2 次元流体運動を
!> 球面調和函数を用いたスペクトル法と MPI 並列化によって数値計算するための
!> Fortran90 関数を提供する.
!>
!> w_mpi_module は実際には以下の下部モジュールからなる．
!> * w_mpi_module_base: 基本変換
!> * w_mpi_module_deriv: 微分計算
!> * w_mpi_module_integral: 積分・平均計算
!> * w_mpi_module_spectrum: スペクトル解析
!> * w_mpi_module_interpolate: 補間計算
!>
!> 内部で ISPACK3 の SYPACK のサブルーチンを呼んでいる.
!> スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!> ついては ISPACK3/SYPACK のマニュアルを参照されたい.
!
! 履歴  2019/03/17  竹広真一
!       2019/10/05  佐々木洋平
!       2020/11/07  竹広真一   セクター計算オプション導入
!
!>
!> ## 関数・変数の名前と型について
!>
!> ### 命名法
!>
!> * 関数名の先頭 (`w_`, `nm_`, `n_`, `xy_`, `xv_`, `x_`, `y_`, `v_`) は,
!>   返す値の形を示している.
!>   *  `w_` :: スペクトルデータ
!>   * `xy_` :: 2 次元格子点データ
!>   * `xv_` :: 2 次元分散格子点データ
!>   * `nm_` :: スペクトルデータの並んだ 3 次元配列(スペクトルデータの並びは
!>              全波数 n, 帯状波数 m で指定される 2 次元配列)
!>   *  `n_` :: スペクトルデータの並んだ 2 次元配列 (スペクトルデータの並びは
!>              全波数 n で指定される 1 次元配列)
!>   *  `x_` :: 経度方向 1 次元格子点データ
!>   *  `y_` :: 緯度方向 1 次元格子点データ
!>   *  `v_ `:: 緯度方向 1 次元分散格子点データ
!>
!> * 関数名の間の文字列
!>   `DLon`, `GradLat`, `GradLat`, `DivLon`, `DivLat`,
!>   `Lapla`, `LaplaInv`, `Jacobian` は, その関数の作用を表している.
!>
!> * 関数名の最後 (`_w_w`, `_w`, `_xy`, `_x`, `_y`) は,
!>   入力変数の形スペクトルデータおよび格子点データであることを示している.
!>   * `_w`   :: スペクトルデータ
!>   * `_w_w` :: 2 つのスペクトルデータ
!>   * `_xy`  :: 2 次元格子点データ
!>   * `_xv`  :: 2 次元分散格子点データ
!>   * `_x`   :: 経度方向 1 次元格子点データ
!>   * `_y`   :: 緯度方向 1 次元格子点データ
!>   * `_v`   :: 緯度方向 1 次元分散格子点データ
!>
!> ### 各データの種類の説明
!>
!> * 以下 `DP  = kind(1.0D)` である
!>
!> * `xy` : 2 次元格子点データ.
!>   * 変数の種類と次元は `real(DP), dimension(im,jm)`
!>   * `im`, `jm` はそれぞれ経度, 緯度座標の格子点数であり, サブルーチン
!>     @ref w_mpi_initial にてあらかじめ設定しておく.
!>
!> * `xv` : 2 次元分散格子点データ.
!>   * 変数の種類と次元は `real(DP), dimension(im,jc)`.
!>   * `im`, `jc` はそれぞれ経度の格子点数および
!>     プロセスで保有する緯度方向の格子点数である.
!>   * `jc` はサブルーチン @ref w_mpi_initial にて設定して
!>     おくことにより設定される public 変数である.
!>
!> * `w` : スペクトルデータ.
!>   * 変数の種類と次元は `real(DP), dimension((nm+1)*(nm+1))`.
!>     `nm` は球面調和函数の最大全波数であり,
!>     サブルーチン @ref w_mpi_initial にてあらかじめ設定しておく.
!>   * スペクトルデータの格納のされ方は
!>     関数  [l_nm](@ref w_mpi_module_base#l_nm),
!>     [nm_l](@ref w_mpi_module_base#nm_l) によって調べることができる.
!>
!> * `nm` : スペクトルデータの並んだ 2 次元配列.
!>   * 変数の種類と次元は `real(DP), dimension(0:nm,-nm:nm)`.
!>     第 1 次元が水平全波数,  第 2 次元が帯状波数を表す.
!>   * `nm` は球面調和函数の最大全波数であり,
!>     サブルーチン w_mpi_initial にてあらかじめ設定しておく.
!>
!> * `n` : スペクトルデータの並んだ 1 次元配列.
!>   * 変数の種類と次元は `real(DP), dimension(0:nm)`.
!>   * 第 1 次元が水平全波数を表す. `nm` は球面調和函数の最大全波数であり,
!>     サブルーチン w_mpi_initial にてあらかじめ設定しておく.
!>
!> * `x`, `y` : 経度, 緯度方向 1 次元格子点データ.
!>   * 変数の種類と次元はそれぞれ `real(DP), dimension(im)`
!>     および `real(DP), dimension(jm)`.
!>
!> * `w_` で始まる関数が返す値はスペクトルデータに同じ.
!>
!> * `xy_` で始まる関数が返す値は 2 次元格子点データに同じ.
!>
!> * `xv_` で始まる関数が返す値は 2 次元分散格子点データに同じ.
!>
!> * `x_`, `y_`, `v_` で始まる関数が返す値は 1 次元格子点データに同じ.
!>
!> * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
!>   微分などを作用させたデータをスペクトル変換したものことである.
!>
!> ## 変数・手続き群の要約
!>
!> ### 初期化・終了処理
!>
!> * @ref w_mpi_initial
!>   : スペクトル変換の格子点数, 波数, 領域の大きさの設定
!> * @ref w_mpi_finalize
!>   : モジュールの終了処理(割り付け配列の解放)をおこなう.
!> * [x_Lon](@ref w_mpi_module_base#x_lon),
!>   [x_Lon_Weight](@ref w_mpi_module_base#x_lon_weight)
!>   : 格子点座標(経度)と重みを格納した 1 次元配列
!> * [y_Lat](@ref w_mpi_module_base#y_lat),
!>   [y_Lat_Weight](@ref w_mpi_module_base#y_lat_weight)
!>   : 格子点座標(緯度)と重みを格納した 1 次元配列
!> * [v_Lat](@ref w_mpi_module_base#v_lat),
!>   [v_Lat_Weight](@ref w_mpi_module_base#v_lat_weight)
!>   : 分散格子点座標(緯度)と重みを格納した 1 次元配列
!> * [xy_Lon](@ref w_mpi_module_base#xy_lon),
!>   [xy_Lat](@ref w_mpi_module_base#xy_lat)
!>   : 格子点データの経度・緯度座標(X,Y), 格子点データ型 2 次元配列
!> * [xv_Lon](@ref w_mpi_module_base#xv_lon),
!>   [xv_Lat](@ref w_mpi_module_base#xv_lat),
!>   : 分散格子点データの経度・緯度座標(X,Y), 格子点データ型 2 次元配列
!>
!> ### 基本変換
!>
!> 格子点データとスペクトルの変換．@ref w_mpi_module_base で定義.
!>
!> * [xv_w](@ref w_mpi_module_base#xv_w)
!>   : スペクトルデータから分散格子データへの変換
!> * [w_xv](@ref w_mpi_module_base#w_xv)
!>   : 分散格子データからスペクトルデータへの変換
!> * [l_nm](@ref w_mpi_module_base#l_nm),
!>   [nm_l](@ref w_mpi_module_base#nm_l)
!>   : スペクトルデータの格納位置と全波数・帯状波数の変換
!> * [w_StreamPotential2VectorMPI]
!>   (@ref w_mpi_module_base#w_streampotential2vectormpi)
!>   : 流線ポテンシャルから速度場計算
!> * [w_Vector2VorDivMPI]
!>   (@ref w_mpi_module_base#w_vector2vordivmpi)
!>   : 速度場から渦度発散を計算
!> * [w_VectorCosLat2VorDivMPI]
!>   (@ref w_mpi_module_base#w_vectorcoslat2vordivmpi)
!>   : 速度場から渦度発散を計算
!>
!> ### 微分
!>
!> 微分演算. @ref w_mpi_module_deriv で定義.
!>
!> * [w_Lapla_w]
!>   (@ref w_mpi_module_deriv#w_lapla_w)
!>   : スペクトルデータにラプラシアンを作用させる
!> * [rn]
!>   (@ref w_mpi_module_deriv#rn)
!>   : スペクトルデータのラプラシアンを計算するための係数.
!> * [w_LaplaInv_w]
!>   (@ref w_mpi_module_deriv#w_laplainv_w)
!>  : スペクトルデータにラプラシアンの逆変換を作用させる
!> * [w_DLon_w]
!>   (@ref w_mpi_module_deriv#w_dlon_w)
!>  : スペクトルデータに経度微分∂/∂λを作用させる
!> * [xv_GradLon_w]
!>   (@ref w_mpi_module_deriv#xv_gradlon_w)
!>  : スペクトルデータに勾配型経度微分
!>  \f$1/\cosφ・∂/∂λ\f$ を作用させる
!> * [xv_GradLat_w]
!>   (@ref w_mpi_module_deriv#xv_gradlat_w)
!>  : スペクトルデータに勾配型緯度微分\f$∂/∂φ\f$を作用させる
!> * [w_DivLon_xv]
!>   (@ref w_mpi_module_deriv#w_divlon_xv)
!>  : 格子データに発散型経度微分
!>  \f$1/\cosφ・∂/∂λ\f$ を作用させる
!> * [w_DivLat_xv]
!>   (@ref w_mpi_module_deriv#w_divlat_xv)
!>   : 格子データに発散型緯度微分
!>   \f$1/\cosφ・∂(g \cosφ)/∂φ\f$ を作用させる
!> * [xv_GradLon_w]
!>   (@ref w_mpi_module_deriv#xv_gradlon_w)
!>   : スペクトルデータに勾配型経度微分
!>   \f$ 1/\cosφ・∂/∂λ\f$ を作用させる
!> * [xv_GradLat_w]
!>   (@ref w_mpi_module_deriv#xv_gradlat_w)
!>   : スペクトルデータに勾配型緯度微分\f$∂/∂φ\f$を作用させる
!> * [w_DivLon_xv]
!>   (@ref w_mpi_module_deriv#w_divlon_xv)
!>   : 格子データに発散型経度微分 \f$1/\cosφ・∂/∂λ\f$を作用させる
!> * [w_DivLat_xv]
!>   (@ref w_mpi_module_deriv#w_divlat_xv)
!>   : 格子データに発散型緯度微分
!>   \f$1/\cosφ・∂(g \cosφ)/∂φ\f$ を作用させる
!> * [w_Div_xv_xv]
!>   (@ref w_mpi_module_deriv#w_div_xv_xv)
!>   : ベクトル成分である 2 つの格子データに発散を作用させる
!> * [w_JacobianMPI_w_w]
!>   (@ref w_mpi_module_deriv#w_jacobianmpi_w_w)
!>   : 2 つのスペクトルデータからヤコビアンを計算する
!>
!> \f$(λ,μ)\f$ 座標(\f$μ=\sinφ)\f$での微分
!>
!>  * [xv_GradLambda_w]
!>   (@ref w_mpi_module_deriv#xv_gradlambda_w)
!>    : スペクトルデータに勾配型経度微分\f$∂/∂λ\f$を作用させる
!>  * [xv_Gradmu_w]
!>   (@ref w_mpi_module_deriv#xv_gradmu_w)
!>    : スペクトルデータに勾配型緯度微分
!>    \f$(1-μ^2)∂/∂μ\f$を作用させる
!>  * [xv_GradLambda_w]
!>   (@ref w_mpi_module_deriv#xv_gradlambda_w)
!>    : スペクトルデータに勾配型経度微分\f$∂/∂λ\f$を作用させる
!>  * [xv_GradMu_w]
!>   (@ref w_mpi_module_deriv#xv_gradmu_w)
!>    : スペクトルデータに勾配型緯度微分
!>    \f$(1-μ^2)∂/∂μ\f$を作用させる
!>  * [w_DivLambda_xv]
!>   (@ref w_mpi_module_deriv#w_divlambda_xv)
!>    : 格子データに発散型経度微分
!>    \f$1/(1-μ^2)・∂/∂λ\f$を作用させる
!>  * [w_DivMu_xv]
!>   (@ref w_mpi_module_deriv#w_divmu_xv)
!>    : 格子データに発散型緯度微分\f$∂/∂μ\f$を作用させる
!>
!> ### 補間
!>
!> 補間. @ref w_mpi_module_interpolate で定義.
!>
!>  * [Interpolate_w]
!>    (@ref w_mpi_module_interpolate#interpolate_w)
!>    : スペクトルデータから任意の点での値を求める.
!>
!> ### 積分・平均
!>
!> @ref w_mpi_module_integral で定義.
!>
!>  * [IntLonLat_xv]
!>    (@ref w_mpi_module_integral#intlonlat_xv)
!>    : 2 次元格子点データの全領域積分
!>  * [avrlonlat_xv]
!>    (@ref w_mpi_module_integral#avrlonlat_xv)
!>    : 2 次元格子点データの全領域平均
!>  * [v_intlon_xv]
!>    (@ref w_mpi_module_integral#v_intlon_xv)
!>    : 2 次元格子点データの経度方向積分
!>  * [v_avrlon_xv]
!>    (@ref w_mpi_module_integral#v_avrlon_xv)
!>    : 2 次元格子点データの経度方向平均
!>  * [intlon_x]
!>    (@ref w_mpi_module_integral#intlon_x)
!>    : 1 次元(x)格子点データの経度方向積分
!>  * [avrlon_x]
!>    (@ref w_mpi_module_integral#avrlon_x)
!>    : 1 次元(x)格子点データの経度方向平均
!>  * [x_intlat_xv]
!>    (@ref w_mpi_module_integral#x_intlat_xv)
!>    : 2 次元格子点データの緯度方向積分
!>  * [x_avrlat_xv]
!>    (@ref w_mpi_module_integral#x_avrlat_xv)
!>    : 2 次元格子点データの緯度方向平均
!>  * [IntLat_v]
!>    (@ref w_mpi_module_integral#intlat_v)
!>    : 1 次元(Y)格子点データの緯度方向積分
!>  * [AvrLat_v]
!>    (@ref w_mpi_module_integral#avrlat_v)
!>    : 1 次元(Y)格子点データの緯度方向平均
!>
!> ### スペクトル解析
!>
!> @ref w_mpi_module_spectrum で定義.
!>
!> * [nm_EnergyFromStreamfunc_w]
!>   (@ref w_mpi_module_spectrum#nm_energyfromstreamfunc_w)
!>   : 流線関数からエネルギースペクトルを計算する
!>   (水平全波数 n, 帯状波数 m 空間)
!> * [n_EnergyFromStreamfunc_w]
!>   (@ref w_mpi_module_spectrum#n_energyfromstreamfunc_w)
!>   : 流線関数からエネルギースペクトルを計算する
!>   (水平全波数 n 空間)
!> * [nm_EnstrophyFromStreamfunc_w]
!>   (@ref w_mpi_module_spectrum#nm_enstrophyfromstreamfunc_w)
!>   : 流線関数からエンストロフィースペクトルを計算する
!>   (水平全波数 n, 帯状波数 m 空間)
!> * [n_EnstrophyFromStreamfunc_w]
!>   (@ref w_mpi_module_spectrum#n_enstrophyfromstreamfunc_w)
!>   : 流線関数からエンストロフィースペクトルを計算する
!>    (水平全波数 n 空間)
!> * [w_spectrum_VMiss]
!>   (@ref w_mpi_module_spectrum#w_spectrum_vmiss)
!>   :  欠損値
module w_mpi_module_mint
  use dc_message, only : MessageNotify
  use w_mpi_module_base_mint, only :    &
    & jc => jl, nc,                &
    & x_Lon, x_Lon_Weight,         &
    & v_Lat, v_Lat_Weight,         &
    & xv_Lon, xv_Lat,              &
    & y_Lat, y_lat_Weight,         &
    & xy_Lon, xy_Lat,              &
    & w_base_mpi_Initial,          &
    & w_base_mpi_Finalize,         &
    & nm_l, l_nm, xv_w, w_xv,      &
    & w_StreamPotential2VectorMPI, &
    & w_Vector2VorDivMPI,          &
    & w_VectorCosLat2VorDivMPI,    &
    & ay_av, av_ay
  use w_mpi_module_deriv_mint, only :    &
    & rn,                           &
    & w_Lapla_w, w_LaplaInv_w,      &
    & w_DLon_w,                     &
    & w_CosLatDLat_w,               &
    ! & w_DLatCosLat_w,               &
    & xv_GradLon_w, xv_GradLat_w,   &
    & w_DivLon_xv, w_DivLat_xv,     &
    & w_Div_xv_xv,                  &
    & w_JacobianMPI_w_w,            &
    & xv_GradLambda_w, xv_GradMu_w, &
    & w_DivLambda_xv, w_DivMu_xv,   &
    & w_deriv_mpi_initial,          &
    & w_deriv_mpi_Finalize
  use w_mpi_module_integral_mint
  use w_mpi_module_interpolate_mint
  use w_mpi_module_spectrum_mint

  implicit none
  private

  public jc                                   ! 緯度方向分散格子点数
  public nc                                   ! スペクトル分散格子点数

  public w_mpi_Initial                        ! 初期化
  public w_mpi_Finalize                       ! 終了処理

  public x_Lon, x_Lon_Weight                  ! 経度分散格子座標・重み
  public v_Lat, v_Lat_Weight                  ! 経度分散格子座標・重み
  public y_Lat, y_Lat_Weight                  ! 経度分散格子座標・重み
  public xv_Lon, xv_Lat                       ! 分散格子座標(im,jc)
  public xy_Lon, xy_Lat                       ! 全格子座標(im,jm)
  public l_nm, nm_l                           ! 波数格納位置
  public xv_w, w_xv                           ! 変換関数
  public w_StreamPotential2VectorMPI          ! 流線ポテンシャルから速度場計算
  public w_Vector2VorDivMPI                   ! 速度場から渦度発散を計算
  public w_VectorCosLat2VorDivMPI             ! 速度場から渦度発散を計算
  public ay_av, av_ay                         ! 分散データ/全データ変換

  public rn                                   ! ラプラシアン演算用配列
  public w_Lapla_w, w_LaplaInv_w              ! ラプラシアンと逆演算
  public w_DLon_w                             ! 経度微分
  public w_CosLatDLat_w                       ! 緯度微分
  ! public w_DLatCosLat_w                       ! 緯度微分  
  public xv_GradLon_w, xv_GradLat_w           ! 勾配型微分
  public w_DivLon_xv, w_DivLat_xv             ! 発散型微分
  public w_Div_xv_xv                          ! 発散型微分
  public w_JacobianMPI_w_w                    ! ヤコビアン
  public xv_GradLambda_w, xv_GradMu_w         ! 勾配型微分(λ,μ座標)
  public w_DivLambda_xv, w_DivMu_xv           ! 発散型微分(λ,μ座標)

  public IntLonLat_xv                         ! 緯度経度積分
  public v_IntLon_xv, IntLon_x                ! 経度積分
  public x_IntLat_xv, IntLat_v                ! 緯度積分
  public AvrLonLat_xv                         ! 緯度経度平均
  public v_AvrLon_xv, AvrLon_x                ! 経度平均
  public x_AvrLat_xv, AvrLat_v                ! 緯度平均

  !> 補間関数
  public Interpolate_w

  public nm_EnergyFromStreamfunc_w            ! エネルギースペクトル
                                              ! (水平全波数 n, 帯状波数 m 空間)
  public n_EnergyFromStreamfunc_w             ! エネルギースペクトル
                                              ! (水平全波数 n 空間)
  public nm_EnstrophyFromStreamfunc_w         ! エンストロフィースペクトル
                                              ! (水平全波数 n, 帯状波数 m 空間)
  public n_EnstrophyFromStreamfunc_w          ! エンストロフィースペクトル
                                              ! (水平全波数 n 空間)
  public w_spectrum_VMiss                     ! 欠損値

contains

  !----------------------------------------------------------
  !> スペクトル変換の格子点数, 波数および OPENMP 使用時の
  !> 最大スレッド数を設定する.
  !>
  !> 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
  !> しなければならない.
  !>
  !> np_in に 1 より大きな値を指定すれば ISPACK の球面調和函数変換
  !> OPENMP 並列計算ルーチンが用いられる. 並列計算を実行するには,
  !> 実行時に環境変数 OMP_NUM_THREADS を np_in 以下の数字に設定する等の
  !> システムに応じた準備が必要となる.
  !>
  !> np_in に 1 より大きな値を指定しなければ並列計算ルーチンは呼ばれない.
  !----------------------------------------------------------
  subroutine w_mpi_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_mpi_base_Initial', &
            'Optional argument "NP_IN" is dummy in w_mpi_module')
    endif

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

    call w_deriv_mpi_initial

    call MessageNotify('M','w_mpi_initial', &
      'w_mpi_module (2020/11/07) is initialized')

  end subroutine w_mpi_initial

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

    call MessageNotify('M','w_mpi_Finalize',&
      'w_mpi_module (2020/11/07) is finalized')

  end subroutine w_mpi_Finalize

end module w_mpi_module_mint