!----------------------------------------------------------------------
!     Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wq_mpi_module_supack テストプログラム
!
!      磁場ポロイダルポテンシャルの境界値問題
!
!履歴  2016/04/01  竹広真一   wq_mpi_module_sjpack_polmagbc_test.f90 を supack 化
!
program wq_mpi_module_supack_polmagbc

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

  implicit none

  integer,parameter  :: im=32, jm=16, km=8   ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: nm=10, lm=15         ! 切断波数の設定(水平, 動径)
  real(8),parameter  :: ra=2.5               ! 球半径

  real(8), allocatable :: xvr_POLMAG(:,:,:)
  real(8), allocatable :: xvr_POLMAG_orig(:,:,:)
  real(8), allocatable :: wq_POLMAG(:,:)
  real(8), allocatable :: wq_POLMAG_orig(:,:)
  real(8), allocatable :: wr_POLMAG(:,:)

  real(8), allocatable :: wr_TopBoundary(:,:)
  real(8), allocatable :: wr_Zero(:,:)

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

  ! 判定誤差設定
  integer, parameter :: check_digits = 12
  integer, parameter :: ignore = -13

  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','wq_mpi_module_supack_polmagbc_test', &
       'wq_mpi_module_supack wq_PolmagBoundary subroutine test')

  call wq_mpi_initial(im,jm,km,nm,lm,ra)

  allocate(xvr_POLMAG(0:im-1,1:jc,km))
  allocate(xvr_POLMAG_orig(0:im-1,1:jc,km))
  allocate(wq_POLMAG(nc,0:lm))
  allocate(wq_POLMAG_orig(nc,0:lm))
  allocate(wr_POLMAG(nc,km))
                 
  allocate(wr_TopBoundary(nc,km))
  allocate(wr_Zero(nc,km))
                 
  allocate(wr_n(nc,km))

  wr_Zero = 0.0D0
  
  !=================== wq_PolmagBoundary =======================
  ! P_10
  xvr_POLMAG = sin(xvr_lat) * sin( pi*xvr_rad/ra )

  ! P_1_1
  !xvr_POLMAG = cos(xvr_lat)*cos(xvr_lon)* sin( pi*(xvr_rad-ra)/ra )
  !xvr_POLMAG = 2*sin(xvr_lat)**2 * sin( pi*(xvr_rad-ra)/ra )

  xvr_POLMAG_orig = xvr_POLMAG
  wq_POLMAG = wq_xvr(xvr_POLMAG)
  wq_POLMAG_orig = wq_POLMAG
  call wq_PolmagBoundary(wq_POLMAG)

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

  do k=1,km
     do n=1,nc
        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 = w_cleanup_w(wr_TopBoundary(:,km)),                   &
       check = wr_Zero(:,km),                                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !=================== wq_PolmagBoundaryGrid =======================
  ! P_10
  !xvr_POLMAG = sin(xvr_lat) * sin( pi*(xvr_rad-ri)/(ro-ri) )

  ! P_1_1
  !xvr_POLMAG = cos(xvr_lat)*cos(xvr_lon)* sin( pi*(xvr_rad-ri)/(ro-ri) )
  xvr_POLMAG = 2*sin(xvr_lat)**2 * sin( pi*(xvr_rad-ra)/ra ) * xvr_Rad

  xvr_POLMAG_orig = xvr_POLMAG
  wr_POLMAG = wr_xvr(xvr_POLMAG)
  call wr_PolmagBoundaryGrid(wr_POLMAG)
  xvr_POLMAG = xvr_wr(wr_POLMAG)
  wq_POLMAG  = wq_wr(wr_POLMAG)

  call AssertEqual(&
       message='wq_PolmagBoundaryGrid (internal value)',             &
       answer = xvr_Polmag(:,:,1:km-1),                              &
       check = xvr_Polmag_orig(:,:,1:km-1),                          &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  do k=1,km
     do n=1,nc
        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 = w_cleanup_w(wr_TopBoundary(:,km)),                   &
       check = wr_Zero(:,km),                                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','wq_mpi_module_supack_polmagbc_test', &
       'wq_mpi_module_supack wq_PolmagBoundary 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 wq_mpi_module_supack_polmagbc