!----------------------------------------------------------------------
!     Copyright (C) 2019--2020 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ut_mpi_module テストプログラム
!
!      磁場ポロイダルポテンシャルの境界値問題
!
!履歴  2019/05/29  竹広真一
!      2020/11/11  竹広真一  セクター計算オプション導入
!
program ut_mpi_module_polmagbc_mint_test

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

  implicit none

!!$  integer,parameter  :: im=32, jm=16, km=16  ! 格子点の設定(経度, 緯度, 動径)
!!$  integer,parameter  :: mint=1               ! 経度方向対称性
  integer,parameter  :: im=16, jm=16, km=16  ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: mint=2               ! 経度方向対称性
  integer,parameter  :: nm=10, lm=16         ! 切断波数の設定(水平, 動径)
  real(8),parameter  :: ri=0.5, ro=1.5       ! 内外半径

  integer :: npv = 2
  integer :: iproc, np, ierr

  real(8), allocatable  :: xvh_POLMAG(:,:,:)
  real(8), allocatable  :: ut_POLMAG(:,:)
  real(8), allocatable  :: ut_POLMAG0(:,:)

  real(8), allocatable  :: uz_POLMAG(:,:)
  real(8), allocatable  :: uz_POLMAG0(:,:)

  real(8), allocatable  :: uz_TopBoundary(:,:)
  real(8), allocatable  :: uz_BottomBoundary(:,:)

  real(8), allocatable  :: uz_Zero(:,:)

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

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

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

  integer :: k, n, nn(2)

  !---------------- MPI スタート ---------------------
  !call get_ENV_INT('NPV',npv)

  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','ut_mpi_module_polmagbc_mint_test', &
       'ut_mpi_module_mint ut_PolmagBoundaries subroutine test')

  call ut_mpi_initial(im,jm,km,nm,lm,ri,ro,npv=npv,mint=mint)

  allocate(xvh_POLMAG(0:im-1,jc,kc))
  allocate(ut_POLMAG(nmc,0:lm))
  allocate(ut_POLMAG0(nmc,0:lm))
  allocate(uz_POLMAG(nmc,0:km))
  allocate(uz_POLMAG0(nmc,0:km))

  allocate(uz_TopBoundary(nmc,0:km))
  allocate(uz_BottomBoundary(nmc,0:km))

  allocate(uz_Zero(nmc,0:km))

  allocate(uz_n(nmc,0:km))


  !=================== ut_PolmagBoundaries =======================
  ! P_10
  xvh_POLMAG = sin(xvh_lat) * sin( pi*(xvh_rad-ri)/(ro-ri) )

  ! P_1_1
  !xvh_POLMAG = cos(xvh_lat)*cos(xvh_lon)* sin( pi*(xvh_rad-ri)/(ro-ri) )
  !xvh_POLMAG = 2*sin(xvh_lat)**2 * sin( pi*(xvh_rad-ri)/(ro-ri) )

  ut_POLMAG = ut_xvh(xvh_POLMAG)
  ut_POLMAG0 = ut_POLMAG
  call ut_PolmagBoundaries(ut_POLMAG)

  do k=0,km
     do n=1,nmc
        nn=nm_l(n)
        uz_n(n,k) = nn(1)
     enddo
  enddo

  uz_Zero = 0.0D0
  uz_TopBoundary = uz_ut(ut_DRad_ut(ut_POLMAG)) &
                     + (uz_n +1)*uz_ut(ut_POLMAG)/uz_RAD
  uz_BottomBoundary = uz_ut(ut_DRad_ut(ut_POLMAG)) &
                     - uz_n * uz_ut(ut_POLMAG)/uz_RAD

  call AssertEqual(&
       message='ut_PolmagBoundaries (internal value)',           &
       answer = ut_POLMAG0(:,0:lm-2),                            &
       check = ut_POLMAG(:,0:lm-2),                              &
       significant_digits = check_digits, ignore_digits = ignore &
       )

  call AssertEqual(&
       message='ut_PolmagBoundaries (Top B.C.)',                 &
       answer = uz_Zero(:,0),                                    &
       check = uz_TopBoundary(:,0),                              &
       significant_digits = check_digits, ignore_digits = ignore &
       )

  call AssertEqual(&
       message='ut_PolmagBoundaries (Bottom B.C.)',              &
       answer = uz_Zero(:,km),                                   &
       check = uz_BottomBoundary(:,km),                          &
       significant_digits = check_digits, ignore_digits = ignore &
       )

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

  ! P_1_1
  !xvh_POLMAG = cos(xvh_lat)*cos(xvh_lon)* sin( pi*(xvh_rad-ri)/(ro-ri) )
  xvh_POLMAG = 2*sin(xvh_lat)**2 * sin( pi*(xvh_rad-ri)/(ro-ri) )

  ut_POLMAG = ut_xvh(xvh_POLMAG)
  uz_POLMAG0 = uz_ut(ut_POLMAG)

  call ut_PolmagBoundariesGrid(ut_POLMAG)
  uz_POLMAG = uz_ut(ut_POLMAG)

  do k=0,km
     do n=1,nmc
        nn=nm_l(n)
        uz_n(n,k) = nn(1)
     enddo
  enddo

  uz_TopBoundary = uz_ut(ut_DRad_ut(ut_POLMAG)) &
                     + (uz_n +1)*uz_ut(ut_POLMAG)/uz_RAD
  uz_BottomBoundary = uz_ut(ut_DRad_ut(ut_POLMAG)) &
                     - uz_n * uz_ut(ut_POLMAG)/uz_RAD

  uz_Zero = 0.0D0
  uz_TopBoundary = uz_ut(ut_DRad_ut(ut_POLMAG)) &
                     + (uz_n +1)*uz_ut(ut_POLMAG)/uz_RAD
  uz_BottomBoundary = uz_ut(ut_DRad_ut(ut_POLMAG)) &
                     - uz_n * uz_ut(ut_POLMAG)/uz_RAD

  call AssertEqual(&
       message='ut_PolmagBoundariesGrid (internal value)',       &
       answer = uz_POLMAG0(:,1:km-1),                            &
       check = uz_POLMAG(:,1:km-1),                              &
       significant_digits = check_digits, ignore_digits = ignore &
       )

  call AssertEqual(&
       message='ut_PolmagBoundarieGrid (Top B.C.)',                  &
       answer = uz_Zero(:,0),                                        &
       check = uz_TopBoundary(:,0),                                  &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_PolmagBoundariesGrid (Bottom B.C.)',              &
       answer = uz_Zero(:,km),                                       &
       check = uz_BottomBoundary(:,km),                              &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','ut_mpi_module_test_polmagbc', &
       'ut_mpi_module_mint ut_PolmagBoundaries subroutine test succeeded!')

  !------ MPIの終了 ------
  call MPI_FINALIZE(IERR)

end program ut_mpi_module_polmagbc_mint_test