! -*- mode: f90; coding: utf-8 -*-
!-------------------------------------------------------------------------
! Copyright (c) 2019 SPMODEL Development Group. All rights reserved.
!-------------------------------------------------------------------------
!> @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 を呼ぶこと．
!>
!
! 履歴  2019/03/21 竹広真一
!       2019/10/09 佐々木洋平
!
module w_mpi_module_integral
  use dc_types, only: DP
  use w_mpi_module_base, only : &
    & im, jc=>jl, x_Lon_Weight, v_Lat_Weight, icom
  use mpi

  implicit none
  private

  public IntLonLat_xv ! 緯度経度積分
  public v_IntLon_xv  ! 経度積分
  public IntLon_x     ! 経度積分
  public x_IntLat_xv  ! 緯度積分
  public IntLat_v     ! 緯度積分
  public AvrLonLat_xv ! 緯度経度平均
  public v_AvrLon_xv  ! 経度平均
  public AvrLon_x     ! 経度平均
  public x_AvrLat_xv  ! 緯度平均
  public AvrLat_v     ! 緯度平均

  integer :: ierr

contains

  !---------------------------------------------------------------------
  !> 2 次元緯度経度格子点データの全領域積分(1 層用).
  !>
  !> 実際には格子点データ各点毎に x_X_Weight, v_V_Weight をかけた
  !> 総和を計算している.
  !>
  function IntLonLat_xv(xv_data)
    !> 2 次元経度緯度格子点データ(0:im-1,jc)
    real(DP), intent(in)   :: xv_data(0:im-1,jc)
    !> 積分値
    real(DP) :: IntLonLat_xv

    IntLonLat_xv = IntLon_x(x_IntLat_xv(xv_data))

  end function IntLonLat_xv

  !---------------------------------------------------------------------
  !> 2 次元緯度経度格子点データの緯度(Y)方向積分(1 層用).
  !>
  !> 実際には格子点データ各点毎に v_Y_Weight をかけた総和を計算している.
  !>
  function x_IntLat_xv(xv_data)
    !> 2 次元経度緯度格子点データ(0:im-1,jc)
    real(DP), intent(in) :: xv_data(0:im-1,jc)
    !> 積分された 1 次元経度(X)格子点データ
    real(DP)             :: x_IntLat_xv(0:im-1)

    real(DP)             :: x_IntLatTmp(0:im-1)
    integer :: j

    x_IntLat_xv = 0.0D0
    do j=1,jc
      x_IntLat_xv = x_IntLat_xv + xv_data(:,j) * v_Lat_weight(j)
    enddo

    x_IntLatTmp=x_IntLat_xv
    CALL MPI_ALLREDUCE(x_IntLatTMP, x_IntLat_xv, int(im), MPI_REAL8, &
      & MPI_SUM, int(icom), IERR)
  end function x_IntLat_xv

  !---------------------------------------------------------------------
  !> 2 次元緯度経度格子点データの経度(X)方向積分(1 層用).
  !>
  !> 実際には格子点データ各点毎に x_X_Weight をかけた総和を計算している.
  !>
  function v_IntLon_xv(xv_data)
    !> 2 次元経度緯度格子点データ(0:im-1,jc)
    real(DP), intent(in) :: xv_data(0:im-1,jc)
    !> 積分された 1 次元緯度(Y)格子点データ
    real(DP)             :: v_IntLon_xv(jc)

    integer :: i

    v_IntLon_xv = 0.0D0
    do i=0,int(im)-1
      v_IntLon_xv = v_IntLon_xv + xv_data(i,:) * x_Lon_weight(i)
    enddo

  end function v_IntLon_xv

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

    IntLon_x = sum(x_data * x_Lon_weight)

  end function IntLon_x

  !---------------------------------------------------------------------
  !> 1 次元緯度(Y)格子点データの Y 方向積分(1 層用).
  !>
  !> 実際には格子点データ各点毎に v_Y_Weight をかけた総和を計算している.
  !>
  function IntLat_v(v_data)
    !> 1 次元緯度(Y)格子点データ
    real(DP), intent(in) :: v_data(jc)
    !> 積分値
    real(DP)             :: IntLat_v

    real(DP)             :: IntLatTmp

    IntLat_v = sum(v_data * v_Lat_weight)
    IntLatTmp=IntLat_v
    CALL MPI_ALLREDUCE(IntLatTMP, IntLat_v, 1, MPI_REAL8, &
      & MPI_SUM, int(icom), IERR)

  end function IntLat_v

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

    AvrLonLat_xv = AvrLon_x(x_AvrLat_xv(xv_data))

  end function AvrLonLat_xv

  !---------------------------------------------------------------------
  !> 2 次元緯度経度格子点データの緯度(Y)方向平均(1 層用).
  !>
  !> 実際には格子点データ各点毎に v_Lat_Weight をかけた総和を計算し,
  !> y_Lat_Weight の総和で割ることで平均している.
  !>
  function x_AvrLat_xv(xv_data)
    !> 2 次元経度緯度格子点データ(0:im-1,jc)
    real(DP), intent(in) :: xv_data(0:im-1,jc)
    !> 平均された 1 次元経度(X)格子点データ
    real(DP)             :: x_AvrLat_xv(im)

    real(DP)             :: Tmp
    real(DP)             :: Lat_Weight_Sum

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

    x_AvrLat_xv = x_IntLat_xv(xv_data)/Lat_Weight_Sum

  end function x_AvrLat_xv

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

    v_AvrLon_xv = v_IntLon_xv(xv_data)/sum(x_Lon_weight)

  end function v_AvrLon_xv

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

    AvrLon_x = IntLon_x(x_data)/sum(x_Lon_weight)

  end function AvrLon_x

  !---------------------------------------------------------------------
  !> 1 次元(Y)格子点データの緯度(Y)方向平均(1 層用).
  !>
  !> 実際には格子点データ各点毎に v_Y_Weight をかけた総和を計算し,
  !> v_Y_Weight の総和で割ることで平均している.
  !>
  function AvrLat_v(v_data)
    !> 1 次元緯度格子点データ
    real(DP), intent(in) :: v_data(jc)
    !> 平均値
    real(DP)             :: AvrLat_v

    real(DP)             :: Tmp
    real(DP)             :: Lat_Weight_Sum

    Tmp = sum(v_Lat_weight)
    CALL MPI_ALLREDUCE(TMP, Lat_Weight_Sum, 1, MPI_REAL8, &
      & MPI_SUM, int(icom), IERR)
    AvrLat_v = IntLat_v(v_data)/Lat_Weight_Sum

  end function AvrLat_v

end module w_mpi_module_integral
