!--
!----------------------------------------------------------------------
!     Copyright (c) 2024 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  wa_zonal_module
!
!   spml/wa_zonal_module モジュールは球面上での帯状的な 2 次元流体運動を
!   球面調和函数を用いたスペクトル法によって数値計算するための
!   Fortran90 関数を提供する. 
!
!   球面上の 1 層モデル用 w_zonal_module モジュールを多層モデル用に
!   拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
!   対する変換が行える.
!
!   wa_zonal_module は実際には基本変換, 微分計算をそれぞれ担っている
!   下部モジュール wa_zonal_module_base, wa_zonal_module_deriv,
!   wa_zonal_module_interpolate, および 1 層用のモジュール 
!   w_zonal_module からなっている.
!
!   内部で ISPACK3 の LXPACK の Fortran サブルーチンを呼んでいる.
!   スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!   ついては ISPACK3/LXPACK のマニュアルを参照されたい.
!
!   wa_zonal_modulde で提供される関数・サブルーチンは 2 次元的
!   流体運動を扱う wa_module モジュールで用いられているものと
!   名前およびインターフェースが共通になるように設計してある. 
!   したがって, wa_module を用いて構成された 2 次元モデルを
!   経度方向に 1 波数成分のモデルへと改造するには次の手順が必要となる.
!
!      * use 文での wa_module の引用を wa_zonal_module に変更する.
!      * 配列の大きさを経度方向格子点数 im -> 1 に,
!        水平波数を (nm+1)**2 -> nm-m+1 に変更する.
!      * DO 文で水平波数に関してループを回しているところを
!        (nm+1)**2 -> nm-m+1 に変更する.
!      * gtool 出力の次元変数変更する.
!
!
!履歴  2024/02/17  竹広真一  wa_wave_deriv_module_svpack より改造
!      2025/12/10  竹広真一  経度座標変数を追加
!      2025/12/15  竹広真一  領域積分・平均関数を追加
!
!++
module wa_zonal_module
  !
  != wa_zonal_module
  !
  !== 概要
  !
  !   spml/wa_zonal_module モジュールは球面上での 2 次元流体運動を
  !   球面調和函数を用いたスペクトル法によって東西 1 波数数値計算するための
  !   Fortran90 関数を提供する. 
  !
  !   球面上の 1 層モデル用 w_zonal_module モジュールを多層モデル用に
  !   拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
  !   対する変換が行える.
  !
  !   wa_zonal_module は実際には基本変換, 微分計算をそれぞれ担っている
  !   下部モジュール wa_zonal_module_base, wa_zonal_module_deriv,
  !   wa_zonal_module_interpolate, および 1 層用のモジュール 
  !   w_zonal_module からなっている.
  !
  !   内部で ISPACK3 の LXPACK の Fortran サブルーチンを呼んでいる.
  !   スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  !   ついては ISPACK3/LXPACK のマニュアルを参照されたい.
  !
  !   wa_zonal_modulde で提供される関数・サブルーチンは 2 次元的
  !   流体運動を扱う wa_module モジュールで用いられているものと
  !   名前およびインターフェースが共通になるように設計してある. 
  !   したがって, wa_module を用いて構成された 2 次元モデルを
  !   経度方向に 1 波数成分のモデルへと改造するには次の手順が必要となる.
  !
  !      * use 文での wa_module の引用を wa_zonal_module に変更する.
  !      * 配列の大きさを経度方向格子点数 im -> 1 に,
  !        水平波数を (nm+1)**2 -> nm-m+1 に変更する.
  !      * DO 文で水平波数に関してループを回しているところを
  !        (nm+1)**2 -> nm-m+1 に変更する.
  !      * gtool 出力の次元変数変更する.
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (wa_, nma_, na_, xya_, ya_, w_, xy_, y_, a_) は, 
  !   返す値の形を示している.
  !   wa_  :: 水平スペクトルデータの並んだ 2 次元配列(スペクトルデータの
  !           並びは SNPACK/ISPACK に従ったもの)
  !   xya_ :: 2 次元格子点データの並んだ 3 次元配列
  !   ya_  :: 緯度方向 1 次元格子点データの並んだ 2 次元配列
  !
  ! * 関数名の間の文字列(DLon, GradLat, GradLat, DivLon, DivLat, Lapla, 
  !   LaplaInv, Jacobian)は, その関数の作用を表している.
  !
  ! * 関数名の最後 (_wa_wa, _wa, _xya, _ya, _w, _xy, _y) は, 
  !   入力変数の形スペクトルデータおよび格子点データであることを示している.
  !   _wa    :: スペクトルデータの並んだ 2 次元配列
  !   _wa_wa :: 2 つのスペクトルデータの並んだ 2 次元配列
  !   _xya   :: 2 次元格子点データの並んだ 3 次元配列
  !   _ya    :: 緯度方向 1 次元格子点データの並んだ 2 次元配列
  !
  !=== 各データの種類の説明
  !
  ! * xya : 2 次元格子点データの並んだ 3 次元配列.
  !   * 変数の種類と次元は real(8), dimension(0:im-1,1:jm,:). 
  !   * im(=2), jm はそれぞれ経度, 緯度座標の格子点数であり, サブルーチン
  !     wa_Initial にてあらかじめ設定しておく. 
  !   * 扱う第 3 次元の大きさの最大値を wa_Initial で設定しておく.
  !
  ! * wa : スペクトルデータの並んだ 2 次元配列.
  !   * 変数の種類と次元は real(8), dimension(2*(nm-m+1),:). 
  !   * nm は球面調和函数の最大全波数, m は東西波数であり,
  !     サブルーチン wa_Initial にてあらかじめ設定しておく. 
  !   * スペクトルデータの格納のされ方は関数 
  !     l_nm, nm_l によって調べることができる. 
  !   * 扱う第 3 次元の大きさの最大値を wa_Initial で設定しておく.
  !
  ! * ya : 緯度方向 1 次元格子点データ.
  !   * 変数の種類と次元は real(8), dimension(1:jm,:).
  !
  ! * wa_ で始まる関数が返す値はスペクトルデータの並んだ 2 次元配列に同じ.
  !
  ! * xya_ で始まる関数が返す値は 2 次元格子点データの並んだ 3 次元配列に
  !   同じ.
  !
  ! * ya_ で始まる関数が返す値は 1 次元格子点データの並んだ 2 次元配列に同じ.
  !   
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したものことである.
  !
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化 
  !
  ! wa_Initial :: スペクトル変換の格子点数, 波数, 領域の大きさ, 
  !            :: および同時に計算するデータの個数の最大値設定
  !
  !==== 終了処理 
  !
  ! wa_Finalize :: モジュールの終了処理(割り付け配列の解放)をおこなう. 
  ! 
  !==== 座標変数
  !
  ! y_Lat        :: 格子点座標(経度座標)を格納した 1 次元配列
  ! y_Lat_Weight :: 重み座標を格納した 1 次元配列
  ! xy_Lat       :: 格子点データの緯度座標(格子点データ型 2 次元配列)
  !
  !==== 基本変換
  !
  ! xya_wa, xy_w :: スペクトルデータから格子データへの変換(多層, 1 層用)
  ! wa_xya, w_xy :: 格子データからスペクトルデータへの変換(多層, 1 層用)
  ! l_nm, nm_l :: スペクトルデータの格納位置と全波数・帯状波数の変換 
  !
  !==== 微分
  !
  ! wa_Lapla_wa, w_Lapla_w        :: スペクトルデータにラプラシアンを
  !                                  作用させる(多層, 1 層用)
  ! rn                            :: スペクトルデータのラプラシアンを
  !                                  計算するための係数
  ! irm                           :: 経度微分演算用配列(未定義)
  ! wa_LaplaInv_wa, w_LaplaInv_w  :: スペクトルデータにラプラシアンの逆変換を
  !                                  作用させる(多層, 1 層用)
  ! wa_DLon_wa, w_DLon_w          :: スペクトルデータに
  !                                  経度微分∂/∂λを作用させる(多層, 1 層用)
  ! xya_GradLon_wa, xy_GradLon_w  :: スペクトルデータに
  !                                  勾配型経度微分 1/cosφ・∂/∂λ を
  !                                  作用させる(多層, 1 層用)
  ! xya_GradLat_wa, xy_GradLat_w  :: スペクトルデータに
  !                                  勾配型緯度微分 ∂/∂φ を
  !                                  作用させる(多層, 1 層用)
  ! wa_DivLon_xya, w_DivLon_xy    :: 格子データに
  !                                  発散型経度微分 1/cosφ・∂/∂λを
  !                                  作用させる(多層, 1 層用)
  ! wa_DivLat_xya, w_DivLat_xy    :: 格子データに
  !                                  発散型緯度微分 1/cosφ・∂(g cosφ)/∂φを
  !                                  作用させる (多層, 1 層用)
  ! wa_Div_xya_xya, w_Div_xy_xy   :: ベクトル成分である 2 つの格子データに
  !                                  発散を作用させる(多層, 1 層用)
  !
  !==== 微分(λ,μ=sinφ 座標)
  !
  ! xya_GradLambda_wa, xy_GradLambda_w  :: スペクトルデータに
  !                                        勾配型経度微分∂/∂λを作用させる
  !                                        (多層, 1 層用)
  ! xya_GradMu_wa, xy_GradMu_w          :: スペクトルデータに
  !                                        勾配型緯度微分(1-μ^2)∂/∂μ
  !                                        を作用させる
  !                                        (多層, 1 層用)
  ! wa_DivLambda_xya,w_DivLambda_xy     :: 格子データに発散型経度微分
  !                                        1 /(1-μ^2)・∂/∂λを作用させる
  !                                        (多層, 1 層用)
  !
  ! wa_DivMu_xya,w_DivMu_xy             :: 格子データに発散型緯度微分
  !                                        ∂/∂μを作用させる(多層, 1 層用)
  !
  !==== 補間
  !
  ! a_Interpolate_wa :: スペクトルデータから任意の点での値を求める.
  !
  !
  use dc_message, only : MessageNotify
  use w_zonal_module
  use wa_zonal_module_base
  use wa_zonal_module_deriv
  use wa_zonal_module_integral
  use wa_zonal_module_interpolate

  private

  public wa_Initial                           ! 初期化
  public wa_Finalize                          ! 終了処理

  public x_Lon, y_Lat                         ! 格子座標
  public x_Lon_weight, y_Lat_Weight           ! 格子座標重み
  public xy_Lon, xy_Lat                       ! 格子座標(im,jm)
  public xy_w, w_xy, l_nm, nm_l               ! 変換関数
  public xya_wa, wa_xya                       ! 変換関数
  
  public rn                                   ! ラプラシアンの係数配列
  public w_Lapla_w, w_LaplaInv_w              ! ラプラシアンと逆演算
  public w_DLon_w                             ! 経度微分
  public xy_GradLon_w, xy_GradLat_w           ! 勾配型微分
  public w_DivLon_xy, w_DivLat_xy             ! 発散型微分
  public w_Div_xy_xy                          ! 発散型微分
  public xy_GradLambda_w, xy_GradMu_w         ! 勾配型微分(λ,μ座標)
  public w_DivLambda_xy, w_DivMu_xy           ! 発散型微分(λ,μ座標)

  public wa_Lapla_wa, wa_LaplaInv_wa          ! ラプラシアンと逆演算
  public wa_DLon_wa                           ! 経度微分
  public xya_GradLon_wa, xya_GradLat_wa       ! 勾配型微分
  public wa_DivLon_xya, wa_DivLat_xya         ! 発散型微分
  public wa_Div_xya_xya                       ! 発散型微分
  public xya_GradLambda_wa, xya_GradMu_wa     ! 勾配型微分(λ,μ座標)
  public wa_DivLambda_xya, wa_DivMu_xya       ! 発散型微分(λ,μ座標)

  public a_IntLonLat_xya                      ! 緯度経度積分
  public ya_IntLon_xya, a_IntLon_xa           ! 経度積分
  public xa_IntLat_xya, a_IntLat_ya           ! 緯度積分
  public a_AvrLonLat_xya                      ! 緯度経度平均
  public ya_AvrLon_xya, a_AvrLon_xa           ! 経度平均
  public xa_AvrLat_xya, a_AvrLat_ya           ! 緯度平均

  public IntLonLat_xy                         ! 緯度経度積分
  public y_IntLon_xy, IntLon_x                ! 経度積分
  public x_IntLat_xy, IntLat_y                ! 緯度積分
  public AvrLonLat_xy                         ! 緯度経度積分
  public y_AvrLon_xy, AvrLon_x                ! 経度積分
  public x_AvrLat_xy, AvrLat_y                ! 緯度積分

  public Interpolate_w, a_Interpolate_wa      ! 補間関数

contains
  subroutine wa_Initial(n_in,i_in,j_in,k_in)
    !
    ! スペクトル変換の格子点数, 波数, 最大データ数(層数)および 
    ! OPENMP 使用時の最大スレッド数を設定する.
    !
    ! 他の関数を呼ぶ前に, 最初にこのサブルーチンを呼んで初期設定を
    ! しなければならない. 
    !
    ! SNPACK 用 wa_initial との互換性のために引数 k_in が用意されているが, 
    ! このパラメターによって同時に処理できる変数の数(層数)を制限される事はない.
    !
    integer,intent(in) :: i_in                ! 格子点数(東西)
    integer,intent(in) :: j_in                ! 格子点数(南北)
    integer,intent(in) :: n_in                ! 切断波数
    integer,intent(in), optional :: k_in      ! 最大データ数(層数)

    call w_Initial(n_in,i_in,j_in)
!!$    call wa_base_Initial(k_in)

    call MessageNotify('M','wa_initial', &
                       'wa_zonal_module (2025/12/10) is initialized')

  end subroutine wa_Initial

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

    call MessageNotify('M','wa_Finalize', &
                       'wa_zonal_module (2025/12/10) is finalized')

  end subroutine wa_Finalize

end module wa_zonal_module
