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

  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=10, lmi=16       ! 切断波数の設定(球殻動径, 球動径)
  real(8),parameter  :: ri=0.5, ro=1.5       ! 内外半径
 
  real(8), allocatable :: xvr_TorVel0(:,:,:)
  real(8), allocatable :: xvr_TorVel(:,:,:)
  real(8), allocatable :: wq_TorVel(:,:)
  real(8), allocatable :: wq_TorVel0(:,:)
  real(8), allocatable :: wr_TorVel(:,:)
  character(len=1), dimension(2),parameter :: BCond=(/'F','R'/)

  real(8), allocatable :: xvr_True(:,:,:)
  real(8), allocatable :: xvr_Zero(:,:,:)

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

  integer :: ll
  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 wq_TorvelBoundary subroutine test')

  call wtq_mpi_Initial(im,jm,kmi,kmo,nm,lmi,lmo,ri,ro)
 
  allocate(xvr_TorVel0(0:im-1,1:jc,kmi))
  allocate(xvr_TorVel(0:im-1,1:jc,kmi))
  allocate(wq_TorVel(nc,0:lmi))
  allocate(wq_TorVel0(nc,0:lmi))
  allocate(wr_TorVel(nc,kmi))
              
  allocate(xvr_True(0:im-1,1:jc,kmi))
  allocate(xvr_Zero(0:im-1,1:jc,kmi))

  xvr_Zero = 0.0D0
  
  do ll=1,2

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

     xvr_TorVel0 = xvr_TorVel
     wq_TorVel = wq_xvr(xvr_TorVel)
     wq_TorVel0 = wq_TorVel  
     wr_TorVel = wr_xvr(xvr_TorVel)


     !---- wq_TorBoundary test
     call wq_TorBoundary(wq_TorVel,cond=BCond(ll),new=.true.)
     xvr_TorVel = xvr_wq(wq_TorVel)

     call AssertEqual(&
          message='wq_TorvelgBoundary (internal value)',                &
          answer = wa_cleanup_wa(wq_Torvel(:,0:lmi-2)),                 &
          check = wa_cleanup_wa(wq_Torvel0(:,0:lmi-2)),                 &
          significant_digits = check_digits, ignore_digits = ignore     &
       )

     ! 上端境界チェック
     if( BCond(ll) == 'F' ) then
        xvr_True = xvr_wq(wq_RadDRad_wq(wq_Torvel)- wq_Torvel)/xvr_Rad
     else
        xvr_True = xvr_TorVel
     endif

     call AssertEqual(&
          message='wq_TorlBoundary (Top boundary)',                     &
          answer = xvr_True(:,:,kmi),                                   &
          check = xvr_Zero(:,:,kmi),                                    &
          significant_digits = check_digits, ignore_digits = ignore     &
          )

     call MessageNotify('M','wq_module_supack_torvelbc_test', &
          'wq_TorBoundarry: '//BCond(ll)//'-Top B.C. test succeeded!')


     !---- wr_TorBoundaryGrid test
     call wr_TorBoundaryGrid(wr_TorVel,cond=BCond(ll),new=.true.)
     xvr_TorVel = xvr_wr(wr_TorVel)
     wq_TorVel = wq_wr(wr_TorVel)

     ! 内部チェック
     xvr_True = xvr_TorVel - xvr_TorVel0

     call AssertEqual(&
          message='wq_TorBoundary (internal value)',                    &
          answer = xvr_Torvel(:,:,1:kmi-1),                             &
          check = xvr_Torvel0(:,:,1:kmi-1),                             &
          significant_digits = check_digits, ignore_digits = ignore     &
       )

     ! 上端境界チェック
     if( BCond(ll) == 'F' ) then
        xvr_True = xvr_wq(wq_RadDRad_wq(wq_Torvel)- wq_Torvel)/xvr_Rad
     else
        xvr_True = xvr_TorVel
     endif

     call AssertEqual(&
          message='wq_TorlBoundaryGrid (Top Booundary)',                &
          answer = xvr_True(:,:,kmi),                                   &
          check = xvr_Zero(:,:,kmi),                                    &
          significant_digits = check_digits, ignore_digits = ignore     &
          )

     call MessageNotify('M','wtq_module_supack_torvelbc_test', &
          'wtq_TorBoundaryGrid: '//BCond(ll)//'-Top B.C. test succeeded!')
          
  enddo

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

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

  call MPI_FINALIZE(IERR)

contains
  
  function wa_cleanup_wa(wa_data)
    real(8), intent(IN) :: wa_data(:,:)
    real(8)             :: wa_cleanup_wa(nc,size(wa_data,2))

    integer l,nmary(2)

    wa_cleanup_wa = wa_data

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

  end function wa_cleanup_wa


end program wtq_mpi_module_supack_torvelbc
