!----------------------------------------------------------------------
!     Copyright (c) 2016 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ut_mpi_module_supack テストプログラム
!
!      トロイダルポロイダル関係微分関数のテスト
!  
!履歴  2016/08/04  竹広真一   wt_mpi_module_supack_test_derivavit2.f90 より改造
!
program ut_mpi_module_svpack_derivative2_test

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use ut_mpi_module_svpack
  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, parameter :: check_digits = 8
  integer, parameter :: ignore = -9

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

  real(8), allocatable :: xyh_DData(:,:,:)
  real(8), allocatable :: xyh_Psi(:,:,:)

  real(8), allocatable :: xyh_VLon(:,:,:)
  real(8), allocatable :: xyh_VLat(:,:,:)
  real(8), allocatable :: xyh_VRad(:,:,:)
  real(8), allocatable :: xyh_RadRot(:,:,:)
  real(8), allocatable :: xyh_RadRotRot(:,:,:)

  integer, parameter :: n=2

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

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

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

  allocate(xyh_DData(0:im-1,1:jm,kc))
  allocate(xyh_Psi(0:im-1,1:jm,kc))

  allocate(xyh_VLon(0:im-1,1:jm,kc))
  allocate(xyh_VLat(0:im-1,1:jm,kc))
  allocate(xyh_VRad(0:im-1,1:jm,kc))
  allocate(xyh_RadRot(0:im-1,1:jm,kc))
  allocate(xyh_RadRotRot(0:im-1,1:jm,kc))

  !========== ut_KxRGrad_ut, xyh_KGrad_ut, ut_QOperator_ut ==========

  ! ----------------- 例 1 --------------------
  xyh_Psi = xyh_rad**n * cos(xyh_lat)*sin(xyh_lon)   ! r**2 Y_1^-1

  xyh_DData = xyh_rad**n * cos(xyh_lat)*cos(xyh_lon)
  call AssertEqual(&
       message='ut_KxRGrad_ut with r^2 Y_1^-1',                      &
       answer = xyh_DData,                                           &
       check = xyh_ut(ut_KxRGrad_ut(ut_xyh(xyh_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! k ・▽ r**n Y_1^1 = (n-1)*r**(n-1)* Y_2^1
  xyh_DData = (n-1)*xyh_rad**(n-1)* cos(xyh_lat)*sin(xyh_lat)*sin(xyh_lon) 
  call AssertEqual(&
       message='xyh_KGrad_ut with r^2 Y_1^-1',                       &
       answer = xyh_DData,                                           &
       check = xyh_KGrad_ut(ut_xyh(xyh_Psi)),                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! Q r**n Y_1^1 = -3*(n-1)*r**(n-1)* Y_2^1
  xyh_DData = - 3*(n-1)*xyh_rad**(n-1)* cos(xyh_lat)*sin(xyh_lat)*sin(xyh_lon) 
  call AssertEqual(&
       message='ut_QOperator_ut r^2 Y_1^-1',                         &
       answer = xyh_DData,                                           &
       check = xyh_ut(ut_QOperatorMPI_ut(ut_xyh(xyh_Psi))),          &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

! ----------------- 例 2 --------------------
  xyh_Psi = cos(xyh_lat)*sin(xyh_lat) * sin(xyh_lon) ! Y_2^1

  xyh_DData = cos(xyh_lat)*sin(xyh_lat) * cos(xyh_lon) ! Y_2^1
  call AssertEqual(&
       message='ut_KxRGrad_ut with Y_2^1',                           &
       answer = xyh_DData,                                           &
       check = xyh_ut(ut_KxRGrad_ut(ut_xyh(xyh_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! k・▽ Y_2^1 = (- 4 Y_3^1 / 15 - Y_1^1 /5) 
  xyh_DData = cos(2*xyh_lat)*cos(xyh_lat)*sin(xyh_lon)/xyh_rad
  call AssertEqual(&
       message='xyh_KGrad_ut with Y_2^1',                            &
       answer = xyh_DData,                                           &
       check = xyh_KGrad_ut(ut_xyh(xyh_Psi)),                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyh_DData = (16*sin(xyh_lat)**2 - 5)*cos(xyh_lat)*sin(xyh_lon)/xyh_rad
  call AssertEqual(&
       message='ut_QOperator_ut Y_2^1',                              &
       answer = xyh_DData,                                           &
       check = xyh_ut(ut_QOperatorMPI_ut(ut_xyh(xyh_Psi))),          &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !========== ut_RadRot_xyh_xyh, ut_RadRotRot_xyh_xyh_xyh ==========

  !---------- 単純流れ v_r=r ----------
  xyh_VLon = 0 ; xyh_VLat = 0 ; xyh_VRad = xyh_Rad
  xyh_RadRot = 0 ; xyh_RadRotRot = 0

  call AssertEqual(&
       message='ut_RadRot_xyh_xyh with v_r=r',                       &
       answer = xyh_RadRot,                                          &
       check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xyh_xyh with v_r=r',                    &
       answer = xyh_RadRotRot,                                       &
       check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !---------- 剛体回転解非線形項(東西流れ) ----------
  xyh_VLon = 0
  xyh_VLat = xyh_Rad*sin(xyh_Lat)*cos(xyh_Lat)
  xyh_VRad = -xyh_Rad*cos(xyh_Lat)**2

  xyh_RadRot = 0 ; xyh_RadRotRot = 0

  call AssertEqual(&
       message='ut_RadRot_xyh_xyh with rigid rotation(E-W)',         &
       answer = xyh_RadRot,                                          &
       check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xyh_xyh with rigid rotation(E-W)',      &
       answer = xyh_RadRotRot,                                       &
       check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !---------- 剛体回転解非線形項(南北流れ) ----------
  xyh_VLon = xyh_Rad*cos(xyh_Lat)*sin(xyh_Lon)*cos(xyh_Lon)
  xyh_VLat = -xyh_Rad*sin(xyh_Lat)*cos(xyh_Lat)*sin(xyh_Lon)**2
  xyh_VRad = -xyh_Rad*(sin(xyh_Lat)**2*sin(xyh_Lon)**2 + cos(xyh_Lon)**2)
  xyh_RadRot = 0 ; xyh_RadRotRot = 0

  call AssertEqual(&
      message='ut_RadRot_xyh_xyh with rigid rotation(N-S)',          &
       answer = xyh_RadRot,                                          &
       check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xyh_xyh with rigid rotation(N-S)',      &
       answer = xyh_RadRotRot,                                       &
       check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !--------- 鉛直渦度を伴うベクトル場 ----------
  xyh_Psi = xyh_Rad**2 * cos(xyh_Lat)*sin(xyh_Lon)   ! r**2 P_1^1

  xyh_VLon =   xyh_gRadLat_ut(ut_xyh(xyh_Psi*xyh_Rad))
  xyh_VLat = - xyh_gRadLon_ut(ut_xyh(xyh_Psi*xyh_Rad))
  xyh_VRad = 0
  xyh_RadRot = 2 * xyh_Psi                           ! r・▽×▽×(ψr) = L_2ψ
  xyh_RadRotRot = 0

  call AssertEqual(&
      message='ut_RadRot_xyh_xyh with r^2 Y_1^-1 toroidal field',    &
       answer = xyh_RadRot,                                          &
       check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xyh_xyh with r^2 Y_1^-1 toroidal field',&
       answer = xyh_RadRotRot,                                       &
       check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !---------- 鉛直速度を伴うベクトル場 -----------
                 ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ
  xyh_VRad = xyh_ut(ut_l2_ut(ut_xyh(xyh_Psi/xyh_Rad)))
  xyh_VLat = xyh_gRadLat_ut(ut_dRad_ut(ut_xyh(xyh_Psi*xyh_Rad)))
  xyh_VLon = xyh_gRadLon_ut(ut_dRad_ut(ut_xyh(xyh_Psi*xyh_Rad)))

  xyh_RadRot = 0
  xyh_RadRotRot = -xyh_ut(ut_l2_ut(ut_lapla_ut(ut_xyh(xyh_Psi))))
                 ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ

  call AssertEqual(&
      message='ut_RadRot_xyh_xyh with non-vortical field',           &
       answer = xyh_RadRot,                                          &
       check = xyh_ut(ut_RadRot_xyh_xyh(xyh_VLon,xyh_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='ut_RadRotRot_xyh_xyh with r^2 Y_1^-1 toroidal field',&
       answer = xyh_RadRotRot,                                       &
       check = xyh_ut(ut_RadRotRot_xyh_xyh_xyh(xyh_VLon,xyh_VLat,xyh_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


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

  !----------------- 剛体回転場 --------------------
  ! 
  ut_Psi = ut_xyh(xyh_Rad * sin(xyh_Lat))
  ut_Phi = 0.0D0

  call ut_Potential2VectorMPI(&
       xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi )

  xyh_DData = xyh_Rad * cos(xyh_Lat)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lon) with rigid rotation field(E-W)', &
       answer = xyh_VLon,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

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

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

  call ut_Potential2RotationMPI(&
       xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi )

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

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

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

  !----------------- 剛体回転場(南北流) --------------------
  ! 
  ut_Psi = ut_xyh(xyh_Rad * cos(xyh_Lat) * sin(xyh_Lon))
  ut_Phi = 0.0D0

  call ut_Potential2VectorMPI(&
       xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi )

  xyh_DData = -xyh_Rad*sin(xyh_Lat)*sin(xyh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lon) with rigid rotation field(N-S)', &
       answer = xyh_VLon,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyh_DData = -xyh_Rad*cos(xyh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lat) with rigid rotation field(N-S)', &
       answer = xyh_VLat,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

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

  call ut_Potential2RotationMPI(&
       xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi )

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

  xyh_DData = -2*sin(xyh_Lon)*sin(xyh_Lat)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lat) with rigid rotation field(N-S)', &
       answer = xyh_VLat,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyh_DData = 2*sin(xyh_Lon)*cos(xyh_Lat)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Rad) with rigid rotation field(N-S)', &
       answer = xyh_VRad,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 渦無し場 --------------------
  ut_Psi = 0.0D0
  ut_Phi = ut_xyh(xyh_Rad * sin(xyh_Lat))

  call ut_Potential2VectorMPI(&
       xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi )

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

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

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

  call ut_Potential2RotationMPI(&
       xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi )

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

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

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

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

  ut_Psi = 0.0D0
  ut_Phi = ut_xyh(xyh_Rad**2 * cos(xyh_Lat)*sin(xyh_Lon))

  call ut_Potential2VectorMPI(&
       xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi )

  xyh_DData = 3 * xyh_Rad * cos(xyh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lon) with poloidal field',   &
       answer = xyh_VLon,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyh_DData = - 3 * xyh_Rad * sin(xyh_Lat) * sin(xyh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Lat) with poloidal field', &
       answer = xyh_VLat,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyh_DData = 2 * xyh_Rad * cos(xyh_Lat) * sin(xyh_Lon)
  call AssertEqual(&
      message='ut_Potential2VectorMPI (Rad) with poloidal field', &
       answer = xyh_VRad,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call ut_Potential2RotationMPI(&
       xyh_VLon,xyh_VLat,xyh_VRad, ut_Psi, ut_Phi )

  xyh_DData = 4*sin(xyh_Lat)*sin(xyh_Lon)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lon) with poloidal field', &
       answer = xyh_VLon,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyh_DData = 4*cos(xyh_Lon)
  call AssertEqual(&
      message='ut_Potential2RotationMPI (Lat) with poloidal field', &
       answer = xyh_VLat,                                            &
       check = xyh_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

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

  call MessageNotify('M','ut_mpi_module_svpack_derivative_test2', &
       'ut_mpi_module_svpack derivative function test #2 succeeded!')

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

  call MPI_FINALIZE(IERR)

end program ut_mpi_module_svpack_derivative2_test
