!----------------------------------------------------------------------
!     Copyright (c) 2008--2012 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wtq_module テストプログラム
!
!      ポロイダルポテンシャルの境界値問題
!
!履歴  2009/12/11  竹広真一   wq_test_polvelbc.f90 より SJPACK 用に改造
!
!メモ
!      2008/07/05  佐々木洋平 1D-1 しか精度保証できないの?
!      2011/09/09  竹広真一   wq_module_polvelbc_test.f90 を MPI 化
!      2011/09/12  竹広真一   wtq_mpi_module 用に改造
!      2012/04/02  竹広真一   wtq_mpi_module_sjpack 用に改造
!      2012/04/03  竹広真一   wtq_module_sjpack 用に改造
!      2012/04/03  竹広真一   wtq_module 用に改造
!
program wtq_module_polvelbc_wq_test

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use wtq_module

  implicit none

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

  real(8), dimension(0:im-1,1:jm,kmi)     :: xyr_Poloidal
  real(8), dimension(0:im-1,1:jm,kmi)     :: xyr_LaplaPol
  real(8), dimension(0:im-1,1:jm,kmi)     :: xyr_True
  real(8), dimension(0:im-1,1:jm,kmi)     :: xyr_Zero=0.0D0
  real(8), dimension((nm+1)**2,0:lmi)     :: wq_Poloidal
  real(8), dimension((nm+1)**2,0:lmi)     :: wq_LaplaPol
  real(8), dimension((nm+1)**2,0:lmi)     :: wq_LaplaPol0
  character(len=1), dimension(2), parameter :: BCond=(/'F','R'/)

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

  integer :: l

  call MessageNotify('M','wtq_module_polvelbc_wq_test', &
       'wtq_module wq_LaplaPol2Pol_wq subroutine test')

  call wtq_Initial(im,jm,kmi,kmo,nm,lmi,lmo,ri,ro)

  do l=1,2

     ! P_10
     xyr_Poloidal = sin(xyr_Lat) * xyr_Rad * ((xyr_Rad-ri)*(xyr_Rad+ri))**3
!!$     xyr_Poloidal = sin(xyr_Lat) * xyr_Rad * ((xyr_Rad-ri)*(xyr_Rad+ri))**2
     wq_LaplaPol0 = wq_Lapla_wq(wq_xyr(xyr_Poloidal))
     !xyr_LaplaPol = sin(xyr_Lat) * sin( pi*(xyr_Rad-ri)/(ro-ri) )
     ! P_1_1
     !xyr_LaplaPol = cos(xyr_Lat)*cos(xyr_Lon)* xyr_Rad * (xyr_Rad-ri) 
     !xyr_LaplaPol = 2*sin(xyr_Lat)**2 * xyr_Rad * (xyr_Rad-ri) 

     !xyr_Poloidal = xyr_wz(wr_LaplaPol2pol_wr(wr_xyr(xyr_LaplaPol),BCond(l)))
     !xyr_Poloidal = xyr_wq(wq_LaplaPol2PolTau_wq(wq_xyr(xyr_LaplaPol),BCond(l)))

     wq_Poloidal  = wq_LaplaPol2Pol_wq(wq_LaplaPol0,BCond(l),new=.true.)
     wq_LaplaPol  = wq_Lapla_wq(wq_Poloidal)
     xyr_Poloidal = xyr_wq(wq_Poloidal)

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

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

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

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

     call AssertEqual(&
          message='wq_LaplaPol2Pol_wq (Top Booundary)',                 &
          answer = xyr_True(:,:,kmi),                                    &
          check = xyr_Zero(:,:,kmi),                                     &
          significant_digits = check_digits, ignore_digits = ignore     &
          )

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

  call MessageNotify('M','wtq_module_polvelbc_wq_test', &
       'wtq_module wq_LaplaPol2Pol_wq subroutine test succeeded!')

end program wtq_module_polvelbc_wq_test