!--
!----------------------------------------------------------------------
!     Copyright (c) 2020 SPMDOEL Development Group
!----------------------------------------------------------------------
!表題  ua_base_mpi_module_sypack
!
!  spml/ua_base_mpi_module_sypack モジュールは球面上での流体運動を
!  球面調和函数を用いたスペクトル法によって数値計算するための
!  モジュール ua_mpi_module_sypack の下部モジュールであり,
!  スペクトル計算の基本的な Fortran90 関数を提供する.
!
!  球面上の 1 層モデル用 u_base_module_sypack モジュールを多層モデル用に
!  並列化拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
!  対する変換が行える.
!
!  内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる.
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK3 の SYPACK のマニュアルを参照されたい.
!
!  このモジュールを使うためには前もって u_mpi_base_initial を呼んで
!  切断波数, 格子点数の設定をしておく必要がある.
!
!
!履歴  2020/07/31  竹広真一  ua_mpi_module よリ改造
!
!++
module wb_mpi_module_base
  !
  != wb_mpi_module_base
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  !  spml/wb_mpi_module_base モジュールは球面上での流体運動を
  !  球面調和函数を用いたスペクトル法によって数値計算するための
  !  モジュール wb_mpi_module の下部モジュールであり,
  !  スペクトル計算の基本的な Fortran90 関数を提供する.
  !
  !  球面上の 1 層モデル用 w_base_module モジュールを多層モデル用に
  !  並列化拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
  !  対する変換が行える.
  !
  !  内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる.
  !  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  !  ついては ISPACK3 の SYPACK のマニュアルを参照されたい.
  !
  !  このモジュールを使うためには前もって w_mpi_base_initial を呼んで
  !  切断波数, 格子点数の設定をしておく必要がある.
  !
  use mpi
  use dc_message, only: MessageNotify
  use w_mpi_module_base,  only : im, jm, jcc=>jl, nn, mm, &
    xv_Lon, xv_Lat, x_Lon, v_Lat, x_Lon_Weight, v_Lat_Weight, &
       xv_w, w_xv, &
       w_StreamPotential2VectorMPI, w_Vector2VorDivMPI, &
       w_VectorCosLat2VorDivMPI, &
       w_base_mpi_initial, w_base_mpi_finalize, ncc=>nc, l_nm, nm_l


  implicit none

  integer              :: km=16         ! 同時に処理する最大データ数(層の数).

  integer              :: nproc, nproc_h, nproc_v
  integer              :: iproc, iproc_h, iproc_v

  integer              :: jc
  integer              :: nc
  
  integer              :: is, ie, ic
  integer, allocatable :: isp(:), iep(:), icp(:)

  integer              :: ks, ke, kc
  integer, allocatable :: ksp(:), kep(:), kcp(:)

  real(8), allocatable :: pv_Lon(:,:), pv_Lat(:,:)
  real(8), allocatable :: p_Lon(:)
  real(8), allocatable :: p_Lon_Weight(:)
  
  integer              :: MPI_COMM_W
  integer              :: MPI_COMM_H
  integer              :: MPI_COMM_V

  logical              :: first_xvb_pva=.true.
  logical              :: first_pva_xvb=.true.
  logical              :: first_aa_ab=.true.
  logical              :: first_xa_pa=.true.

  private

  public wb_base_mpi_Initial                   ! 初期化サブルーチン
  public wb_base_mpi_Finalize                  ! 終了処理

  public MPI_COMM_W, MPI_COMM_H, MPI_COMM_V
  public nproc, nproc_h, nproc_v
  public iproc, iproc_h, iproc_v

  public jc, nc
  public is, ie, ic
  public km
  public ks, ke, kc

  public l_nm, nm_l
  public kp_kl

  public xv_Lon, xv_Lat
  public x_Lon, v_Lat
  public x_Lon_Weight, v_Lat_Weight
  public pv_Lon, pv_Lat
  public p_Lon
  public p_Lon_Weight

  public wb_pva, pva_wb                        ! 変換関数
  public wb_xvb, xvb_wb                        ! 変換関数
  public xvb_pva, pva_xvb                      ! 変換関数
  public pv_w, w_pv                            ! 変換関数
  public xv_w, w_xv                            ! 変換関数

  public xa_pa, pa_xa                          ! 変換関数
  public aa_ab, ab_aa                          ! 変換関数
  public wb_StreamPotential2VectorMPI          ! 流線ポテンシャルから速度場計算
  public wb_Vector2VorDivMPI                   ! 速度場から渦度発散を計算
  public wb_VectorCosLat2VorDivMPI             ! 速度場から渦度発散を計算

  save jc  
  save nc
  save km                                      ! 最大データ数(層数)を記憶
  save nproc, iproc, nproc_h, iproc_h, nproc_v, iproc_v
  save icp, isp, iep
  save ic, is, ie
  save kcp, ksp, kep
  save kc, ks, ke
  save pv_Lon, pv_Lat
  save p_Lon
  save p_Lon_Weight
  save MPI_COMM_W, MPI_COMM_H, MPI_COMM_V
  save first_xvb_pva, first_pva_xvb
  save first_aa_ab, first_xa_pa

  contains
  !--------------- 初期化 -----------------
    subroutine wb_base_mpi_initial(n_in,i_in,j_in,k_in,npv,mpi_comm)
      !
      ! ua_mpi_module_base 初期化
      !
      integer,intent(in) :: n_in                !(in) 切断波数
      integer,intent(in) :: i_in                !(in) 格子点数(東西)
      integer,intent(in) :: j_in                !(in) 格子点数(南北)
      integer,intent(in) :: k_in                !(in) 最大データ数(層数)
      integer,intent(in),optional :: npv        !(in) 層方向分割数
      integer,intent(in),optional :: mpi_comm   !(in) 全体コミュニケーター

      integer :: ierr, iv

      km = k_in

      if ( present(npv) ) then
         nproc_v = npv
      else
         nproc_v = 1
      endif

      if ( present(mpi_comm) ) then
         MPI_COMM_W = mpi_comm
      else
         MPI_COMM_W = MPI_COMM_WORLD
      endif

      CALL MPI_COMM_SIZE(MPI_COMM_W,NPROC,IERR)
      CALL MPI_COMM_RANK(MPI_COMM_W,IPROC,IERR)

      NPROC_H = NPROC/NPROC_V
      IPROC_V = IPROC/NPROC_H
      IPROC_H = MOD(IPROC,NPROC_H)

      CALL MPI_COMM_SPLIT(MPI_COMM_W, iproc_v, iproc_h, MPI_COMM_H, IERR )
      CALL MPI_COMM_SPLIT(MPI_COMM_W, iproc_h, iproc_v, MPI_COMM_V, IERR )

      call w_base_mpi_initial(n_in,i_in,j_in,mpi_comm_h)
      jc = jcc
      nc = ncc      

      allocate(isp(0:NPROC_V-1))
      allocate(iep(0:NPROC_V-1))
      allocate(icp(0:NPROC_V-1))
      call mpidivide(int(im),nproc_v,isp,iep,icp)
      isp = isp-1 ; iep = iep-1
      is = isp(iproc_v) ; ie = iep(iproc_v) ; ic = icp(iproc_v)

      IF(iproc.EQ.0)THEN
         WRITE(6,*)'   RANK  KSTART  KEND    IC '
         WRITE(6,*)'------------------------------'
         DO IV=0,NPROC_V-1
            WRITE(6,'(4i7)') IV, ISP(IV), IEP(IV), ICP(IV)
         ENDDO
         WRITE(6,*)''
      END IF
      
      allocate(ksp(0:NPROC_V-1))
      allocate(kep(0:NPROC_V-1))
      allocate(kcp(0:NPROC_V-1))
      call mpidivide(km,nproc_v,ksp,kep,kcp)
      ks = ksp(iproc_v) ; ke = kep(iproc_v) ; kc = kcp(iproc_v)

      IF(iproc.EQ.0)THEN
         WRITE(6,*)'   RANK  KSTART  KEND    KC '
         WRITE(6,*)'------------------------------'
         DO IV=0,NPROC_V-1
            WRITE(6,'(4i7)') IV, KSP(IV), KEP(IV), KCP(IV)
         ENDDO
         WRITE(6,*)''
      END IF

      allocate(pv_Lon(ic,jc),pv_Lat(ic,jc))
      allocate(p_Lon(ic),p_Lon_Weight(ic))

      pv_Lon = xv_Lon(is:ie,:)
      pv_Lat = xv_Lat(is:ie,:)      
      p_Lon  = x_Lon(is:ie)
      p_Lon_Weight  = x_Lon_Weight(is:ie)
      
      call MessageNotify('M','wb_base_mpi_initial', &
          'wb_mpi_module_base (2020/07/31) is initialized')

    end subroutine wb_base_mpi_initial

    subroutine mpidivide(mm,np,msp,mep,mcp)
      !
      ! 鉛直分割情報
      !
      integer, intent(IN)  :: mm              ! 全格子点数
      integer, intent(IN)  :: np              ! 分割プロセス数
      integer, intent(OUT) :: msp(0:np-1)     ! 分割開始インデックス
      integer, intent(OUT) :: mep(0:np-1)     ! 分割終了インデックス
      integer, intent(OUT) :: mcp(0:np-1)     ! 分割インデックス数

      integer :: i, mint, mmod

      mint = mm/np
      mmod = mod(mm,np)
      mcp  = mint
      do i = 0,mmod-1
         mcp(i) = mcp(i)+1
      enddo

      msp(0)=1 ; mep(0) = msp(0) + mcp(0) -1
      do i = 1,np-1
         msp(i) = msp(i-1) + mcp(i-1)
         mep(i) = msp(i) + mcp(i) - 1
      end do

    end subroutine mpidivide

    !--------------- 基本関数 -----------------
    function kp_kl(kl)
      integer, intent(IN) :: kl             ! 層番号(全体)
      integer             :: kp_kl          ! 層番号(分割プロセス中)

      if ( ks .le. kl .AND. kl .le. ke ) then
         kp_kl = kl - ks + 1
      else
         kp_kl = 0
      endif

    end function kp_kl


    !--------------- データ入れ換え -----------------
    function xvb_pva(pva)
      !
      ! 格子点データ(経度・緯度分割)から格子点データ(緯度・層分割)へ変換する.
      !
      real(8), dimension(ic,jc,km), intent(in)       :: pva
      !(in) 2 次元球面調和函数(水平再分割) (nc,:)
      real(8), dimension(0:im-1,jc,kc)               :: xvb_pva
      !(out) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), allocatable :: pvb_recv(:,:,:)

      integer :: IRQSTS
      integer :: IST( MPI_STATUS_SIZE )

      integer :: ierr, iv
      integer, parameter   :: itag=110
      integer, allocatable :: snd_rank_v(:)
      integer, allocatable :: rcv_rank_v(:)

      save snd_rank_v, rcv_rank_v !, rcv_idx

      if ( first_xvb_pva ) then
         first_xvb_pva = .false.
         if ( allocated(snd_rank_v) ) deallocate(snd_rank_v)
         if ( allocated(rcv_rank_v) ) deallocate(rcv_rank_v)
         allocate(snd_rank_v(1:nproc_v-1))
         allocate(rcv_rank_v(1:nproc_v-1))

         do iv=1,nproc_v-1
            rcv_rank_v(iv) = mod(iproc_v+iv,nproc_v)
            snd_rank_v(iv) = mod(iproc_v-iv+nproc_v,nproc_v)
         enddo

      endif

      call MPI_Barrier(MPI_COMM_V,IERR)
      
      if ( jc > 0 ) then
         do iv=1,nproc_v-1
            call MPI_ISEND(pva(1,1,ksp(rcv_rank_v(iv))),ic*jc*kcp(rcv_rank_v(iv)), &
                 MPI_REAL8, rcv_rank_v(iv),&
                 itag,MPI_COMM_V,IRQSTS,IERR)
            allocate(pvb_recv(icp(snd_rank_v(iv)),jc,kc))
            call MPI_RECV(pvb_recv,kc*jc*icp(snd_rank_v(iv)), MPI_REAL8, &
                 snd_rank_v(iv), itag, MPI_COMM_V, IST, IERR)
            xvb_pva(isp(snd_rank_v(iv)):iep(snd_rank_v(iv)),:,:) = pvb_recv
            deallocate(pvb_recv)
            CALL MPI_WAIT(IRQSTS,IST,IERR)
         enddo
      end if

      xvb_pva(is:ie,:,:) = pva(:,:,ks:ke)

    end function xvb_pva

    function pva_xvb(xvb)
      !
      ! 格子点データ(経度・緯度分割)から格子点データ(緯度・層分割)へ変換する.
      !
      real(8), dimension(0:im-1,jc,kc), intent(in)   :: xvb
      !(in) 2 次元球面調和函数(水平再分割) (nc,:)
      real(8), dimension(ic,jc,km)                   :: pva_xvb
      !(out) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), allocatable :: pvb_send(:,:,:)

      integer :: IRQSTS
      integer :: IST( MPI_STATUS_SIZE )

      integer :: ierr, iv
      integer, parameter :: itag=111
      integer, allocatable ::snd_rank_v(:)
      integer, allocatable ::rcv_rank_v(:)
      integer, allocatable ::rcv_idx(:)

      save snd_rank_v, rcv_rank_v, rcv_idx

      if ( first_pva_xvb ) then
         first_pva_xvb = .false.
         if ( allocated(snd_rank_v) ) deallocate(snd_rank_v)
         if ( allocated(rcv_rank_v) ) deallocate(rcv_rank_v)
         if ( allocated(rcv_idx) ) deallocate(rcv_idx)
         allocate(snd_rank_v(1:nproc_v-1))
         allocate(rcv_rank_v(1:nproc_v-1))
         allocate(rcv_idx(0:nproc_v-1))

         do iv=1,nproc_v-1
            rcv_rank_v(iv) = mod(iproc_v+iv,nproc_v)
            snd_rank_v(iv) = mod(iproc_v-iv+nproc_v,nproc_v)
         enddo

         rcv_idx(0) = 1
         do iv=1,nproc_v-1
            rcv_idx(iv) = rcv_idx(iv-1)+kcp(iv-1)
         enddo
      endif

      call MPI_Barrier(MPI_COMM_V,IERR)

      if ( jc > 0 ) then
         do iv=1,nproc_v-1
            allocate(pvb_send(icp(rcv_rank_v(iv)),jc,kc))
            pvb_send = xvb(isp(rcv_rank_v(iv)):iep(rcv_rank_v(iv)),:,:)
            call MPI_ISEND(pvb_send, icp(rcv_rank_v(iv))*jc*kc,&
                 MPI_REAL8, rcv_rank_v(iv),&
                 itag,MPI_COMM_V,IRQSTS,IERR)
            deallocate(pvb_send)
            call MPI_RECV(pva_xvb(1,1,rcv_idx(snd_rank_v(iv))),&
                 ic*jc*kcp(snd_rank_v(iv)), MPI_REAL8,    &
                 snd_rank_v(iv), itag, MPI_COMM_V, IST, IERR)
            CALL MPI_WAIT(IRQSTS,IST,IERR)
         enddo
      endif
      
      pva_xvb(:,:,ks:ke) = xvb(is:ie,:,:)

    end function pva_xvb

    function xa_pa(pa)
      !
      ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
      !
      real(8), dimension(:,:), intent(in)       :: pa
      !(in) 2 次元格子点データ(経度分割) (ic,:)
      real(8), dimension(0:im-1,size(pa,2))     :: xa_pa
      !(out) 2 次元格子点データ

      real(8), dimension(size(pa,2),ic)      :: ap
      real(8), dimension(size(pa,2),0:im-1)  :: ax

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

      integer :: ierr, ip_v
      integer, parameter :: itag=120
      ! integer, allocatable ::snd_rank(:)
      integer, allocatable ::snd_rank_v(:)
      ! integer, allocatable ::rcv_rank(:)
      integer, allocatable ::rcv_rank_v(:)
      integer, allocatable ::rcv_idx(:)

      save snd_rank_v, rcv_rank_v, rcv_idx !, snd_rank, rcv_rank

      if ( first_xa_pa ) then
         first_xa_pa = .false.
         if ( allocated(snd_rank_v) ) deallocate(snd_rank_v)
         if ( allocated(rcv_rank_v) ) deallocate(rcv_rank_v)
         if ( allocated(rcv_idx)    ) deallocate(rcv_idx)
         allocate(snd_rank_v(1:nproc_v-1))
         allocate(rcv_rank_v(1:nproc_v-1))
         allocate(rcv_idx(0:nproc_v-1))

         do ip_v=1,nproc_v-1
            rcv_rank_v(ip_v) = mod(iproc_v+ip_v,nproc_v)
            snd_rank_v(ip_v) = mod(iproc_v-ip_v+nproc_v,nproc_v)
         enddo

         rcv_idx(0) = 0
         do ip_v=1,nproc_v-1
            rcv_idx(ip_v) = rcv_idx(ip_v-1)+icp(ip_v-1)
         enddo
      endif

      ap = transpose(pa)
      a_size = size(pa,2)

      call MPI_Barrier(MPI_COMM_V,IERR)

      if ( a_size > 0 ) then
         do ip_v=1,nproc_v-1
            call MPI_ISEND(ap, a_size*ic, MPI_REAL8, rcv_rank_v(ip_v),&
                 itag,MPI_COMM_V,IRQSTS,IERR)
            call MPI_RECV(ax(1,rcv_idx(snd_rank_v(ip_v))),  &
                 a_size*icp(snd_rank_v(ip_v)), MPI_REAL8,   &
                 snd_rank_v(ip_v), itag, MPI_COMM_V, IST, IERR)
            CALL MPI_WAIT(IRQSTS,IST,IERR)
         enddo
      endif
      
      ax(:,is:ie) = ap

      xa_pa = transpose(ax)

    end function xa_pa

    function pa_xa(xa)
      !
      ! 水平格子データを再分割する
      !
      real(8), dimension(0:,:), intent(in) :: xa
      !(in) 2 次元格子点数データ (0:im-1,:)
      real(8), dimension(ic,size(xa,2))  :: pa_xa
      !(out) 2 次元格子点データ(経度分割)

      pa_xa = xa(is:ie,:)

    end function pa_xa

    function aa_ab(ab)
      !
      ! 水平スペクトル・動径格子点データを集める.
      !
      real(8), dimension(:,:), intent(in) :: ab
      !(in)  動径格子点データ(分割)
      real(8), dimension(size(ab,1),km)   :: aa_ab
      !(out) 動径格子点データ

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

      integer :: ierr, ip_v
      integer, parameter :: itag=111
      integer, allocatable ::snd_rank_v(:)
      integer, allocatable ::rcv_rank_v(:)
      integer, allocatable ::rcv_idx(:)

      save snd_rank_v, rcv_rank_v, rcv_idx

      if ( first_aa_ab ) then
         first_aa_ab = .false.
         if ( allocated(snd_rank_v) ) deallocate(snd_rank_v)
         if ( allocated(rcv_rank_v) ) deallocate(rcv_rank_v)
         if ( allocated(rcv_idx) )    deallocate(rcv_idx)
         allocate(snd_rank_v(1:nproc_v-1))
         allocate(rcv_rank_v(1:nproc_v-1))
         allocate(rcv_idx(0:nproc_v-1))

         do ip_v=1,nproc_v-1
            rcv_rank_v(ip_v) = mod(iproc_v+ip_v,nproc_v)
            snd_rank_v(ip_v) = mod(iproc_v-ip_v+nproc_v,nproc_v)
         enddo

         rcv_idx(0) = 1
         do ip_v=1,nproc_v-1
            rcv_idx(ip_v) = rcv_idx(ip_v-1)+kcp(ip_v-1)
         enddo
      endif

      a_size = size(ab,1)

      call MPI_Barrier(MPI_COMM_V,IERR)
      do ip_v=1,nproc_v-1
         call MPI_ISEND(ab, a_size*kc, MPI_REAL8, rcv_rank_v(ip_v),&
              itag,MPI_COMM_V,IRQSTS,IERR)
         call MPI_RECV(aa_ab(1,rcv_idx(snd_rank_v(ip_v))),&
              a_size*kcp(snd_rank_v(ip_v)), MPI_REAL8,    &
              snd_rank_v(ip_v), itag, MPI_COMM_V, IST, IERR)
         CALL MPI_WAIT(IRQSTS,IST,IERR)
      enddo

      aa_ab(:,rcv_idx(iproc_v):rcv_idx(iproc_v)+kcp(iproc_v)-1) = ab

    end function aa_ab

    function ab_aa(aa)
      !
      ! 水平スペクトル・動径格子点データを分割する.
      !
      real(8), dimension(:,:), intent(in) :: aa
      !(in)  動径格子点データ(:,0:km)
      real(8), dimension(size(aa,1),kc)    :: ab_aa
      !(out) 動径格子点データ(分割)

      ab_aa = aa(:,ks:ke)

    end function ab_aa

    !--------------- 基本変換 -----------------

    function pva_wb(wb_data,ipow,iflag)    ! 球面調和関数スペクトル -> 格子点
      !
      ! スペクトルデータから格子データへ変換する(多層用).
      !
      real(8), intent(in)   :: wb_data(nc,kc)
      !(in) スペクトルデータ
      !
      real(8)               :: pva_wb(1:ic,1:jc,km)
      !(out) 格子点データ
      !
      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 ipval, ifval

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

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

      pva_wb = pva_xvb(xvb_wb(wb_data,ipow=ipval,iflag=ifval))

    end function pva_wb

    function wb_pva(pva_data,ipow,iflag)    ! 球面調和関数スペクトル -> 格子点
      !
      ! スペクトルデータから格子データへ変換する(多層用).
      !
      real(8), intent(in)   :: pva_data(1:ic,1:jc,km)
      !(in) スペクトルデータ
      !
      real(8)               :: wb_pva(nc,kc)
      !(out) 格子点データ
      !
      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 ipval, ifval

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

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

      wb_pva = wb_xvb(xvb_pva(pva_data),ipow=ipval,iflag=ifval)

    end function wb_pva

    function w_pv(pv_data,ipow,iflag) ! 格子点 -> 球面調和関数スペクトル
      !
      ! 格子データからスペクトルデータへ(正)変換する(多層用).
      !
      real(8), intent(in)   :: pv_data(ic,jc)
      !(in) 格子点データ

      real(8)               :: w_pv(nc)
      !(out) スペクトルデータ

      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 ipval, ifval

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

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

      w_pv = w_xv(xa_pa(pv_data),ipow=ipval,iflag=ifval)

    end function w_pv

    function pv_w(w_data,ipow,iflag)    ! 球面調和関数スペクトル -> 格子点
      !
      ! スペクトルデータから格子データへ変換する(多層用).
      !
      real(8), intent(in)   :: w_data(nc)
      !(in) スペクトルデータ
      !
      real(8)               :: pv_w(1:ic,1:jc)
      !(out) 格子点データ
      !
      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 ipval, ifval

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

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

      pv_w = pa_xa(xv_w(w_data,ipow=ipval,iflag=ifval))

    end function pv_w

    function xvb_wb(wb_data,ipow,iflag)    ! 球面調和関数スペクトル -> 格子点
      !
      ! スペクトルデータから格子データへ変換する(多層用).
      !
      real(8), intent(in)   :: wb_data(nc,kc)
      !(in) スペクトルデータ
      !
      real(8)               :: xvb_wb(0:im-1,1:jc,kc)
      !(out) 格子点データ
      !
      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 ipval, ifval

      integer :: k

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

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

      do k=1,kc
         xvb_wb(:,:,k) = xv_w(wb_data(:,k),iflag=ifval,ipow=ipval)
      enddo

    end function xvb_wb

    function wb_xvb(xvb_data,ipow,iflag) ! 格子点 -> 球面調和関数スペクトル
      !
      ! 格子データからスペクトルデータへ(正)変換する(多層用).
      !
      real(8), intent(in)   :: xvb_data(0:im-1,jc,kc)
      !(in) 格子点データ

      real(8)               :: wb_xvb(nc,kc)
      !(out) スペクトルデータ

      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 ipval, ifval

      integer k

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

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

      do k=1,kc
         wb_xvb(:,k) = w_xv(xvb_data(:,:,k),iflag=ifval,ipow=ipval)
      enddo

    end function wb_xvb

    !----------- 速度, 渦度・発散, 流線・ポテンシャル計算 -------------

    subroutine wb_StreamPotential2VectorMPI(wb_Psi, wb_Chi, pva_U, pva_V)
      !
      ! 流線・ポテンシャル(スペクトルデータ)から速度場(格子データ)に
      ! (逆)変換する(多層用)
      !
      ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
      !
      !   u cosφ =      ∂χ/∂λ - cosφ∂ψ/∂φ,
      !   v cosφ = cosφ∂χ/∂φ +      ∂ψ/∂λ
      !
      real(8), intent(in)    :: wb_Psi(nc,kc)
      !(in) 流線関数

      real(8), intent(in)    :: wb_Chi(nc,kc)
      !(in) 速度ポテンシャル((nm+1)*(nm+1),:)

      real(8), intent(out)   :: pva_U(ic,jc,km)
      !(out) 速度経度成分

      real(8), intent(out)   :: pva_V(ic,jc,km)
      !(out) 速度緯度成分

      real(8) :: xvb_U(0:im-1,jc,kc)
      !(out) 速度経度成分

      real(8) :: xvb_V(0:im-1,jc,kc)
      !(out) 速度緯度成分
      integer :: k

      do k=1,kc
         call w_StreamPotential2VectorMPI &
              ( wb_Psi(:,k), wb_Chi(:,k), xvb_U(:,:,k), xvb_V(:,:,k) )
      enddo

      pva_U = pva_xvb(xvb_U)
      pva_V = pva_xvb(xvb_V)      
      
    end subroutine wb_StreamPotential2VectorMPI

    subroutine wb_Vector2VorDivMPI(pva_U, pva_V, wb_Vor, wb_Div)
      !
      ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に
      ! (正)変換する(多層用)
      !
      ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
      !
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      real(8), intent(in)   :: pva_U(ic,jc,km)
      !(in) 速度経度成分

      real(8), intent(in)   :: pva_V(ic,jc,km)
      !(in) 速度緯度成分

      real(8), intent(out)   :: wb_Vor(nc,kc)
      !(out) 渦度
      real(8), intent(out)   :: wb_Div(nc,kc)
      !(out) 発散

      real(8) :: xvb_U(0:im-1,jc,kc)
      !速度経度成分

      real(8) :: xvb_V(0:im-1,jc,kc)
      !速度緯度成分

      integer :: k

      xvb_U = xvb_pva(pva_U)
      xvb_V = xvb_pva(pva_V)

      do k=1,kc
         call w_Vector2VorDivMPI &
               (xvb_U(:,:,k), xvb_V(:,:,k), wb_Vor(:,k), wb_Div(:,k))
      enddo

    end subroutine wb_Vector2VorDivMPI

    subroutine wb_VectorCosLat2VorDivMPI(pva_UCosLat, pva_VCosLat, wb_Vor, wb_Div)
      !
      ! 速度場(格子データ)から渦度・発散(スペクトルデータ)に
      ! (正)変換する(多層用)
      !
      ! スペクトル変換を用いず微分を計算するために, 変換回数が 2 回ですむ.
      !
      !   ζ = 1/cosφ∂v/∂λ - 1/cosφ ∂(u cosφ)/∂φ
      !    D = 1/cosφ∂u/∂λ + 1/cosφ ∂(v cosφ)/∂φ
      !
      real(8), intent(in)  :: pva_UCosLat(ic,jc,km)
      !(in) 速度経度成分 * cos(lat)

      real(8), intent(in)  :: pva_VCosLat(ic,jc,km)
      !(in) 速度緯度成分 * cos(lat)

      real(8), intent(out) :: wb_Vor(nc,kc)
      !(out) 流線関数
      real(8), intent(out) :: wb_Div(nc,kc)
      !(out) 速度ポテンシャル

      real(8) :: xvb_UCosLat(0:im-1,jc,kc)
      !速度経度成分 * cos(lat)

      real(8) :: xvb_VCosLat(0:im-1,jc,kc)
      !速度緯度成分 * cos(lat)

      integer :: k

      xvb_UCosLat = xvb_pva(pva_UCosLat)
      xvb_VCosLat = xvb_pva(pva_VCosLat)      
      
      do k=1,kc
         call w_VectorCosLat2VorDivMPI &
               (xvb_UCosLat(:,:,k), xvb_VCosLat(:,:,k), wb_Vor(:,k), wb_Div(:,k))
      enddo

    end subroutine wb_VectorCosLat2VorDivMPI

    !--------------- 終了処理 -----------------
    subroutine wb_base_mpi_Finalize
      !
      !
      ! このサブルーチンを単独で用いるのでなく,
      ! 上位サブルーチン wa_Finalize を使用すること.
      !
      deallocate(isp)
      deallocate(iep)
      deallocate(icp)
      deallocate(ksp)
      deallocate(kep)
      deallocate(kcp)

      deallocate(pv_Lon,pv_Lat)
      deallocate(p_Lon,p_Lon_Weight)

      call w_base_mpi_Finalize

      first_xvb_pva=.true.
      first_pva_xvb=.true.
      first_aa_ab=.true.
      first_xa_pa=.true.

      call MessageNotify('M',     &
           & 'wb_base_mpi_finalize', &
           & 'wb_mpi_module_base is finalized (2020/07/31) ')

    end subroutine wb_base_mpi_Finalize

end module wb_mpi_module_base
