!----------------------------------------------------------------------
!     Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ua_base_mpi_module_sypack テストプログラム :: ベクトル用サブルーチンのテスト
!
!履歴  2019/04/30  竹広真一
!
program ua_base_mpi_module_sypack_vector_test

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use ua_mpi_module

  use mpi
  implicit none

  integer, parameter :: im=128, jm=64, nm=42, km=6, npv=2

  real(8), allocatable ::  xvb_U(:,:,:)              ! 速度経度成分
  real(8), allocatable ::  xvb_V(:,:,:)              ! 速度緯度成分

  real(8), allocatable ::  xvb_UCosLat(:,:,:)        ! 速度経度成分
  real(8), allocatable ::  xvb_VCosLat(:,:,:)        ! 速度緯度成分

  real(8), allocatable ::  xvb_Uans(:,:,:)           ! 速度経度成分
  real(8), allocatable ::  xvb_Vans(:,:,:)           ! 速度緯度成分

  real(8), allocatable ::  xvb_Psi(:,:,:)            ! 流線関数
  real(8), allocatable ::  xvb_Chi(:,:,:)            ! 速度ポテンシャル

  real(8), allocatable ::  ua_Vor(:,:)               ! 渦度
  real(8), allocatable ::  ua_Div(:,:)               ! 発散

  real(8), allocatable ::  xvb_Vorans(:,:,:)         ! 渦度
  real(8), allocatable ::  xvb_Divans(:,:,:)         ! 発散


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

  integer :: ip, np, ierr
  integer :: k

  call MessageNotify('M','ua_base_mpi_module_vector_sypack_test', &
                         'ua_base_mpi_module_sypack subroutine tests')

 !---------------- MPI スタート ---------------------
  call MPI_INIT(IERR)
  call MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)

  call ua_mpi_Initial( nm, im, jm, km, npv=npv )

  allocate(xvb_U(0:im-1,jc,kc))
  allocate(xvb_V(0:im-1,jc,kc))
  allocate(xvb_UCosLat(0:im-1,jc,kc))
  allocate(xvb_VCosLat(0:im-1,jc,kc))
  allocate(xvb_Uans(0:im-1,jc,kc))
  allocate(xvb_Vans(0:im-1,jc,kc))
  allocate(xvb_Psi(0:im-1,jc,kc))
  allocate(xvb_Chi(0:im-1,jc,kc))
  allocate(xvb_Vorans(0:im-1,jc,kc))
  allocate(xvb_Divans(0:im-1,jc,kc))

  allocate(ua_Vor(nmc,km))
  allocate(ua_Div(nmc,km))

  xvb_Psi    = 0.0D0
  xvb_Chi    = 0.0D0
  xvb_Uans   = 0.0D0
  xvb_Vans   = 0.0D0
  xvb_Vorans = 0.0D0
  xvb_Divans = 0.0D0

  if ( kp_kl(1) /= 0 ) then
     xvb_Psi(:,:,kp_kl(1)) = cos(xv_Lat)*sin(xv_Lon)       ! Y_1^{-1}
     xvb_Chi(:,:,kp_kl(1)) = 0.0D0

     xvb_Uans(:,:,kp_kl(1)) = sin(xv_Lat)*sin(xv_Lon)
     xvb_Vans(:,:,kp_kl(1)) = cos(xv_Lon)

     xvb_Vorans(:,:,kp_kl(1)) = -2*cos(xv_Lat)*sin(xv_Lon)
     xvb_Divans(:,:,kp_kl(1)) = 0.0D0
  endif

  if ( kp_kl(2) /= 0 ) then
     xvb_Psi(:,:,kp_kl(2)) = 0.0D0
     xvb_Chi(:,:,kp_kl(2)) = cos(xv_Lat)*sin(xv_Lon)       ! Y_1^{-1}

     xvb_Uans(:,:,kp_kl(2)) = cos(xv_Lon)
     xvb_Vans(:,:,kp_kl(2)) = - sin(xv_Lat)*sin(xv_Lon)

     xvb_Vorans(:,:,kp_kl(2)) = 0.0D0
     xvb_Divans(:,:,kp_kl(2)) = -2*cos(xv_Lat)*sin(xv_Lon)
  endif

  if ( kp_kl(3) /= 0 ) then
     xvb_Psi(:,:,kp_kl(3)) = sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon)       ! Y_2^1
     xvb_Chi(:,:,kp_kl(3)) = 0.0D0

     xvb_Uans(:,:,kp_kl(3)) = - cos(2*xv_Lat)*cos(xv_Lon)
     xvb_Vans(:,:,kp_kl(3)) = - sin(xv_Lat)*sin(xv_Lon)

     xvb_Vorans(:,:,kp_kl(3)) = -6*sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon) ! Y_2^1
     xvb_Divans(:,:,kp_kl(3)) = 0.0D0
  endif

  if ( kp_kl(4) /= 0 ) then
     xvb_Psi(:,:,kp_kl(4)) = 0.0D0
     xvb_Chi(:,:,kp_kl(4)) = sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon)       ! Y_2^1

     xvb_Uans(:,:,kp_kl(4)) = - sin(xv_Lat)*sin(xv_Lon)
     xvb_Vans(:,:,kp_kl(4)) =   cos(2*xv_Lat)*cos(xv_Lon)

     xvb_Vorans(:,:,kp_kl(4)) = 0.0D0
     xvb_Divans(:,:,kp_kl(4)) = -6*sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon) ! Y_2^1
  endif

  if ( kp_kl(5) /= 0 ) then
     xvb_Psi(:,:,kp_kl(5)) = sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon)       ! Y_2^1
     xvb_Chi(:,:,kp_kl(5)) = cos(xv_Lat)*sin(xv_Lon)                     ! Y_1^{-1}

     xvb_Uans(:,:,kp_kl(5)) = - cos(2*xv_Lat)*cos(xv_Lon) + cos(xv_Lon)
     xvb_Vans(:,:,kp_kl(5)) = - sin(xv_Lat)*sin(xv_Lon) - sin(xv_Lat)*sin(xv_Lon)

     xvb_Vorans(:,:,kp_kl(5)) = -6*sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon) ! Y_2^1
     xvb_Divans(:,:,kp_kl(5)) = -2*cos(xv_Lat)*sin(xv_Lon)            ! Y_1^{-1}
  endif

  if ( kp_kl(6) /= 0 ) then
     xvb_Psi(:,:,kp_kl(6)) = cos(xv_Lat)*sin(xv_Lon)                ! Y_1^{-1}
     xvb_Chi(:,:,kp_kl(6)) = sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon)  ! Y_2^1

     xvb_Uans(:,:,kp_kl(6)) = 0.0D0
     xvb_Vans(:,:,kp_kl(6)) = cos(2*xv_Lat)*cos(xv_Lon) + cos(xv_Lon)

     xvb_Vorans(:,:,kp_kl(6)) = -2*cos(xv_Lat)*sin(xv_Lon)            ! Y_1^{-1}
     xvb_Divans(:,:,kp_kl(6)) = -6*sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon) ! Y_2^1
  endif

  call ua_StreamPotential2VectorMPI &
       ( ua_xvb(xvb_Psi), ua_xvb(xvb_Chi), xvb_U, xvb_V )

  call AssertEqual(&
    message='Test of ua_StreamPotential2VectorMPI(U)',            &
    answer = xvb_Uans,                                            &
    check  = xvb_U,                                               &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  call AssertEqual(&
    message='Test of ua_StreamPotential2VectorMPI(V)',            &
    answer = xvb_Vans,                                            &
    check  = xvb_V,                                               &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call ua_Vector2VorDivMPI( xvb_U, xvb_V, ua_Vor, ua_Div )

  call AssertEqual(&
    message='Test of ua_Vector2VorDivMPI(Vor)',                   &
    answer = xvb_Vorans,                                          &
    check  = xvb_ua(ua_Vor),                                      &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  call AssertEqual(&
    message='Test of ua_Vector2VorDivMPI(Div)',                   &
    answer = xvb_Divans,                                          &
    check  = xvb_ua(ua_Div),                                      &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  do k=1,kc
     xvb_UCosLat(:,:,k) = xvb_U(:,:,k) * cos(xv_Lat)
     xvb_VCosLat(:,:,k) = xvb_V(:,:,k) * cos(xv_Lat)
  enddo

  call ua_VectorCosLat2VorDivMPI( xvb_UCosLat, xvb_VCosLat, ua_Vor, ua_Div )

  call AssertEqual(&
    message='Test of ua_VectorCosLat2VorDivMPI(Vor)',             &
    answer = xvb_Vorans,                                          &
    check  = xvb_ua(ua_Vor),                                      &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  call AssertEqual(&
    message='Test of ua_VectorCosLat2VorDivMPI(Div)',             &
    answer = xvb_Divans,                                          &
    check  = xvb_ua(ua_Div),                                      &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call MessageNotify('M','ua_base_mpi_module_vector sypack_test', &
                         'ua_base_mpi_module_sypack subroutine tests succeeded!')

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

  call MPI_FINALIZE(IERR)

end program ua_base_mpi_module_sypack_vector_test
