!----------------------------------------------------------------------
!     Copyright (c) 2012 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wtq_module_sjpack テストプログラム
!
!      磁場ポロイダルポテンシャルの境界値問題
!
!履歴  2012/04/02  竹広真一   wtq_mpi_module_polmagbc_test.f90 を sjpack 化
!      2012/04/03  竹広真一   wtq_module_sjpack_polmagbc_test.f90 に改造
!
program wtq_module_sjpack_polmagbc_wq

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

  implicit none

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

  real(8), dimension(0:im-1,1:jm,kmi)       :: xyr_POLMAG
  real(8), dimension(0:im-1,1:jm,kmi)       :: xyr_POLMAG_orig
  real(8), dimension((nm+1)*(nm+1),0:lmi)   :: wq_POLMAG
  real(8), dimension((nm+1)*(nm+1),0:lmi)   :: wq_POLMAG_orig
  real(8), dimension((nm+1)*(nm+1),kmi)     :: wr_POLMAG

  real(8), dimension((nm+1)*(nm+1),kmi)     :: wr_TopBoundary
  real(8), dimension((nm+1)*(nm+1),kmi)     :: wr_Zero = 0.0D0

  real(8), dimension((nm+1)*(nm+1),kmi)     :: wr_n   ! 全波数

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

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

  integer :: k, n, nn(2)

  call MessageNotify('M','wtq_module_sjpack_polmagbc_wq_test', &
       'wtq_module_sjpack wq_PolmagBoundary subroutine test')

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

  !=================== wq_PolmagBoundary =======================
  ! P_10
  xyr_POLMAG = sin(xyr_lat) * sin( pi*xyr_rad/ri )

  ! P_1_1
  !xyr_POLMAG = cos(xyr_lat)*cos(xyr_lon)* sin( pi*(xyr_rad-ri)/ri )
  !xyr_POLMAG = 2*sin(xyr_lat)**2 * sin( pi*(xyr_rad-ri)/ri )

  xyr_POLMAG_orig = xyr_POLMAG
  wq_POLMAG = wq_xyr(xyr_POLMAG)
  wq_POLMAG_orig = wq_POLMAG
  call wq_PolmagBoundary(wq_POLMAG)

  call AssertEqual(&
       message='wq_PolmagBoundary (internal value)',                 &
       answer = wq_Polmag(:,0:lmi-1),                                 &
       check = wq_Polmag_orig(:,0:lmi-1),                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  do k=1,kmi
     do n=1,(nm+1)**2
        nn=nm_l(n)
        wr_n(n,k) = nn(1)
     enddo
  enddo

  wr_TopBoundary = wr_wq(wq_RadDRad_wq(wq_POLMAG))/wr_RAD &
                     + (wr_n +1)*wr_wq(wq_POLMAG)/wr_RAD

  call AssertEqual(&
       message='wq_PolmagBoundary (Top B.C.)',                       &
       answer = wr_TopBoundary(:,kmi),                                &
       check = wr_Zero(:,kmi),                                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !=================== wq_PolmagBoundaryGrid =======================
  ! P_10
  !xyr_POLMAG = sin(xyr_lat) * sin( pi*(xyr_rad-ri)/(ro-ri) )

  ! P_1_1
  !xyr_POLMAG = cos(xyr_lat)*cos(xyr_lon)* sin( pi*(xyr_rad-ri)/(ro-ri) )
  xyr_POLMAG = 2*sin(xyr_lat)**2 * sin( pi*(xyr_rad-ri)/ri ) * xyr_Rad

  xyr_POLMAG_orig = xyr_POLMAG
  wr_POLMAG = wr_xyr(xyr_POLMAG)
  call wr_PolmagBoundaryGrid(wr_POLMAG)
  xyr_POLMAG = xyr_wr(wr_POLMAG)
  wq_POLMAG  = wq_wr(wr_POLMAG)

  call AssertEqual(&
       message='wq_PolmagBoundaryGrid (internal value)',             &
       answer = xyr_Polmag(:,:,1:kmi-1),                              &
       check = xyr_Polmag_orig(:,:,1:kmi-1),                          &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  do k=1,kmi
     do n=1,(nm+1)**2
        nn=nm_l(n)
        wr_n(n,k) = nn(1)
     enddo
  enddo

  wr_TopBoundary = wr_wq(wq_RadDRad_wq(wq_POLMAG))/wr_RAD &
                      + (wr_n +1)*wr_POLMAG/wr_RAD
  call AssertEqual(&
       message='wq_PolmagBoundaryGrid (Top B.C.)',                   &
       answer = wr_TopBoundary(:,kmi),                                &
       check = wr_Zero(:,kmi),                                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','wtq_module_sjpack_polmagbc_wq_test', &
       'wtq_module_sjpack wq_PolmagBoundary subroutine test succeeded!')

end program wtq_module_sjpack_polmagbc_wq