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

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use wt_mpi_module_supack
  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  :: xvz_POLMAG(:,:,:)
  real(8), allocatable  :: wt_POLMAG(:,:)

  real(8), allocatable  :: wz_TopBoundary(:,:)
  real(8), allocatable  :: wz_BottomBoundary(:,:)

  real(8), allocatable  :: wz_Zero(:,:)

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

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

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

  integer :: k, n, 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','wt_mpi_module_supack_polmagbc_test', &
       'wt_mpi_module_supack wt_PolmagBoundaries subroutine test')

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

  allocate(xvz_POLMAG(0:im-1,jc,0:km))
  allocate(wt_POLMAG(nc,0:lm))

  allocate(wz_TopBoundary(nc,0:km))
  allocate(wz_BottomBoundary(nc,0:km))

  allocate(wz_Zero(nc,0:km))

  allocate(wz_n(nc,0:km))


  !=================== wt_PolmagBoundaries =======================
  ! P_10
  xvz_POLMAG = sin(xvz_lat) * sin( pi*(xvz_rad-ri)/(ro-ri) )

  ! P_1_1
  !xvz_POLMAG = cos(xvz_lat)*cos(xvz_lon)* sin( pi*(xvz_rad-ri)/(ro-ri) )
  !xvz_POLMAG = 2*sin(xvz_lat)**2 * sin( pi*(xvz_rad-ri)/(ro-ri) )

  wt_POLMAG = wt_xvz(xvz_POLMAG)
  call wt_PolmagBoundaries(wt_POLMAG)
  xvz_POLMAG = xvz_wt(wt_POLMAG)

  do k=0,km
     do n=1,nc
        nn=nm_l(n)
        wz_n(n,k) = nn(1)
     enddo
  enddo

  wz_Zero = 0.0D0
  wz_TopBoundary = wz_wt(wt_DRad_wt(wt_POLMAG)) &
                     + (wz_n +1)*wz_wt(wt_POLMAG)/wz_RAD
  wz_BottomBoundary = wz_wt(wt_DRad_wt(wt_POLMAG)) &
                     - wz_n * wz_wt(wt_POLMAG)/wz_RAD 
  call AssertEqual(&
       message='wt_PolmagBoundaries (Top B.C.)',                     &
       answer = wz_Zero(:,0),                                        &
       check = w_cleanup_w(wz_TopBoundary(:,0)),                     &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  call AssertEqual(&
       message='wt_PolmagBoundaries (Bottom B.C.)',                  &
       answer = wz_Zero(:,km),                                       &
       check = w_cleanup_w(wz_BottomBoundary(:,km)),                 &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !=================== wt_PolmagBoundariesGrid =======================
  ! P_10
  !xyz_POLMAG = sin(xyz_lat) * sin( pi*(xyz_rad-ri)/(ro-ri) )

  ! P_1_1
  !xvz_POLMAG = cos(xvz_lat)*cos(xvz_lon)* sin( pi*(xvz_rad-ri)/(ro-ri) )
  xvz_POLMAG = 2*sin(xvz_lat)**2 * sin( pi*(xvz_rad-ri)/(ro-ri) )

  wt_POLMAG = wt_xvz(xvz_POLMAG)
  call wt_PolmagBoundariesGrid(wt_POLMAG)
  xvz_POLMAG = xvz_wt(wt_POLMAG)

  do k=0,km
     do n=1,nc
        nn=nm_l(n)
        wz_n(n,k) = nn(1)
     enddo
  enddo

  wz_TopBoundary = wz_wt(wt_DRad_wt(wt_POLMAG)) &
                     + (wz_n +1)*wz_wt(wt_POLMAG)/wz_RAD
  wz_BottomBoundary = wz_wt(wt_DRad_wt(wt_POLMAG)) &
                     - wz_n * wz_wt(wt_POLMAG)/wz_RAD 

  wz_Zero = 0.0D0
  wz_TopBoundary = wz_wt(wt_DRad_wt(wt_POLMAG)) &
                     + (wz_n +1)*wz_wt(wt_POLMAG)/wz_RAD
  wz_BottomBoundary = wz_wt(wt_DRad_wt(wt_POLMAG)) &
                     - wz_n * wz_wt(wt_POLMAG)/wz_RAD 
  call AssertEqual(&
       message='wt_PolmagBoundariesTau (Top B.C.)',                  &
       answer = wz_Zero(:,0),                                        &
       check = w_cleanup_w(wz_TopBoundary(:,0)),                     &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  call AssertEqual(&
       message='wt_PolmagBoundariesTau (Bottom B.C.)',               &
       answer = wz_Zero(:,km),                                       &
       check = w_cleanup_w(wz_BottomBoundary(:,km)),                 &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','wt_mpi_module_supack_polmagbc_test', &
       'wt_mpi_module_supack wt_PolmagBoundaries subroutine test succeeded!')

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

  call MPI_FINALIZE(IERR)

contains

  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 wt_mpi_module_supack_polmagbc_test