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

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

  private

  public km                                    ! 層数
  public ua_base_mpi_Initial                   ! 初期化サブルーチン
  public ua_base_mpi_Finalize                  ! 終了処理

  public nc
  public kc
  public k1, k2
  
  public l_nm, nm_l
  public kp_kl
  
  public xyb_ua, ua_xyb                        ! 変換関数
  public wb_ua, ua_wb                          ! 変換関数
  public ua_StreamPotential2VectorMPI          ! 流線ポテンシャルから速度場計算
  public ua_Vector2VorDivMPI                   ! 速度場から渦度発散を計算
  public ua_VectorCosLat2VorDivMPI             ! 速度場から渦度発散を計算

  save km                                      ! 最大データ数(層数)を記憶
  save np, ip
  save kc, k1, k2
 
  contains
  !--------------- 初期化 -----------------
    subroutine ua_base_mpi_initial(n_in,i_in,j_in,k_in)
      ! 
      ! 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 :: ierr
      
      km = k_in

      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)

      call set_layerdiv(ip,k1,k2)
      kc = k2-k1+1                              ! 分散多層データの大きさ

      call w_base_mpi_initial(n_in,i_in,j_in)

      call MessageNotify('M','ua_base_mpi_initial', &
          'ua_base_mpi_module_svpack (2016/08/15) 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)/np+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,np-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
    
    !--------------- 基本関数 -----------------

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

      integer :: ipk

      call layer2procinfo(kl,ipk,kp_kl)
      if ( ip /= ipk ) kp_kl = 0

    end function kp_kl
    
    !--------------- データ入れ換え -----------------
    function wb_ua(ua_data)
      real(8), intent(IN) :: ua_data(nc,km)
      real(8)             :: wb_ua((mm+1)**2,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 xyb_ua(ua_data,ipow,iflag)    ! 球面調和関数スペクトル -> 格子点
      !
      ! スペクトルデータから格子データへ変換する(多層用).
      !
      real(8), intent(in)   :: ua_data(nc,km)
      !(in) スペクトルデータ
      !
      real(8)               :: xyb_ua(0:im-1,1:jm,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((mm+1)**2,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
        xyb_ua(:,:,k) = xy_w(wb_data(:,k),iflag=ifval,ipow=ipval)
      enddo

    end function xyb_ua

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

      real(8)               :: ua_xyb(nc,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((mm+1)**2,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_xy(xyb_data(:,:,k),iflag=ifval,ipow=ipval)
      enddo

      ua_xyb = ua_wb(wb_data)
      
    end function ua_xyb

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

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

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

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

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

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

      real(8) :: wb_Chi((mm+1)**2,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_StreamPotential2Vector &
              ( wb_Psi(:,k), wb_Chi(:,k), xyb_U(:,:,k), xyb_V(:,:,k) )
      enddo

    end subroutine ua_StreamPotential2VectorMPI

    subroutine ua_Vector2VorDivMPI(xyb_U, xyb_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)   :: xyb_U(0:im-1,jm,kc)
      !(in) 速度経度成分(0:im-1,1:jm,:)

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

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

      real(8)  :: wb_Vor((mm+1)**2,kc)  ! 渦度
      real(8)  :: wb_Div((mm+1)**2,kc)  ! 発散

      integer :: k

      do k=1,kc
         call w_Vector2VorDiv &
               (xyb_U(:,:,k), xyb_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(xyb_UCosLat, xyb_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)   :: xyb_UCosLat(0:im-1,jm,kc)
      !(in) 速度経度成分 * cos(lat) (0:im-1,1:jm,:)

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

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

      real(8)  :: wb_Vor((mm+1)**2,kc)  ! 渦度
      real(8)  :: wb_Div((mm+1)**2,kc)  ! 発散

      integer :: k

      do k=1,kc
         call w_VectorCosLat2VorDiv &
               (xyb_UCosLat(:,:,k), xyb_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 を使用すること.
      !
      call MessageNotify('M','ua_base_mpi_finalize',&
           'ua_base_mpi_module_svpack is finalized (dummy, 2016/08/15) ')

    end subroutine ua_base_mpi_Finalize

end module ua_base_mpi_module_svpack
