!----------------------------------------------------------------------
! Copyright (c) 2019-2020 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!表題  ua_integral_mpi_module_mint
!
!  spml/ua_integral_mpi_module_mint モジュールは球面上での流体運動を
!  球面調和函数と MPI 並列ライブラリを用いたスペクトル法によって
!  数値計算するための モジュール ua_mpi_module_sypack の下部モジュールであり,
!  積分・平均計算のための Fortran90 関数を提供する.
!
!  球面上の 1 層モデル用 u_integral_module_sypack モジュールを多層モデル用に
!  拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
!  対する変換が行える.
!
!  内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる.
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK3/SYPACK のマニュアルを参照されたい.
!
!
!履歴  2019/05/01  竹広真一
!      2020/08/04  竹広真一  pva 変数の積分を追加
!
module ua_mpi_module_integral_mint
  !
  ! ua_mpi_module_integral_mint
  !
  !  spml/ua_mpi_module_integral_mint モジュールは球面上での流体運動を
  !  球面調和函数と MPI 並列ライブラリを用いたスペクトル法によって
  !  数値計算するための モジュール ua_mpi_module の下部モジュールであり,
  !  積分・平均計算のための Fortran90 関数を提供する.
  !
  !  球面上の 1 層モデル用 u_integral_module モジュールを多層モデル用に
  !  拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
  !  対する変換が行える.
  !
  !  内部で ISPACK3 の SYPACK の Fortran77 サブルーチンを呼んでいる.
  !  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  !  ついては ISPACK3/SYPACK のマニュアルを参照されたい.
  !
  !
  use mpi
  use w_mpi_module_base_mint, only : im
  use w_mpi_module_mint, only : IntLat_v, AvrLat_v
  use ua_mpi_module_base_mint, only : MPI_COMM_H, MPI_COMM_V, &
       ic, p_Lon_Weight, x_Lon_Weight, jc, v_Lat_Weight

  implicit none
  integer :: ierr

  private

  public a_IntLonLat_pva, IntLonLat_pv        ! 緯度経度積分
  public va_IntLon_pva, a_IntLon_pa, IntLon_p ! 経度積分
  public pa_IntLat_pva, a_IntLat_va           ! 緯度積分
  public a_AvrLonLat_pva, AvrLonLat_pv        ! 緯度経度平均
  public va_AvrLon_pva, a_AvrLon_pa, AvrLon_p ! 経度平均
  public pa_AvrLat_pva, a_AvrLat_va           ! 緯度平均

  public b_IntLonLat_xvb                      ! 緯度経度積分
  public vb_IntLon_xvb, b_IntLon_xb           ! 経度積分
  public xb_IntLat_xvb, b_IntLat_vb           ! 緯度積分
  public b_AvrLonLat_xvb                      ! 緯度経度平均
  public vb_AvrLon_xvb, b_AvrLon_xb           ! 経度平均
  public xb_AvrLat_xvb, b_AvrLat_vb           ! 緯度平均

contains

  !--------------- 積分計算 (pva) -----------------
  function a_IntLonLat_pva(pva_data)
    !
    ! 2 次元緯度経度格子点データの全領域積分(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
    ! 総和を計算している.
    !
    real(8), intent(in)   :: pva_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:ic-1,jc,*)

    real(8) :: a_IntLonLat_pva(size(pva_data,3))
    !(out) 積分されたデータの並び(*)

    a_IntLonLat_pva = a_IntLon_pa(pa_IntLat_pva(pva_data))
  end function a_IntLonLat_pva

  function IntLonLat_pv(pv_data)
    !
    ! 2 次元緯度経度格子点データの全領域積分(一層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
    ! 総和を計算している.
    !
    real(8), intent(in)   :: pv_data(:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:ic-1,jc)

    real(8) :: IntLonLat_pv
    !(out) 積分されたデータの並び(*)

    IntLonLat_pv = IntLat_v(a_IntLon_pa(pv_Data))
  end function IntLonLat_pv

  function pa_IntLat_pva(pva_data)
    !
    ! 2 次元緯度経度格子点データの緯度(Y)方向積分(多層用).
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: pva_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8)             :: pa_IntLat_pva(0:ic-1,size(pva_data,3))
    !(out) 積分された 1 次元経度(X)格子点データの並び

    real(8)             :: pa_IntLatTmp(0:ic-1,size(pva_data,3))
    Integer :: j

    pa_IntLat_pva = 0
    do j=1,jc
       pa_IntLat_pva = pa_IntLat_pva + pva_data(:,j,:) * v_Lat_Weight(j)
    enddo

    pa_IntLatTmp=pa_IntLat_pva
    CALL MPI_ALLREDUCE(pa_IntLatTMP,pa_IntLat_pva,ic*size(pva_data,3),&
         MPI_REAL8, MPI_SUM,MPI_COMM_H,IERR)

  end function pa_IntLat_pva

  function va_IntLon_pva(pva_data)
    !
    ! 2 次元緯度経度格子点データの経度(X)方向積分(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: pva_data(0:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:ic-1,jc,*)

    real(8)             :: va_IntLon_pva(jc,size(pva_data,3))
    !(out) 積分された 1 次元緯度(Y)格子点データの並び

    real(8)             :: va_IntLonTmp(jc,size(pva_data,3))
    integer :: i

    va_IntLon_pva = 0
    do i=0,ic-1
       va_IntLon_pva = va_IntLon_pva + pva_data(i,:,:) * p_Lon_Weight(i)
    enddo

    va_IntLonTmp=va_IntLon_pva
    CALL MPI_ALLREDUCE(va_IntLonTMP,va_IntLon_pva,jc*size(pva_data,3),&
         MPI_REAL8, MPI_SUM, MPI_COMM_V, IERR)

  end function va_IntLon_pva

  function a_IntLon_pa(pa_data)          ! 経度積分
    !
    ! 1 次元経度(X)格子点データの X 方向積分(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: pa_data(0:,:)
    !(in)  1 次元経度(X)格子点データの並び(0:ic-1,*)
    real(8)             :: a_IntLon_pa(size(pa_data,2))
    !(out) 積分値の並び(*)
    integer :: i

    real(8)             :: a_IntLonTmp(size(pa_data,2))

    a_IntLon_pa = 0.0D0
    do i=0,ic-1
       a_IntLon_pa = a_IntLon_pa + pa_data(i,:) * p_Lon_Weight(i)
    enddo

    a_IntLonTmp = a_IntLon_pa
    CALL MPI_ALLREDUCE(a_IntLonTMP,a_IntLon_pa,size(pa_data,2),&
         MPI_REAL8, MPI_SUM, MPI_COMM_V, IERR)

  end function a_IntLon_pa

  function IntLon_p(p_data)          ! 経度積分
    !
    ! 1 次元経度(X)格子点データの X 方向積分(一層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: p_data(:)
    !(in)  1 次元経度(X)格子点データの並び(0:ic-1,*)
    real(8)             :: IntLon_p
    !(out) 積分値の並び(*)

    real(8)             :: IntLonTmp

    IntLonTmp= sum(p_data*p_Lon_Weight)

    CALL MPI_ALLREDUCE(IntLonTMP,IntLon_p,1,&
         MPI_REAL8, MPI_SUM, MPI_COMM_V, IERR)

  end function IntLon_p

  function a_IntLat_va(va_data)
    !
    ! 1 次元緯度(Y)格子点データの Y 方向積分(多層用).
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: va_data(:,:)
    !(in)  1 次元緯度(Y)格子点データの並び(jc,*)

    real(8)             :: a_IntLat_va(size(va_data,2))
    !(out) 積分値の並び(*)

    real(8)             :: a_IntLatTmp(size(va_data,2))
    integer :: j

    a_IntLat_va = 0
    do j=1,jc
       a_IntLat_va = a_IntLat_va + va_data(j,:) * v_Lat_Weight(j)
    enddo

    a_IntLatTmp=a_IntLat_va
    CALL MPI_ALLREDUCE(a_IntLatTMP,a_IntLat_va,size(va_data,2),MPI_REAL8, &
         MPI_SUM,MPI_COMM_H,IERR)

  end function a_IntLat_va

  !--------------- 平均計算 (pva) -----------------
  function a_AvrLonLat_pva(pva_data)
    !
    ! 2 次元緯度経度格子点データの全領域平均(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
    ! 総和を計算し, x_X_Weight*y_Y_Weight の総和で割ることで平均している.
    !
    real(8), intent(in)   :: pva_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8) :: a_AvrLonLat_pva(size(pva_data,3))
    !(out) 平均値の並び(*)

    a_AvrLonLat_pva = a_AvrLon_pa(pa_AvrLat_pva(pva_data))
  end function a_AvrLonLat_pva

  function AvrLonLat_pv(pv_data)
    !
    ! 2 次元緯度経度格子点データの全領域平均(一層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
    ! 総和を計算し, x_X_Weight*y_Y_Weight の総和で割ることで平均している.
    !
    real(8), intent(in)   :: pv_data(:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:ic-1,jc,*)

    real(8) :: AvrLonLat_pv
    !(out) 平均値の並び(*)

    AvrLonLat_pv = AvrLat_v(a_AvrLon_pa(pv_data))
  end function AvrLonLat_pv

  function pa_AvrLat_pva(pva_data)
    !
    ! 2 次元緯度経度格子点データの緯度(Y)方向平均(多層用).
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算し,
    ! y_Y_Weight の総和で割ることで平均している.
    !
    real(8), intent(in) :: pva_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8)             :: pa_AvrLat_pva(0:ic-1,size(pva_data,3))
    !(out) 平均された 1 次元経度(X)格子点データの並び(0:im-1,*)

    real(8)             :: Tmp
    real(8)             :: Lat_Weight_Sum

    Tmp = sum(v_Lat_weight)
    CALL MPI_ALLREDUCE(TMP,Lat_Weight_Sum,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_H,IERR)

    pa_AvrLat_pva = pa_IntLat_pva(pva_data)/Lat_Weight_Sum

  end function pa_AvrLat_pva

  function va_AvrLon_pva(pva_data)
    !
    ! 2 次元緯度経度格子点データの経度(X)方向平均(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し,
    ! x_X_Weight の総和で割ることで平均している.
    !
    real(8), intent(in) :: pva_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8)             :: va_AvrLon_pva(jc,size(pva_data,3))
    !(out) 平均された 1 次元緯度(Y)格子点の並び(jc,*)

    real(8)             :: Tmp
    real(8)             :: Lon_Weight_Sum

    Tmp = sum(p_Lon_weight)
    CALL MPI_ALLREDUCE(TMP,Lon_Weight_Sum,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_V,IERR)

    va_AvrLon_pva = va_IntLon_pva(pva_data)/Lon_Weight_Sum

  end function va_AvrLon_pva

  function a_AvrLon_pa(pa_data)          ! 経度平均
    !
    ! 1 次元(X)格子点データの経度(X)方向平均(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し,
    ! x_X_Weight の総和で割ることで平均している.
    !
    real(8), intent(in) :: pa_data(0:,:)
    !(in)  1 次元経度(X)格子点データの並び(0:im-1,*)

    real(8)             :: a_AvrLon_pa(size(pa_data,2))
    !(out) 平均値の並び(*)

    real(8)             :: Tmp
    real(8)             :: Lon_Weight_Sum

    Tmp = sum(p_Lon_weight)
    CALL MPI_ALLREDUCE(TMP,Lon_Weight_Sum,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_V,IERR)

    a_AvrLon_pa = a_IntLon_pa(pa_data)/Lon_Weight_Sum

  end function a_AvrLon_pa

  function AvrLon_p(p_data)          ! 経度積分
    !
    ! 1 次元経度(X)格子点データの X 方向積分(一層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: p_data(:)
    !(in)  1 次元経度(X)格子点データの並び(0:im-1,*)
    real(8)             :: AvrLon_p
    !(out) 積分値の並び(*)

    real(8)             :: Tmp
    real(8)             :: Lon_Weight_Sum

    Tmp = sum(p_Lon_weight)
    CALL MPI_ALLREDUCE(TMP,Lon_Weight_Sum,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_V,IERR)

    AvrLon_p = IntLon_p(p_data)/Lon_Weight_Sum

  end function AvrLon_p

  function a_AvrLat_va(va_data)
    !
    ! 1 次元(Y)格子点データの緯度(Y)方向平均(多層用).
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算し,
    ! y_Y_Weight の総和で割ることで平均している.
    !
    real(8), intent(in) :: va_data(:,:)
    !(in)  1 次元緯度格子点データの並び(jc,*)

    real(8)             :: a_AvrLat_va(size(va_data,2))
    !(out) 平均値の並び(*)

    real(8)             :: Tmp
    real(8)             :: Lat_Weight_Sum

    Tmp = sum(v_Lat_weight)
    CALL MPI_ALLREDUCE(TMP,Lat_Weight_Sum,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_H,IERR)

    a_AvrLat_va = a_IntLat_va(va_data)/Lat_Weight_Sum

  end function a_AvrLat_va

  !--------------- 積分計算 (xvb) -----------------
  function b_IntLonLat_xvb(xvb_data)
    !
    ! 2 次元緯度経度格子点データの全領域積分(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
    ! 総和を計算している.
    !
    real(8), intent(in)   :: xvb_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8) :: b_IntLonLat_xvb(size(xvb_data,3))
    !(out) 積分されたデータの並び(*)

    b_IntLonLat_xvb = b_IntLon_xb(xb_IntLat_xvb(xvb_data))
  end function b_IntLonLat_xvb

  function xb_IntLat_xvb(xvb_data)
    !
    ! 2 次元緯度経度格子点データの緯度(Y)方向積分(多層用).
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: xvb_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8)             :: xb_IntLat_xvb(0:im-1,size(xvb_data,3))
    !(out) 積分された 1 次元経度(X)格子点データの並び

    real(8)             :: xb_IntLatTmp(0:im-1,size(xvb_data,3))
    Integer :: j

    xb_IntLat_xvb = 0
    do j=1,jc
       xb_IntLat_xvb = xb_IntLat_xvb + xvb_data(:,j,:) * v_Lat_Weight(j)
    enddo

    xb_IntLatTmp=xb_IntLat_xvb
    CALL MPI_ALLREDUCE(xb_IntLatTMP,xb_IntLat_xvb,int(im)*size(xvb_data,3),&
         MPI_REAL8, MPI_SUM,MPI_COMM_H,IERR)

  end function xb_IntLat_xvb

  function vb_IntLon_xvb(xvb_data)
    !
    ! 2 次元緯度経度格子点データの経度(X)方向積分(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: xvb_data(0:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8)             :: vb_IntLon_xvb(jc,size(xvb_data,3))
    !(out) 積分された 1 次元緯度(Y)格子点データの並び

    integer :: i

    vb_IntLon_xvb = 0
    do i=0,int(im)-1
       vb_IntLon_xvb = vb_IntLon_xvb + xvb_data(i,:,:) * x_Lon_Weight(i)
    enddo

  end function vb_IntLon_xvb

  function b_IntLon_xb(xb_data)          ! 経度積分
    !
    ! 1 次元経度(X)格子点データの X 方向積分(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: xb_data(0:,:)
    !(in)  1 次元経度(X)格子点データの並び(0:im-1,*)
    real(8)             :: b_IntLon_xb(size(xb_data,2))
    !(out) 積分値の並び(*)
    integer :: i

    b_IntLon_xb = 0.0D0
    do i=0,int(im)-1
       b_IntLon_xb = b_IntLon_xb + xb_data(i,:) * x_Lon_Weight(i)
    enddo

  end function b_IntLon_xb

  function b_IntLat_vb(vb_data)
    !
    ! 1 次元緯度(Y)格子点データの Y 方向積分(多層用).
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している.
    !
    real(8), intent(in) :: vb_data(:,:)
    !(in)  1 次元緯度(Y)格子点データの並び(jc,*)

    real(8)             :: b_IntLat_vb(size(vb_data,2))
    !(out) 積分値の並び(*)

    real(8)             :: b_IntLatTmp(size(vb_data,2))
    integer :: j

    b_IntLat_vb = 0
    do j=1,jc
       b_IntLat_vb = b_IntLat_vb + vb_data(j,:) * v_Lat_Weight(j)
    enddo

    b_IntLatTmp=b_IntLat_vb
    CALL MPI_ALLREDUCE(b_IntLatTMP,b_IntLat_vb,size(vb_data,2),MPI_REAL8, &
         MPI_SUM,MPI_COMM_H,IERR)

  end function b_IntLat_vb

  !--------------- 平均計算 (xvb) -----------------
  function b_AvrLonLat_xvb(xvb_data)
    !
    ! 2 次元緯度経度格子点データの全領域平均(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight, y_Y_Weight をかけた
    ! 総和を計算し, x_X_Weight*y_Y_Weight の総和で割ることで平均している.
    !
    real(8), intent(in)   :: xvb_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8) :: b_AvrLonLat_xvb(size(xvb_data,3))
    !(out) 平均値の並び(*)

    b_AvrLonLat_xvb = b_AvrLon_xb(xb_AvrLat_xvb(xvb_data))
  end function b_AvrLonLat_xvb

  function xb_AvrLat_xvb(xvb_data)
    !
    ! 2 次元緯度経度格子点データの緯度(Y)方向平均(多層用).
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算し,
    ! y_Y_Weight の総和で割ることで平均している.
    !
    real(8), intent(in) :: xvb_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8)             :: xb_AvrLat_xvb(0:im-1,size(xvb_data,3))
    !(out) 平均された 1 次元経度(X)格子点データの並び(0:im-1,*)

    real(8)             :: Tmp
    real(8)             :: Lat_Weight_Sum

    Tmp = sum(v_Lat_weight)
    CALL MPI_ALLREDUCE(TMP,Lat_Weight_Sum,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_H,IERR)

    xb_AvrLat_xvb = xb_IntLat_xvb(xvb_data)/Lat_Weight_Sum

  end function xb_AvrLat_xvb

  function vb_AvrLon_xvb(xvb_data)
    !
    ! 2 次元緯度経度格子点データの経度(X)方向平均(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し,
    ! x_X_Weight の総和で割ることで平均している.
    !
    real(8), intent(in) :: xvb_data(:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,jc,*)

    real(8)             :: vb_AvrLon_xvb(jc,size(xvb_data,3))
    !(out) 平均された 1 次元緯度(Y)格子点の並び(jc,*)

    vb_AvrLon_xvb = vb_IntLon_xvb(xvb_data)/sum(x_Lon_Weight)

  end function vb_AvrLon_xvb

  function b_AvrLon_xb(xb_data)          ! 経度平均
    !
    ! 1 次元(X)格子点データの経度(X)方向平均(多層用).
    !
    ! 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算し,
    ! x_X_Weight の総和で割ることで平均している.
    !
    real(8), intent(in) :: xb_data(0:,:)
    !(in)  1 次元経度(X)格子点データの並び(0:im-1,*)

    real(8)             :: b_AvrLon_xb(size(xb_data,2))
    !(out) 平均値の並び(*)

    b_AvrLon_xb = b_IntLon_xb(xb_data)/sum(x_Lon_Weight)

  end function b_AvrLon_xb

  function b_AvrLat_vb(vb_data)
    !
    ! 1 次元(Y)格子点データの緯度(Y)方向平均(多層用).
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算し,
    ! y_Y_Weight の総和で割ることで平均している.
    !
    real(8), intent(in) :: vb_data(:,:)
    !(in)  1 次元緯度格子点データの並び(jc,*)

    real(8)             :: b_AvrLat_vb(size(vb_data,2))
    !(out) 平均値の並び(*)

    real(8)             :: Tmp
    real(8)             :: Lat_Weight_Sum

    Tmp = sum(v_Lat_weight)
    CALL MPI_ALLREDUCE(TMP,Lat_Weight_Sum,1,MPI_REAL8, &
         MPI_SUM,MPI_COMM_H,IERR)

    b_AvrLat_vb = b_IntLat_vb(vb_data)/Lat_Weight_Sum

  end function b_AvrLat_vb

end module ua_mpi_module_integral_mint