!----------------------------------------------------------------------
!     Copyright (c) 2019-2020 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ut_mpi_module テストプログラム
!
!      微分のテスト
!
!履歴  2019/05/28  竹広真一
!      2020/11/11  竹広真一  セクター計算オプション導入
!
program ut_mpi_module_derivative_mint_test1

  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         ! 切断波数の設定(水平, 動径)
  real(8),parameter  :: ri=0.5, ro=1.5       ! 内外半径

  integer :: npv=2

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

  real(8), allocatable              :: uz_Data(:,:)
  real(8), allocatable              :: uz_DData(:,:)

  real(8), allocatable              :: xvh_VLon(:,:,:)
  real(8), allocatable              :: xvh_VLat(:,:,:)
  real(8), allocatable              :: xvh_VRad(:,:,:)
  real(8), allocatable              :: xvh_Div(:,:,:)
  real(8), allocatable              :: xvh_Psi(:,:,:)

  ! 判定誤差設定
  integer, parameter :: check_digits = 10
  integer, parameter :: ignore = -11

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

  integer :: n
  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_test1', &
       'ut_mpi_module_mint derivative function test #1')

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

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

  allocate(uz_Data(nmc,0:km))
  allocate(uz_DData(nmc,0:km))

  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_Div(0:im-1,1:jc,kc))
  allocate(xvh_Psi(0:im-1,1:jc,kc))

  uz_Data = 0.0D0 ; uz_Data = uz_ut(ut_xvh(xvh_ut(uz_ut(uz_Data))))
  !==================== ut_DivRad_ut ====================
  do n=nmin,nmax
     write(6,*) 'n=',n

     uz_Data = uz_Rad**n
     uz_DData = (n+2)*uz_Rad**(n-1)

     call AssertEqual(&
       message='ut_Divrad_ut with r^n',                             &
       answer = uz_DData,                                           &
       check = uz_ut(ut_DivRad_ut(ut_uz(uz_Data))),                 &
       significant_digits = check_digits, ignore_digits = ignore    &
       )

!!$     xvh_Data = xvh_Rad**n
!!$     xvh_DData = (n+2)*xvh_Rad**(n-1)
!!$
!!$     call AssertEqual(&
!!$       message='ut_Divrat_ut with r^n',                              &
!!$       answer = xvh_DData,                                           &
!!$       check = xvh_ut(ut_DivRad_ut(ut_xvh(xvh_Data))),               &
!!$       significant_digits = check_digits, ignore_digits = ignore     &
!!$       )
  enddo

  !========== xvz_GradLon_ut, xvz_GradLat_ut, wt_Div_xvz_xvz_xvz ==========

  ! ------ r^2 Y_1^0 Toridal field -----
  n = 2
  write(6,*) 'n=',n
  xvh_Psi = xvh_Rad**n * sin(xvh_Lat)                       ! r^2 Y_1^0
  xvh_VRad = 0

  xvh_VLon =  0.0D0
  call AssertEqual(&
       message='xvh_GradLon_ut with r^2 Y_1^0',                      &
       answer = xvh_VLon,                                            &
       check = xvh_GradLon_ut(ut_xvh(xvh_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_VLat = xvh_Rad**(n-1)*cos(xvh_Lat)
  call AssertEqual(&
       message='xvh_GradLat_ut with r^2 Y_1^0',                      &
       answer = xvh_VLat,                                            &
       check = xvh_GradLat_ut(ut_xvh(xvh_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_Div = - 2 * xvh_Psi/xvh_Rad**2
  call AssertEqual(&
       message='ut_Div_xvh_xvh_xvh_ut with Toroiral r^2 Y_1^0',       &
       answer = xvh_Div,                                              &
       check = xvh_ut(ut_Div_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)),&
       significant_digits = check_digits, ignore_digits = ignore      &
       )

  ! ------ r^2 Y_2^-2 Toridal field -----
  n = 2
  xvh_VRad = 0
  xvh_Psi = xvh_Rad**n * cos(xvh_Lat)**2 * sin(2*xvh_Lon) ! r^2 Y_2^-2

  xvh_VLon =  xvh_Rad**(n-1)*cos(xvh_Lat) * 2*cos(2*xvh_Lon)
  call AssertEqual(&
       message='xvh_GradLon_ut with r^2 Y_2^-2',                     &
       answer = xvh_VLon,                                            &
       check =  xvh_GradLon_ut(ut_xvh(xvh_Psi)),                     &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_VLat = -xvh_Rad**(n-1)*2*cos(xvh_Lat)*sin(xvh_Lat)*sin(2*xvh_Lon)
  call AssertEqual(&
       message='xvh_GradLat_ut with r^2 Y_2^-2',                     &
       answer = xvh_VLat,                                            &
       check = xvh_GradLat_ut(ut_xvh(xvh_Psi)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_Div = - 6* xvh_Psi/xvh_Rad**2
  call AssertEqual(&
       message='ut_Div_xvh_xvh_xvh_ut with Toroiral r2 Y_1^-2',      &
       answer = xvh_Div,                                             &
       check = xvh_ut(ut_Div_xvh_xvh_xvh(xvh_VLon,xvh_VLat,xvh_VRad)),&
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !==================== ut_DivLon_xyz, ut_DivLat_xyz ====================
  n=2
  write(6,*)'n =',n

  xvh_VLon  = xvh_rad**n * cos(xvh_Lat)**3*sin(2*xvh_Lon)
  xvh_DData = xvh_rad**(n-1)*cos(xvh_Lat)**2* 2*cos(2*xvh_Lon)
  call AssertEqual(&
       message='ut_DivLon_xvh with r^n Y_2^-2 cos^2 φ',             &
       answer = xvh_DData,                                           &
       check = xvh_ut(ut_DivLon_xvh(xvh_VLon)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvh_VLat   = xvh_rad**n * cos(xvh_Lat)**3*sin(2*xvh_Lon)
  xvh_DData = -4*xvh_rad**(n-1)*cos(xvh_Lat)**2*sin(xvh_Lat)*sin(2*xvh_Lon)
  call AssertEqual(&
       message='ut_DivLat_xvh with r^n Y_2^-2 cos φ',               &
       answer = xvh_DData,                                           &
       check = xvh_ut(ut_DivLat_xvh(xvh_VLat)),                      &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','ut_mpi_module_derivative_test1', &
       'ut_mpi_module_mint derivative function test #1 succeeded!')

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

  call MPI_FINALIZE(IERR)

end program ut_mpi_module_derivative_mint_test1
