!----------------------------------------------------------------------
!     Copyright (c) 2018--2020 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ut_mpi_module テストプログラム
!
!      トロイダルポロイダル関係微分関数のテスト
!
!履歴  2018/08/09  竹広真一
!      2020/11/11  竹広真一  セクター計算オプション導入
!      2021/11/15  竹広真一  ut_RadRot_RadRotRot テスト追加
!
program ut_mpi_module_derivative_mint_test2

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

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

  integer :: npv=2

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

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

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

  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))

  !========== 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     &
       )

  call ut_RadRot_RadRotRot( &
       xvh_VLon, xvh_VLat, xvh_VRad, ut_Psi, ut_Phi )
  
  call AssertEqual(&
       message='ut_RadRot_RadRotRot_xvh (RadRot) with v_r=r',        &
       answer = xvh_RadRot,                                          &
       check = xvh_ut(ut_Psi),                                       &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  call AssertEqual(&
       message='ut_RadRot_RadRotRot_xvh (RadRotRot) with v_r=r',     &
       answer = xvh_RadRotRot,                                       &
       check = xvh_ut(ut_Phi),                                       &
       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     &
       )

  call ut_RadRot_RadRotRot( &
       xvh_VLon, xvh_VLat, xvh_VRad, ut_Psi, ut_Phi )
  
  call AssertEqual(&
       message='ut_RadRot_RadRotRot_xvh (RadRot) with rigid rotation(E-W)',        &
       answer = xvh_RadRot,                                          &
       check = xvh_ut(ut_Psi),                                       &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  call AssertEqual(&
       message='ut_RadRot_RadRotRot_xvh (RadRotRot) with rigid rotation(E-W)',     &
       answer = xvh_RadRotRot,                                       &
       check = xvh_ut(ut_Phi),                                       &
       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**3 * cos(xvh_Lat)**2*sin(2*xvh_Lon)   ! r**2 P_2^-2

  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 = 6 * xvh_Psi                           ! r・▽×▽×(ψr) = L_2ψ
  xvh_RadRotRot = 0

  call AssertEqual(&
      message='ut_RadRot_xvh_xvh with r^2 Y_2^-2 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_2^-2 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     &
       )

  call ut_RadRot_RadRotRot( &
       xvh_VLon, xvh_VLat, xvh_VRad, ut_Psi, ut_Phi )
  
  call AssertEqual(&
       message='ut_RadRot_RadRotRot_xvh (RadRot) with r^2 Y_2^-2 toroidal field',&
       answer = xvh_RadRot,                                          &
       check = xvh_ut(ut_Psi),                                       &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  call AssertEqual(&
       message='ut_RadRot_RadRotRot_xvh (RadRotRot) with r^2 Y_2^-2 toroidal field',&
       answer = xvh_RadRotRot,                                       &
       check = xvh_ut(ut_Phi),                                       &
       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_2^-2 Poloidal 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     &
       )

  call ut_RadRot_RadRotRot( &
       xvh_VLon, xvh_VLat, xvh_VRad, ut_Psi, ut_Phi )
  
  call AssertEqual(&
       message='ut_RadRot_RadRotRot_xvh (RadRot) with r^2 Y_2^-2 Poloidal field',&
       answer = xvh_RadRot,                                          &
       check = xvh_ut(ut_Psi),                                       &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  call AssertEqual(&
       message='ut_RadRot_RadRotRot_xvh (RadRotRot) with r^2 Y_2^-2 Poloidal field',&
       answer = xvh_RadRotRot,                                       &
       check = xvh_ut(ut_Phi),                                       &
       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**3 * cos(xvh_Lat)**2*cos(2*xvh_Lon))

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

  xvh_DData = 4*xvh_Rad**2 * cos(xvh_Lat) *(-2)*sin(2*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 = 4*xvh_Rad**2 * (-2)*cos(xvh_Lat)*sin(xvh_Lat)*cos(2*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 = 6 * xvh_Rad**2 * cos(xvh_Lat)**2 * cos(2*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 = 6* xvh_Rad * 2*cos(xvh_Lat)*sin(xvh_Lat)*cos(2*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 = 6* xvh_Rad * cos(xvh_Lat)*(-2)*sin(2*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_mint_test2', &
       'ut_mpi_module_mint derivative function test #2 succeeded!')

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

end program ut_mpi_module_derivative_mint_test2
