!----------------------------------------------------------------------
!     Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wq_mpi_module_supack テストプログラム
!
!      ポロイダルポテンシャルの境界値問題
!
!履歴  2016/04/01  竹広真一   wq_mpi_module_sjpack_polvelbc_test.f90 を supack 化
!
!
program wq_mpi_module_supack_polvelbc

  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=8   ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: nm=10, lm=15         ! 切断波数の設定(水平, 動径)
  real(8),parameter  :: ra=0.5D0             ! 球半径

  real(8), allocatable :: xvr_Poloidal(:,:,:)
  real(8), allocatable :: xvr_LaplaPol(:,:,:)
  real(8), allocatable :: xvr_True(:,:,:)
  real(8), allocatable :: xvr_Zero(:,:,:)
  real(8), allocatable :: wq_Poloidal(:,:)
  real(8), allocatable :: wq_LaplaPol(:,:)
  real(8), allocatable :: wq_LaplaPol0(:,:)
  character(len=1), dimension(2), parameter :: BCond=(/'F','R'/)

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

  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','wq_mpi_module_supack_polvelbc_test', &
       'wq_mpi_module_supack wq_LaplaPol2Pol_wq subroutine test')

  call wq_mpi_initial(im,jm,km,nm,lm,ra)

  allocate(xvr_Poloidal(0:im-1,1:jc,km))
  allocate(xvr_LaplaPol(0:im-1,1:jc,km))
  allocate(xvr_True(0:im-1,1:jc,km))
  allocate(xvr_Zero(0:im-1,1:jc,km))
  allocate(wq_Poloidal(nc,0:lm))
  allocate(wq_LaplaPol(nc,0:lm))
  allocate(wq_LaplaPol0(nc,0:lm))

  xvr_Zero = 0.0D0

  
  do l=1,2
     ! P_10
     xvr_Poloidal = sin(xvr_Lat) * xvr_Rad * ((xvr_Rad-ra)*(xvr_Rad+ra))**3
!!$     xvr_Poloidal = sin(xvr_Lat) * xvr_Rad * ((xvr_Rad-ra)*(xvr_Rad+ra))**2
     wq_LaplaPol0 = wq_Lapla_wq(wq_xvr(xvr_Poloidal))
     !xvr_LaplaPol = sin(xvr_Lat) * sin( pi*(xvr_Rad-ri)/(ro-ri) )
     ! P_1_1
     !xvr_LaplaPol = cos(xvr_Lat)*cos(xvr_Lon)* xvr_Rad * (xvr_Rad-ra) 
     !xvr_LaplaPol = 2*sin(xvr_Lat)**2 * xvr_Rad * (xvr_Rad-ra) 

     !xvr_Poloidal = xvr_wz(wr_LaplaPol2pol_wr(wr_xvr(xvr_LaplaPol),BCond(l)))
     !xvr_Poloidal = xvr_wq(wq_LaplaPol2PolTau_wq(wq_xvr(xvr_LaplaPol),BCond(l)))

     wq_Poloidal  = wq_LaplaPol2Pol_wq(wq_LaplaPol0,BCond(l),new=.true.)
     wq_LaplaPol  = wq_Lapla_wq(wq_Poloidal)
     xvr_Poloidal = xvr_wq(wq_Poloidal)

     !---------------- 内部チェック -----------------------
     call AssertEqual(&
          message='wq_Laplapol2Pol_wq (internal value)',                &
          answer = wq_cleanup_wq(wq_LaplaPol(:,0:lm-4)),                &
          check = wq_cleanup_wq(wq_LaplaPol0(:,0:lm-4)),                &
          significant_digits = check_digits, ignore_digits = ignore     &
       )

     !--------- 上端境界チェック ----------

     !----- Φ=0 at the top ---------
     call AssertEqual(&
          message='wq_LaplaPol2Pol_wq (Φ=0 Top Booundary)',            &
          answer = xvr_Poloidal(:,:,km),                                &
          check = xvr_Zero(:,:,km),                                     &
          significant_digits = check_digits, ignore_digits = ignore     &
          )

     !----- dΦ/dr=0, d^2Φ/dr^2 at the top ---------
     if( BCond(l)(1:1) == 'F' ) then
        xvr_True = xvr_wq(wq_RadDRad_wq(wq_RadDRad_wq(wq_Poloidal))&
                                       -wq_RadDRad_wq(wq_Poloidal))/xvr_Rad**2
     else
        xvr_True = xvr_wq(wq_RadDRad_wq(wq_Poloidal))/xvr_Rad
     endif

     call AssertEqual(&
          message='wq_LaplaPol2Pol_wq (Top Booundary)',                 &
          answer = xvr_True(:,:,km),                                    &
          check = xvr_Zero(:,:,km),                                     &
          significant_digits = check_digits, ignore_digits = ignore     &
          )

     call MessageNotify('M','wq_module_supack_polvelbc_test', &
          'wq_LaplaPol2Pol_wq: '//BCond(l)//'-Top B.C. test succeeded!')
          
  end do

  call MessageNotify('M','wq_mpi_module_supack_polvelbc_test', &
       'wq_mpi_module_supack wq_LaplaPol2Pol_wq subroutine test succeeded!')

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

  call MPI_FINALIZE(IERR)

contains

  function wq_cleanup_wq(wq_data)
    real(8), intent(IN) :: wq_data(:,:)
    real(8)             :: wq_cleanup_wq(size(wq_data,1),size(wq_data,2))

    integer l,nmary(2)

    wq_cleanup_wq = wq_data

    do l=1,nc
       nmary = nm_l(l)
       if ( nmary(1) .lt. 0 ) wq_cleanup_wq(l,:) = 0.0D0
    end do

  end function wq_cleanup_wq

end program wq_mpi_module_supack_polvelbc