!----------------------------------------------------------------------
!     Copyright (c) 2019 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ut_mpi_module テストプログラム
!
!      トロイダルポロイダル関係微分関数のテスト
!
!履歴  2018/08/09  竹広真一
!
program ut_mpi_module_derivative_test2

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use ut_mpi_module
  use mpi
  implicit none

  integer,parameter  :: im=32, jm=16, km=16  ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: nm=10, lm=16         ! 切断波数の設定(水平, 動径)
  real(8),parameter  :: ri=0.5, ro=1.5       ! 内外半径

  integer :: npv=2

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

  real(8), allocatable :: ut_Psi(:,:)
  real(8), allocatable :: ut_Phi(:,:)

  real(8), allocatable :: xvh_DData(:,:,:)
  real(8), allocatable :: xvh_Psi(:,:,:)

  real(8), allocatable :: xvh_VLon(:,:,:)
  real(8), allocatable :: xvh_VLat(:,:,:)
  real(8), allocatable :: xvh_VRad(:,:,:)
  real(8), allocatable :: xvh_RadRot(:,:,:)
  real(8), allocatable :: xvh_RadRotRot(:,:,:)

  integer, parameter :: n=2
  integer :: iproc, np, ierr

  !---------------- MPI スタート ---------------------
  !call get_ENV_INT('NPV',npv)

  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','ut_mpi_module_derivative_test2', &
       'ut_mpi_module derivative function test #2')

  call ut_mpi_initial(im,jm,km,nm,lm,ri,ro,npv=npv)

  allocate(ut_Psi(nmc,0:lm))
  allocate(ut_Phi(nmc,0:lm))

  allocate(xvh_DData(0:im-1,1:jc,kc))
  allocate(xvh_Psi(0:im-1,1:jc,kc))

  allocate(xvh_VLon(0:im-1,1:jc,kc))
  allocate(xvh_VLat(0:im-1,1:jc,kc))
  allocate(xvh_VRad(0:im-1,1:jc,kc))
  allocate(xvh_RadRot(0:im-1,1:jc,kc))
  allocate(xvh_RadRotRot(0:im-1,1:jc,kc))

!!$  !========== wt_KxRGrad_wt, xvz_KGrad_wt, wt_QOperator_wt ==========
!!$
!!$  ! ----------------- 例 1 --------------------
!!$  xvz_Psi = xvz_rad**n * cos(xvz_lat)*sin(xvz_lon)   ! r**2 Y_1^-1
!!$
!!$  xvz_DData = xvz_rad**n * cos(xvz_lat)*cos(xvz_lon)
!!$  call AssertEqual(&
!!$       message='wt_KxRGrad_wt with r^2 Y_1^-1',                      &
!!$       answer = xvz_DData,                                           &
!!$       check = xvz_wt(wt_KxRGrad_wt(wt_xvz(xvz_Psi))),               &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
!!$  ! k ・▽ r**n Y_1^1 = (n-1)*r**(n-1)* Y_2^1
!!$  xvz_DData = (n-1)*xvz_rad**(n-1)* cos(xvz_lat)*sin(xvz_lat)*sin(xvz_lon)
!!$  call AssertEqual(&
!!$       message='xvz_KGrad_wt with r^2 Y_1^-1',                       &
!!$       answer = xvz_DData,                                           &
!!$       check = xvz_KGrad_wt(wt_xvz(xvz_Psi)),                        &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
!!$  ! Q r**n Y_1^1 = -3*(n-1)*r**(n-1)* Y_2^1
!!$  xvz_DData = - 3*(n-1)*xvz_rad**(n-1)* cos(xvz_lat)*sin(xvz_lat)*sin(xvz_lon)
!!$  call AssertEqual(&
!!$       message='wt_QOperator_wt r^2 Y_1^-1',                         &
!!$       answer = xvz_DData,                                           &
!!$       check = xvz_wt(wt_QOperatorMPI_wt(wt_xvz(xvz_Psi))),          &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
!!$! ----------------- 例 2 --------------------
!!$  xvz_Psi = cos(xvz_lat)*sin(xvz_lat) * sin(xvz_lon) ! Y_2^1
!!$
!!$  xvz_DData = cos(xvz_lat)*sin(xvz_lat) * cos(xvz_lon) ! Y_2^1
!!$  call AssertEqual(&
!!$       message='wt_KxRGrad_wt with Y_2^1',                           &
!!$       answer = xvz_DData,                                           &
!!$       check = xvz_wt(wt_KxRGrad_wt(wt_xvz(xvz_Psi))),               &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
!!$  ! k・▽ Y_2^1 = (- 4 Y_3^1 / 15 - Y_1^1 /5)
!!$  xvz_DData = cos(2*xvz_lat)*cos(xvz_lat)*sin(xvz_lon)/xvz_rad
!!$  call AssertEqual(&
!!$       message='xvz_KGrad_wt with Y_2^1',                            &
!!$       answer = xvz_DData,                                           &
!!$       check = xvz_KGrad_wt(wt_xvz(xvz_Psi)),                        &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
!!$  xvz_DData = (16*sin(xvz_lat)**2 - 5)*cos(xvz_lat)*sin(xvz_lon)/xvz_rad
!!$  call AssertEqual(&
!!$       message='wt_QOperator_wt Y_2^1',                              &
!!$       answer = xvz_DData,                                           &
!!$       check = xvz_wt(wt_QOperatorMPI_wt(wt_xvz(xvz_Psi))),          &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
!!$
!!$
  !========== ut_RadRot_xvh_xvh, ut_RadRotRot_xvh_xvh_xvh ==========

  !---------- 単純流れ v_r=r ----------
  xvh_VLon = 0 ; xvh_VLat = 0 ; xvh_VRad = xvh_Rad
  xvh_RadRot = 0 ; xvh_RadRotRot = 0

  call AssertEqual(&
       message='ut_RadRot_xvh_xvh with v_r=r',                       &
       answer = xvh_RadRot,                                          &
       check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xvh_xvh with v_r=r',                    &
       answer = xvh_RadRotRot,                                       &
       check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !---------- 剛体回転解非線形項(東西流れ) ----------
  xvh_VLon = 0
  xvh_VLat = xvh_Rad*sin(xvh_Lat)*cos(xvh_Lat)
  xvh_VRad = -xvh_Rad*cos(xvh_Lat)**2

  xvh_RadRot = 0 ; xvh_RadRotRot = 0

  call AssertEqual(&
       message='ut_RadRot_xvh_xvh with rigid rotation(E-W)',         &
       answer = xvh_RadRot,                                          &
       check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xvh_xvh with rigid rotation(E-W)',      &
       answer = xvh_RadRotRot,                                       &
       check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !---------- 剛体回転解非線形項(南北流れ) ----------
  xvh_VLon = xvh_Rad*cos(xvh_Lat)*sin(xvh_Lon)*cos(xvh_Lon)
  xvh_VLat = -xvh_Rad*sin(xvh_Lat)*cos(xvh_Lat)*sin(xvh_Lon)**2
  xvh_VRad = -xvh_Rad*(sin(xvh_Lat)**2*sin(xvh_Lon)**2 + cos(xvh_Lon)**2)
  xvh_RadRot = 0 ; xvh_RadRotRot = 0

  call AssertEqual(&
      message='ut_RadRot_xvh_xvh with rigid rotation(N-S)',          &
       answer = xvh_RadRot,                                          &
       check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xvh_xvh with rigid rotation(N-S)',      &
       answer = xvh_RadRotRot,                                       &
       check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !--------- 鉛直渦度を伴うベクトル場 ----------
  xvh_Psi = xvh_Rad**2 * cos(xvh_Lat)*sin(xvh_Lon)   ! r**2 P_1^1

  xvh_VLon =   xvh_gRadLat_ut(ut_xvh(xvh_Psi*xvh_Rad))
  xvh_VLat = - xvh_gRadLon_ut(ut_xvh(xvh_Psi*xvh_Rad))
  xvh_VRad = 0
  xvh_RadRot = 2 * xvh_Psi                           ! r・▽×▽×(ψr) = L_2ψ
  xvh_RadRotRot = 0

  call AssertEqual(&
      message='ut_RadRot_xvh_xvh with r^2 Y_1^-1 toroidal field',    &
       answer = xvh_RadRot,                                          &
       check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xvh_xvh with r^2 Y_1^-1 toroidal field',&
       answer = xvh_RadRotRot,                                       &
       check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !---------- 鉛直速度を伴うベクトル場 -----------
                 ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ
  xvh_VRad = xvh_ut(ut_l2_ut(ut_xvh(xvh_Psi/xvh_Rad)))
  xvh_VLat = xvh_gRadLat_ut(ut_dRad_ut(ut_xvh(xvh_Psi*xvh_Rad)))
  xvh_VLon = xvh_gRadLon_ut(ut_dRad_ut(ut_xvh(xvh_Psi*xvh_Rad)))

  xvh_RadRot = 0
  xvh_RadRotRot = -xvh_ut(ut_l2_ut(ut_lapla_ut(ut_xvh(xvh_Psi))))
                 ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ

  call AssertEqual(&
      message='ut_RadRot_xvh_xvh with non-vortical field',           &
       answer = xvh_RadRot,                                          &
       check = xvh_ut(ut_RadRot_xvh_xvh(xvh_VLon,xvh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xvh_xvh with r^2 Y_1^-1 toroidal field',&
       answer = xvh_RadRotRot,                                       &
       check = xvh_ut(ut_RadRotRot_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !========== ut_Potential2Vector, ut_Potential2Rotation ==========

  !----------------- 剛体回転場 --------------------
  !
  ut_Psi = ut_xvh(xvh_Rad * sin(xvh_Lat))
  ut_Phi = 0.0D0

  call ut_Potential2VectorMPI(&
       xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi )

  xvh_DData = xvh_Rad * cos(xvh_Lat)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lon) with rigid rotation field(E-W)', &
       answer = xvh_VLon,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 0.0D0
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lat) with rigid rotation field(E-W)', &
       answer = xvh_VLat,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 0.0D0
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Rad) with rigid rotation field(E-W)', &
       answer = xvh_VRad,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call ut_Potential2RotationMPI(&
       xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi )

  xvh_DData = 0.0D0
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lon) with rigid rotation field(E-W)', &
       answer = xvh_VLon,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 2*cos(xvh_LAT)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lat) with rigid rotation field(E-W)', &
       answer = xvh_VLat,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 2*sin(xvh_LAT)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Rad) with rigid rotation field(E-W)', &
       answer = xvh_VRad,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !----------------- 剛体回転場(南北流) --------------------
  !
  ut_Psi = ut_xvh(xvh_Rad * cos(xvh_Lat) * sin(xvh_Lon))
  ut_Phi = 0.0D0

  call ut_Potential2VectorMPI(&
       xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi )

  xvh_DData = -xvh_Rad*sin(xvh_Lat)*sin(xvh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lon) with rigid rotation field(N-S)', &
       answer = xvh_VLon,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = -xvh_Rad*cos(xvh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lat) with rigid rotation field(N-S)', &
       answer = xvh_VLat,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 0.0D0
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Rad) with rigid rotation field(N-S)', &
       answer = xvh_VRad,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call ut_Potential2RotationMPI(&
       xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi )

  xvh_DData = 2*cos(xvh_Lon)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lon) with rigid rotation field(N-S)', &
       answer = xvh_VLon,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = -2*sin(xvh_Lon)*sin(xvh_Lat)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lat) with rigid rotation field(N-S)', &
       answer = xvh_VLat,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 2*sin(xvh_Lon)*cos(xvh_Lat)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Rad) with rigid rotation field(N-S)', &
       answer = xvh_VRad,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 渦無し場 --------------------
  ut_Psi = 0.0D0
  ut_Phi = ut_xvh(xvh_Rad * sin(xvh_Lat))

  call ut_Potential2VectorMPI(&
       xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi )

  xvh_DData = 0
  call AssertEqual(&
       message='ut_Potential2VectorMPI (Lon) with non-vortical field',&
       answer = xvh_VLon,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 2 * cos(xvh_Lat)
  call AssertEqual(&
       message='ut_Potential2VectorMPI (Lat) with non-vortical field',&
       answer = xvh_VLat,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 2 * sin(xvh_Lat)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Rad) with non-vortical field', &
       answer = xvh_VRad,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call ut_Potential2RotationMPI(&
       xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi )

  xvh_DData = 0.0D0
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lon) with non-vortical field', &
       answer = xvh_VLon,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 0.0D0
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lat) with non-vortical field', &
       answer = xvh_VLat,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 0.0D0
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Rad) with non-vortical field', &
       answer = xvh_VRad,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

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

  ut_Psi = 0.0D0
  ut_Phi = ut_xvh(xvh_Rad**2 * cos(xvh_Lat)*sin(xvh_Lon))

  call ut_Potential2VectorMPI(&
       xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi )

  xvh_DData = 3 * xvh_Rad * cos(xvh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lon) with poloidal field',    &
       answer = xvh_VLon,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = - 3 * xvh_Rad * sin(xvh_Lat) * sin(xvh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lat) with poloidal field',    &
       answer = xvh_VLat,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 2 * xvh_Rad * cos(xvh_Lat) * sin(xvh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Rad) with poloidal field',    &
       answer = xvh_VRad,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call ut_Potential2RotationMPI(&
       xvh_VLon,xvh_VLat,xvh_VRad, ut_Psi, ut_Phi )

  xvh_DData = 4*sin(xvh_Lat)*sin(xvh_Lon)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lon) with poloidal field',  &
       answer = xvh_VLon,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 4*cos(xvh_Lon)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lat) with poloidal field',  &
       answer = xvh_VLat,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_DData = 0.0D0
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Rad) with poloidal field',  &
       answer = xvh_VRad,                                            &
       check = xvh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','ut_mpi_module_derivative_test2', &
       'ut_mpi_module derivative function test #2 succeeded!')

  !------ MPIの終了 ------
  call MPI_FINALIZE(IERR)

end program ut_mpi_module_derivative_test2
