!----------------------------------------------------------------------
!     Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ut_mpi_module_supack テストプログラム
!
!      ポロイダルポテンシャルの境界値問題
!
!履歴  2016/03/27  竹広真一   wt_module_nmomp_test_polvelbc.f90 より改造
!
program ut_mpi_module_svpack_polvelbc_test

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use ut_mpi_module_svpack
  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 :: xyh_Poloidal(:,:,:)
  real(8), allocatable :: xyh_LaplaPol(:,:,:)
  real(8), allocatable :: xyh_LaplaPol1(:,:,:)
  real(8), allocatable :: xyh_True(:,:,:)
  real(8), allocatable :: xyh_Zero(:,:,:)
  character(len=2), dimension(4), parameter :: BCond=(/'FF','FR','RF','RR'/)

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

  real(8), parameter  :: pi=3.1415926535897932385D0

  integer :: l

  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)

  call MessageNotify('M','ut_mpi_module_svpack_polvelbc_test', &
       'ut_mpi_module_svpack ut_LaplaPol2polGrid_ut function tests')

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

  allocate(xyh_Poloidal(0:im-1,jm,kc))
  allocate(xyh_LaplaPol(0:im-1,jm,kc))
  allocate(xyh_LaplaPol1(0:im-1,jm,kc))
  allocate(xyh_True(0:im-1,jm,kc))
  allocate(xyh_Zero(0:im-1,jm,kc))

  do l=1,4

     ! P_10
     xyh_Poloidal = sin(xyh_Lat) * sin( pi*(xyh_Rad-ri)/(ro-ri) )
     xyh_LaplaPol = xyh_ut(ut_Lapla_ut(ut_xyh(xyh_Poloidal)))
     !xyh_LaplaPol = sin(xyh_Lat) * sin( pi*(xyh_Rad-ri)/(ro-ri) )
     ! P_1_1
     !xyh_LaplaPol = cos(xyh_Lat)*cos(xyh_Lon)* sin( pi*(xyh_Rad-ri)/(ro-ri) )
     !xyh_LaplaPol = 2*sin(xyh_Lat)**2 * sin( pi*(xyh_Rad-ri)/(ro-ri) )

     !xyh_Poloidal = xyh_wz(wz_LaplaPol2pol_wz(wz_xyh(xyh_LaplaPol),BCond(l)))
     !xyh_Poloidal = xyh_ut(ut_LaplaPol2PolTau_ut(ut_xyh(xyh_LaplaPol),BCond(l)))

     xyh_Poloidal = xyh_ut(ut_LaplaPol2PolGrid_ut(ut_xyh(xyh_LaplaPol),BCond(l),new=.true.))

     xyh_LaplaPol1 = xyh_ut(ut_Lapla_ut(ut_xyh(xyh_Poloidal)))
     xyh_Zero = 0.0D0

     !---------------- 内部チェック -----------------------
     if ( kh_kz(0) /= 0 .and. kh_kz(km) /= 0 )then
        call AssertEqual(&
             message='ut_TorBoundariesGrid_ut (intenal value)',            &
             answer = xyh_LaplaPol1(:,:,3:kc-2),                           &
             check = xyh_LaplaPol(:,:,3:kc-2),                             &
             significant_digits = check_digits, ignore_digits = ignore     &
             )
     else if ( kh_kz(0) == 0 .and. kh_kz(km) /= 0 ) then
        call AssertEqual(&
             message='ut_TorBoundariesGrid_ut (intenal value)',            &
             answer = xyh_LaplaPol1(:,:,1:kc-2),                           &
             check = xyh_LaplaPol(:,:,1:kc-2),                             &
             significant_digits = check_digits, ignore_digits = ignore     &
             )
     else if ( kh_kz(0) /= 0 .and. kh_kz(km) == 0 ) then
        call AssertEqual(&
             message='ut_TorBoundariesGrid_ut (intenal value)',            &
             answer = xyh_LaplaPol1(:,:,kh_kz(0)+2:kc),                    &
             check = xyh_LaplaPol(:,:,kh_kz(0)+2:kc),                      &
             significant_digits = check_digits, ignore_digits = ignore     &
             )
     else 
        call AssertEqual(&
             message='ut_TorBoundariesGrid_ut (intenal value)',            &
             answer = xyh_LaplaPol1,                                       &
             check = xyh_LaplaPol,                                         &
             significant_digits = check_digits, ignore_digits = ignore     &
             )
     endif
     !--------- 上端境界チェック ----------

     !----- Φ=0 at the top ---------
     if ( kh_kz(0) /= 0 )then
        call AssertEqual(&
             message='ut_LaplaPol2polGrid_ut (Top B.C. Φ=0)',             &
             answer = xyh_Zero(:,:,kh_kz(0)),                              &
             check = xyh_Poloidal(:,:,kh_kz(0)),                           &
             significant_digits = check_digits, ignore_digits = ignore     &
             )
     endif
        
     !----- dΦ/dr=0, d^2Φ/dr^2 at the top ---------
     if( BCond(l)(1:1) == 'F' ) then
        xyh_True = xyh_ut(ut_DRad_ut(ut_DRad_ut(ut_xyh(xyh_Poloidal))))
     else
        xyh_True = xyh_ut(ut_DRad_ut(ut_xyh(xyh_Poloidal)))
     endif

     if ( kh_kz(0) /= 0 )then
        call AssertEqual(&
             message='ut_LaplaPol2polGrid_ut ('//BCond(l)//'-Top B.C)',    &
             answer = xyh_Zero(:,:,kh_kz(0)),                              &
             check = xyh_True(:,:,kh_kz(0)),                               &
             significant_digits = check_digits, ignore_digits = ignore     &
             )
     endif
        
     !--------- 下端境界チェック ----------

     !----- Φ=0 at the bottom ---------
     if ( kh_kz(km) /= 0 )then
        call AssertEqual(&
             message='ut_LaplaPol2polGrid_ut (Bottom B.C. Φ=0)',          &
             answer = xyh_Zero(:,:,kh_kz(km)),                             &
             check = xyh_Poloidal(:,:,kh_kz(km)),                          &
             significant_digits = check_digits, ignore_digits = ignore     &
             )
     endif

     !----- dΦ/dr=0, d^2Φ/dr^2 at the bottom ---------
     if( BCond(l)(2:2) == 'F' ) then
        xyh_True = xyh_ut(ut_DRad_ut(ut_DRad_ut(ut_xyh(xyh_Poloidal))))
     else
        xyh_True = xyh_ut(ut_DRad_ut(ut_xyh(xyh_Poloidal)))
     endif

     if ( kh_kz(km) /= 0 )then
        call AssertEqual(&
             message='ut_LaplaPol2polGrid_ut ('//BCond(l)//'-Bottom B.C)', &
             answer = xyh_Zero(:,:,kh_kz(km)),                             &
             check = xyh_True(:,:,kh_kz(km)),                              &
             significant_digits = check_digits, ignore_digits = ignore     &
             )
     endif
  end do

  call MessageNotify('M','ut_mpi_module_svpack_polvelbc_test', &
       'ut_mpi_module_svpack ut_LaplaPol2polGrid_ut function tests suceeded!')

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

  call MPI_FINALIZE(IERR)

end program ut_mpi_module_svpack_polvelbc_test