!----------------------------------------------------------------------
!     Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wq_mpi_module_supack テストプログラム
!
!      wq_KxRGrad_wq, xyr_KGrad_wq, wq_QOperator_wq のテスト
!      wr_RadRot_xyr_xyr, wr_RadRotRot_xyr_xyr_xyr のテスト
!      wq_Potential2Vector, wq_Potential2Rotation のテスト
!
!履歴  2012/04/02  竹広真一  wq_mpi_module_deriv_test2.f90 を supack 用に改造
!
program wq_mpi_module_supack_deriv2

  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=1.5               ! 球半径

  real(8), allocatable :: wq_Torvel(:,:)  ! トロイダルポテンシャル
  real(8), allocatable :: wq_Polvel(:,:)  ! ポロイダルポテンシャル

  real(8), allocatable :: xvr_Psi(:,:,:)
  real(8), allocatable :: xvr_KxRGrad(:,:,:)
  real(8), allocatable :: xvr_KGrad(:,:,:)
  real(8), allocatable :: xvr_QOperator(:,:,:)

  real(8), allocatable :: xvr_VLon(:,:,:)
  real(8), allocatable :: xvr_VLat(:,:,:)
  real(8), allocatable :: xvr_VRad(:,:,:)
  real(8), allocatable :: xvr_V0Lon(:,:,:) ! 速度(正解, 経度)
  real(8), allocatable :: xvr_V0Lat(:,:,:)  ! 速度(正解, 緯度)
  real(8), allocatable :: xvr_V0rad(:,:,:)  ! 速度(正解,動径)
  real(8), allocatable :: xvr_RadRot(:,:,:)
  real(8), allocatable :: xvr_RadRotRot(:,:,:)
  real(8), allocatable :: xvr_Torvel(:,:,:)
  real(8), allocatable :: xvr_Polvel(:,:,:)

  integer, parameter :: n=3

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

  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_derivative_test2', &
       'wq_mpi_module_supack derivative function test #1')

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

  allocate(wq_Torvel(nc,0:lm))
  allocate(wq_Polvel(nc,0:lm))

  allocate(xvr_Psi(0:im-1,1:jc,km))
  allocate(xvr_KxRGrad(0:im-1,1:jc,km))
  allocate(xvr_KGrad(0:im-1,1:jc,km))
  allocate(xvr_QOperator(0:im-1,1:jc,km))

  allocate(xvr_VLon(0:im-1,1:jc,km))
  allocate(xvr_VLat(0:im-1,1:jc,km))
  allocate(xvr_VRad(0:im-1,1:jc,km))
  allocate(xvr_V0Lon(0:im-1,1:jc,km)) ! 速度(正解, 経度)
  allocate(xvr_V0Lat(0:im-1,1:jc,km))  ! 速度(正解, 緯度)
  allocate(xvr_V0rad(0:im-1,1:jc,km))  ! 速度(正解,動径)
  allocate(xvr_RadRot(0:im-1,1:jc,km))
  allocate(xvr_RadRotRot(0:im-1,1:jc,km))
  allocate(xvr_Torvel(0:im-1,1:jc,km))
  allocate(xvr_Polvel(0:im-1,1:jc,km))

  !========== wq_KxRGrad_wq, xyr_KGrad_wq, wq_QOperator_wq =========
  call MessageNotify('M','wq_mpi_module_supack_derivative_test2', &
       'Test for wq_KxRGrad_wq, xyr_KGrad_wq, wq_QOperator_wq')

  xvr_Psi = xvr_rad**n * cos(xvr_lat)*sin(xvr_lon)       ! r**n Y_1^1
  xvr_KxRGrad = xvr_Rad**n * cos(xvr_lat)*cos(xvr_lon)   ! r**n Y_1^-1

  ! k ・▽ r**n Y_1^1 = (n-1)*r**(n-1)* Y_2^1
  xvr_KGrad = (n-1)*xvr_rad**(n-1)* cos(xvr_lat)*sin(xvr_lat)*sin(xvr_lon) 

  ! Q r**n Y_1^1 = -3*(n-1)*r**(n-1)* Y_2^1
  xvr_QOperator = - 3*(n-1)*xvr_rad**(n-1)* &
                   cos(xvr_lat)*sin(xvr_lat)*sin(xvr_lon) 

  call AssertEqual(&
       message='k x r grad with r**n Y_1^1',                &
       answer = xvr_KxRgrad,                                         &
       check = xvr_wq(wq_KxRGrad_wq(wq_xvr(xvr_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='k grad with r**n Y_1^1',                    &
       answer = xvr_KGrad,                                           &
       check = xvr_KGrad_wq(wq_xvr(xvr_Psi)),                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Q operator with r**n Y_1^1',                &
       answer = xvr_QOperator,                                       &
       check = xvr_wq(wq_QOperatorMPI_wq(wq_xvr(xvr_Psi))),          &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvr_Psi = 3*cos(xvr_lat)*sin(xvr_lat) * sin(xvr_lon) * xvr_Rad**2 ! Y_2^1
  xvr_KxRGrad = 3*cos(xvr_lat)*sin(xvr_lat) * cos(xvr_lon) * xvr_Rad**2 ! Y_2^-1

  ! k・▽r^2 Y_2^1 = 3r Y_1^1 
  xvr_KGrad = 3*xvr_Rad*cos(2*xvr_Lat)*cos(xvr_Lat)*sin(xvr_Lon) &
             +3*xvr_Rad*sin(2*xvr_Lat)*sin(xvr_Lat)*sin(xvr_Lon)

  xvr_QOperator = -9*cos(xvr_lat)*sin(xvr_lon)*xvr_rad


  call AssertEqual(&
       message='k x r grad with r**2 Y_2^1',                &
       answer = xvr_KxRgrad,                                         &
       check = xvr_wq(wq_KxRGrad_wq(wq_xvr(xvr_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='k grad with r**2 Y_2^1',                    &
       answer = xvr_KGrad,                                           &
       check = xvr_KGrad_wq(wq_xvr(xvr_Psi)),                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Q operator with r**2 Y_2^1',                &
       answer = xvr_QOperator,                                       &
       check = xvr_wq(wq_QOperatorMPI_wq(wq_xvr(xvr_Psi))),          &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !=========== wr_RadRot_xvr_xvr, wr_RadRotRot_xvr_xvr_xvr ============
  call MessageNotify('M','wq_module_supack_derivative_test2', &
       'Test for wr_RadRot_xvr_xvr, wr_RadRotRot_xvr_xvr_xvr')

  xvr_VLon = 0.0D0 ; xvr_VLat = 0.0D0 ; xvr_VRad = xvr_Rad
  xvr_RadRot = 0.0D0 ; xvr_RadRotRot = 0.0D0

  call AssertEqual(&
       message='r Rot v with simple example (v_r=r)',       &
       answer = xvr_RadRot,                                          &
       check = xvr_wr(wr_RadRot_xvr_xvr(xvr_VLon,xvr_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='r Rot Rot v with simple example (v_r=r)',   &
       answer = xvr_RadRotRot,                                       &
       check = xvr_wr(wr_RadRotRot_xvr_xvr_xvr(xvr_VLon,xvr_VLat,xvr_VRad)),&
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  ! 剛体回転解非線形項(東西流れ)
  xvr_VLon = 0.0D0
  xvr_VLat = xvr_Rad*sin(xvr_Lat)*cos(xvr_Lat)
  xvr_VRad = -xvr_Rad*cos(xvr_Lat)**2
  xvr_RadRot = 0.0D0 ; xvr_RadRotRot = 0.0D0

  call AssertEqual(&
       message='r Rot v with rigid rotation(E-W flow)',     &
       answer = xvr_RadRot,                                          &
       check = xvr_wr(wr_RadRot_xvr_xvr(xvr_VLon,xvr_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='r Rot Rot v with rigid rotation(E-W flow)', &
       answer = xvr_RadRotRot,                                       &
       check = xvr_wr(wr_RadRotRot_xvr_xvr_xvr(xvr_VLon,xvr_VLat,xvr_VRad)),&
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! 剛体回転解非線形項(南北流れ)
  xvr_VLon = xvr_Rad*cos(xvr_Lat)*sin(xvr_Lon)*cos(xvr_Lon)
  xvr_VLat = -xvr_Rad*sin(xvr_Lat)*cos(xvr_Lat)*sin(xvr_Lon)**2
  xvr_VRad = -xvr_Rad*(sin(xvr_Lat)**2*sin(xvr_Lon)**2 + cos(xvr_Lon)**2)
  xvr_RadRot = 0.0D0 ; xvr_RadRotRot = 0.0D0

  call AssertEqual(&
       message='r Rot v with rigid rotation(N-S flow)',     &
       answer = xvr_RadRot,                                          &
       check = xvr_wr(wr_RadRot_xvr_xvr(xvr_VLon,xvr_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='r Rot Rot v with rigid rotation(N-S flow)', &
       answer = xvr_RadRotRot,                                       &
       check = xvr_wr(wr_RadRotRot_xvr_xvr_xvr(xvr_VLon,xvr_VLat,xvr_VRad)),&
       significant_digits = check_digits, ignore_digits = ignore     &
       )


! 鉛直渦度を伴うベクトル場
  xvr_Polvel = xvr_Rad**2 * cos(xvr_Lat)*sin(xvr_Lon)   ! r**2 P_1^1
  !xvr_Psi = xvr_Rad**3 * cos(xvr_Lat)*sin(xvr_Lat)*sin(xvr_Lon)   ! r**3 P_2^1

  xvr_VLon =   xvr_GradLat_wq(wq_xvr(xvr_Polvel*xvr_Rad))
  xvr_VLat = - xvr_GradLon_wq(wq_xvr(xvr_Polvel*xvr_Rad))
  xvr_VRad = 0.0D0
  xvr_RadRot = 2 * xvr_Polvel                       ! r・▽×▽×(ψr) = L_2ψ
  !xvr_RadRot = 6 * xvr_Polvel                      ! r・▽×▽×(ψr) = L_2ψ
  xvr_RadRotRot = 0.0D0

  call AssertEqual(&
       message='r Rot v with vortical field',               &
       answer = xvr_RadRot,                                          &
       check = xvr_wr(wr_RadRot_xvr_xvr(xvr_VLon,xvr_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='r Rot Rot v with vortical field',           &
       answer = xvr_RadRotRot,                                       &
       check = xvr_wr(wr_RadRotRot_xvr_xvr_xvr(xvr_VLon,xvr_VLat,xvr_VRad)),&
       significant_digits = check_digits, ignore_digits = ignore     &
       )


! 鉛直速度を伴うベクトル場
  xvr_Torvel = 0.0D0
  xvr_Polvel = xvr_Rad**3 * cos(xvr_Lat)*sin(xvr_Lon)   ! r**2 P_1^1
  call wq_Potential2VectorMPI(&
         xvr_VLon,xvr_VLat,xvr_VRad, wq_xvr(xvr_Torvel),wq_xvr(xvr_Polvel) )

  xvr_RadRot = 0.0D0
  xvr_RadRotRot = -xvr_wq(wq_L2_wq(wq_Lapla_wq(wq_xvr(xvr_Polvel))))
                 ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ

  call AssertEqual(&
       message='r Rot v with non-vortical field',           &
       answer = xvr_RadRot,                                          &
       check = xvr_wr(wr_RadRot_xvr_xvr(xvr_VLon,xvr_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='r Rot Rot v with non-vortical field',       &
       answer = xvr_RadRotRot,                                       &
       check = xvr_wr(wr_RadRotRot_xvr_xvr_xvr(xvr_VLon,xvr_VLat,xvr_VRad)),&
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !==================== wq_Potential2Vector ====================
  call MessageNotify('M','wq_module_supack_derivative_test2', &
       'Test for wq_Potential2Vector, wq_Potential2VectorMPI')

  ! ----------------- 例 1 --------------------
  ! 剛体回転場
  wq_Torvel = wq_xvr(xvr_Rad * sin(xvr_Lat))
  wq_Polvel = 0.0D0

  xvr_V0lon = xvr_Rad * cos(xvr_Lat)
  xvr_V0lat = 0.0D0
  xvr_V0Rad = 0.0D0

  xvr_V0lon = xvr_Rad * cos(xvr_Lat)
  xvr_V0lat = 0.0D0
  xvr_V0Rad = 0.0D0

  call wq_Potential2VectorMPI(&
       xvr_VLon,xvr_VLat,xvr_VRad, wq_Torvel, wq_Polvel )

  call AssertEqual(&
       message='Potential2VectorMPI(VLon) with rigid rotation(EW)',&
       answer = xvr_VLon,                                            &
       check = xvr_V0Lon,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2VectorMPI(VLat) with rigid rotation(EW)',&
       answer = xvr_VLat,                                            &
       check = xvr_V0Lat,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2VectorMPI(VRad) with rigid rotation(EW)',&
       answer = xvr_VRad,                                            &
       check = xvr_V0Rad,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 例 2 --------------------
  ! 剛体回転場(南北流)
  wq_Torvel = wq_xvr(xvr_Rad * cos(xvr_Lat) * sin(xvr_Lon))
  wq_Polvel = 0.0D0

  xvr_V0Lon = -xvr_Rad*sin(xvr_Lat)*sin(xvr_Lon)
  xvr_V0Lat = -xvr_Rad*cos(xvr_Lon)
  xvr_V0Rad = 0.0D0

  xvr_V0Lon = -xvr_Rad*sin(xvr_Lat)*sin(xvr_Lon)
  xvr_V0Lat = -xvr_Rad*cos(xvr_Lon)
  xvr_V0Rad = 0.0D0

  call wq_Potential2VectorMPI(&
       xvr_VLon,xvr_VLat,xvr_VRad, wq_Torvel, wq_Polvel )

  call AssertEqual(&
       message='Potential2VectorMPI(VLon) with rigid rotation(NS)',&
       answer = xvr_VLon,                                            &
       check = xvr_V0Lon,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2VectorMPI(VLat) with rigid rotation(NS)',&
       answer = xvr_VLat,                                            &
       check = xvr_V0Lat,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2VectorMPI(VRad) with rigid rotation(NS)',&
       answer = xvr_VRad,                                            &
       check = xvr_V0Rad,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 例 3 --------------------
  ! 渦無し場

  wq_Torvel = 0.0D0
  wq_Polvel = wq_xvr(xvr_Rad**3 * sin(xvr_Lat))

  xvr_V0lon = 0.0D0
  xvr_V0lat = 4 * xvr_Rad**2 * cos(xvr_Lat)
  xvr_V0rad = 2 * xvr_Rad**2 * sin(xvr_Lat)

  xvr_V0lon = 0.0D0
  xvr_V0lat = 4 * xvr_Rad**2 * cos(xvr_Lat)
  xvr_V0rad = 2 * xvr_Rad**2 * sin(xvr_Lat)

  call wq_Potential2VectorMPI(&
       xvr_VLon,xvr_VLat,xvr_VRad, wq_Torvel, wq_Polvel )

  call AssertEqual(&
       message='Potential2VectorMPI(VLon) with no rotation',   &
       answer = xvr_VLon,                                            &
       check = xvr_V0Lon,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2VectorMPI(VLat) with no rotation',   &
       answer = xvr_VLat,                                            &
       check = xvr_V0Lat,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2VectorMPI(VRad) with no rotation',   &
       answer = xvr_VRad,                                            &
       check = xvr_V0Rad,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 例 4 --------------------
  ! ポロイダル速度場

  wq_Torvel = 0.0D0
  wq_Polvel = wq_xvr(xvr_Rad**5 * cos(xvr_Lat)*sin(xvr_Lon))

  xvr_V0lon = 6 * xvr_Rad**4 * cos(xvr_Lon)
  xvr_V0lat = - 6 * xvr_Rad**4 * sin(xvr_Lat) * sin(xvr_Lon)
  xvr_V0rad = 2 * xvr_Rad**4 * cos(xvr_Lat) * sin(xvr_Lon)

  call wq_Potential2VectorMPI(&
       xvr_VLon,xvr_VLat,xvr_VRad, wq_Torvel, wq_Polvel )

  call AssertEqual(&
       message='Potential2VectorMPI(VLat) with poloidal field',&
       answer = xvr_VLon,                                            &
       check = xvr_V0Lon,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2VectorMPI(VLat) with poloidal field',&
       answer = xvr_VLat,                                            &
       check = xvr_V0Lat,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2VectorMPI(VRad) with poloidal field',&
       answer = xvr_VRad,                                            &
       check = xvr_V0Rad,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !==================== wq_Potential2Rotation ====================
  call MessageNotify('M','wq_module_supack_derivative_test2', &
       'Test for wq_Potential2Rotation, wq_Potential2RotationMPI')

  ! ----------------- 例 1 --------------------
  ! 剛体回転場
  wq_Torvel = wq_xvr(xvr_Rad * sin(xvr_Lat))
  wq_Polvel = 0.0D0

  ! xvr_Vlon = xvr_Rad * cos(xvr_Lat)'
  ! xvr_Vlat = 0.0'
  ! xvr_VRad = 0.0'

  xvr_V0lon = 0.0D0
  xvr_V0lat = 2*cos(xvr_LAT)
  xvr_V0rad = 2*sin(xvr_LAT)

  call wq_Potential2RotationMPI(&
         xvr_VLon,xvr_VLat,xvr_VRad, wq_Torvel, wq_Polvel )

  call AssertEqual(&
       message='Potential2RotationMPI(VLon) with rigid rotation(EW)',&
       answer = xvr_VLon,                                            &
       check = xvr_V0Lon,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2RotationMPI(VLat) with rigid rotation(EW)',&
       answer = xvr_VLat,                                            &
       check = xvr_V0Lat,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2RotationMPI(VRad) with rigid rotation(EW)',&
       answer = xvr_VRad,                                            &
       check = xvr_V0Rad,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

! ----------------- 例 2 --------------------
  ! 剛体回転場(南北流)
  wq_Torvel = wq_xvr(xvr_Rad * cos(xvr_Lat) * sin(xvr_Lon))
  wq_Polvel = 0.0D0

  !  xvr_VLon = -xvr_Rad*sin(xvr_Lat)*sin(xvr_Lon)
  !  xvr_VLat = -xvr_Rad*cos(xvr_Lon)
  !  xvr_VRad = 0.0

  xvr_V0Lon = 2*cos(xvr_Lon)
  xvr_V0Lat = -2*sin(xvr_Lon)*sin(xvr_Lat)
  xvr_V0Rad = 2*sin(xvr_Lon)*cos(xvr_Lat)

  call wq_Potential2RotationMPI(&
       xvr_VLon,xvr_VLat,xvr_VRad, wq_Torvel, wq_Polvel )

  call AssertEqual(&
       message='Potential2RotationMPI(VLon) with rigid rotation(NS)',&
       answer = xvr_VLon,                                            &
       check = xvr_V0Lon,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2RotationMPI(VLat) with rigid rotation(NS)',&
       answer = xvr_VLat,                                            &
       check = xvr_V0Lat,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2RotationMPI(VRad) with rigid rotation(NS)',&
       answer = xvr_VRad,                                            &
       check = xvr_V0Rad,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 例 3 --------------------
  ! 渦無し場

  wq_Torvel = 0.0D0
  wq_Polvel = wq_xvr(xvr_Rad * sin(xvr_Lat))

  ! xvr_Vlon = 0
  ! xvr_Vlat = 2 * cos(xvr_Lat)
  ! xvr_Vrad = 2 * sin(xvr_Lat)

  xvr_V0lon = 0.0D0
  xvr_V0lat = 0.0D0
  xvr_V0rad = 0.0D0

  call wq_Potential2RotationMPI(&
       xvr_VLon,xvr_VLat,xvr_VRad, wq_Torvel, wq_Polvel )

  call AssertEqual(&
       message='Potential2RotationMPI(VLon) with no rotation',   &
       answer = xvr_VLon,                                            &
       check = xvr_V0Lon,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2RotationMPI(VLat) with no rotation',   &
       answer = xvr_VLat,                                            &
       check = xvr_V0Lat,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2RotationMPI(VRad) with no rotation',   &
       answer = xvr_VRad,                                            &
       check = xvr_V0Rad,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 例 4 --------------------
  ! ポロイダル速度場

  wq_Torvel = 0.0D0
  wq_Polvel = wq_xvr(xvr_Rad **3 * cos(xvr_Lat)*sin(xvr_Lon))

  ! xvr_Vlon = 4 * xvr_Rad**2 * cos(xvr_Lon)
  ! xvr_Vlat = -4 * xvr_Rad**2 * sin(xvr_Lat) * sin(xvr_Lon)
  ! xvr_Vrad = 2 * xvr_Rad**2  * cos(xvr_Lat) * sin(xvr_Lon)

  xvr_V0lon = 10*xvr_Rad*sin(xvr_Lat)*sin(xvr_Lon)
  xvr_V0lat = 10*xvr_Rad*cos(xvr_Lon)
  xvr_V0rad = 0.0D0

  call wq_Potential2RotationMPI(&
       xvr_VLon,xvr_VLat,xvr_VRad, wq_Torvel, wq_Polvel )

  call AssertEqual(&
       message='Potential2RotationMPI(VLon) with poloidal field',&
       answer = xvr_VLon,                                            &
       check = xvr_V0Lon,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2RotationMPI(VLat) with poloidal field',&
       answer = xvr_VLat,                                            &
       check = xvr_V0Lat,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='Potential2RotationMPI(VRad) with poloidal field',&
       answer = xvr_VRad,                                            &
       check = xvr_V0Rad,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  call MessageNotify('M','wq_mpi_module_supack_derivative_test2', &
       'wq_mpi_module_supack derivative function test #2 succeeded!')

  call MPI_FINALIZE(IERR)

end program wq_mpi_module_supack_deriv2
