!----------------------------------------------------------------------
!     Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wq_mpi_module_supack テストプログラム :: スペクトル関数のテスト
!
!履歴  2016/04/01  竹広真一   wq_mpi_module_sjpack_spectrum_test.f90 を supack 化
!
program wq_mpi_module_supack_spectrum

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

  integer,parameter  :: im=32, jm=16, km=16  ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: nm=10, lm=31         ! 切断波数の設定(水平, 動径)
  real(8),parameter  :: ro=1.5               ! 内外半径

  real(8), allocatable ::  xvr_Torvel(:,:,:)
  real(8), allocatable ::  wq_Torvel(:,:)

  real(8), allocatable ::  xvr_Polvel(:,:,:)
  real(8), allocatable ::  wq_Polvel(:,:)

  real(8), allocatable ::  xvr_Vlon(:,:,:)
  real(8), allocatable ::  xvr_Vlat(:,:,:)
  real(8), allocatable ::  xvr_Vrad(:,:,:)

  real(8), dimension(0:nm,-nm:nm,km) ::  nmr_EkTor
  real(8), dimension(0:nm,-nm:nm,km) ::  nmr_EkPol
  real(8), dimension(0:nm,km)        ::  nr_EkTor
  real(8), dimension(0:nm,km)        ::  nr_EkPol
 
  real(8), dimension(km)             ::  r_Data

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

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

 !---------------- MPI スタート ---------------------
  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 wq_mpi_Initial(im,jm,km,nm,lm,ro)

  call MessageNotify('M','wq_mpi_module_supack_spectrum_test',&
                         'wq_mpi_module_supack spectrum functions tests.') 


  allocate(xvr_Torvel(0:im-1,1:jc,km))
  allocate(wq_Torvel(nc,0:lm))
            
  allocate(xvr_Polvel(0:im-1,1:jc,km))
  allocate(wq_Polvel(nc,0:lm))
            
  allocate(xvr_Vlon(0:im-1,1:jc,km))
  allocate(xvr_Vlat(0:im-1,1:jc,km))
  allocate(xvr_Vrad(0:im-1,1:jc,km))

  ! 2Y_1^0 * r
  xvr_Torvel = 2*sqrt(3.0D0)*sin(xvr_Lat) * xvr_Rad

  ! (2Y_2^0 +Y_2^-2) * r^2
  xvr_Polvel =( 2*sqrt(5.0D0)*(3.0/2*sin(xvr_Lat)**2-1/2.0) &
                - sqrt(5.0D0/24)*3.0*cos(xvr_Lat)**2*sin(2*xvr_Lon) ) &
              * xvr_Rad**2

  wq_Torvel = wq_xvr(xvr_Torvel)
  wq_Polvel = wq_xvr(xvr_Polvel)

  !--------- Energy spectrum ---------
  nmr_EkTor=wq_VMiss
  nmr_EkPol=wq_VMiss

  do n=0,nm
     do m=-n,n
        nmr_EkTor(n,m,:) = 0.0D0
        nmr_EkPol(n,m,:) = 0.0D0
     enddo
  enddo

  nmr_EkTor(1,0,:) = 4.0D0 * r_Rad**2 * (4*pi)*r_Rad**2

  nmr_EkPol(2,0,:) = 12.0D0 *(9*r_Rad**4 + 6.0D0*r_Rad**4) * (4*pi)

  nmr_EkPol(2,-2,:)= 2*3*(1.0D0/2)**2/2 &
       *(9*r_Rad**4 + 2*3* r_Rad**4 ) * (4*pi)
  nmr_EkPol(2,2,:) = nmr_EkPol(2,-2,:)

  call AssertEqual(&
    message='nmr_ToroidalEnergySpectrum_wq',                      &
    answer = nmr_ToroidalEnergySpectrum_wq(wq_Torvel),            &
    check =  nmr_EkTor,                                           &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  call AssertEqual(&
    message='nmr_PoloidalEnergySpectrum_wq',                      &
    answer = nmr_PoloidalEnergySpectrum_wq(wq_Polvel),            &
    check =  nmr_EkPol,                                           &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  nr_EkTor = 0.0D0
  nr_EkTor(1,:) = 4.0D0 * r_Rad**2 * (4*pi)*r_Rad**2 

  nr_EkPol = 0.0D0
  nr_EkPol(2,:) = (12.0D0+2*3*(1.0D0/2)**2)*(9*r_Rad**4+2*3*r_Rad**4)*(4*pi)

  call AssertEqual(&
    message='nr_ToroidalEnergySpectrum_wq',                       &
    answer = nr_ToroidalEnergySpectrum_wq(wq_Torvel),             &
    check =  nr_EkTor,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  call AssertEqual(&
    message='nr_PoloidalEnergySpectrum_wq',                       &
    answer = nr_PoloidalEnergySpectrum_wq(wq_Polvel),             &
    check =  nr_EkPol,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call wq_Potential2VectorMPI(xvr_Vlon,xvr_Vlat,xvr_Vrad,wq_Torvel,wq_Polvel)

  nmr_EkTor = nmr_ToroidalEnergySpectrum_wq(wq_Torvel)
  nmr_EkPol = nmr_PoloidalEnergySpectrum_wq(wq_Polvel)

  r_Data = 0.0D0
  do n=0,nm
     do m=-n,n
        r_Data = r_Data + nmr_EkTor(n,m,:) + nmr_EkPol(n,m,:)
     enddo
  enddo

  call AssertEqual(&
    message='Total energy by nmr-spectrum',                        &
    answer = IntLonLatRad_xvr(xvr_VLon**2+xvr_VLat**2+xvr_VRad**2)/2.0D0, &
    check =  IntRad_r(r_Data/r_Rad**2),                           &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  nr_EkTor = nr_ToroidalEnergySpectrum_wq(wq_Torvel)
  nr_EkPol = nr_PoloidalEnergySpectrum_wq(wq_Polvel)

  call AssertEqual(&
    message='Total energy by nr-spectrum',                        &
    answer = IntLonLatRad_xvr(xvr_VLon**2+xvr_VLat**2+xvr_VRad**2)/2.0D0, &
    check =   IntRad_r(sum(nr_EkTor,1)/r_Rad**2) &
            + IntRad_r(sum(nr_EkPol,1)/r_Rad**2), &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call MessageNotify('M','wq_mpi_module_supack_spectrum_test',&
                         'wq_mpi_module_supack spectrum functions tests suceeded!') 

 !------ MPIの終了 ------

  call MPI_FINALIZE(IERR)

end program wq_mpi_module_supack_spectrum