!----------------------------------------------------------------------
!     Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wsc_mpi_module_supack テストプログラム :: スペクトル関数のテスト
!
!履歴  2017/05/31  竹広真一 
!
program wsc_mpi_module_supack_spectrum_test

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

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

  real(8), allocatable ::  xvz_Torvel(:,:,:)
  real(8), allocatable ::  wc_Torvel(:,:)

  real(8), allocatable ::  xvz_Polvel(:,:,:)
  real(8), allocatable ::  ws_Polvel(:,:)

  real(8), allocatable ::  xvz_Vlon(:,:,:)
  real(8), allocatable ::  xvz_Vlat(:,:,:)
  real(8), allocatable ::  xvz_Vrad(:,:,:)

  real(8), allocatable ::  xvz_Xi(:,:,:)  

  real(8), dimension(0:nm,-nm:nm,0:km) ::  nmz_EkTor
  real(8), dimension(0:nm,-nm:nm,0:km) ::  nmz_EkPol
  real(8), dimension(0:nm,0:km)        ::  nz_EkTor
  real(8), dimension(0:nm,0:km)        ::  nz_EkPol

  real(8), dimension(0:km)             ::  z_Data
  real(8), dimension(0:km)             ::  z_Xi  

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

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

  PI = 4.0D0*Atan(1.0D0)

 !---------------- MPI スタート ---------------------
  call MPI_INIT(IERR)
  call MPI_COMM_RANK(MPI_COMM_WORLD,IPROC,IERR)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)

  call wsc_mpi_Initial(im,jm,km,nm,lm,ri,ro)

  call MessageNotify('M','wsc_mpi_module_supack_spectrum_test',&
                         'wsc_mpi_module_supack spectrum functions tests.') 


  allocate(xvz_Torvel(0:im-1,jc,0:km))
  allocate(wc_Torvel(nc,0:lm))

  allocate(xvz_Polvel(0:im-1,jc,0:km))
  allocate(ws_Polvel(nc,lm))

  allocate(xvz_Vlon(0:im-1,jc,0:km))
  allocate(xvz_Vlat(0:im-1,jc,0:km))
  allocate(xvz_Vrad(0:im-1,jc,0:km))

  allocate(xvz_Xi(0:im-1,jc,0:km))

  xvz_xi = (xvz_Rad-ri)/(ro-ri)
  z_xi = xvz_xi(0,1,:)
  
  ! (2Y_1^0+3Y_2^1) 
  xvz_Torvel =( 2*sqrt(3.0D0)*sin(xvz_Lat) &
              + 3*sqrt(5.0D0/6)*3.0*sin(xvz_Lat)*cos(xvz_Lat) * cos(xvz_Lon) )

  ! (3*Y_1^0 + 2Y_1^1 + 2Y_2^0 +Y_2^-2) * sin(z)
  xvz_Polvel =( 3*sqrt(3.0D0)*sin(xvz_Lat) & 
                + 2* sqrt(3.0D0/2)*cos(xvz_Lat)*cos(xvz_Lon) &
                + 2*sqrt(5.0D0)*(3.0/2*sin(xvz_Lat)**2-1/2.0) &
                - sqrt(5.0D0/24)*3.0*cos(xvz_Lat)**2*sin(2*xvz_Lon) ) &
              * sin(PI*xvz_xi)

  wc_Torvel = wc_xvz(xvz_Torvel)
  ws_Polvel = ws_xvz(xvz_Polvel)

  !--------- Energy spectrum ---------
  nmz_EkTor=wsc_VMiss
  nmz_EkPol=wsc_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 * (4*pi)*Ro**2
  nmz_EkTor(2,1,:) = 2*3*(3.0D0/2)**2/2 * (4*pi)*Ro**2 
  nmz_EkTor(2,-1,:)= nmz_EkTor(2,1,:)

  nmz_EkPol(1,0,:) = 9.0D0 * ((PI*Ro*cos(PI*z_xi))**2 + 2.0D0*sin(PI*z_xi)**2)*(4*pi)
  nmz_EkPol(1,1,:) = 1*2*(2/2)**2/2 &
       * ((PI*Ro*cos(PI*z_xi))**2 + 1*2*sin(PI*z_xi)**2) * (4*pi)
  nmz_EkPol(1,-1,:)= nmz_EkPol(1,1,:)
  nmz_EkPol(2,0,:) = 12.0D0 *((PI*Ro*cos(PI*z_xi))**2 + 6.0D0*sin(PI*z_xi)**2) * (4*pi)
  nmz_EkPol(2,-2,:)= 2*3*(1.0D0/2)**2/2 &
       *((PI*Ro*cos(PI*z_xi))**2 + 2*3*sin(PI*z_xi)**2) * (4*pi)
  nmz_EkPol(2,2,:) = nmz_EkPol(2,-2,:)

  call AssertEqual(&
    message='nmz_ToroidalEnergySpectrum_wc',                      &
    answer = nmz_ToroidalEnergySpectrum_wc(wc_Torvel),            &
    check =  nmz_EkTor,                                           &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  call AssertEqual(&
    message='nmz_PoloidalEnergySpectrum_ws',                      &
    answer = nmz_PoloidalEnergySpectrum_ws(ws_Polvel),            &
    check =  nmz_EkPol,                                           &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  nz_EkTor = 0.0D0
  nz_EkTor(1,:) = 4.0D0 * (4*pi)*Ro**2
  nz_EkTor(2,:) = 2 * 3 *(3.0D0/2)**2 * (4*pi)*Ro**2

  nz_EkPol = 0.0D0
  nz_EkPol(1,:) = (9.0D0+1*2*(2/2)**2/2*2 )*((PI*Ro*cos(PI*z_xi))**2+2*sin(PI*z_xi)**2)*(4*pi)
  nz_EkPol(2,:) = (12.0D0+2*3*(1.0D0/2)**2)*((PI*Ro*cos(PI*z_xi))**2+2*3*sin(PI*z_xi)**2)*(4*pi)

  call AssertEqual(&
    message='nz_ToroidalEnergySpectrum_wc',                       &
    answer = nz_ToroidalEnergySpectrum_wc(wc_Torvel),             &
    check =  nz_EkTor,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  call AssertEqual(&
    message='nz_PoloidalEnergySpectrum_ws',                       &
    answer = nz_PoloidalEnergySpectrum_ws(ws_Polvel),             &
    check =  nz_EkPol,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call wsc_Potential2VectorMPI(xvz_Vlon,xvz_Vlat,xvz_Vrad,wc_Torvel,ws_Polvel)

  nmz_EkTor = nmz_ToroidalEnergySpectrum_wc(wc_Torvel)
  nmz_EkPol = nmz_PoloidalEnergySpectrum_ws(ws_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_xvz(xvz_VLon**2+xvz_VLat**2+xvz_VRad**2)/2.0D0, &
    check =  IntRad_z(z_Data/Ro**2),                              &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  nz_EkTor = nz_ToroidalEnergySpectrum_wc(wc_Torvel)
  nz_EkPol = nz_PoloidalEnergySpectrum_ws(ws_Polvel)

  call AssertEqual(&
    message='Total energy by nz-spectrum',                        &
    answer = IntLonLatRad_xvz(xvz_VLon**2+xvz_VLat**2+xvz_VRad**2)/2.0D0, &
    check =   IntRad_z(sum(nz_EkTor,1)/Ro**2) &
            + IntRad_z(sum(nz_EkPol,1)/Ro**2), &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call MessageNotify('M','wsc_mpi_module_supack_spectrum_test',&
                         'wsc_mpi_module_supack spectrum functions tests suceeded!') 

  CALL MPI_FINALIZE(IERR)

end program wsc_mpi_module_supack_spectrum_test
