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

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

  implicit none

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

  character(len=2), dimension(4),parameter :: BCond=(/'FF','FR','RF','RR'/)


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

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

  integer :: l

  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_torvelbc_test', &
       'wtq_mpi_module_supack wt_TorgBoundariesGrid subroutine test')

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

  allocate(xvz_TorVel0(0:im-1,jc,0:kmo))
  allocate(xvz_TorVel(0:im-1,jc,0:kmo))
  allocate(wt_TorVel(nc,0:lmo))
  allocate(wt_TorVel0(nc,0:lmo))
  allocate(xvz_Zero(0:im-1,jc,0:kmo))
  allocate(xvz_True(0:im-1,jc,0:kmo))

  xvz_Zero = 0.0D0

  do l=1,4

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

     xvz_TorVel0 = xvz_TorVel
     wt_TorVel = wt_xvz(xvz_TorVel)
     wt_TorVel0 = wt_TorVel
     
     !---- wt_TorBoundaries test
     call wt_TorBoundaries(wt_TorVel,cond=BCond(l),new=.true.)
     xvz_TorVel = xvz_wt(wt_TorVel)

     call AssertEqual(&
          message='wt_TorvelBoundaries (internal value)',                &
          answer = ua_cleanup_ua(wt_Torvel(:,0:lmi-2)),                 &
          check = ua_cleanup_ua(wt_Torvel0(:,0:lmi-2)),                 &
          significant_digits = check_digits, ignore_digits = ignore     &
       )

     ! 上端境界チェック
     if( BCond(l)(1:1) == 'F' ) then
        xvz_True = xvz_wt(wt_DRad_wt(wt_xvz(xvz_TorVel/xvz_Rad)))
     else
        xvz_True = xvz_TorVel
     endif

     call AssertEqual(&
          message='wt_TorvelBoundaries (Top boundary)',                 &
          answer = xvz_True(:,:,0),                                     &
          check = xvz_Zero(:,:,0),                                      &
          significant_digits = check_digits, ignore_digits = ignore     &
          )


     ! 下端境界チェック
     if( BCond(l)(2:2) == 'F' ) then
        xvz_True = xvz_wt(wt_DRad_wt(wt_xvz(xvz_TorVel/xvz_Rad)))
     else
        xvz_True = xvz_TorVel
     endif

     call AssertEqual(&
          message='wt_TorvelBoundaries_wt ('//BCond(l)//'-Bottom B.C)', &
          answer = xvz_Zero(:,:,kmo),                                   &
          check = xvz_True(:,:,kmo),                                    &
          significant_digits = check_digits, ignore_digits = ignore     &
          )
     call MessageNotify('M','wtq_module_supack_torvelbc_wt_test', &
          'wt_TorBoundarry: '//BCond(l)//'-Top B.C. test succeeded!')


     !---- ur_TorBoundaryGrid test
     wt_TorVel = wt_xvz(xvz_TorVel0)
     call wt_TorBoundariesGrid(wt_TorVel,cond=BCond(l),new=.true.)
     xvz_TorVel = xvz_wt(wt_TorVel)

     ! 内部チェック
     call AssertEqual(&
          message='wt_TorBoundariesGrid_wt (intenal value)',            &
          answer = xvz_Torvel(:,:,1:kmo-1),                             &
          check = xvz_Torvel0(:,:,1:kmo-1),                             &
          significant_digits = check_digits, ignore_digits = ignore     &
          )

     ! 上端境界チェック
     if( BCond(l)(1:1) == 'F' ) then
        xvz_True = xvz_wt(wt_DRad_wt(wt_xvz(xvz_TorVel/xvz_Rad)))
     else
        xvz_True = xvz_TorVel
     endif

     call AssertEqual(&
          message='wt_TorvelBoundariesGrid_wt ('//BCond(l)//'-Top B.C)',&
          answer = xvz_Zero(:,:,0),                                     &
          check = xvz_True(:,:,0),                                      &
          significant_digits = check_digits, ignore_digits = ignore     &
          )

     ! 下端境界チェック
     if( BCond(l)(2:2) == 'F' ) then
        xvz_True = xvz_wt(wt_DRad_wt(wt_xvz(xvz_TorVel/xvz_Rad)))
     else
        xvz_True = xvz_TorVel
     endif

     call AssertEqual(&
          message='wt_TorvelBoundariesGrid_wt ('//BCond(l)//'-Bottom B.C)',&
          answer = xvz_Zero(:,:,kmo),                                   &
          check = xvz_True(:,:,kmo),                                    &
          significant_digits = check_digits, ignore_digits = ignore     &
          )
  enddo

  call MessageNotify('M','wtq_mpi_module_supack_torvelbc_test', &
       'wtq_mpi_module_supack  wt_TorvelBoundariesGrid subroutine test succeeded!')

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

  call MPI_FINALIZE(IERR)

contains
  
  function ua_cleanup_ua(ua_data)
    real(8), intent(IN) :: ua_data(:,:)
    real(8)             :: ua_cleanup_ua(nc,size(ua_data,2))

    integer l,nmary(2)

    ua_cleanup_ua = ua_data

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

  end function ua_cleanup_ua

end program wtq_mpi_module_supack_torvelbc_test
