!--
!----------------------------------------------------------------------
!     Copyright (c) 2019 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 を呼んで
!  切断波数, 格子点数の設定をしておく必要がある. 
!
!
!履歴  2016/08/15  竹広真一  wa_base_module_sypack より多層 MPI 対応改造
!
!++
module ua_base_mpi_module_sypack
  !
  != ua_base_mpi_module_sypack
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  !  spml/ua_base_mpi_module_sypack モジュールは球面上での流体運動を
  !  球面調和函数を用いたスペクトル法によって数値計算するための 
  !  モジュール ua_mpi_module_sypack の下部モジュールであり, 
  !  スペクトル計算の基本的な Fortran90 関数を提供する. 
  !
  !  球面上の 1 層モデル用 w_base_module_sypack モジュールを多層モデル用に
  !  並列化拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
  !  対する変換が行える.
  !
  !  内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる. 
  !  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  !  ついては ISPACK3 の SYPACK のマニュアルを参照されたい.
  !
  !  このモジュールを使うためには前もって w_mpi_base_initial を呼んで
  !  切断波数, 格子点数の設定をしておく必要がある. 
  !
  use dc_message
  use w_base_mpi_module_sypack,  only : im, jm, jc=>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, nc, l_nm, nm_l
  use mpi
  
  implicit none

  integer              :: km=16         ! 同時に処理する最大データ数(層の数).
  
  integer              :: nproc, nproc_h, nproc_v
  integer              :: iproc, iproc_h, iproc_v

  integer              :: ks, ke, kc
  integer, allocatable :: ksp(:), kep(:), kcp(:)
  
  integer              :: nms, nme, nmc  
  integer, allocatable :: nmsp(:), nmep(:), nmcp(:)

  integer              :: MPI_COMM_W
  integer              :: MPI_COMM_H
  integer              :: MPI_COMM_V

  logical              :: first_wb_ua=.true.
  logical              :: first_ua_wb=.true.
  logical              :: first_wa_ua=.true.
  logical              :: first_aa_ab=.true.

  private

  !public km                                    ! 層数
  public ua_base_mpi_Initial                   ! 初期化サブルーチン
  public ua_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 km
  public ks, ke, kc
  public nms, nme, nmc
  
  public l_nm, nm_l
  public lu_nm, nm_lu
  public kp_kl
  
  public xv_Lon, xv_Lat
  public x_Lon, v_Lat
  public x_Lon_Weight, v_Lat_Weight
  public xvb_ua, ua_xvb                        ! 変換関数
  public wb_ua, ua_wb                          ! 変換関数
  public wa_ua, ua_wa                          ! 変換関数
  public aa_ab, ab_aa                          ! 変換関数    
  public ua_StreamPotential2VectorMPI          ! 流線ポテンシャルから速度場計算
  public ua_Vector2VorDivMPI                   ! 速度場から渦度発散を計算
  public ua_VectorCosLat2VorDivMPI             ! 速度場から渦度発散を計算

  save km                                      ! 最大データ数(層数)を記憶
  save nproc, iproc, nproc_h, iproc_h, nproc_v, iproc_v
  save kc, ks, ke
  save MPI_COMM_W, MPI_COMM_H, MPI_COMM_V
  save first_wb_ua, first_ua_wb, first_wa_ua, first_aa_ab
  
  contains
  !--------------- 初期化 -----------------
    subroutine ua_base_mpi_initial(n_in,i_in,j_in,k_in,npv,mpi_comm)
      ! 
      ! ua_base_mpi_module 初期化
      !
      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 MPI_COMM_RANK( MPIA_COMM_H, MYID_H_TMP , NERR )
!!$  CALL MPI_COMM_SIZE( MPIA_COMM_H, NPROC_H_TMP, NERR )
!!$
!!$  CALL MPI_COMM_RANK( MPIA_COMM_V, MYID_V_TMP , NERR )
!!$  CALL MPI_COMM_SIZE( MPIA_COMM_V, NPROC_V_TMP, NERR )

      
!!$      call set_layerdiv(npv,iproc_v,ks,ke)
!!$      kc = ke-ks+1                              ! 分散多層データの大きさ

      call w_base_mpi_initial(n_in,i_in,j_in,mpi_comm_h)

      allocate(ksp(0:NPROC_V-1))
      allocate(kep(0:NPROC_V-1))
      allocate(kcp(0:NPROC_V-1))
      call mpidevide(km,nproc_v,ksp,kep,kcp)
      ksp = ksp ; kep = kep
      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(nmsp(0:NPROC_V-1))
      allocate(nmep(0:NPROC_V-1))
      allocate(nmcp(0:NPROC_V-1))
      call mpidevide(nc,nproc_v,nmsp,nmep,nmcp)
      nms = nmsp(iproc_v) ; nme = nmep(iproc_v) ; nmc = nmcp(iproc_v)

      IF(iproc.EQ.0)THEN
         WRITE(6,*)'   RANK  NSTART  NEND    NC '
         WRITE(6,*)'------------------------------'
         DO IV=0,NPROC_V-1
            WRITE(6,'(4i7)') IV, NMSP(IV), NMEP(IV), NMCP(IV)
         ENDDO
         WRITE(6,*)''
      ENDIF

      call MessageNotify('M','ua_base_mpi_initial', &
          'ua_base_mpi_module_sypack (2019/04/26) is initialized')

    end subroutine ua_base_mpi_initial

!!$    subroutine set_layerdiv(ip,kk1,kk2)
!!$      integer, intent(IN)  :: ip             ! プロセス番号
!!$      
!!$      integer, intent(OUT) :: kk1            ! 担当する層開始番号
!!$      integer, intent(OUT) :: kk2            ! 担当する層終了番号
!!$
!!$      integer :: ki
!!$      
!!$      ki=(km-1)/npv+1
!!$      if(ki*ip.lt.km) then
!!$        kk1=1+ki*ip
!!$        kk2=min(km,ki*(ip+1))
!!$      else              ! 当該プロセスでは層データを担当しない場合
!!$        kk1=0
!!$        kk2=-1
!!$      end if
!!$
!!$    end subroutine set_layerdiv
!!$
!!$    subroutine layer2procinfo(k,ipk,kp)
!!$      integer, intent(IN)  :: k             ! 層番号
!!$      integer, intent(OUT) :: ipk           ! 担当するプロセス番号
!!$      integer, intent(OUT) :: kp            ! プロセス内での分割後層番号
!!$
!!$      integer :: iip, kk1, kk2
!!$
!!$      do iip=0,npv-1
!!$         call set_layerdiv(iip,kk1,kk2)
!!$         if ( kk1 .le. k .AND. k .le. kk2 ) then
!!$            ipk = iip
!!$            kp = k - kk1 + 1
!!$            goto 999
!!$         end if
!!$      end do
!!$
!!$      call MessageNotify('E','layer2porcinfo','Process not found.')
!!$      
!!$999   continue
!!$      
!!$    end subroutine layer2procinfo
    
    subroutine mpidevide(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 mpidevide

    !--------------- 基本関数 -----------------
    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 lu_nm(n,m)
      integer, intent(IN) :: n              ! 全波数
      integer, intent(IN) :: m              ! 帯状波数
      integer             :: lu_nm          ! 位置情報(分割プロセス中)

      integer :: ll                         ! 位置情報(全体)

      ll = l_nm(n,m)
      if ( nms .le. ll .AND. ll .le. nme ) then
         lu_nm = ll-nms+1
      else
         lu_nm = 0
      endif

    end function lu_nm
    
    function nm_lu(l)
      integer, intent(IN) :: l              ! 位置情報(分割プロセス中)
      integer             :: nm_lu(2)        ! 全波数, 帯状波数

      nm_lu = nm_l(l+nms-1)
      
    end function nm_lu
    
    !--------------- データ入れ換え -----------------
    function wb_ua(ua)
      !
      ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
      !
      real(8), dimension(nmc,km), intent(in)       :: ua
      !(in) 2 次元球面調和函数(水平再分割) (nc,:)
      real(8), dimension(nc,kc)                    :: wb_ua
      !(out) 2 次元球面調和函数チェビシェフスペクトルデータ

      real(8), allocatable :: ua_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(:)
      integer, allocatable ::rcv_idx(:)     

      save snd_rank_v, rcv_rank_v, rcv_idx

      if ( first_wb_ua ) then
         first_wb_ua = .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)+nmcp(iv-1)
         enddo
      endif

      call MPI_Barrier(MPI_COMM_V,IERR)
      do iv=1,nproc_v-1
         call MPI_ISEND(ua(1,ksp(rcv_rank_v(iv))),nmc*kcp(rcv_rank_v(iv)), &
              MPI_REAL8, rcv_rank_v(iv),&
              itag,MPI_COMM_V,IRQSTS,IERR)
         allocate(ua_recv(nmcp(snd_rank_v(iv)),kc))
         call MPI_RECV(ua_recv,kc*nmcp(snd_rank_v(iv)), MPI_REAL8,&
              snd_rank_v(iv), itag, MPI_COMM_V, IST, IERR)
         wb_ua(nmsp(snd_rank_v(iv)):nmep(snd_rank_v(iv)),:) = ua_recv
         deallocate(ua_recv)
         CALL MPI_WAIT(IRQSTS,IST,IERR)
      enddo

      wb_ua(nms:nme,:) = ua(:,ks:ke)

    end function wb_ua

    function ua_wb(wb)
      !
      ! 水平スペクトル・層分割データを集める.
      !
      real(8), dimension(nc,kc), intent(in) :: wb
      !(in)  動径格子点データ(分割)
      real(8), dimension(nmc,km) :: ua_wb
      !(out) 動径格子点データ

      real(8), allocatable :: wb_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_ua_wb ) then
         first_ua_wb = .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)
      do iv=1,nproc_v-1
         allocate(wb_send(nmcp(rcv_rank_v(iv)),kc))
         wb_send = wb(nmsp(rcv_rank_v(iv)):nmep(rcv_rank_v(iv)),:)
         call MPI_ISEND(wb_send, nmcp(rcv_rank_v(iv))*kc,&
              MPI_REAL8, rcv_rank_v(iv),&
              itag,MPI_COMM_V,IRQSTS,IERR)
         deallocate(wb_send)
         call MPI_RECV(ua_wb(1,rcv_idx(snd_rank_v(iv))),&
              nmc*kcp(snd_rank_v(iv)), MPI_REAL8,    &
              snd_rank_v(iv), itag, MPI_COMM_V, IST, IERR)
         CALL MPI_WAIT(IRQSTS,IST,IERR)
      enddo

      ua_wb(:,ks:ke) = wb(nms:nme,:)

    end function ua_wb

   function wa_ua(ua)
     !
     ! 水平スペクトル・動径格子点データからスペクトルデータへ(正)変換する.
     !
     real(8), dimension(:,:), intent(in)       :: ua
     !(in) 2 次元球面調和函数(水平再分割) (nc,:)
     real(8), dimension(nc,size(ua,2))   :: wa_ua
     !(out) 2 次元球面調和函数チェビシェフスペクトルデータ

     real(8), dimension(size(ua,2),nmc)  :: au
     real(8), dimension(size(ua,2),nc)   :: aw

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

     integer :: ierr, ip_v
     integer, parameter :: itag=110
     ! 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_wa_ua ) then
        first_wa_ua = .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)+nmcp(ip_v-1)
        enddo
     endif

     au = transpose(ua)
     a_size = size(ua,2)
     
     call MPI_Barrier(MPI_COMM_V,IERR)
     do ip_v=1,nproc_v-1
        call MPI_ISEND(au, a_size*nmc, MPI_REAL8, rcv_rank_v(ip_v),&
             itag,MPI_COMM_V,IRQSTS,IERR)
        call MPI_RECV(aw(1,rcv_idx(snd_rank_v(ip_v))),  &
             a_size*nmcp(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

     aw(:,rcv_idx(iproc_v):rcv_idx(iproc_v)+nmcp(iproc_v)-1) = au

     wa_ua = transpose(aw)

   end function wa_ua
   
   function ua_wa(wa)
     !
     ! 水平スペクトルデータを再分割する
     !
     real(8), dimension(:,:), intent(in) :: wa
     !(in) 2 次元球面調和函数データ (nc,:)
     real(8), dimension(nmc,size(wa,2))  :: ua_wa
     !(out) 2 次元球面調和函数データ(水平再分割)

     ua_wa = wa(nms:nme,:)

   end function ua_wa

    
!!$    function wb_ua(ua_data)
!!$      real(8), intent(IN) :: ua_data(nmc,km)
!!$      real(8)             :: wb_ua(nc,kc)
!!$
!!$      real(8) :: databuf((MM/NP+1)*(2*(NN+1)-MM/NP*NP)*NP*KM)
!!$      integer :: nb, kk1, kk2, kp, m, ipk, k, kk, ns, ipdest, l, la, ierr
!!$
!!$      NB=(MM/NP+1)*(2*(NN+1)-MM/NP*NP)
!!$ 
!!$      do ipk=0,np-1
!!$         call set_layerdiv(ipk,kk1,kk2)
!!$         kp = kk2-kk1+1
!!$
!!$         CALL MPI_GATHER(ua_data(:,kk1),NB*kp,MPI_REAL8, &
!!$                         databuf,NB*kp,MPI_REAL8,&
!!$                         IPK,MPI_COMM_WORLD,IERR)
!!$      enddo
!!$
!!$!$omp parallel do default(shared) private(m,kk,ns,ipdest,l,la)      
!!$      do k=1,kc
!!$         CALL SUCOPY(nn+1,databuf(nb*(k-1)+1),wb_ua(1,k))
!!$         DO M=1,MM
!!$            KK=M/NP
!!$            NS=KK*(2*(NN+1)-(KK-1)*NP)+1
!!$            CALL SUNM2L(NN,M,M,IPDEST,L)
!!$            CALL SVNM2L(NN,M,M,LA)
!!$            CALL SUCOPY(2*(NN+1-M),databuf(ns+nb*kc*ipdest+nb*(k-1)),wb_ua(la,k))   
!!$         enddo
!!$      END DO
!!$!$omp end parallel do      
!!$      
!!$    end function wb_ua
!!$
!!$    function ua_wb(wb_data)
!!$      real(8), intent(IN) :: wb_data((mm+1)**2,kc)
!!$      real(8)             :: ua_wb(nc,km)
!!$
!!$      real(8) :: databuf( (mm/np+1)*(2*(nn+1)-mm/np*np)*np*km )
!!$      integer :: ipk, kk1, kk2, kp, k, m, kk, NS, NB, IPDEST, L, LA, IERR
!!$
!!$      NB=(MM/NP+1)*(2*(NN+1)-MM/NP*NP)
!!$      
!!$!$omp parallel do default(shared) private(m,kk,ns,ipdest,l,la)      
!!$      do k=1,kc
!!$         CALL SUCOPY(nn+1,wb_data(1,k),databuf(nb*(k-1)+1))
!!$         DO M=1,MM
!!$            KK=M/NP
!!$            NS=KK*(2*(NN+1)-(KK-1)*NP)+1
!!$            CALL SUNM2L(NN,M,M,IPDEST,L)
!!$            CALL SVNM2L(NN,M,M,LA)
!!$            CALL SUCOPY(2*(NN+1-M),wb_data(la,k),databuf(ns+nb*kc*ipdest+nb*(k-1)))
!!$         enddo
!!$      END DO
!!$!$omp end parallel do
!!$      
!!$      do ipk=0,np-1
!!$         call set_layerdiv(ipk,kk1,kk2)
!!$         kp = kk2-kk1+1
!!$         if ( kp == 0 ) then
!!$            CALL MPI_SCATTER(databuf,NB*KP,MPI_REAL8,&
!!$                             ua_wb(1,1),NB*KP,MPI_REAL8, &
!!$                             ipk,MPI_COMM_WORLD,IERR)
!!$         else
!!$            CALL MPI_SCATTER(databuf,NB*KP,MPI_REAL8,&
!!$                             ua_wb(1,kk1),NB*KP,MPI_REAL8, &
!!$                             ipk,MPI_COMM_WORLD,IERR)
!!$         endif
!!$      enddo
!!$
!!$    end function ua_wb

    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) = 0
         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 xvb_ua(ua_data,ipow,iflag)    ! 球面調和関数スペクトル -> 格子点
      !
      ! スペクトルデータから格子データへ変換する(多層用).
      !
      real(8), intent(in)   :: ua_data(nmc,km)
      !(in) スペクトルデータ
      !
      real(8)               :: xvb_ua(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

      real(8) :: wb_data(nc,kc)
      integer :: k

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

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

      wb_data = wb_ua(ua_data)

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

    end function xvb_ua

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

      real(8)               :: ua_xvb(nmc,km)
      !(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

      real(8) :: wb_data(nc,kc)
      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_data(:,k) = w_xv(xvb_data(:,:,k),iflag=ifval,ipow=ipval)
      enddo

      ua_xvb = ua_wb(wb_data)
      
    end function ua_xvb

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

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

      real(8), intent(in)    :: ua_Chi(nmc,km)
      !(in) 速度ポテンシャル((nm+1)*(nm+1),:)

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

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

      real(8) :: wb_Psi(nc,kc)
      !(in) 流線関数((nm+1)*(nm+1),:)

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

      integer :: k

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

    end subroutine ua_StreamPotential2VectorMPI

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

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

      real(8), intent(out)   :: ua_Vor(nmc,km)
      !(out) 渦度
      real(8), intent(out)   :: ua_Div(nmc,km)
      !(out) 発散

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

      integer :: k

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

      ua_Vor = ua_wb(wb_Vor)
      ua_Div = ua_wb(wb_Div)      
      
    end subroutine ua_Vector2VorDivMPI

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

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

      real(8), intent(out)   :: ua_Vor(nmc,km)
      !(out) 流線関数
      real(8), intent(out)   :: ua_Div(nmc,km)
      !(out) 速度ポテンシャル

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

      integer :: k

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

      ua_Vor = ua_wb(wb_Vor)
      ua_Div = ua_wb(wb_Div)      

    end subroutine ua_VectorCosLat2VorDivMPI

  !--------------- 終了処理 -----------------
    subroutine ua_base_mpi_Finalize
      ! 
      !
      ! このサブルーチンを単独で用いるのでなく, 
      ! 上位サブルーチン wa_Finalize を使用すること.
      !
      deallocate(ksp)
      deallocate(kep)
      deallocate(kcp)
      deallocate(nmsp)
      deallocate(nmep)
      deallocate(nmcp)
      
      call w_base_mpi_Finalize

      first_wb_ua=.true.
      first_ua_wb=.true.
      first_wa_ua=.true.
      first_aa_ab=.true.
      
      call MessageNotify('M','ua_base_mpi_finalize',&
           'ua_base_mpi_module_sypack is finalized (2019/05/29) ')

    end subroutine ua_base_mpi_Finalize

end module ua_base_mpi_module_sypack
