!----------------------------------------------------------------------
!     Copyright (c) 2019--2020 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ut_mpi_module テストプログラム :: スペクトル関数のテスト
!
!履歴  2019/05/29  竹広真一
!      2020/11/11  竹広真一  セクター計算オプション導入
!
program ut_mpi_module_spectrum_mint_test

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use ut_mpi_module_mint
  use mpi
  implicit none

!!$  integer,parameter  :: im=32, jm=16, km=16  ! 格子点の設定(経度, 緯度, 動径)
!!$  integer,parameter  :: mint=1               ! 経度方向対称性
  integer,parameter  :: im=16, jm=16, km=16  ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: mint=2               ! 経度方向対称性
  integer,parameter  :: nm=10, lm=16         ! 切断波数の設定(水平, 動径)
  real(8),parameter  :: ri=0.5, ro=1.5       ! 内外半径

  integer :: npv=2

  real(8), allocatable ::  xvh_Torvel(:,:,:)
  real(8), allocatable ::  ut_Torvel(:,:)

  real(8), allocatable ::  xvh_Polvel(:,:,:)
  real(8), allocatable ::  ut_Polvel(:,:)

  real(8), allocatable ::  xvh_Vlon(:,:,:)
  real(8), allocatable ::  xvh_Vlat(:,:,:)
  real(8), allocatable ::  xvh_Vrad(:,:,:)

  real(8), allocatable ::  nmz_EkTor(:,:,:)
  real(8), allocatable ::  nmz_EkPol(:,:,:)

  real(8), allocatable ::  nz_EkTor(:,:)
  real(8), allocatable ::  nz_EkPol(:,:)

  real(8), allocatable ::  z_Data(:)

  ! 判定誤差設定
  integer, parameter :: check_digits = 10
  integer, parameter :: ignore = -11

  real(8) :: PI
  integer :: n,m
  integer :: iproc, np, ierr

  !---------------- MPI スタート ---------------------
  !call get_ENV_INT('NPV',npv)

  call MPI_INIT(IERR)
  call MPI_COMM_RANK(MPI_COMM_WORLD,IPROC,IERR)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)

  PI = 4.0D0*Atan(1.0D0)

  call MessageNotify('M','ut_mpi_module_spectrum_test',&
                         'ut_mpi_module_mint spectrum functions tests.')

  call ut_mpi_Initial(im,jm,km,nm,lm,ri,ro,npv=npv,mint=mint)


  allocate(xvh_Torvel(0:im-1,1:jc,kc))
  allocate(ut_Torvel(nmc,0:lm))
  allocate(xvh_Polvel(0:im-1,1:jc,kc))
  allocate(ut_Polvel(nmc,0:lm))

  allocate(xvh_Vlon(0:im-1,1:jc,kc))
  allocate(xvh_Vlat(0:im-1,1:jc,kc))
  allocate(xvh_Vrad(0:im-1,1:jc,kc))

  allocate(nmz_EkTor(0:nm,-nm:nm,0:km))
  allocate(nmz_EkPol(0:nm,-nm:nm,0:km))

  allocate(nz_EkTor(0:nm,0:km))
  allocate(nz_EkPol(0:nm,0:km))

  allocate(z_Data(0:km))

  ut_Torvel = 0.0D0
  ut_Torvel = ut_xvh(xvh_ut(ut_Torvel))

  ! 2Y_1^0  * r
  xvh_Torvel = 2*sqrt(3.0D0)*sin(xvh_Lat) * xvh_Rad

  ! (3*Y_1^0 + 2Y_2^0 +Y_2^-2) * r^2
  xvh_Polvel =( 3*sqrt(3.0D0)*sin(xvh_Lat) &
                + 2*sqrt(5.0D0)*(3.0/2*sin(xvh_Lat)**2-1/2.0) &
                - sqrt(5.0D0/24)*3.0*cos(xvh_Lat)**2*sin(2*xvh_Lon) ) &
              * xvh_Rad**2

  ut_Torvel = ut_xvh(xvh_Torvel)
  ut_Polvel = ut_xvh(xvh_Polvel)

  !--------- Energy spectrum ---------
  nmz_EkTor=ut_VMiss
  nmz_EkPol=ut_VMiss

  do n=0,nm
     do m=-n,n
        nmz_EkTor(n,m,:) = 0.0D0
        nmz_EkPol(n,m,:) = 0.0D0
     enddo
  enddo

  nmz_EkTor(1,0,:) = 4.0D0 * z_Rad**2 * (4*pi)*z_Rad**2

  nmz_EkPol(1,0,:)  = 9.0D0 * (9*z_Rad**4 + 2.0D0*z_Rad**4) * (4*pi)
  nmz_EkPol(2,0,:)  = 12.0D0 * (9*z_Rad**4 + 6.0D0*z_Rad**4) * (4*pi)
  nmz_EkPol(2,-2,:) = 2*3*(1.0D0/2)**2/2 * (9*z_Rad**4 + 2*3* z_Rad**4 ) * (4*pi)
  nmz_EkPol(2,2,:)  = nmz_EkPol(2,-2,:)

  call AssertEqual(&
    message='nmz_ToroidalEnergySpectrum_ut',                      &
    answer = nmz_ToroidalEnergySpectrum_ut(ut_Torvel),            &
    check =  nmz_EkTor,                                           &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  call AssertEqual(&
    message='nmz_PoloidalEnergySpectrum_ut',                      &
    answer = nmz_PoloidalEnergySpectrum_ut(ut_Polvel),            &
    check =  nmz_EkPol,                                           &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  nz_EkTor = 0.0D0
  nz_EkTor(1,:) = 4.0D0  * z_Rad**2 * (4*pi)*z_Rad**2

  nz_EkPol = 0.0D0
  nz_EkPol(1,:) = (9.0D0)*(9*z_Rad**4+2*z_Rad**4)*(4*pi)
  nz_EkPol(2,:) = (12.0D0+2*3*(1.0D0/2)**2)*(9*z_Rad**4+2*3*z_Rad**4)*(4*pi)

  call AssertEqual(&
    message='nz_ToroidalEnergySpectrum_ut',                       &
    answer = nz_ToroidalEnergySpectrum_ut(ut_Torvel),             &
    check =  nz_EkTor,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  call AssertEqual(&
    message='nz_PoloidalEnergySpectrum_ut',                       &
    answer = nz_PoloidalEnergySpectrum_ut(ut_Polvel),             &
    check =  nz_EkPol,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call ut_Potential2VectorMPI(xvh_Vlon,xvh_Vlat,xvh_Vrad,ut_Torvel,ut_Polvel)

  nmz_EkTor = nmz_ToroidalEnergySpectrum_ut(ut_Torvel)
  nmz_EkPol = nmz_PoloidalEnergySpectrum_ut(ut_Polvel)

  z_Data = 0.0D0
  do n=0,nm
     do m=-n,n
        z_Data = z_Data + nmz_EkTor(n,m,:) + nmz_EkPol(n,m,:)
     enddo
  enddo

  call AssertEqual(&
    message='Total energy by nmz-spectrum',                        &
    answer = IntLonLatRad_xvh(xvh_VLon**2+xvh_VLat**2+xvh_VRad**2)/2.0D0*mint, &
    check =  IntRad_z(z_Data/z_Rad**2),                            &
    significant_digits = check_digits, ignore_digits = ignore      &
    )

  nz_EkTor = nz_ToroidalEnergySpectrum_ut(ut_Torvel)
  nz_EkPol = nz_PoloidalEnergySpectrum_ut(ut_Polvel)

  call AssertEqual(&
    message='Total energy by nz-spectrum',                        &
    answer = IntLonLatRad_xvh(xvh_VLon**2+xvh_VLat**2+xvh_VRad**2)/2.0D0*mint, &
    check =   IntRad_z(sum(nz_EkTor,1)/z_Rad**2) &
            + IntRad_z(sum(nz_EkPol,1)/z_Rad**2), &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call MessageNotify('M','ut_mpi_module_spectrum_mint_test',&
                         'ut_mpi_module_mint spectrum functions tests suceeded!')

  !------ MPIの終了 ------
  call MPI_FINALIZE(IERR)

end program ut_mpi_module_spectrum_mint_test