!----------------------------------------------------------------------
!     Copyright (c) 2012 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wq_mpi_module_sjpack テストプログラム
!
!      wr_RotRad_wq のテスト (1/r d/dr r)
!      wr_DivRad_wq のテスト (1/r^2 d/dr r^2)
!      wq_Rad2Inv_wq のテスト (1/r^2)
!      wq_Lapla_wq のテスト (1/r^2 d/dr r^2)
!      xyr_GradLon_wq のテスト
!      xyr_GradLat_wq のテスト
!      wq_Div_xyr_xyr_xyr のテスト
!      wr_DivLon_xyr のテスト (1/r cosφ d/dλ)
!      wr_DivLat_xyr のテスト (1/r cosφ d/dφ cosφ)
!      wr_DivRad_xyr のテスト (1/r^2 d/dr r^2)
!      xyr_RotLon_wq_wq のテスト (1/r ∂Vrad/∂φ-1/r ∂(r Vlat)/∂r)
!      xyr_RotLat_wq_wq のテスト (1/r ∂(r Vlon)/∂r - 1/rcosφ・∂Vrad/∂λ)
!      wr_RotRad_xyr_xyr のテスト (1/rcosφ・∂Vlat/∂λ - 1/rcosφ・∂(Vlon cosφ)/∂φ)
!
!      xvr_GradLon_wq のテスト
!      xvr_GradLat_wq のテスト
!      wq_Div_xvr_xvr_xvr のテスト
!      wr_DivLon_xvr のテスト (1/r cosφ d/dλ)
!      wr_DivLat_xvr のテスト (1/r cosφ d/dφ cosφ)
!      wr_DivRad_xvr のテスト (1/r^2 d/dr r^2)
!      xvr_RotLon_wq_wq のテスト (1/r ∂Vrad/∂φ-1/r ∂(r Vlat)/∂r)
!      xvr_RotLat_wq_wq のテスト (1/r ∂(r Vlon)/∂r - 1/rcosφ・∂Vrad/∂λ)
!      wr_RotRad_xvr_xvr のテスト (1/rcosφ・∂Vlat/∂λ - 1/rcosφ・∂(Vlon cosφ)/∂φ)
!
!履歴  2012/04/02  竹広真一  wq_mpi_module_deriv_test1.f90 を sjpack 用に改造
!
program wq_mpi_module_sjpack_deriv1

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use wq_mpi_module_sjpack
  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), dimension(0:im-1,1:jm,km)     :: xyr_Data
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_Data1

  real(8), dimension(0:im-1,1:jm,km)     :: xyr_VLon
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_VLat
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_VRad
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_Psi
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_GradLon
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_GradLat
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_Div
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_DivLon
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_DivLat
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_DivRad
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_RotLon
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_RotLat
  real(8), dimension(0:im-1,1:jm,km)     :: xyr_RotRad

  real(8), allocatable :: xvr_Data(:,:,:)
  real(8), allocatable :: xvr_Data1(:,:,:)

  real(8), allocatable :: xvr_VLon(:,:,:)
  real(8), allocatable :: xvr_VLat(:,:,:)
  real(8), allocatable :: xvr_VRad(:,:,:)
  real(8), allocatable :: xvr_Psi(:,:,:)
  real(8), allocatable :: xvr_GradLon(:,:,:)
  real(8), allocatable :: xvr_GradLat(:,:,:)
  real(8), allocatable :: xvr_Div(:,:,:)
  real(8), allocatable :: xvr_DivLon(:,:,:)
  real(8), allocatable :: xvr_DivLat(:,:,:)
  real(8), allocatable :: xvr_DivRad(:,:,:)
  real(8), allocatable :: xvr_RotLon(:,:,:)
  real(8), allocatable :: xvr_RotLat(:,:,:)
  real(8), allocatable :: xvr_RotRad(:,:,:)

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

  integer, parameter :: nmin=1, nmax=10

  integer :: iproc, np, ierr
  integer :: n

 !---------------- 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_sjpack_derivative_test1', &
       'wq_mpi_module_sjpack derivative function test #1')

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

  allocate(xvr_Data(0:im-1,1:jc,km))
  allocate(xvr_Data1(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_Psi(0:im-1,1:jc,km))
  allocate(xvr_GradLon(0:im-1,1:jc,km))
  allocate(xvr_GradLat(0:im-1,1:jc,km))
  allocate(xvr_Div(0:im-1,1:jc,km))
  allocate(xvr_DivLon(0:im-1,1:jc,km))
  allocate(xvr_DivLat(0:im-1,1:jc,km))
  allocate(xvr_DivRad(0:im-1,1:jc,km))
  allocate(xvr_RotLon(0:im-1,1:jc,km))
  allocate(xvr_RotLat(0:im-1,1:jc,km))
  allocate(xvr_RotRad(0:im-1,1:jc,km))

  !==================== wr_RotRad_wq ====================

  write( 6,* ) 'Test for wr_RotDRad_wq (even mode)'
  do n=nmin+1,nmax,2
     write(6,*) 'n=',n
     xyr_Data = xyr_Rad**n
     xyr_Data1 = (n+1)*xyr_Rad**(n-1)

     xyr_Data = xyr_wr(wr_RotDRad_wq(wq_xyr(xyr_Data)))

     call AssertEqual(&
       message='wq_RotRad_wq with r^n (even n)',                     &
       answer = xyr_Data1,                                           &
       check = xyr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

     write(6,*) 'n=',n
     xvr_Data = xvr_Rad**n
     xvr_Data1 = (n+1)*xvr_Rad**(n-1)

     xvr_Data = xvr_wr(wr_RotDRad_wq(wq_xvr(xvr_Data)))

     call AssertEqual(&
       message='wq_RotRad_wq with r^n (even n)',                     &
       answer = xvr_Data1,                                           &
       check = xvr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  enddo

  do n=nmin,nmax,2
     xyr_Data = xyr_Rad**n*cos(xyr_Lat)*sin(xyr_Lon)
     xyr_Data1 = (n+1)*xyr_Rad**(n-1)*cos(xyr_Lat)*sin(xyr_Lon)
     xyr_Data = xyr_wr(wr_RotDRad_wq(wq_xyr(xyr_Data)))

     call AssertEqual(&
       message='wq_RotRad_wq with r^n Y_1^-1 (odd n)',               &
       answer = xyr_Data1,                                           &
       check = xyr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

     xvr_Data = xvr_Rad**n*cos(xvr_Lat)*sin(xvr_Lon)
     xvr_Data1 = (n+1)*xvr_Rad**(n-1)*cos(xvr_Lat)*sin(xvr_Lon)
     xvr_Data = xvr_wr(wr_RotDRad_wq(wq_xvr(xvr_Data)))

     call AssertEqual(&
       message='wq_RotRad_wq with r^n Y_1^-1 (odd n)',               &
       answer = xvr_Data1,                                           &
       check = xvr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  enddo

  !==================== wr_DivRad_wq ====================
  do n=nmin,nmax,2
     write(6,*) 'n=',n
     xyr_Data = xyr_Rad**n * cos(xyr_Lat)*sin(xyr_Lon)
     xyr_Data1 = (n+2)*xyr_Rad**(n-1)*cos(xyr_Lat)*sin(xyr_Lon)

     xyr_Data = xyr_wr(wr_DivRad_wq(wq_xyr(xyr_Data)))

     call AssertEqual(&
       message='wr_DivRad_wq with r^n Y_1^-1 (odd n)',               &
       answer = xyr_Data1,                                           &
       check = xyr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

     write(6,*) 'n=',n
     xvr_Data = xvr_Rad**n * cos(xvr_Lat)*sin(xvr_Lon)
     xvr_Data1 = (n+2)*xvr_Rad**(n-1)*cos(xvr_Lat)*sin(xvr_Lon)

     xvr_Data = xvr_wr(wr_DivRad_wq(wq_xvr(xvr_Data)))

     call AssertEqual(&
       message='wr_DivRad_wq with r^n Y_1^-1 (odd n)',               &
       answer = xvr_Data1,                                           &
       check = xvr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  enddo

  do n=nmin+1,nmax,2
     write(6,*) 'n=',n
     xyr_Data = xyr_Rad**n
     xyr_Data1 = (n+2)*xyr_Rad**(n-1)

     xyr_Data = xyr_wr(wr_DivRad_wq(wq_xyr(xyr_Data)))

     call AssertEqual(&
       message='wr_DivRad_wq with r^n (even n)',                     &
       answer = xyr_Data1,                                           &
       check = xyr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

     write(6,*) 'n=',n
     xvr_Data = xvr_Rad**n
     xvr_Data1 = (n+2)*xvr_Rad**(n-1)

     xvr_Data = xvr_wr(wr_DivRad_wq(wq_xvr(xvr_Data)))

     call AssertEqual(&
       message='wr_DivRad_wq with r^n (even n)',                     &
       answer = xvr_Data1,                                           &
       check = xvr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  enddo

  !==================== wq_Rad2Inv_wq ====================
  do n=nmin+2,nmax,2
     write(6,*) 'n=',n
     xyr_Data = xyr_Rad**n * cos(xyr_Lat)*sin(xyr_Lon)
     xyr_Data1 = xyr_Data/xyr_Rad**2

     xyr_Data = xyr_wq(wq_Rad2Inv_wq(wq_xyr(xyr_Data)))

     call AssertEqual(&
       message='wq_Rad2Inv_wq with r^n Y_1^-1 (odd n)',              &
       answer = xyr_Data1,                                           &
       check = xyr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

     write(6,*) 'n=',n
     xvr_Data = xvr_Rad**n * cos(xvr_Lat)*sin(xvr_Lon)
     xvr_Data1 = xvr_Data/xvr_Rad**2

     xvr_Data = xvr_wq(wq_Rad2Inv_wq(wq_xvr(xvr_Data)))

     call AssertEqual(&
       message='wq_Rad2Inv_wq with r^n Y_1^-1 (odd n)',              &
       answer = xvr_Data1,                                           &
       check = xvr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  enddo

  do n=nmin+1,nmax,2
     write(6,*) 'n=',n
     xyr_Data = xyr_Rad**n
     xyr_Data1 = xyr_Data/xyr_Rad**2

     xyr_Data = xyr_wq(wq_Rad2Inv_wq(wq_xyr(xyr_Data)))

     call AssertEqual(&
       message='wq_Rad2Inv_wq with r^n (even n)',                    &
       answer = xyr_Data1,                                           &
       check = xyr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

     write(6,*) 'n=',n
     xvr_Data = xvr_Rad**n
     xvr_Data1 = xvr_Data/xvr_Rad**2

     xvr_Data = xvr_wq(wq_Rad2Inv_wq(wq_xvr(xvr_Data)))

     call AssertEqual(&
       message='wq_Rad2Inv_wq with r^n (even n)',                    &
       answer = xvr_Data1,                                           &
       check = xvr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  enddo

  !==================== wq_Lapla_wq ====================
  do n=nmin,nmax,2
     write(6,*) 'n=',n
     xyr_Data = xyr_Rad**n * cos(xyr_Lat)*sin(xyr_Lon)
     xyr_Data1 = (n*(n+1)-2)*xyr_Data/xyr_Rad**2

     xyr_Data = xyr_wq(wq_Lapla_wq(wq_xyr(xyr_Data)))

     call AssertEqual(&
       message='wq_Lapla_wq with r^n Y_1^-1 (odd n)',                &
       answer = xyr_Data1,                                           &
       check = xyr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

     write(6,*) 'n=',n
     xvr_Data = xvr_Rad**n * cos(xvr_Lat)*sin(xvr_Lon)
     xvr_Data1 = (n*(n+1)-2)*xvr_Data/xvr_Rad**2

     xvr_Data = xvr_wq(wq_Lapla_wq(wq_xvr(xvr_Data)))

     call AssertEqual(&
       message='wq_Lapla_wq with r^n Y_1^-1 (odd n)',                &
       answer = xvr_Data1,                                           &
       check = xvr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  enddo

  write( 6,* ) 'Test for wq_Lapla_wq(even mode)'
  do n=nmin+1,nmax,2
     write(6,*) 'n=',n
     xyr_Data = xyr_Rad**n
     xyr_Data1 = n*(n+1)*xyr_Data/xyr_Rad**2

     xyr_Data = xyr_wq(wq_Lapla_wq(wq_xyr(xyr_Data)))

     call AssertEqual(&
       message='wq_Lapla_wq with r^n (even n)',                      &
       answer = xyr_Data1,                                           &
       check = xyr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

     write(6,*) 'n=',n
     xvr_Data = xvr_Rad**n
     xvr_Data1 = n*(n+1)*xvr_Data/xvr_Rad**2

     xvr_Data = xvr_wq(wq_Lapla_wq(wq_xvr(xvr_Data)))

     call AssertEqual(&
       message='wq_Lapla_wq with r^n (even n)',                      &
       answer = xvr_Data1,                                           &
       check = xvr_Data,                                             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  enddo

  !========== xyr_GradLon_wq, xyr_GradLat_wq, wr_Div_xyr_xyr_xyr ==========
  n=2

  ! ------ r^3 Y_1^1 Poloidal field -----
  xyr_VRad = 0.0D0
  xyr_Psi = xyr_Rad**(n+1) * cos(xyr_Lat)*sin(xyr_Lon)   ! r**3 P_1^1

  xyr_GradLon =  xyr_Rad**n * cos(xyr_Lon)
  xyr_GradLat = -xyr_Rad**n * sin(xyr_Lat)*sin(xyr_Lon)

  xyr_Div = - 2* xyr_Psi/xyr_Rad**2

  call AssertEqual(&
       message='xyr_GradLon_wq with r^3 Y_1^1',                      &
       answer = xyr_GradLon,                                         &
       check = xyr_GradLon_wq(wq_xyr(xyr_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='xyr_GradLat_wq with r^3 Y_1^1',                      &
       answer = xyr_GradLat,                                         &
       check = xyr_GradLat_wq(wq_xyr(xyr_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wr_Div_xyr_xyr_xyr with r^3 Y_1^1',                  &
       answer = xyr_Div,                                             &
       check = xyr_wr(wr_Div_xyr_xyr_xyr(xyr_GradLon,xyr_GradLat,xyr_Vrad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvr_VRad = 0.0D0
  xvr_Psi = xvr_Rad**(n+1) * cos(xvr_Lat)*sin(xvr_Lon)   ! r**3 P_1^1

  xvr_GradLon =  xvr_Rad**n * cos(xvr_Lon)
  xvr_GradLat = -xvr_Rad**n * sin(xvr_Lat)*sin(xvr_Lon)

  xvr_Div = - 2* xvr_Psi/xvr_Rad**2

  call AssertEqual(&
       message='xvr_GradLon_wq with r^3 Y_1^1',                      &
       answer = xvr_GradLon,                                         &
       check = xvr_GradLon_wq(wq_xvr(xvr_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='xvr_GradLat_wq with r^3 Y_1^1',                      &
       answer = xvr_GradLat,                                         &
       check = xvr_GradLat_wq(wq_xvr(xvr_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wr_Div_xvr_xvr_xvr with r^3 Y_1^1',                  &
       answer = xvr_Div,                                             &
       check = xvr_wr(wr_Div_xvr_xvr_xvr(xvr_GradLon,xvr_GradLat,xvr_Vrad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ------ r^2 Y_2^1 Poloidal field -----
  xyr_VRad = 0.0D0
  xyr_Psi = xyr_Rad**n * cos(xyr_Lat)*sin(xyr_Lat) * sin(xyr_Lon) ! P_2^1

  xyr_GradLon =  xyr_Rad**(n-1)*sin(xyr_Lat)*cos(xyr_Lon)
  xyr_GradLat =  xyr_Rad**(n-1)*cos(2*xyr_Lat)*sin(xyr_Lon)

  xyr_Div = - 6* xyr_Psi/xyr_Rad**2

  call AssertEqual(&
       message='xyr_GradLon_wq with r^2 Y_2^1',                      &
       answer = xyr_GradLon,                                         &
       check = xyr_GradLon_wq(wq_xyr(xyr_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='xyr_GradLat_wq with r^2 Y_2^1',                      &
       answer = xyr_GradLat,                                         &
       check = xyr_GradLat_wq(wq_xyr(xyr_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wr_Div_xyr_xyr_xyr with r^2 Y_2^1',                  &
       answer = xyr_Div,                                             &
       check = xyr_wr(wr_Div_xyr_xyr_xyr(xyr_GradLon,xyr_GradLat,xyr_Vrad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvr_VRad = 0.0D0
  xvr_Psi = xvr_Rad**n * cos(xvr_Lat)*sin(xvr_Lat) * sin(xvr_Lon) ! P_2^1

  xvr_GradLon =  xvr_Rad**(n-1)*sin(xvr_Lat)*cos(xvr_Lon)
  xvr_GradLat =  xvr_Rad**(n-1)*cos(2*xvr_Lat)*sin(xvr_Lon)

  xvr_Div = - 6* xvr_Psi/xvr_Rad**2

  call AssertEqual(&
       message='xvr_GradLon_wq with r^2 Y_2^1',                      &
       answer = xvr_GradLon,                                         &
       check = xvr_GradLon_wq(wq_xvr(xvr_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='xvr_GradLat_wq with r^2 Y_2^1',                      &
       answer = xvr_GradLat,                                         &
       check = xvr_GradLat_wq(wq_xvr(xvr_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wr_Div_xvr_xvr_xvr with r^2 Y_2^1',                  &
       answer = xvr_Div,                                             &
       check = xvr_wr(wr_Div_xvr_xvr_xvr(xvr_GradLon,xvr_GradLat,xvr_Vrad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !========== wr_DivLon_xyr, wr_DivLat_xyr, wr_DivRad_xyr ==========

  xyr_VLon   = xyr_rad**n * cos(xyr_Lat)**2*sin(xyr_Lon)
  xyr_DivLon = xyr_rad**(n-1)*cos(xyr_Lat)*cos(xyr_Lon)

  xyr_VLat   = xyr_rad**n * cos(xyr_Lat)**2*sin(xyr_Lon)
  xyr_DivLat = -3*xyr_rad**(n-1)*cos(xyr_Lat)*sin(xyr_Lat)*sin(xyr_Lon)

  xyr_VRad = xyr_Rad**n * cos(xyr_Lat)*sin(xyr_Lon)
  xyr_DivRad = (n+2)*xyr_Rad**(n-1)*cos(xyr_Lat)*sin(xyr_Lon)

  call AssertEqual(&
       message='wr_DivLon_xyr',                                      &
       answer = xyr_DivLon,                                          &
       check = xyr_wr(wr_DivLon_xyr(xyr_VLon)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wr_DivLat_xyr',                                      &
       answer = xyr_DivLat,                                          &
       check = xyr_wr(wr_DivLat_xyr(xyr_VLat)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wr_DivRad_xyr',                                      &
       answer = xyr_DivRad,                                          &
       check = xyr_wr(wr_DivRad_xyr(xyr_VRad)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvr_VLon   = xvr_rad**n * cos(xvr_Lat)**2*sin(xvr_Lon)
  xvr_DivLon = xvr_rad**(n-1)*cos(xvr_Lat)*cos(xvr_Lon)

  xvr_VLat   = xvr_rad**n * cos(xvr_Lat)**2*sin(xvr_Lon)
  xvr_DivLat = -3*xvr_rad**(n-1)*cos(xvr_Lat)*sin(xvr_Lat)*sin(xvr_Lon)

  xvr_VRad = xvr_Rad**n * cos(xvr_Lat)*sin(xvr_Lon)
  xvr_DivRad = (n+2)*xvr_Rad**(n-1)*cos(xvr_Lat)*sin(xvr_Lon)

  call AssertEqual(&
       message='wr_DivLon_xvr',                                      &
       answer = xvr_DivLon,                                          &
       check = xvr_wr(wr_DivLon_xvr(xvr_VLon)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wr_DivLat_xyr',                                      &
       answer = xvr_DivLat,                                          &
       check = xvr_wr(wr_DivLat_xvr(xvr_VLat)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wr_DivRad_xvr',                                      &
       answer = xvr_DivRad,                                          &
       check = xvr_wr(wr_DivRad_xvr(xvr_VRad)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !========== xyr_RotLon_wq_wq, xyr_DivLat_wq_wq, wr_DivRad_xyr_xyr ==========
  n=3

  xyr_VLon = xyr_rad**n * cos(xyr_Lat)*sin(xyr_Lon)
  xyr_VLat = xyr_rad**n * cos(xyr_Lat)*sin(xyr_Lon)
  xyr_VRad = xyr_Rad**n * cos(xyr_Lat)*sin(xyr_Lon)

  xyr_RotLon = - xyr_rad**(n-1) * sin(xyr_Lat) * sin(xyr_Lon)  &
               - (n+1)*xyr_rad**(n-1)*cos(xyr_Lat)*sin(xyr_Lon)

  xyr_RotLat =   (n+1)*xyr_rad**(n-1)*cos(xyr_Lat)*sin(xyr_Lon) &
               - xyr_rad**(n-1) * cos(xyr_Lon) 

  call AssertEqual(&
       message='xyr_RotLon_wq_wq',                                   &
       answer = xyr_RotLon,                                          &
       check = xyr_RotLon_wq_wq(wq_xyr(xyr_VRad),wq_xyr(xyr_VLat)),  &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='xyr_RotLat_wq_wq',                                   &
       answer = xyr_RotLat,                                          &
       check = xyr_RotLat_wq_wq(wq_xyr(xyr_VLon),wq_xyr(xyr_VRad)),  &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyr_VLon = xyr_rad**n * cos(xyr_Lat) 
  xyr_VLat = xyr_rad**n * cos(xyr_Lat)**2 * sin(xyr_Lon)

  xyr_RotRad =   xyr_rad**(n-1) * cos(xyr_Lat) * cos(xyr_Lon) &
               + 2 * xyr_rad**(n-1) * sin(xyr_Lat) 

  call AssertEqual(&
       message='wr_RotRad_xyr_xyr',                                  &
       answer = xyr_RotRad,                                          &
       check = xyr_wr(wr_RotRad_xyr_xyr(xyr_VLat,xyr_VLon)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvr_VLon = xvr_rad**n * cos(xvr_Lat)*sin(xvr_Lon)
  xvr_VLat = xvr_rad**n * cos(xvr_Lat)*sin(xvr_Lon)
  xvr_VRad = xvr_Rad**n * cos(xvr_Lat)*sin(xvr_Lon)

  xvr_RotLon = - xvr_rad**(n-1) * sin(xvr_Lat) * sin(xvr_Lon)  &
               - (n+1)*xvr_rad**(n-1)*cos(xvr_Lat)*sin(xvr_Lon)

  xvr_RotLat =   (n+1)*xvr_rad**(n-1)*cos(xvr_Lat)*sin(xvr_Lon) &
               - xvr_rad**(n-1) * cos(xvr_Lon) 

  call AssertEqual(&
       message='xvr_RotLon_wq_wq',                                   &
       answer = xvr_RotLon,                                          &
       check = xvr_RotLon_wq_wq(wq_xvr(xvr_VRad),wq_xvr(xvr_VLat)),  &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='xvr_RotLat_wq_wq',                                   &
       answer = xvr_RotLat,                                          &
       check = xvr_RotLat_wq_wq(wq_xvr(xvr_VLon),wq_xvr(xvr_VRad)),  &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvr_VLon = xvr_rad**n * cos(xvr_Lat) 
  xvr_VLat = xvr_rad**n * cos(xvr_Lat)**2 * sin(xvr_Lon)

  xvr_RotRad =   xvr_rad**(n-1) * cos(xvr_Lat) * cos(xvr_Lon) &
               + 2 * xvr_rad**(n-1) * sin(xvr_Lat) 

  call AssertEqual(&
       message='wr_RotRad_xvr_xvr',                                  &
       answer = xvr_RotRad,                                          &
       check = xvr_wr(wr_RotRad_xvr_xvr(xvr_VLat,xvr_VLon)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','wq_mpi_module_sjpack_derivative_test1', &
       'wq_mpi_module_sjpack derivative function test #1 succeeded!')

  call MPI_FINALIZE(IERR)

end program wq_mpi_module_sjpack_deriv1
