!--
!----------------------------------------------------------------------
!     Copyright (c) 2024-2025 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  w_zonal_module_base
!
!  spml/w_zonal_module_base モジュールは球面上での 2 次元流体運動を
!  球面調和函数を用いたスペクトル法によって数値計算するためのモジュール 
!  w_zonal_module の下部モジュールであり, スペクトル計算の基本的な 
!  Fortran90 関数を提供する.
!
!  内部で ISPACK3 の LXPACK(SXPACK の Fortran サブルーチンを呼んでいる. 
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算
!  法については ISPACK3/SXPACK のマニュアルを参照されたい.
!
!== 履歴
!
!      2024/02/15  竹広真一  w_wave_module_base.90 より改造
!      2025/11/27  竹広真一  w_base_initial interface を変更
!      2025/12/10  竹広真一  経度座標変数を追加
!
!      制限
!         ・変換する格子点データ, スペクトルデータの配列の大きさは決めうち
!         ・波数切断の仕方は三角波数切断に決めうち. 
!
!++
module w_zonal_module_base
  !
  != w_zonal_module_base
  !
  !== 概要.
  !
  !  spml/w_zonal_module_base モジュールは球面上での 2 次元流体運動を
  !  球面調和函数を用いたスペクトル法によって数値計算するためのモジュール 
  !  w_zonal_module の下部モジュールであり, スペクトル計算の基本的な 
  !  Fortran90 関数を提供する.
  !
  !  内部で ISPACK3 の LXPACK(SXPACK の Fortran サブルーチンを呼んでいる. 
  !  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算
  !  法については ISPACK3/SXPACK のマニュアルを参照されたい.
  !
  !
  use dc_message, only : MessageNotify
  implicit none

  integer(8)            :: im=1             ! 東西波数成分
  integer(8)            :: jm=32            ! 格子点の設定(南北)
  integer(8)            :: nm=21            ! 計算する最大の全波数の設定
  integer(8)            :: nn=22            ! 切断波数(全波数)の設定
  integer(8)            :: mm=21            ! 切断波数(東西波数)の設定
  integer(8)            :: ig=1             ! 緯度格子点種類
  
  real(8),    allocatable  :: pz(:,:)       ! 変換用配列
  real(8),    allocatable  :: rm(:)         ! 変換用配列
  real(8),    allocatable  :: cm(:)         ! 作業配列

  real(8), allocatable  :: x_Lon(:)         ! 経度経度
  real(8), allocatable  :: x_Lon_Weight(:)  ! 座標重み
  real(8), allocatable  :: y_Lat(:)         ! 緯度経度
  real(8), allocatable  :: y_Lat_Weight(:)  ! 座標重み
  real(8), allocatable  :: xy_Lon(:,:), xy_Lat(:,:)

  logical               :: w_zonal_base_initialize=.false.   ! 初期化フラッグ

  real(8), parameter    :: pi=3.1415926535897932385D0

  private

  public im, jm, nn, mm, nm                   ! 格子点数, 切断波数
  public pz, rm, cm                           ! 変換用作業配列

  public w_base_Initial                       ! 初期化サブルーチン
  public w_base_Finalize                      ! 終了サブルーチン

  public x_Lon, y_Lat                         ! 格子座標
  public x_Lon_weight, y_Lat_Weight           ! 格子座標重み
  public xy_Lon, xy_Lat                       ! 格子座標(im,jm)
  public l_nm, nm_l                           ! 波数格納位置
  public xy_w, w_xy                           ! 変換関数
  public w_DLon_w                             ! 経度微分

  interface l_nm
     module procedure l_nm_array00
     module procedure l_nm_array01
     module procedure l_nm_array10
     module procedure l_nm_array11
     module procedure l_nm_darray00
     module procedure l_nm_darray01
     module procedure l_nm_darray10
     module procedure l_nm_darray11
  end interface

  interface nm_l
     module procedure nm_l_int
     module procedure nm_l_array
     module procedure nm_l_dint
     module procedure nm_l_darray
  end interface

  save im, jm, nm, mm, nn, ig                 ! 格子点数, 切断波数を記憶
  save pz, rm                                 ! 変換用配列を記憶
  save cm                                     ! 変換用配列の大きさ  
  save w_zonal_base_initialize                ! 初期化フラグ

contains
  !--------------- 初期化 -----------------
  subroutine w_base_Initial(n_in,i_in,j_in)
    !
    ! スペクトル変換の格子点数, 波数を設定する.
    !
    ! 実際の使用には上位サブルーチン w_zonal_Initial を用いること.
    !
    integer,intent(in) :: i_in              !(in) 格子点数(東西), 2の巾乗(<=2048)
    integer,intent(in) :: j_in              !(in) 格子点数(南北), 4 の倍数
    integer,intent(in) :: n_in              !(in) 切断全波数

    integer :: i, j

    im = i_in   ; jm = j_in
    nn = n_in   ; nm = n_in+1 ;  mm = n_in    ! default は三角波数切断

    allocate(x_Lon(0:im-1))                   ! 格子点座標格納配列(経度)
    allocate(x_Lon_Weight(0:im-1))
    allocate(y_Lat(1:jm))
    allocate(y_Lat_Weight(1:jm))              ! 格子点座標格納配列
    allocate(xy_Lon(0:im-1,1:jm))             ! 格子点座標格納配列
    allocate(xy_Lat(0:im-1,1:jm))             ! 格子点座標格納配列

    allocate(pz(jm/2,5))                      ! 変換用配列
    allocate(rm(nm/2*3+nm+1))                 ! 変換用配列
    allocate(cm(2*nn+1))                      ! 変換用配列      


    do i=0,im-1
       x_Lon(i)  = 2*pi/im*i                  ! 経度座標
       x_Lon_Weight(i) = 2*pi/im              ! 経度座標重み
    enddo

    call lxinig(jm,pz,ig)
    call lxinir(nm,0_8,rm)
    call lxinic(nn,0_8,cm)

    do j=1,jm/2
       y_Lat(jm/2+j)   =  asin(pz(j,1))        ! 緯度座標
       y_Lat(jm/2-j+1) = -asin(pz(j,1))        ! 緯度座標
       y_Lat_Weight(jm/2+j)   = 2*pz(j,2)      ! 緯度重み(Gauss grid)
       y_Lat_Weight(jm/2-j+1) = 2*pz(j,2)      ! 緯度重み(Gauss grid)
    enddo

    do j=1,jm
       xy_Lon(:,j) = x_Lon
    enddo

    do i=0,im-1
       xy_Lat(i,:) = y_Lat
    enddo

    w_zonal_base_initialize = .true.

    call MessageNotify('M','w_base_initial',&
         'w_zonal_module_base (2025/12/10) is initialized')

  end subroutine w_base_Initial

  !--------------- 基本関数 -----------------

  function l_nm_array00(n_in,m_in)
    !
    ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
    ! 
    ! 引数 n,m がともに整数値の場合, 整数値を返す. 
    !
    integer               :: l_nm_array00
    !(out) スペクトルデータの格納位置 

    integer, intent(in)   :: n_in     !(in) 全波数
    integer, intent(in)   :: m_in     !(in) 帯状波数           

    if ( m_in /= 0 ) then
       call MessageNotify('E','l_nm_array00', &
            'longitudinal wavenumber m must be zero')
    else
       l_nm_array00 = n_in-m_in+1
    endif

  end function l_nm_array00

  function l_nm_darray00(n_in,m_in)
    !
    ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
    ! 
    ! 引数 n,m がともに整数値の場合, 整数値を返す. 
    !
    integer               :: l_nm_darray00   
    !(out) スペクトルデータの格納位置 

    integer(8), intent(in)   :: n_in     !(in) 全波数
    integer(8), intent(in)   :: m_in     !(in) 帯状波数           

    if ( m_in /= 0 ) then
       call MessageNotify('E','l_nm_darray00', &
            'longitudinal wavenumber m must be zero')
    else
    l_nm_darray00 = n_in+m_in+1
    endif

  end function l_nm_darray00

  function l_nm_array01(n,marray)           ! スペクトルデータの格納位置 
    !
    ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
    ! 
    ! 第 1 引数 n が整数, 第 2 引数 marray が整数 1 次元配列の場合, 
    ! marray と同じ大きさの 1 次元整数配列を返す. 
    !
    integer, intent(in)  :: n               !(in) 全波数
    integer, intent(in)  :: marray(:)       !(in) 帯状波数
    integer              :: l_nm_array01(size(marray))
    !(out) スペクトルデータ位置

    integer              :: i 

    do i=1, size(marray)
       l_nm_array01(i) = l_nm_array00(n,marray(i))
    enddo
  end function l_nm_array01

  function l_nm_darray01(n,marray)           ! スペクトルデータの格納位置 
    !
    ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
    ! 
    ! 第 1 引数 n が整数, 第 2 引数 marray が整数 1 次元配列の場合, 
    ! marray と同じ大きさの 1 次元整数配列を返す. 
    !
    integer(8), intent(in)  :: n               !(in) 全波数
    integer(8), intent(in)  :: marray(:)       !(in) 帯状波数
    integer                 :: l_nm_darray01(size(marray))
    !(out) スペクトルデータ位置

    integer              :: i 

    do i=1, size(marray)
       l_nm_darray01(i) = l_nm_darray00(n,marray(i))
    enddo
  end function l_nm_darray01

  function l_nm_array10(narray,m_in)
    !
    ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
    ! 
    ! 第 1 引数 narray が整数 1 次元配列, 第 2 引数  m が整数の場合, 
    ! narray と同じ大きさの 1 次元整数配列を返す. 
    !
    integer, intent(in)  :: narray(:)           !(in) 全波数  
    integer, intent(in)  :: m_in                !(in) 帯状波数
    integer              :: l_nm_array10(size(narray))
    !(out) スペクトルデータ位置

    integer              :: i 

    do i=1, size(narray)
       l_nm_array10(i) = l_nm_array00(narray(i),m_in)
    enddo
  end function l_nm_array10

  function l_nm_darray10(narray,m_in)
    !
    ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
    ! 
    ! 第 1 引数 narray が整数 1 次元配列, 第 2 引数  m が整数の場合, 
    ! narray と同じ大きさの 1 次元整数配列を返す. 
    !
    integer(8), intent(in)  :: narray(:)           !(in) 全波数  
    integer(8), intent(in)  :: m_in                !(in) 帯状波数
    integer                 :: l_nm_darray10(size(narray))
    !(out) スペクトルデータ位置

    integer              :: i 

    do i=1, size(narray)
       l_nm_darray10(i) = l_nm_darray00(narray(i),m_in)
    enddo
  end function l_nm_darray10

  function l_nm_array11(narray,marray)
    !
    ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
    ! 
    ! 第 1,2 引数 narray, marray がともに整数 1 次元配列の場合, 
    ! narray, marray と同じ大きさの 1 次元整数配列を返す. 
    ! narray, marray は同じ大きさでなければならない. 
    !
    integer, intent(in)  :: narray(:)          !(in) 全波数  
    integer, intent(in)  :: marray(:)          !(in) 帯状波数
    integer              :: l_nm_array11(size(narray))
    !(out) スペクトルデータ位置

    integer              :: i 

    if ( size(narray) .ne. size(marray) ) then
       call MessageNotify('E','l_nm_array11',&
            'dimensions of input arrays  n and m are different.')
    endif

    do i=1, size(narray)
       l_nm_array11(i) = l_nm_array00(narray(i),marray(i))
    enddo
  end function l_nm_array11

  function l_nm_darray11(narray,marray)
    !
    ! 全波数(n)と東西波数(m)からそのスペクトルデータの格納位置を返す.
    ! 
    ! 第 1,2 引数 narray, marray がともに整数 1 次元配列の場合, 
    ! narray, marray と同じ大きさの 1 次元整数配列を返す. 
    ! narray, marray は同じ大きさでなければならない. 
    !
    integer(8), intent(in)  :: narray(:)          !(in) 全波数  
    integer(8), intent(in)  :: marray(:)          !(in) 帯状波数
    integer                 :: l_nm_darray11(size(narray))
    !(out) スペクトルデータ位置

    integer              :: i 

    if ( size(narray) .ne. size(marray) ) then
       call MessageNotify('E','l_nm_darray11',&
            'dimensions of input arrays  n and m are different.')
    endif

    do i=1, size(narray)
       l_nm_darray11(i) = l_nm_darray00(narray(i),marray(i))
    enddo
  end function l_nm_darray11

  function nm_l_int(l)
    ! 
    ! スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す.
    !
    ! 引数 l が整数値の場合, 対応する全波数と帯状波数を
    ! 長さ 2 の 1 次元整数値を返す. 
    ! nm_l(1) が全波数, nm_l(2) が帯状波数である. 
    !
    integer               :: nm_l_int(2)  !(out) 全波数, 帯状波数
    integer, intent(in)   :: l            !(in) スペクトルデータの格納位置

    nm_l_int(1) = l-1
    nm_l_int(2) = 0

  end function nm_l_int

  function nm_l_dint(l)
    ! 
    ! スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す.
    !
    ! 引数 l が整数値の場合, 対応する全波数と帯状波数を
    ! 長さ 2 の 1 次元整数値を返す. 
    ! nm_l(1) が全波数, nm_l(2) が帯状波数である. 
    !
    integer                  :: nm_l_dint(2)  !(out) 全波数, 帯状波数
    integer(8), intent(in)   :: l             !(in) スペクトルデータの格納位置

    nm_l_dint(1) = l-1
    nm_l_dint(2) = 0

  end function nm_l_dint

  function nm_l_array(larray)
    ! 
    ! スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す.
    !
    ! 引数 larray が整数 1 次元配列の場合, 
    ! larray に対応する n, m を格納した 2 次元整数配列を返す. 
    ! nm_l_array(:,1) が全波数, nm_l_array(:,2) が帯状波数である. 
    !
    integer, intent(in)  :: larray(:)
    !(out) 全波数, 帯状波数

    integer              :: nm_l_array(size(larray),2)
    !(in) スペクトルデータの格納位置

    integer              :: i

    do i=1, size(larray)
       nm_l_array(i,:) = nm_l_int(larray(i))
    enddo
  end function nm_l_array

  function nm_l_darray(larray)
    ! 
    ! スペクトルデータの格納位置(l)から全波数(n)と東西波数(m)を返す.
    !
    ! 引数 larray が整数 1 次元配列の場合, 
    ! larray に対応する n, m を格納した 2 次元整数配列を返す. 
    ! nm_l_array(:,1) が全波数, nm_l_array(:,2) が帯状波数である. 
    !
    integer(8), intent(in)  :: larray(:)
    !(out) 全波数, 帯状波数

    integer                 :: nm_l_darray(size(larray),2)
    !(in) スペクトルデータの格納位置

    integer              :: i

    do i=1, size(larray)
       nm_l_darray(i,:) = nm_l_dint(larray(i))
    enddo
  end function nm_l_darray

  !--------------- 基本変換 -----------------
  function w_DLon_w(w_data)
    !
    ! 入力スペクトルデータに経度微分∂/∂λを作用する. 
    !
    ! スペクトルデータの経度微分とは, 対応する格子点データに
    ! 経度微分を作用させたデータのスペクトル変換のことである. 
    !
    real(8)              :: w_DLon_w(nn+1)
    !(out) 入力スペクトルデータのラプラシアン

    real(8), intent(in)  :: w_data(nn+1)
    !(in) 入力スペクトルデータ

    integer :: n, n1, n2

    w_DLon_w = 0.0D0

  end function w_DLon_w

  function xy_w(w_data,ipow,iflag)
    !
    ! スペクトルデータから格子データへ変換する(1 層用).
    !
    real(8)               :: xy_w(0:im-1,1:jm)
    !(out) 格子点データ

    real(8), intent(in)   :: w_data(nn+1)
    !(in) スペクトルデータ

    integer, intent(in), optional  :: ipow      
    !(in) 作用させる 1/cosφ の次数. 省略時は 0. 

    integer, intent(in), optional  :: iflag
    !(in) 変換の種類
    !    0 : 通常の正変換
    !   -1 : 経度微分を作用させた逆変換
    !    1 : 緯度微分 cosφ・∂/∂φ を作用させた逆変換
    !    2 : sinφを作用させた逆変換
    !    省略時は 0.
    !
    integer, parameter  :: ipow_default  = 0
    integer, parameter  :: iflag_default = 0

    integer(8) ipval, ifval

    real(8)             :: w_Ydata(nn+2)
    ! 作業用スペクトルデータ(LXCSWY 出力用)

    if ( .not. w_zonal_base_initialize ) then
       call MessageNotify('E','xy_w',&
            'w_zonal_base_module_svpack not initialize yet.')
    endif

    if (present(ipow)) then
       ipval = ipow
    else
       ipval = ipow_default
    endif

    if (present(iflag)) then
       ifval = iflag
    else
       ifval = iflag_default
    endif

    if ( ifval==0 ) then
       call lxtszg(nm,nn,jm,w_Data,xy_w,pz,rm,ipval)
    else if( ifval==-1 ) then
       call lxtszg(nm,nn,jm,w_DLon_w(w_Data),xy_w,pz,rm,ipval)
    else if( ifval==1 ) then
       call lxcszy(nn,w_data,w_Ydata,cm)
       call lxtszg(nm,nm,jm,w_YData,xy_w,pz,rm,ipval)
    else if( ifval==2 ) then
       call lxtszg(nm,nn,jm,w_data,xy_w,pz,rm,ipval)
       xy_w = xy_w * sin(xy_Lat)
    else
       call MessageNotify('E','xy_w','invalid value of iflag')
    endif

  end function xy_w

  function w_xy(xy_data,ipow,iflag)
    !
    ! 格子データからスペクトルデータへ(正)変換する(1 層用).
    !
    real(8)               :: w_xy(nn+1)
    !(out) スペクトルデータ

    real(8), intent(in)   :: xy_data(0:im-1,1:jm)
    !(in) 格子点データ

    integer, intent(in), optional  :: ipow
    !(in) 変換時に同時に作用させる 1/cosφ の次数. 省略時は 0.

    integer, intent(in), optional  :: iflag
    ! 変換の種類
    !    0 : 通常の正変換
    !   -1 : 経度微分を作用させた正変換 
    !    1 : 緯度微分 1/cosφ・∂(f cos^2φ)/∂φ を作用させた正変換
    !    2 : sinφを作用させた正変換
    !  省略時は 0.


    integer, parameter  :: ipow_default  = 0    ! スイッチデフォルト値
    integer, parameter  :: iflag_default = 0    ! スイッチデフォルト値

    integer(8) ipval, ifval

    real(8)             :: w_Ydata(2*(nn+2))
    ! 作業用スペクトルデータ(LXCYWS 出力用)

    if ( .not. w_zonal_base_initialize ) then
       call MessageNotify('E','w_xy',&
            'w_zonal_module_base not initialize yet.')
    endif

    if (present(ipow)) then
       ipval = ipow
    else
       ipval = ipow_default
    endif

    if (present(iflag)) then
       ifval = iflag
    else
       ifval = iflag_default
    endif

    if ( ifval == 0 ) then
       call lxtgzs(nm,nn,jm,w_xy,xy_data,pz,rm,ipval)
    else if ( ifval == -1 ) then
       call lxtgzs(nm,nn,jm,w_xy,xy_data,pz,rm,ipval)
       w_xy = w_Dlon_w(w_xy)
    else if ( ifval == 1 ) then
       call lxtgzs(nm,nm,jm,w_Ydata,xy_data,pz,rm,ipval)
       call lxcyzs(nn,w_Ydata,w_xy,cm)
    else if ( ifval == 2 ) then
       call lxtgzs(nm,nn,jm,w_xy,xy_data*sin(xy_Lat),pz,rm,ipval)
    end if

  end function w_xy

  !--------------- 終了処理 -----------------
  subroutine w_base_Finalize
    !
    ! モジュールの終了処理(割り付け配列の解放)をおこなう. 
    !
    ! 実際の使用には上位サブルーチン w_Finalize を用いること.
    !
    if ( .not. w_zonal_base_initialize ) then
       call MessageNotify('W','w_base_Finalize',&
            'w_zonal_module_base not initialized yet')
       return
    endif

    deallocate(pz)                  ! 変換用配列
    deallocate(rm)                  ! 変換用配列
    deallocate(cm)                  ! 変換用作業配列      

    deallocate(x_Lon)
    deallocate(x_Lon_Weight)       ! 格子点座標格納配列
    deallocate(xy_Lon)             ! 格子点座標格納配列
    deallocate(y_Lat)
    deallocate(y_Lat_Weight)       ! 格子点座標格納配列
    deallocate(xy_Lat)             ! 格子点座標格納配列

    w_zonal_base_initialize = .false.

    call MessageNotify('M','w_base_Finalize',&
         'w_zonal_module_base (2025/12/10) is finalized')

  end subroutine w_base_Finalize

end module w_zonal_module_base