!--
!----------------------------------------------------------------------
!     Copyright (c) 2025 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!表題  wa_zonal_module_integral
!
!  spml/wa_zonal_module_integral モジュールは球面上での帯状的な流体運動を
!  球面調和函数を用いたスペクトル法によって数値計算するためのモジュール
!  wa_zonal_module の下部モジュールであり, 積分・平均計算のための 
!  Fortran90 関数を提供する.
!
!  球面上の 1 層モデル用 w_zonal_module_integral モジュールを多層モデル用に
!  拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
!  対する変換が行える.
!
!  内部で ISPACK3 の LXPACK の Fortran77 サブルーチンを呼んでいる.
!  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
!  ついては ISPACK3/LXPACK のマニュアルを参照されたい.
!
!
!履歴  2025/12/25  竹広真一
!
!++
module wa_zonal_module_integral
  !
  != wa_zonal_integral_module
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  !  spml/wa_zonal_integral_module モジュールは球面上での帯状的な流体運動を
  !  球面調和函数を用いたスペクトル法によって数値計算するための
  !  モジュール wa_zonal_module の下部モジュールであり,
  !  積分・平均計算のための Fortran90 関数を提供する.
  !
  !  球面上の 1 層モデル用 w_integral_module モジュールを多層モデル用に
  !  拡張したものであり, 同時に複数個のスペクトルデータ, 格子点データに
  !  対する変換が行える.
  !
  !  内部で ISPACK3 の LXPACK の Fortran77 サブルーチンを呼んでいる.
  !  スペクトルデータおよび格子点データの格納方法や変換の詳しい計算法に
  !  ついては ISPACK3/LXPACK のマニュアルを参照されたい.
  !
  !
  use dc_types, only: DP
  use w_zonal_module_base, only : im, jm, x_Lon_Weight, y_Lat_Weight
  implicit none

  private

  public a_IntLonLat_xya                      ! 緯度経度積分
  public ya_IntLon_xya, a_IntLon_xa           ! 経度積分
  public xa_IntLat_xya, a_IntLat_ya           ! 緯度積分
  public a_AvrLonLat_xya                      ! 緯度経度平均
  public ya_AvrLon_xya, a_AvrLon_xa           ! 経度平均
  public xa_AvrLat_xya, a_AvrLat_ya           ! 緯度平均

contains

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

    a_IntLonLat_xya = a_IntLon_xa(xa_IntLat_xya(xya_data))

  end function a_IntLonLat_xya

  function xa_IntLat_xya(xya_data)
    !
    ! 2 次元緯度経度格子点データの緯度(Y)方向積分(多層用).
    !
    ! 実際には格子点データ各点毎に y_Y_Weight をかけた総和を計算している.
    !
    real(DP), intent(in) :: xya_data(0:,:,:)
    !(in) 2 次元経度緯度格子点データの並び(0:im-1,1:jm,*)
    real(DP)             :: xa_IntLat_xya(0:im-1,size(xya_data,3))
    !(out) 積分された 1 次元経度(X)格子点データの並び
    integer :: j
    integer :: i_jm

    i_jm = int(jm)
    xa_IntLat_xya = 0.0D0
    do j=1,i_jm
      xa_IntLat_xya = xa_IntLat_xya + xya_data(:,j,:) * y_Lat_Weight(j)
    enddo

  end function xa_IntLat_xya

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

    real(DP)             :: ya_IntLon_xya(1:jm,size(xya_data,3))
    !(out) 積分された 1 次元緯度(Y)格子点データの並び

    integer :: i
    integer :: i_im

    i_im = int(im)

    ya_IntLon_xya = 0.0D0
    do i=0,i_im-1
      ya_IntLon_xya = ya_IntLon_xya + xya_data(i,:,:) * x_Lon_Weight(i)
    enddo

  end function ya_IntLon_xya

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

    real(DP)             :: a_IntLat_ya(size(ya_data,2))
    !(out) 積分値の並び(*)

    integer :: i_jm
    integer :: j

    i_jm = int(jm)

    a_IntLat_ya = 0.0D0
    do j=1,i_jm
      a_IntLat_ya = a_IntLat_ya + ya_data(j,:) * y_Lat_Weight(j)
    enddo

  end function a_IntLat_ya

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

    i_im = int(im)
    a_IntLon_xa = 0.0D0
    do i=0,i_im-1
      a_IntLon_xa = a_IntLon_xa + xa_data(i,:) * x_Lon_Weight(i)
    enddo

  end function a_IntLon_xa

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

    real(DP) :: a_AvrLonLat_xya(size(xya_data,3))
    !(out) 平均値の並び(*)

    a_AvrLonLat_xya = a_AvrLon_xa(xa_AvrLat_xya(xya_data))

  end function a_AvrLonLat_xya

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

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

    xa_AvrLat_xya = xa_IntLat_xya(xya_data)/sum(y_Lat_Weight)

  end function xa_AvrLat_xya

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

    real(DP)             :: ya_AvrLon_xya(1:jm,size(xya_data,3))
    !(out) 平均された 1 次元緯度(Y)格子点の並び(1:jm,*)

    ya_AvrLon_xya = ya_IntLon_xya(xya_data)/sum(x_Lon_Weight)

  end function ya_AvrLon_xya

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

    real(DP)             :: a_AvrLat_ya(size(ya_data,2))
    !(out) 平均値の並び(*)

    a_AvrLat_ya = a_IntLat_ya(ya_data)/sum(y_Lat_Weight)

  end function a_AvrLat_ya

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

    real(DP)             :: a_AvrLon_xa(size(xa_data,2))
    !(out) 平均値の並び(*)

    a_AvrLon_xa = a_IntLon_xa(xa_data)/sum(x_Lon_Weight)

  end function a_AvrLon_xa

end module wa_zonal_module_integral