!----------------------------------------------------------------------
!     Copyright (c) 2017 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wsc_mpi_module_supack テストプログラム
!
!      トロイダルポロイダル関係微分関数のテスト
!  
!履歴  2017/05/31  竹広真一   wt_module_nmomp_test_derivavit2.f90 より改造
!
program wsc_mpi_module_supack_derivative2_test

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use wsc_mpi_module_supack
  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 :: wc_Psi(:,:)
  real(8), allocatable :: ws_Phi(:,:)

  real(8), allocatable :: xvz_DData(:,:,:)
  real(8), allocatable :: xvz_Psi(:,:,:)

  real(8), allocatable :: xvz_VLon(:,:,:)
  real(8), allocatable :: xvz_VLat(:,:,:)
  real(8), allocatable :: xvz_VRad(:,:,:)
  real(8), allocatable :: xvz_RadRot(:,:,:)
  real(8), allocatable :: xvz_RadRotRot(:,:,:)

  real(8), allocatable :: xvz_Xi(:,:,:)  

  integer, parameter :: n=2

  integer :: iproc, np, ierr

  real(8), parameter  :: pi=3.1415926535897932385D0

  !---------------- 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','wsc_mpi_module_supack_derivative2_test', &
       'wsc_mpi_module_supack derivative function test #2')

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

  allocate(wc_Psi(nc,0:lm))
  allocate(ws_Phi(nc,lm))

  allocate(xvz_DData(0:im-1,1:jc,0:km))
  allocate(xvz_Psi(0:im-1,1:jc,0:km))

  allocate(xvz_VLon(0:im-1,1:jc,0:km))
  allocate(xvz_VLat(0:im-1,1:jc,0:km))
  allocate(xvz_VRad(0:im-1,1:jc,0:km))
  allocate(xvz_RadRot(0:im-1,1:jc,0:km))
  allocate(xvz_RadRotRot(0:im-1,1:jc,0:km))
  allocate(xvz_Xi(0:im-1,1:jc,0:km))

  xvz_xi = (xvz_Rad - ri)/(ro-ri)

  !========== wc_KxRGrad_wc, ws_KxRGrad_ws ==========

  ! ----------------- 例 1 --------------------
  xvz_Psi = sin(Pi*xvz_xi) * cos(xvz_lat)*sin(xvz_lon)   ! Y_1^-1 sin(z)

  xvz_DData = sin(Pi*xvz_xi) *  cos(xvz_lat)*cos(xvz_lon)
  call AssertEqual(&
       message='ws_KxRGrad_ws with Y_1^-1 sin(z)',                   &
       answer = xvz_DData,                                           &
       check = xvz_ws(ws_KxRGrad_ws(ws_xvz(xvz_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 例 2 --------------------
  xvz_Psi = cos(xvz_lat)*sin(xvz_lat) * sin(xvz_lon) * cos(PI*xvz_xi)  ! Y_2^1 cos(z)

  xvz_DData = cos(xvz_lat)*sin(xvz_lat) * cos(xvz_lon) * cos(PI*xvz_xi) ! Y_2^1
  call AssertEqual(&
       message='wc_KxRGrad_wc with Y_2^1 cos(z)',                    &
       answer = xvz_DData,                                           &
       check = xvz_wc(wc_KxRGrad_wc(wc_xvz(xvz_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !========== wsc_Potential2VectorMPI, wsc_Potential2RotationMPI, wc_RadRot, ws_RadRotRot ==========

  !----------------- 剛体回転場 --------------------
  ! 
  wc_Psi = wc_xvz(sin(xvz_Lat))
  ws_Phi = 0.0D0

  call wsc_Potential2VectorMPI(&
       xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi )

  xvz_DData = cos(xvz_Lat)
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Lon) with rigid rotation field(E-W)', &
       answer = xvz_VLon,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Lat) with rigid rotation field(E-W)', &
       answer = xvz_VLat,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Rad) with rigid rotation field(E-W)', &
       answer = xvz_VRad,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvz_RadRot = xvz_wc(wc_L2_wc(wc_Psi))
  call AssertEqual(&
       message='wc_RadRot_xvz_xvz with with rigid rotation field(E-W)', &
       answer = xvz_RadRot,                                          &
       check = xvz_wc(wc_RadRot_xvz_xvz(xvz_VLon,xvz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvz_RadRotRot = 0.0D0
  call AssertEqual(&
       message='ws_RadRotRot_xvz_xvz with rigid rotation field(E-W)', &
       answer = xvz_RadRotRot,                                       &
       check = xvz_ws(ws_RadRotRot_xvz_xvz_xvz(xvz_VLon,xvz_VLat,xvz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wsc_Potential2RotationMPI(&
       xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi )

  xvz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Lon) with rigid rotation field(E-W)', &
       answer = xvz_VLon,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Lat) with rigid rotation field(E-W)', &
       answer = xvz_VLat,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = xvz_RadRot/Ro
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Rad) with rigid rotation field(E-W)', &
       answer = xvz_VRad,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !----------------- 剛体回転場(南北流) --------------------
  ! 
  wc_Psi = wc_xvz(cos(xvz_Lat) * sin(xvz_Lon))
  ws_Phi = 0.0D0

  call wsc_Potential2VectorMPI(&
       xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi )

  xvz_DData = -sin(xvz_Lat)*sin(xvz_Lon)
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Lon) with rigid rotation field(N-S)', &
       answer = xvz_VLon,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = -cos(xvz_Lon)
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Lat) with rigid rotation field(N-S)', &
       answer = xvz_VLat,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Rad) with rigid rotation field(N-S)', &
       answer = xvz_VRad,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvz_RadRot = xvz_wc(wc_L2_wc(wc_Psi))
  call AssertEqual(&
       message='wc_RadRot_xvz_xvz with with rigid rotation field(N-S)', &
       answer = xvz_RadRot,                                          &
       check = xvz_wc(wc_RadRot_xvz_xvz(xvz_VLon,xvz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvz_RadRotRot = 0.0D0
  call AssertEqual(&
       message='ws_RadRotRot_xvz_xvz with rigid rotation field(N-S)', &
       answer = xvz_RadRotRot,                                       &
       check = xvz_ws(ws_RadRotRot_xvz_xvz_xvz(xvz_VLon,xvz_VLat,xvz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wsc_Potential2RotationMPI(&
       xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi )

  xvz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Lon) with rigid rotation field(N-S)', &
       answer = xvz_VLon,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Lat) with rigid rotation field(N-S)', &
       answer = xvz_VLat,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = xvz_RadRot/Ro
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Rad) with rigid rotation field(N-S)', &
       answer = xvz_VRad,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 渦無し場 --------------------
  wc_Psi = 0.0D0
  ws_Phi = ws_xvz(sin(PI*xvz_xi) * sin(xvz_Lat))

  call wsc_Potential2VectorMPI(&
       xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi )

  xvz_DData = 0
  call AssertEqual(&
       message='wsc_Potential2VectorMPI (Lon) with non-vortical field', &
       answer = xvz_VLon,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = PI * cos(xvz_Lat) * cos(PI*xvz_xi) 
  call AssertEqual(&
       message='wsc_Potential2VectorMPI (Lat) with non-vortical field', &
       answer = xvz_VLat,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = 2 * sin(xvz_Lat) * sin(PI*xvz_xi)/Ro
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Rad) with non-vortical field',  &
       answer = xvz_VRad,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvz_RadRot = 0.0D0
  call AssertEqual(&
       message='wc_RadRot_xvz_xvz with non-vortical field',          &
       answer = xvz_RadRot,                                          &
       check = xvz_wc(wc_RadRot_xvz_xvz(xvz_VLon,xvz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvz_RadRotRot = - xvz_ws(ws_L2_ws(ws_Lapla_ws((ws_Phi))))
  call AssertEqual(&
       message='ws_RadRotRot_xvz_xvz with non-vortical field',       &
       answer = xvz_RadRotRot,                                       &
       check = xvz_ws(ws_RadRotRot_xvz_xvz_xvz(xvz_VLon,xvz_VLat,xvz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wsc_Potential2RotationMPI(&
       xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi )

  xvz_DData = (2/Ro**2 + PI**2) * cos(xvz_Lat) * sin(PI*xvz_xi)
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Lon) with non-vortical field', &
       answer = xvz_VLon,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = 0.0D0
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Lat) with non-vortical field', &
       answer = xvz_VLat,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = xvz_RadRot/Ro
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Rad) with non-vortical field', &
       answer = xvz_VRad,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

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

  wc_Psi = 0.0D0
  ws_Phi = ws_xvz(sin(2*PI*xvz_xi) * cos(xvz_Lat)*sin(xvz_Lon))

  call wsc_Potential2VectorMPI(&
       xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi )

  xvz_DData = 2*PI * cos(2*PI*xvz_xi) * cos(xvz_Lon)
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Lon) with poloidal field',   &
       answer = xvz_VLon,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = - 2*PI* cos(2*PI*xvz_xi) * sin(xvz_Lat)*sin(xvz_Lon)
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Lat) with poloidal field', &
       answer = xvz_VLat,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = 2*sin(2*PI*xvz_xi) * cos(xvz_Lat)*sin(xvz_Lon)/Ro
  call AssertEqual(&
      message='wsc_Potential2VectorMPI (Rad) with poloidal field', &
       answer = xvz_VRad,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvz_RadRot = 0.0D0
  call AssertEqual(&
       message='wc_RadRot_xvz_xvz with with poloidal field',         &
       answer = xvz_RadRot,                                          &
       check = xvz_wc(wc_RadRot_xvz_xvz(xvz_VLon,xvz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xvz_RadRotRot = - xvz_ws(ws_L2_ws(ws_Lapla_ws((ws_Phi))))
  call AssertEqual(&
       message='ws_RadRotRot_xvz_xvz with with poloidal field',      &
       answer = xvz_RadRotRot,                                       &
       check = xvz_ws(ws_RadRotRot_xvz_xvz_xvz(xvz_VLon,xvz_VLat,xvz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wsc_Potential2RotationMPI(&
       xvz_VLon,xvz_VLat,xvz_VRad, wc_Psi, ws_Phi )

  xvz_DData = -(2/Ro**2 + 4*PI**2)*sin(2*PI*xvz_xi)*sin(xvz_Lat)*sin(xvz_Lon)
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Lon) with poloidal field',    &
       answer = xvz_VLon,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = -(2/Ro**2+4*PI**2) * sin(2*PI*xvz_xi) * cos(xvz_Lon)
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Lat) with poloidal field', &
       answer = xvz_VLat,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xvz_DData = xvz_RadRot/Ro
  call AssertEqual(&
      message='wsc_Potential2RotationMPI (Rad) with poloidal field',    &
       answer = xvz_VRad,                                            &
       check = xvz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','wsc_mpi_module_supack_derivative_test2', &
       'wsc_mpi_module_supack derivative function test #2 succeeded!')

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

  call MPI_FINALIZE(IERR)

end program wsc_mpi_module_supack_derivative2_test
