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

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use mpi
  use wtq_mpi_module_supack

  implicit none

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

  real(8), allocatable :: xvz_POLMAG(:,:,:)
  real(8), allocatable :: wz_POLMAG(:,:)
  real(8), allocatable :: xvr_POLMAG(:,:,:)
  real(8), allocatable :: wr_POLMAG(:,:)
  real(8), allocatable :: xvz_DPOLDR(:,:,:)
  real(8), allocatable :: xvr_DPOLDR(:,:,:)

  real(8), allocatable :: wz_TopBoundary(:,:)
  real(8), allocatable :: wz_n(:,:)   ! 全波数
  real(8), allocatable :: w_Null(:)

  real(8), parameter  :: pi=3.1415926535897932385D0
  ! 判定誤差設定
  integer, parameter :: check_digits = 11
  integer, parameter :: ignore = -12

  integer :: n,k, nn(2)
  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','wtq_mpi_module_supack_polmagbcgrid_test', &
       'wtq_mpi_module_supack  wtq_TormagBoundariesGrid subroutine test')

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

  allocate(xvz_POLMAG(im,jc,0:kmo))
  allocate(wz_POLMAG(nc,0:kmo))
  allocate(xvr_POLMAG(im,jc,kmi))
  allocate(wr_POLMAG(nc,kmi))
  allocate(xvz_DPOLDR(im,jc,0:kmo))
  allocate(xvr_DPOLDR(im,jc,kmi))              
  allocate(wz_TopBoundary(nc,0:kmo))

  allocate(wz_n(nc,0:kmo))
  allocate(w_Null(nc))

  w_Null = 0.0D0
  
  do k=0,kmo
     do n=1,nc
        nn=nm_l(n)
        wz_n(n,k) = nn(1)
     enddo
  enddo

  ! P_10

  xvz_POLMAG = sin(xvz_lat) * cos( pi*(xvz_rad-ri)/(ro-ri) )
  xvr_POLMAG = sin(xvr_lat) * cos( pi*(xvr_rad-ri)/ri ) * xvr_Rad

  wz_POLMAG = wz_xvz(xvz_POLMAG)
  wr_POLMAG = wr_xvr(xvr_POLMAG)
  call wtq_PolmagBoundariesGrid(wz_POLMAG,wr_POLMAG)
  xvz_POLMAG = xvz_wz(wz_POLMAG)
  xvr_POLMAG = xvr_wr(wr_POLMAG)
  xvz_DPOLDR = xvz_wt(wt_DRad_wt(wt_wz(wz_POLMAG)))
  xvr_DPOLDR = xvr_wq(wq_RadDRad_wq(wq_wr(wr_POLMAG)))/xvr_Rad
  wz_TopBoundary = wz_wt(wt_DRad_wt(wt_wz(wz_POLMAG))) &
                     + (wz_n +1)*wz_POLMAG/wz_RAD

  call checkresults
  ! P_1_1
  xvz_POLMAG = cos(xvz_lat)*cos(xvz_lon)* cos( pi*(xvz_rad-ri)/(ro-ri) )
  xvr_POLMAG = cos(xvr_lat)*cos(xvr_lon)* cos( pi*(xvr_rad-ri)/ri )*xvr_Rad

  wz_POLMAG = wz_xvz(xvz_POLMAG)
  wr_POLMAG = wr_xvr(xvr_POLMAG)
  call wtq_PolmagBoundariesGrid(wz_POLMAG,wr_POLMAG)
  xvz_POLMAG = xvz_wz(wz_POLMAG)
  xvr_POLMAG = xvr_wr(wr_POLMAG)

  xvz_DPOLDR = xvz_wt(wt_DRad_wt(wt_wz(wz_POLMAG)))
  xvr_DPOLDR = xvr_wq(wq_RadDRad_wq(wq_wr(wr_POLMAG)))/xvr_Rad
  wz_TopBoundary = wz_wt(wt_DRad_wt(wt_wz(wz_POLMAG))) &
                     + (wz_n +1)*wz_POLMAG/wz_RAD

  call checkresults

  call MessageNotify('M','wtq_test_polmagbcgrid', &
       'wtq_module_supack  wtq_PolmagBoundariesGrid subroutine test succeded')

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

  call MPI_FINALIZE(IERR)

contains
  subroutine checkresults

    call AssertEqual(&
         message='Top B.C',                                            &
         answer = w_cleanup_w(wz_TopBoundary(:,0)),                    &
         check = w_Null,                                               &
         significant_digits = check_digits, ignore_digits = ignore     &
         )

    call AssertEqual(&
         message='Inner B.C (value conituity) ',                       &
         answer = xvz_POLMAG(:,:,kmo),                                 &
         check = xvr_POLMAG(:,:,kmi),                                  &
         significant_digits = check_digits, ignore_digits = ignore     &
         )

    call AssertEqual(&
         message='Inner B.C (derivative conituity) ',                  &
         answer = xvz_DPOLDR(:,:,kmo),                                 &
         check = xvr_DPOLDR(:,:,kmi),                                  &
         significant_digits = check_digits, ignore_digits = ignore     &
         )
  end subroutine checkresults

  function w_cleanup_w(w_data)
    real(8), intent(IN) :: w_data(nc)
    real(8)             :: w_cleanup_w(nc)

    integer l,nmary(2)

    w_cleanup_w = w_data

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

  end function w_cleanup_w

end program wtq_mpi_supack_polmagbcgrid_test