! -*- mode: f90; coding: utf-8 -*-
!-------------------------------------------------------------------------
! Copyright (c) 2019 SPMODEL Development Group. All rights reserved.
!
! 履歴  2019/03/25 竹広真一
!       2019/10/05 佐々木洋平
!       2020/08/12 竹広真一
!
!-------------------------------------------------------------------------
!> @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_mpi_module の下部モジュール: 基本変換
!>
!> @warning
!> 本モジュールを直接呼ぶ事は想定されていないことに注意されたい．
!> ユーザは w_mpi_module を呼ぶこと．
!>
module w_mpi_module_base
  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
  use mpi
  implicit none
  private

  public im, jm, nn, mm, nm                   ! 格子点数, 切断波数
  public it, t, p, r, jc                      ! 変換用作業配列
  public j1, j2                               ! 分散格子データの大きさ
  public jl                                   ! 分散スペクトルデータの大きさ
  public nc
  public np, ip
  public icom                                 ! コミュニケーター

  public w_base_mpi_Initial                   ! 初期化サブルーチン
  public w_base_mpi_Finalize                  ! 終了サブルーチン

  public x_Lon                                ! 格子座標
  public x_Lon_Weight                         ! 格子座標重み
  public v_Lat                                ! 格子座標
  public v_Lat_Weight                         ! 格子座標重み
  public y_Lat                                ! 格子座標
  public y_Lat_Weight                         ! 格子座標重み
  public xv_Lon, xv_Lat                       ! 格子座標(im,jl)
  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

  !> 格子点の設定(東西)
  integer(8)              :: im=64
  !> 格子点の設定(南北)
  integer(8)              :: jm=32
  !> 切断波数(東西波数)の設定
  integer(8)              :: mm
  !> 切断波数(全波数)の設定
  integer(8)              :: nn
  !> 計算する最大の全波数の設定
  integer(8)              :: nm
  !> 経度座標格子点
  real(DP),   allocatable :: x_Lon(:)
  !> 経度座標格子点重み
  real(DP),   allocatable :: x_Lon_Weight(:)
  !> 緯度座標分散格子点
  real(DP),   allocatable :: v_Lat(:)
  !> 緯度座標分散格子点重み
  real(DP),   allocatable :: v_Lat_Weight(:)
  !> 緯度座標格子点
  real(DP),   allocatable :: y_Lat(:)
  !> 緯度座標格子点重み
  real(DP),   allocatable :: y_Lat_Weight(:)
  !> 分散座標経度格子点, 2次元配列
  real(DP),   allocatable :: xv_Lon(:,:)
  !> 分散座標緯度格子点, 2次元配列
  real(DP),   allocatable :: xv_Lat(:,:)
  !> 座標経度格子点, 2次元配列
  real(DP),   allocatable :: xy_Lon(:,:)
  !> 座標緯度格子点, 2次元配列
  real(DP),   allocatable :: xy_Lat(:,:)
  !> コミュニケータ
  integer(8)              :: icom
  !> 変換用配列
  integer(8), allocatable :: it(:)
  !> 変換用配列
  integer(8), allocatable :: jc(:)
  !> 変換用配列
  real(DP),   allocatable :: t(:)
  !> 変換用配列
  real(DP),   allocatable :: r(:)
  !> 変換用配列
  real(DP),   pointer     :: p(:,:)
  !> 分散配置用変数
  integer(8)              :: j1
  !> 分散配置用変数
  integer(8)              :: j2
  !> 分散格子データの大きさ
  integer                 :: jl
  !> 分散スペクトルデータの大きさ
  integer                 :: nc
  !> MPI プロセス数
  integer                 :: np
  !> MPI プロセス番号
  integer                 :: ip

  !> グリッド種類
  integer(8)              :: ig=1
  !> 変換用配列
  real(DP),   pointer     :: g(:,:)
  !> 変換用配列
  real(DP),   pointer     :: w(:)
  !> 変換用配列
  real(DP),   pointer     :: g1(:,:)
  !> 変換用配列
  real(DP),   pointer     :: g2(:,:)
  !> 変換用配列
  real(DP),   pointer     :: ww(:)
  !> 変換用配列
  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(:)
  !> 分散配置用変数
  integer,    allocatable :: j1p(:)
  !> 分散配置用変数
  integer,    allocatable :: j2p(:)
  !> 分散格子データの大きさ
  integer,    allocatable :: jlp(:)
  !> 分散スペクトルデータの大きさ
  integer,    allocatable :: ncp(:)
  !> 分散/全データ変換用配列
  integer,    allocatable :: snd_rank(:)
  !> 分散/全データ変換用配列
  integer,    allocatable :: rcv_rank(:)
  !> ベクトル長
  integer(8)              :: jv
  !> 初期化フラグ
  logical                 :: w_base_mpi_initialize=.false.

  !---------------------------------------------------------------------
  !> 波数 (n,m) よりスペクトルデータの格納位置 l を返す
  !>
  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

  !---------------------------------------------------------------------
  !> スペクトルデータの位置(l) より，対応する波数 (n,m) を返す
  !>
  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

  save im, jm, nm, mm, nn                     ! 格子点数, 切断波数を記憶
  save it, t, p, r, jc                        ! 変換用配列を記憶
  save icom
  save jv
  save j1, j2
  save jl, nc
  save j1p, j2p
  save jlp, ncp
  save np, ip
  save snd_rank, rcv_rank
  save c                                      ! 変換用配列
  save w, g, pw, pg
  save ww, g1, g2, pww, pg1, pg2
  save w_base_mpi_initialize                  ! 初期化フラグ

contains
  !---------------------------------------------------------------------
  !> スペクトル変換の格子点数, 波数および OPENMP 使用時の
  !> 最大スレッド数を設定する.
  !>
  !> 実際の使用には上位サブルーチン w_mpi_initial() を用いること.
  !>
  subroutine w_base_mpi_Initial(n_in,i_in,j_in,mpi_comm)
    !> 格子点数(東西), 2の巾乗(<=2048)
    integer,intent(in) :: i_in
    !> 格子点数(南北), 4 の倍数
    integer,intent(in) :: j_in
    !> 切断全波数
    integer,intent(in) :: n_in
    !> コミュニケーター
    integer,intent(in), optional :: mpi_comm

    integer :: i, j
    integer :: ierr

    ! 三角波数切断
    im = i_in
    jm = j_in
    nn = n_in
    nm = n_in+1
    mm = n_in
    if ( present(mpi_comm) ) then
      icom = mpi_comm
    else
      icom = MPI_COMM_WORLD
    endif

    CALL MPI_COMM_SIZE(int(ICOM),NP,IERR)
    CALL MPI_COMM_RANK(int(ICOM),IP,IERR)

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

    ! 変換用作業配列
    allocate(c((mm/np+1)*(2*(nn+1)-mm/np*np)))

    call syini1(mm,nm,im,it,t,r,icom)
    call syini2(mm,nm,jm,ig,p,r,jc,icom)
    call syinic(mm,nn,c,icom)

    ! 変換用変数割付
    call syqrjv(jm,jv)

    call mxallc(pg,im*((jm/jv-1)/np+1)*jv)
    call c_f_pointer(pg, g, [im,((jm/jv-1)/np+1)*jv])
    call mxallc(pw,2*jv*((jm/jv-1)/np+1)*(mm/np+1)*np*2)
    call c_f_pointer(pw, w, [2*jv*((jm/jv-1)/np+1)*(mm/np+1)*np*2])

    call mxallc(pg1,im*((jm/jv-1)/np+1)*jv)
    call c_f_pointer(pg1, g1, [im,((jm/jv-1)/np+1)*jv])
    call mxallc(pg2,im*((jm/jv-1)/np+1)*jv)
    call c_f_pointer(pg2, g2, [im,((jm/jv-1)/np+1)*jv])
    call mxallc(pww,2*jv*((jm/jv-1)/np+1)*(mm/np+1)*np*2*2)
    call c_f_pointer(pww, ww, [2*jv*((jm/jv-1)/np+1)*(mm/np+1)*np*2*2])


    ! 各プロセスの担当緯度, 波数情報
    call syqrnj(jm,jv,j1,j2,icom)
    ! 分散格子データの大きさ
    jl = int(j2-j1+1)
    ! 分散スペクトルデータの大きさ
    nc = int((mm/np+1)*(2*(nn+1)-mm/np*np))

    allocate(j1p(0:np-1),j2p(0:np-1))
    allocate(jlp(0:np-1))
    allocate(ncp(0:np-1))

    call MPI_ALLGATHER( jl, 1, MPI_INTEGER, &
      jlp(0), 1, MPI_INTEGER, int(icom), ierr)
    call MPI_ALLGATHER( int(j1), 1, MPI_INTEGER, &
      j1p(0), 1, MPI_INTEGER, int(icom), ierr)
    call MPI_ALLGATHER( int(j2), 1, MPI_INTEGER, &
      j2p(0), 1, MPI_INTEGER, int(icom), ierr)
!!$    call MPI_ALLGATHER( j1, 1, MPI_INTEGER, &
!!$      j1p(0), 1, MPI_INTEGER, int(icom), ierr)
!!$    call MPI_ALLGATHER( j2, 1, MPI_INTEGER, &
!!$      j2p(0), 1, MPI_INTEGER, int(icom), ierr)
    call MPI_ALLGATHER( nc, 1, MPI_INTEGER, &
      ncp(0), 1, MPI_INTEGER, int(icom), ierr)

    allocate(snd_rank(1:np-1))
    allocate(rcv_rank(1:np-1))

    do i=1,np-1
      rcv_rank(i) = mod(ip+i,np)
      snd_rank(i) = mod(ip+np-i,np)
    enddo

    ! 座標変数設定
    allocate(x_Lon(0:im-1))                       ! 格子点座標格納配列(経度)
    allocate(x_Lon_Weight(0:im-1))
    allocate(v_Lat(jl),v_Lat_Weight(jl))          ! 格子点座標格納配列
    allocate(y_Lat(jm),y_Lat_Weight(jm))          ! 格子点座標格納配列
    allocate(xv_Lon(0:im-1,jl),xv_Lat(0:im-1,jl)) ! 格子点座標格納配列
    allocate(xy_Lon(0:im-1,jm),xy_Lat(0:im-1,jm)) ! 格子点座標格納配列

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

    do j=int(j1),int(j2)
      if ( j .gt. jm/2 ) then
        v_Lat(j-j1+1)        =  asin(p(j-jm/2,1))    ! 緯度座標
        v_Lat_Weight(j-j1+1) = 2*p(j -jm/2,2)        ! 緯度重み(Gauss grid)
      else
        v_Lat(j-j1+1)        = - asin(p(jm/2-j+1,1)) ! 緯度座標
        v_Lat_Weight(j-j1+1) = 2*p(jm/2-j+1,2)       ! 緯度重み(Gauss grid)
      end if
    end do

    do i=0,int(im)-1
      xv_Lat(i,:) = v_Lat
    enddo
    xy_Lat = ay_av(xv_Lat)
    y_Lat = xy_Lat(0,:)
    y_Lat_Weight = reshape(ay_av(reshape(v_Lat_Weight, (/1,jl/))),(/jm/))

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

    w_base_mpi_initialize = .true.

    call MessageNotify('M','w_base_mpi_initial', &
      'w_base_mpi_module (2019/10/05) is initialized')

  end subroutine w_base_mpi_Initial

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

    integer, parameter  :: ipow_default  = 0
    integer, parameter  :: iflag_default = 0

    integer(8) :: ipval, ifval

    ! 作業用スペクトルデータ(SYCS2X 出力用)
    real(DP)             :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np))
    ! 作業用スペクトルデータ(SYCS2Y 出力用)
    real(DP)             :: w_Ydata((mm/np+1)*(2*(nm+1)-mm/np*np))

    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('E','xv_w',&
        'w_base_mpi_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 syts2g(mm,nm,nn,im,jm,jv,w_data,g,it,t,p,r,jc,w,ipval,icom)
      xv_w = g(:,1:jl)
    else if( ifval == -1 ) then
      call sycs2x(mm,nn,w_data,w_Xdata,icom)
      call syts2g(mm,nm,nn,im,jm,jv,w_Xdata,g,it,t,p,r,jc,w,ipval,icom)
      xv_w = g(:,1:jl)
    else if( ifval == 1 ) then
      call sycs2y(mm,nn,w_data,w_Ydata,c,icom)
      call syts2g(mm,nm,nm,im,jm,jv,w_Ydata,g,it,t,p,r,jc,w,ipval,icom)
      xv_w = g(:,1:jl)
    else if( ifval == 2 ) then
      call syts2g(mm,nm,nn,im,jm,jv,w_data,g,it,t,p,r,jc,w,ipval,icom)
      xv_w = g(:,1:jl)*sin(xv_Lat)
    else
      call MessageNotify('E','xv_w','invalid value of iflag')
    endif

  end function xv_w

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


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

    integer(8) :: ipval, ifval

    ! 作業用スペクトルデータ(SYCS2X 出力用)
    real(DP)             :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np))
    ! 作業用スペクトルデータ(SYCS2Y 出力用)
    real(DP)             :: w_Ydata((mm/np+1)*(2*(nm+1)-mm/np*np))

    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('E','xv_w',&
        'w_base_mpi_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 = 0.0D0 ; g(:,1:jl) = xv_data
      call sytg2s(mm,nm,nn,im,jm,jv,w_xv,g,it,t,p,r,jc,w,ipval,icom)
    else if ( ifval == -1 ) then
      g = 0.0D0 ; g(:,1:jl) = xv_data
      call sytg2s(mm,nm,nn,im,jm,jv,w_Xdata,g,it,t,p,r,jc,w,ipval,icom)
      call sycs2x(mm,nn,w_Xdata,w_xv,icom)
    else if ( ifval == 1 ) then
      g = 0.0D0 ; g(:,1:jl) = xv_data
      call sytg2s(mm,nm,nm,im,jm,jv,w_Ydata,g,it,t,p,r,jc,w,ipval,icom)
      call sycy2s(mm,nn,w_Ydata,w_xv,c,icom)
    else if ( ifval == 2 ) then
      g = 0.0D0 ; g(:,1:jl) = xv_data*sin(xv_Lat)
      call sytg2s(mm,nm,nn,im,jm,jv,w_xv,g,it,t,p,r,jc,w,ipval,icom)
    end if

  end function w_xv


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

    real(DP)             :: w_Rdata1((mm/np+1)*(2*(nm+1)-mm/np*np))
    real(DP)             :: w_Rdata2((mm/np+1)*(2*(nm+1)-mm/np*np))
    ! 作業用スペクトルデータ(SYTS2G 出力用)
    real(DP)             :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np))
    ! 作業用スペクトルデータ(SYCS2X 出力用)
    real(DP)             :: w_Ydata((mm/np+1)*(2*(nm+1)-mm/np*np))
    ! 作業用スペクトルデータ(SYCS2Y 出力用)

    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('E','w_StreamPotential2VectorMPI',&
        'w_base_module not initialize yet.')
    endif
    !
    !   u cosφ = ∂χ/∂λ - cosφ∂ψ/∂φ
    !
    call sycs2x(mm,nn,w_Chi,w_Xdata,icom)
    call sycs2y(mm,nn,w_Psi,w_Ydata,c,icom)
    call sycrpk(mm,nn,nm,w_Xdata,w_Rdata1,icom)
    w_Rdata1 = w_Rdata1 - w_Ydata
    !
    !   v cosφ = cosφ∂χ/∂φ + ∂ψ/∂λ
    !
    call sycs2y(mm,nn,w_Chi,w_Ydata,c,icom)
    call sycs2x(mm,nn,w_Psi,w_Xdata,icom)
    call sycrpk(mm,nn,nm,w_Xdata,w_Rdata2,icom)
    w_Rdata2= w_Rdata2 + w_Ydata
    !
    call syts2v(mm,nm,nm,im,jm,jv, &
      & w_Rdata1,w_Rdata2,g1,g2,it,t,p,r,jc,ww,1_8,icom)
    ! u
    xv_U = g1(:,1:jl)
    ! v
    xv_V = g2(:,1:jl)

  end subroutine w_StreamPotential2VectorMPI

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

    real(DP)             :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np))
    ! 作業用スペクトルデータ(SUCS2X 出力用)
    real(DP)             :: w_Ydata1((mm/np+1)*(2*(nm+1)-mm/np*np))
    real(DP)             :: w_Ydata2((mm/np+1)*(2*(nm+1)-mm/np*np))
    ! 作業用スペクトルデータ(SUCS2Y 出力用)
    real(DP)             :: w_Data1((mm/np+1)*(2*(nn+1)-mm/np*np))
    real(DP)             :: w_Data2((mm/np+1)*(2*(nn+1)-mm/np*np))

    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('E','w_Vector2VorDivMPI',&
        'w_base_module not initialize yet.')
    endif
    !
    ! 1/cosφ∂u/∂λ, 1/cosφ ∂(u cosφ)/∂φ
    !
    g1 = 0.0D0 ; g1(:,1:jl) = xv_U
    g2 = 0.0D0 ; g2(:,1:jl) = xv_V
    call sytv2s(mm,nm,nn+1,im,jm,jv,&
      & w_Ydata1,w_Ydata2,g1,g2,it,t,p,r,jc,ww,1_8,icom)

    call sycrpk(mm,nm,nn,w_Ydata1,w_Xdata,icom)
    call sycs2x(mm,nn,w_Xdata,w_Div,icom)
    call sycy2s(mm,nn,w_Ydata2,w_Data2,c,icom)
    !
    ! 1/cosφ∂v/∂λ, 1/cosφ ∂(v cosφ)/∂φ
    !
    call sycrpk(mm,nm,nn,w_Ydata2,w_Xdata,icom)
    call sycs2x(mm,nn,w_Xdata,w_Vor,icom)
    call sycy2s(mm,nn,w_Ydata1,w_Data1,c,icom)
    !
    !  渦度・発散の計算
    !
    !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ
    !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
    !
    w_Vor = w_Vor - w_Data1
    w_Div = w_Div + w_Data2

  end subroutine w_Vector2VorDivMPI

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

    real(DP)             :: w_Xdata((mm/np+1)*(2*(nn+1)-mm/np*np))
    ! 作業用スペクトルデータ(SUCS2X 出力用)
    real(DP)             :: w_Ydata1((mm/np+1)*(2*(nm+1)-mm/np*np))
    real(DP)             :: w_Ydata2((mm/np+1)*(2*(nm+1)-mm/np*np))
    ! 作業用スペクトルデータ(SUCS2Y 出力用)
    real(DP)             :: w_Data1((mm/np+1)*(2*(nn+1)-mm/np*np))
    real(DP)             :: w_Data2((mm/np+1)*(2*(nn+1)-mm/np*np))

    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('E','w_VectorCosLat2VorDivMPI',&
        'w_base_module not initialize yet.')
    endif
    !
    ! 1/cosφ∂u/∂λ, 1/cosφ ∂(u cosφ)/∂φ の計算
    !
    g1 = 0.0D0 ; g1(:,1:jl) = xv_UCosLat
    g2 = 0.0D0 ; g2(:,1:jl) = xv_VCosLat
    call sytv2s(mm,nm,nn+1,im,jm,jv,w_Ydata1,w_Ydata2,g1,g2,it,t,p,r,jc,ww,2_8,icom)

    call sycrpk(mm,nm,nn,w_Ydata1,w_Xdata,icom)
    call sycs2x(mm,nn,w_Xdata,w_Div,icom)
    call sycy2s(mm,nn,w_Ydata2,w_Data2,c,icom)
    !
    ! 1/cosφ∂v/∂λ, 1/cosφ ∂(v cosφ)/∂φ の計算
    !
    call sycrpk(mm,nm,nn,w_Ydata2,w_Xdata,icom)
    call sycs2x(mm,nn,w_Xdata,w_Vor,icom)
    call sycy2s(mm,nn,w_Ydata1,w_Data1,c,icom)
    !
    !  渦度・発散の計算
    !
    !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ
    !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
    !
    w_Vor = w_Vor - w_Data1
    w_Div = w_Div + w_Data2

  end subroutine w_VectorCosLat2VorDivMPI

  !---------------------------------------------------------------------
  !> 分割緯度データを集める.
  !>
  function ay_av(av)
    !> 緯度格子点データ(分割) (:,jl)
    real(DP), dimension(:,:), intent(in) :: av
    !> 緯度格子点データ
    real(DP), dimension(size(av,1),jm)   :: ay_av

    integer :: a_size
    integer :: IRQSTS
    integer :: IST( MPI_STATUS_SIZE )

    integer :: ierr, i
    integer, parameter :: itag=103

    a_size = size(av,1)

    call MPI_Barrier(int(icom),IERR)
    do i=1,np-1
      call MPI_ISEND(av, a_size*jl, MPI_REAL8, rcv_rank(i),&
        itag,int(icom),IRQSTS,IERR)
      if ( j1p(snd_rank(i)) .ge. 1 ) then
        call MPI_RECV(ay_av(1,j1p(snd_rank(i))),&
          a_size*jlp(snd_rank(i)), MPI_REAL8,  &
          snd_rank(i), itag, int(icom), IST, IERR)
      endif
      CALL MPI_WAIT(IRQSTS,IST,IERR)
    enddo

    ay_av(:,j1:j2) = av

  end function ay_av

  !---------------------------------------------------------------------
  !> 緯度データを分割する
  !>
  function av_ay(ay)
    !> 緯度格子点データ(分割) (:,jc)
    real(DP), dimension(:,:), intent(in) :: ay
    !> 緯度格子点データ
    real(DP), dimension(size(ay,1),jl) :: av_ay

    av_ay = ay(:,j1:j2)

  end function av_ay


  !---------------------------------------------------------------------
  !> モジュールの終了処理(割り付け配列の解放)をおこなう.
  !>
  !> 実際の使用には上位サブルーチン w_mpi_finalize() を用いること.
  !>
  subroutine w_base_mpi_Finalize
    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('W','w_base_mpi_Finalize',&
        'w_base_mpi_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(j1p,j2p)
    deallocate(jlp)
    deallocate(ncp)
    deallocate(snd_rank)
    deallocate(rcv_rank)
    deallocate(x_Lon)              ! 格子点座標格納配列(経度)
    deallocate(x_Lon_weight)
    deallocate(xv_Lon)
    deallocate(v_Lat)
    deallocate(v_Lat_Weight)       ! 格子点座標格納配列
    deallocate(xv_Lat)             ! 格子点座標格納配列
    deallocate(y_Lat,y_Lat_Weight) ! 格子点座標格納配列
    deallocate(xy_Lon,xy_Lat)      ! 格子点座標格納配列

    w_base_mpi_initialize = .false.

    call MessageNotify('M','w_base_mpi_Finalize',&
      'w_base_mpi_module (2019/10/05) is finalized')

  end subroutine w_base_mpi_Finalize

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

    integer(8) :: ip_w
    integer(8) :: l, n, m

    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('E','l_nm_array00',&
        'w_base_mpi_module not initialize yet. Use synm2l routine in ISPACK3 directly.')
    end if
    n = n_in
    m = m_in
    call synm2l(nn,n,m,ip_w,l,icom)
    if ( ip_w /= ip ) then
      l_nm_array00 = 0
    else
      l_nm_array00 = int(l)
    endif
  end function l_nm_array00

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

    integer(8) :: ip_w
    integer(8) :: l, n, m

    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('E','l_nm_array00',&
        'w_base_mpi_module not initialize yet. Use synm2l routine in ISPACK3 directly.')
    end if
    n = n_in
    m = m_in
    call synm2l(nn,n,m,ip_w,l,icom)
    if ( ip_w /= ip ) then
      l_nm_darray00 = 0
    else
      l_nm_darray00 = int(l)
    endif
  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, intent(in)   :: l
    !> 全波数, 帯状波数
    integer               :: nm_l_int(2)

    integer(8)            :: n,m          ! 全波数, 帯状波数

    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('E','nm_l_int', &
        'w_base_mpi_module not initialize yet. Use syl2nm routine in ISPACK3 directly.')
    endif
    call syl2nm(mm,nn,int(l,kind=8),n,m,icom)
    nm_l_int = (/int(n),int(m)/)

  end function nm_l_int

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

    integer(8)            :: n,m          ! 全波数, 帯状波数

    if ( .not. w_base_mpi_initialize ) then
      call MessageNotify('E','nm_l_int', &
        'w_base_mpi_module not initialize yet. Use syl2nm routine in ISPACK3 directly.')
    endif
    call syl2nm(mm,nn,l,n,m,icom)
    nm_l_dint = (/int(n),int(m)/)

  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_mpi_module_base
