! -*- mode: f90; coding: utf-8 -*-
!-------------------------------------------------------------------------
! Copyright (c) 2019-2020 SPMODEL Development Group. All rights reserved.
!
! 履歴  2019/03/13  竹広真一
!       2019/10/04  佐々木洋平
!       2020/08/31  佐々木洋平 pointer 版に変更
!       2020/11/06  竹広真一   セクター計算オプション導入
!-------------------------------------------------------------------------
!> @copyright Copyright (C) SPMODEL Development Group, 2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @brief
!> @ref w_module の下部モジュール: 基本変換
!>
!> @warning
!> 本モジュールを直接呼ぶ事は想定されていないことに注意されたい．
!> ユーザは @ref w_module を呼ぶこと．
!>
module w_module_base_mint
  use dc_types, only: DP
  use dc_message, only: MessageNotify
  use spml_utils, only: pi
  use iso_c_binding, only: c_ptr, c_f_pointer

  implicit none
  private
  public im
  public jm
  public mm
  public nm
  public nn
  public mi
  public nc
  public it
  public t
  public p
  public r
  public jc
  public w_base_Initial
  public w_base_Finalize
  public x_Lon
  public y_Lat
  public x_Lon_Weight
  public y_Lat_Weight
  public xy_Lon
  public xy_Lat
  public l_nm
  public nm_l
  public xy_w
  public w_xy
  public w_StreamPotential2Vector
  public w_Vector2VorDiv
  public w_VectorCosLat2VorDiv

  !> 緯度座標
  real(DP), allocatable    :: x_Lon(:)
  !> 緯度座標重み
  real(DP), allocatable    :: x_Lon_Weight(:)
  !> 経度座標
  real(DP), allocatable    :: y_Lat(:)
  !> 経度座標重み
  real(DP), allocatable    :: y_Lat_Weight(:)
  !> 緯度座標(2次元配列)
  real(DP), allocatable    :: xy_Lon(:,:)
  !> 経度座標(2次元配列)
  real(DP), allocatable    :: xy_Lat(:,:)
  !> 格子点の設定(東西)
  integer(8)               :: im=64
  !> 格子点の設定(南北)
  integer(8)               :: jm=32
  !> m の切断波数
  integer(8)               :: mm
  !> n の切断波数
  integer(8)               :: nn
  !> n の最大値
  integer(8)               :: nm
  !> 経度対称性
  integer(8)               :: mi
  !> スペクトルデータの大きさ
  integer                  :: nc
  !> 変換用配列
  integer(8), allocatable  :: it(:)
  !> 変換用配列
  integer(8), allocatable  :: jc(:)
  !> 変換用配列
  real(DP), allocatable    :: t(:)
  !> 変換用配列
  real(DP), allocatable    :: r(:)
  !> 変換用配列
  real(DP), pointer        :: p(:,:)

  real(DP), pointer        :: w(:)
  real(DP), pointer        :: g(:)
  real(DP), pointer        :: g1(:)
  real(DP), pointer        :: g2(:)
  real(DP), pointer        :: ww(:)
  integer(8)               :: ig=1
  type(c_ptr)              :: pp
  type(c_ptr)              :: pg
  type(c_ptr)              :: pw
  type(c_ptr)              :: pww
  type(c_ptr)              :: pg1
  type(c_ptr)              :: pg2
  real(DP), allocatable    :: c(:)
  logical                  :: w_base_initialize=.false.

  save im, jm, mm, nn, nm, mi, nc
  save it, t, p, r, jc, c
  save w_base_initialize

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

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

contains

  !----------------------------------------------------------
  !> スペクトル変換の格子点数, 波数および OPENMP 使用時の
  !> 最大スレッド数を設定する.
  !>
  !> @warning
  !> 実際の使用には上位サブルーチン
  !> [w_Initial](@ref w_module#w_initial) を用いること.
  !>
  subroutine w_base_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

    integer :: i, j
    im = i_in
    jm = j_in
    nn = n_in
    mm = nn
    nm = nn+1

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

    if ( present(mint) )then
       mi = mint
    else
       mi = 1_8
    end if

    allocate(t(im*3/2))
    allocate(it(im/2))
    allocate(r(5*(mm/mi+1)*(2*nm-mm/mi*mi)/4+mm/mi+1))
    allocate(jc((mm/mi+1)*(2*nm-mm/mi*mi)/16+mm/mi+1))
    call mxallc(pp,jm/2*(5+2*(mm/mi+1)))
    call c_f_pointer(pp,p,[jm/2,5+2*(mm/mi+1)])

    call mxallc(pg,jm*im)
    call mxallc(pw,jm*im)
    call mxallc(pg1,jm*im)
    call mxallc(pg2,jm*im)
    call mxallc(pww,jm*im*2)
    call c_f_pointer(pg, g, [jm*im])
    call c_f_pointer(pw, w, [jm*im])
    call c_f_pointer(pg1, g1, [jm*im])
    call c_f_pointer(pg2, g2, [jm*im])
    call c_f_pointer(pww, ww, [jm*im]*2)

    allocate(c((mm/mi+1)*(2*(nn+1)-mm/mi*mi)))
    allocate(x_Lon(0:im-1))
    allocate(x_Lon_Weight(0:im-1))
    allocate(xy_Lon(0:im-1,1:jm))
    allocate(y_Lat(1:jm))
    allocate(y_Lat_Weight(1:jm))
    allocate(xy_Lat(0:im-1,1:jm))

    call sxinm1(mm,nm,im,it,t,r,mi)
    call sxinm2(mm,nm,jm,ig,p,r,jc,mi)
    call sxinmc(mm,nn,c,mi)

    ! 分散スペクトルデータの大きさ
    nc = int((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    do i=0,int(im)-1
      x_Lon(i)        = 2*pi/im*i/mi
      x_Lon_Weight(i) = 2*pi/im/mi
    enddo

    do j=1,int(jm)/2
      y_Lat(jm/2+j)          =   asin(p(j,1))
      y_Lat(jm/2-j+1)        = - asin(p(j,1))
      y_Lat_Weight(jm/2+j)   =   2*p(j,2)
      y_Lat_Weight(jm/2-j+1) =   2*p(j,2)
    enddo

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

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

    w_base_initialize = .true.

    call MessageNotify('M','w_base_initial',&
      'w_base_module_mint (2020/11/18) is initialized')

  end subroutine w_base_Initial

  !----------------------------------------------------------
  !> スペクトルデータから格子データへ変換する(1 層用).
  !>
  function xy_w(w_data,ipow,iflag)
    !> 格子点データ
    real(DP)               :: xy_w(0:im-1,1:jm)
    !> スペクトルデータ
    real(DP), intent(in)   :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> 作用させる \f$1/\cosφ\f$ の次数. 省略時は 0.
    integer, intent(in), optional  :: ipow
    !> 変換の種類. 省略時は 0
    !> *   0 : 通常の逆変換(default)
    !> *  -1 : 経度微分を作用させた逆変換
    !> *   1 : 緯度微分 \f$\cosφ \cdot ∂/∂ φ\f$
    !>         を作用させた逆変換
    !> *   2 : \f$\sinφ\f$を作用させた逆変換
    integer, intent(in), optional  :: iflag
    !
    integer, parameter  :: ipow_default  = 0
    integer, parameter  :: iflag_default = 0

    real(DP)             :: w_Xdata((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    real(DP)             :: w_Ydata((mm/mi+1)*(2*(nn+2)-mm/mi*mi))

    integer(8) ipval, ifval

    if ( .not. w_base_initialize ) then
      call MessageNotify('E','xy_w',&
        'w_base_module 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 sxtsmg(mm,nm,nn,im,jm,w_data,g,it,t,p,r,jc,w,ipval,mi)
      xy_w = reshape(g,(/im,jm/))
    else if( ifval == -1 ) then
      call sxcsmx(mm,nn,w_data,w_Xdata,mi)
      call sxtsmg(mm,nm,nn,im,jm,w_Xdata,g,it,t,p,r,jc,w,ipval,mi)
      xy_w = reshape(g,(/im,jm/))
    else if( ifval == 1 ) then
      call sxcsmy(mm,nn,w_data,w_Ydata,c,mi)
      call sxtsmg(mm,nm,nn+1_8,im,jm,w_Ydata,g,it,t,p,r,jc,w,ipval,mi)
      xy_w = reshape(g,(/im,jm/))
    else if( ifval == 2 ) then
      call sxtsmg(mm,nm,nn,im,jm,w_data,g,it,t,p,r,jc,w,ipval,mi)
      xy_w = reshape(g,(/im,jm/)) * sin(xy_Lat)
    else
      call MessageNotify('E','xy_w','invalid value of iflag')
    endif

  end function xy_w

  !----------------------------------------------------------
  !> 格子データからスペクトルデータへ(正)変換する(1 層用).
  !>
  function w_xy(xy_data,ipow,iflag)
    !> 格子点データ
    real(DP), intent(in)   :: xy_data(0:im-1,1:jm)
    !> 変換時に同時に作用させる \f$1/\cosφ\f$ の次数. 省略時は 0.
    integer, intent(in), optional  :: ipow
    !> 変換の種類. 省略時は 0
    !> *   0 : 通常の正変換(default)
    !> *  -1 : 経度微分を作用させた正変換
    !> *   1 : 緯度微分 \f$\cosφ \cdot ∂/∂φ\f$
    !>         を作用させた正変換
    !> *   2 : \f$\sinφ\f$を作用させた正変換
    integer, intent(in), optional  :: iflag
    !> スペクトルデータ
    real(DP)               :: w_xy((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    integer, parameter  :: ipow_default  = 0
    integer, parameter  :: iflag_default = 0
    integer(8) ipval, ifval
    real(DP)             :: w_Xdata((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    real(DP)             :: w_Ydata((mm/mi+1)*(2*(nn+2)-mm/mi*mi))

    if ( .not. w_base_initialize ) then
      call MessageNotify('E','w_xy',&
        'w_base_module 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
      g = reshape(xy_data,(/im*jm/))
      call sxtgms(mm,nm,nn,im,jm,w_xy,g,it,t,p,r,jc,w,ipval,mi)
    else if ( ifval == -1 ) then
      g = reshape(xy_data,(/im*jm/))
      call sxtgms(mm,nm,nn,im,jm,w_Xdata,g,it,t,p,r,jc,w,ipval,mi)
      call sxcsmx(mm,nn,w_Xdata,w_xy,mi)
    else if ( ifval == 1 ) then
      g = reshape(xy_data,(/im*jm/))
      call sxtgms(mm,nm,nn+1,im,jm,w_Ydata,g,it,t,p,r,jc,w,ipval,mi)
      call sxcyms(mm,nn,w_Ydata,w_xy,c,mi)
    else if ( ifval == 2 ) then
      g = reshape(xy_data*sin(xy_Lat),(/im*jm/))
      call sxtgms(mm,nm,nn,im,jm,w_xy,g,it,t,p,r,jc,w,ipval,mi)
    else
      call MessageNotify('E','w_xy','invalid value of iflag')
    endif

  end function w_xy

  !----------------------------------------------------------
  !> 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に
  !> (逆)変換する(1 層用)
  !>
  !> スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
  !> \f{align*}{
  !>   u \cosφ & =       ∂χ/∂λ - \cosφ∂ψ/∂φ, \\
  !>   v \cosφ & = \cosφ∂χ/∂φ +       ∂ψ/∂λ
  !> \f}
  !>
  subroutine w_StreamPotential2Vector(w_Psi, w_Chi, xy_U, xy_V)
    !> 流線関数
    real(DP), intent(in)   :: w_Psi((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> 速度ポテンシャル
    real(DP), intent(in)   :: w_Chi((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> 速度経度成分
    real(DP), intent(out)   :: xy_U(0:im-1,1:jm)
    !> 速度緯度成分
    real(DP), intent(out)   :: xy_V(0:im-1,1:jm)

    real(DP)             :: w_Xdata((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    real(DP)             :: w_Ydata((mm/mi+1)*(2*(nn+2)-mm/mi*mi))
    real(DP)             :: w_RData1((mm/mi+1)*(2*(nn+2)-mm/mi*mi))
    real(DP)             :: w_RData2((mm/mi+1)*(2*(nn+2)-mm/mi*mi))

    if ( .not. w_base_initialize ) then
      call MessageNotify('E','w_StreamPotential2Vector',&
        'w_base_module not initialize yet.')
    endif
    !
    ! $ u \cosφ = ∂χ/∂λ - \cosφ∂ψ/∂φ $
    !
    call sxcsmx(mm,nn,w_Chi,w_Xdata,mi)
    call sxcsmy(mm,nn,w_Psi,w_Ydata,c,mi)
    call sxcmpk(mm,nn,nm,w_Xdata,w_Rdata1,mi)
    w_Rdata1 = w_Rdata1 - w_Ydata
    !
    ! $ v \cosφ = \cosφ∂χ/∂φ + ∂ψ/∂λ $
    !
    call sxcsmy(mm,nn,w_Chi,w_Ydata,c,mi)
    call sxcsmx(mm,nn,w_Psi,w_Xdata,mi)
    call sxcmpk(mm,nn,nm,w_Xdata,w_Rdata2,mi)
    w_Rdata2= w_Rdata2 + w_Ydata
    !
    ! u,v
    !
    call sxtsmv(mm,nm,nm,im,jm,w_Rdata1,w_Rdata2,g1,g2,it,t,p,r,jc,ww,1_8,mi)
    xy_U = reshape(g1,(/im,jm/))
    xy_V = reshape(g2,(/im,jm/))

  end subroutine w_StreamPotential2Vector

  !----------------------------------------------------------
  !> 速度場(格子データ)から渦度・発散(スペクトルデータ)に
  !> (正)変換する(1 層用)
  !>
  !> スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
  !> \f{align}{
  !>   ζ = 1/\cosφ∂v/∂λ - 1/\cosφ ∂(u \cosφ)/∂φ, \\
  !>    D = 1/\cosφ∂u/∂λ + 1/\cosφ ∂(v \cosφ)/∂φ
  !> \f}
  !>
  subroutine w_Vector2VorDiv(xy_U, xy_V, w_Vor, w_Div)
    !> 速度経度成分
    real(DP), intent(in)   :: xy_U(0:im-1,1:jm)
    !> 速度緯度成分
    real(DP), intent(in)   :: xy_V(0:im-1,1:jm)
    !> 流線関数
    real(DP), intent(out)   :: w_Vor((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> 速度ポテンシャル
    real(DP), intent(out)   :: w_Div((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    real(DP)             :: w_Xdata((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    real(DP)             :: w_Ydata1((mm/mi+1)*(2*(nn+2)-mm/mi*mi))
    real(DP)             :: w_Ydata2((mm/mi+1)*(2*(nn+2)-mm/mi*mi))
    real(DP)             :: w_Data1((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    real(DP)             :: w_Data2((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    if ( .not. w_base_initialize ) then
      call MessageNotify('E','w_Vector2VorDiv',&
        'w_base_module not initialize yet.')
    endif

    !
    ! $ 1/\cosφ∂u/∂λ, 1/\cosφ ∂(u \cosφ)/∂φ$
    !
    g1 = reshape(xy_U,(/im*jm/))
    g2 = reshape(xy_V,(/im*jm/))
    call sxtvms(mm,nm,nn+1,im,jm,w_Ydata1,w_Ydata2,g1,g2,it,t,p,r,jc,ww,1_8,mi)

    call sxcmpk(mm,nm,nn,w_Ydata1,w_Xdata,mi)
    call sxcsmx(mm,nn,w_Xdata,w_Div,mi)
    call sxcyms(mm,nn,w_Ydata2,w_Data2,c,mi)
    !
    ! $ 1/\cosφ∂v/∂λ, 1/\cosφ ∂(v \cosφ)/∂φ $
    !
    call sxcmpk(mm,nm,nn,w_Ydata2,w_Xdata,mi)
    call sxcsmx(mm,nn,w_Xdata,w_Vor,mi)
    call sxcyms(mm,nn,w_Ydata1,w_Data1,c,mi)
    !
    ! Vorticity:
    ! $ ζ = 1/\cosφ∂v/∂λ - 1/\cosφ ∂(u cosφ)/∂φ $
    !
    w_Vor = w_Vor - w_Data1
    !
    ! Divergence:
    ! $  D = 1/\cosφ∂u/∂λ + 1/\cosφ ∂(v \cosφ)/∂φ $
    !
    w_Div = w_Div + w_Data2

  end subroutine w_Vector2VorDiv

  !----------------------------------------------------------
  !> 速度場(格子データ)から渦度・発散(スペクトルデータ)に
  !> (正)変換する(1 層用)
  !>
  !> スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
  !> \f{align*}{
  !>   ζ = 1/\cosφ∂v/∂λ - 1/\cosφ ∂(u \cosφ)/∂φ,\\
  !>    D = 1/\cosφ∂u/∂λ + 1/\cosφ ∂(v \cosφ)/∂φ
  !> \f}
  !>
  subroutine w_VectorCosLat2VorDiv(xy_UCosLat, xy_VCosLat, w_Vor, w_Div)
    !> 速度経度成分 * cos(lat)
    real(DP), intent(in)   :: xy_UCosLat(0:im-1,1:jm)
    !> 速度緯度成分 * cos(lat)
    real(DP), intent(in)   :: xy_VCosLat(0:im-1,1:jm)
    !> 流線関数
    real(DP), intent(out)   :: w_Vor((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> 速度ポテンシャル
    real(DP), intent(out)   :: w_Div((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    real(DP)             :: w_Xdata((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    real(DP)             :: w_Ydata1((mm/mi+1)*(2*(nn+2)-mm/mi*mi))
    real(DP)             :: w_Ydata2((mm/mi+1)*(2*(nn+2)-mm/mi*mi))
    real(DP)             :: w_Data1((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    real(DP)             :: w_Data2((mm/mi+1)*(2*(nn+1)-mm/mi*mi))

    if ( .not. w_base_initialize ) then
      call MessageNotify('E','w_VectorCosLat2VorDiv',&
        'w_base_module not initialize yet.')
    endif
    !
    ! $1/\cosφ∂u/∂λ, 1/\cosφ∂(u cosφ)/∂φ$
    !
    g1 = reshape(xy_UCosLat,(/im*jm/))
    g2 = reshape(xy_VCosLat,(/im*jm/))
    call sxtvms(mm,nm,nn+1,im,jm,w_Ydata1,w_Ydata2,g1,g2,it,t,p,r,jc,ww,2_8,mi)

    call sxcmpk(mm,nm,nn,w_Ydata1,w_Xdata,mi)
    call sxcsmx(mm,nn,w_Xdata,w_Div,mi)
    call sxcyms(mm,nn,w_Ydata2,w_Data2,c,mi)
    !
    ! $1/\cosφ∂v/∂λ, 1/\cosφ∂(v cosφ)/∂φ$
    !
    call sxcmpk(mm,nm,nn,w_Ydata2,w_Xdata,mi)
    call sxcsmx(mm,nn,w_Xdata,w_Vor,mi)
    call sxcyms(mm,nn,w_Ydata1,w_Data1,c,mi)
    !
    ! Vorticity:
    !   $ζ = 1/\cosφ∂v/∂λ - 1/\cosφ∂(u \cosφ)/∂φ$
    !
    w_Vor = w_Vor - w_Data1
    !
    ! Divergence:
    !   $D = 1/\cosφ∂u/∂λ + 1/\cosφ∂(v \cosφ)/∂φ$
    !
    w_Div = w_Div + w_Data2

  end subroutine w_VectorCosLat2VorDiv

  !----------------------------------------------------------
  !> 終了処理(割り付け配列の解放)をおこなう
  !>
  subroutine w_base_Finalize
    if ( .not. w_base_initialize ) then
      call MessageNotify('W','w_base_Finalize',&
        'w_base_module not initialized yet')
      return
    endif

    deallocate(it)
    deallocate(t)
    deallocate(r)
    deallocate(jc)
    deallocate(c)

    call mxfree(pp)
    call mxfree(pg)
    call mxfree(pg1)
    call mxfree(pg2)
    call mxfree(pw)
    call mxfree(pww)

    deallocate(x_Lon)
    deallocate(x_Lon_weight)
    deallocate(xy_Lon)
    deallocate(y_Lat)
    deallocate(y_Lat_Weight)
    deallocate(xy_Lat)
    w_base_initialize = .false.

    call MessageNotify('M','w_base_Finalize',&
      'w_base_module_mint (2020/11/06) is finalized')

  end subroutine w_base_Finalize

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

    integer(8)   :: nd     ! 全波数
    integer(8)   :: md     ! 帯状波数
    integer(8)   :: ld     ! スペクトルデータの格納位置

    nd = n
    md = m
    call sxnmml(nn,nd,md,ld,mi)
    l_nm_array00 = int(ld)

  end function l_nm_array00

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

    integer(8)   :: nd     ! 全波数
    integer(8)   :: md     ! 帯状波数
    integer(8)   :: ld     ! スペクトルデータの格納位置

    nd = n
    md = m
    call sxnmml(nn,nd,md,ld,mi)
    l_nm_darray00 = int(ld)

  end function l_nm_darray00

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

    integer              :: i

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

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

    integer              :: i

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

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

    integer              :: i

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

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

    integer              :: i

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


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

    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

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

    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_darray11(i) = l_nm_darray00(narray(i),marray(i))
    enddo
  end function l_nm_darray11

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

    integer(8)            :: ld
    integer(8)            :: nd
    integer(8)            :: md

    ld = l
    call sxlmnm(mm,nn,ld,nd,md,mi)
    nm_l_int = int((/nd,md/))

  end function nm_l_int

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

    integer(8)             :: ld
    integer(8)             :: nd
    integer(8)             :: md

    ld = l
    call sxlmnm(mm,nn,ld,nd,md,mi)
    nm_l_dint = int((/nd,md/))

  end function nm_l_dint

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

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

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

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

end module w_module_base_mint