!----------------------------------------------------------------------
!     Copyright (c) 2019 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  ua_deriv_mpi_module_sypack テストプログラム :: 微分関数のテスト
!
!履歴  2019/04/30  竹広真一
!履歴  2020/01/25  竹広真一
!
program ua_deriv_mpi_module_sypack_test

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use ua_mpi_module
!!$  use ua_base_mpi_module_sypack, only : ua_base_mpi_Initial, kp_kl, &
!!$       xvb_ua, ua_xvb, jc, kc, xv_Lon, xv_Lat
!!$  use ua_deriv_mpi_module_sypack
  use mpi
  implicit none

  integer, parameter :: im=64, jm=32, nm=21, km=2, npv=2

  real(8), allocatable  ::  xvb_data1(:,:,:)              ! 元の関数
  real(8), allocatable  ::  xvb_data2(:,:,:)              ! 元の関数
  real(8), allocatable  ::  xvb_ddata(:,:,:)              ! 微分の正解
  real(8), allocatable  ::  xv_data1(:,:)                 ! 元の関数
  real(8), allocatable  ::  xv_data2(:,:)                 ! 元の関数
  real(8), allocatable  ::  xv_ddata(:,:)                 ! 微分の正解
  real(8), allocatable  ::  xv_mu(:,:)                    ! μ=sinφ

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

  integer :: ip, np, ierr


 !---------------- MPI スタート ---------------------
  call MPI_INIT(IERR)
  call MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)

  call MessageNotify('M','ua_deriv_mpi_module_sypack_test', &
                         'ua_deriv_mpi_module_sypack function tests')

!!$  call ua_base_mpi_Initial( nm, im, jm, km, npv=npv )
!!$  call ua_deriv_mpi_Initial
  call ua_mpi_Initial( nm, im, jm, km, npv=npv )

  allocate(xvb_data1(0:im-1,jc,kc))
  allocate(xvb_data2(0:im-1,jc,kc))
  allocate(xvb_ddata(0:im-1,jc,kc))
  allocate(xv_data1(0:im-1,jc))
  allocate(xv_data2(0:im-1,jc))
  allocate(xv_ddata(0:im-1,jc))
  allocate(xv_mu(0:im-1,jc))

  !============== 微分計算 (lon,lat座標系用) のテスト ==============

  !---- Y_1^-1 Y_2^1 のテスト ----
  if ( kp_kl(1) /= 0 ) &
       xvb_data1(:,:,kp_kl(1)) = -cos(xv_Lat)*sin(xv_Lon)             ! Y_1^{-1}
  if ( kp_kl(2) /= 0 ) &
       xvb_data1(:,:,kp_kl(2))  = sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon) ! Y_2^1

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = 2*cos(xv_Lat)*sin(xv_Lon)             ! ua_Lapla_ua
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = -6*sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon) ! ua_Lapla_ua
  call AssertEqual(&
    message='ua_Lapla_ua with Y_1^-1 and Y^2_1',                  &
    answer = xvb_ua(ua_Lapla_ua(ua_xvb(xvb_data1))),              &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

 
  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = 1.0/2.0*cos(xv_Lat)*sin(xv_Lon)      ! ua_LaplaInv_ua
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = -1.0D0/6.0*sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon)
  call AssertEqual(&
    message='ua_LaplaInv_ua with Y_1^-1 and Y^2_1',               &
    answer = xvb_ua(ua_LaplaInv_ua(ua_xvb(xvb_data1))),           &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = -cos(xv_Lat)*cos(xv_Lon)         ! ua_DLon_ua
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = -sin(xv_Lat)*cos(xv_Lat) * sin(xv_Lon)
  call AssertEqual(&
    message='ua_DLon_ua with Y_1^-1 and Y^2_1',                   &
    answer = xvb_ua(ua_DLon_ua(ua_xvb(xvb_data1))),               &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = -cos(xv_Lon)                     ! xvb_GradLon_ua
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = -sin(xv_Lat) * sin(xv_Lon)
  call AssertEqual(&
    message='xvb_GradLon_ua with Y_1^-1 and Y^2_1',               &
    answer = xvb_GradLon_ua(ua_xvb(xvb_data1)),                   &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = sin(xv_Lat)*sin(xv_Lon)         ! xvb_GradLat_ua
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = cos(2*xv_Lat) * cos(xv_Lon)
  call AssertEqual(&
    message='xvb_GradLat_ua with Y_1^-1 and Y^2_1',               &
    answer = xvb_GradLat_ua(ua_xvb(xvb_data1)),                   &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xv_data1 = -cos(xv_Lat)*sin(xv_Lon)                 ! Y_1^{-1}
  xv_ddata = 2*cos(xv_Lat)*sin(xv_Lon)                ! ua_Lapla_ua
  call AssertEqual(&
    message='u_Lapla_u with Y_1^-1',                              &
    answer = xv_u(u_Lapla_u(u_xv(xv_data1))),                     &
    check = xv_ddata,                                             &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  xv_ddata = 1.0/2.0*cos(xv_Lat)*sin(xv_Lon)          ! ua_LaplaInv_ua
  call AssertEqual(&
    message='u_LaplaInv_u with Y_1^-1',                           &
    answer = xv_u(u_LaplaInv_u(u_xv(xv_data1))),                  &
    check = xv_ddata,                                             &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  xv_ddata = sin(xv_Lat)*sin(xv_Lon)                  ! xvb_GradLat_ua
  call AssertEqual(&
    message='xv_GradLat_u with Y_1^-1',                           &
    answer = xv_GradLat_u(u_xv(xv_data1)),                        &
    check = xv_ddata,                                             &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xv_data1 = sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon)    ! Y_2^1
  xv_ddata = -6*sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon) ! ua_Lapla_ua
  call AssertEqual(&
    message='u_Lapla_u with Y_2^1',                               &
    answer = xv_u(u_Lapla_u(u_xv(xv_data1))),                     &
    check = xv_ddata,                                             &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  xv_ddata = -1.0D0/6.0*sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon) ! ua_LaplaInv_ua
  call AssertEqual(&
    message='u_LaplaInv_u with Y_2^1',                            &
    answer = xv_u(u_LaplaInv_u(u_xv(xv_data1))),                  &
    check = xv_ddata,                                             &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  xv_ddata = cos(2*xv_Lat) * cos(xv_Lon)                      ! xvb_GradLat_ua
  call AssertEqual(&
    message='xv_GradLat_u with Y_2^1',                            &
    answer = xv_GradLat_u(u_xv(xv_data1)),                        &
    check = xv_ddata,                                             &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  
  !---- Y_1^-1 cosφ, Y_2^1 cosφ のテスト ----
  if ( kp_kl(1) /= 0 ) &
       xvb_data1(:,:,kp_kl(1)) = -cos(xv_Lat)**2*sin(xv_Lon)              ! Y_1^-1 cosφ
  if ( kp_kl(2) /= 0 ) &
       xvb_data1(:,:,kp_kl(2)) = sin(xv_Lat)*cos(xv_Lat)**2 * cos(xv_Lon) ! Y_2^1 cosφ

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = -cos(xv_Lat)*cos(xv_Lon)                ! ua_DivLon_xv
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = -sin(xv_Lat)*cos(xv_Lat)*sin(xv_Lon)    ! ua_DivLon_xv
  call AssertEqual(&
    message='ua_DivLon_xvb with Y_1^-1 cosφ and Y^2_1 cosφ',    &
    answer = xvb_ua(ua_DivLon_xvb(xvb_data1)),                    &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = 3*sin(xv_Lat)*cos(xv_Lat)*sin(xv_Lon)        !ua_DivLat_ua
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = cos(xv_Lat)*(1-4*sin(xv_Lat)**2)*cos(xv_Lon) !ua_DivLat_ua
  call AssertEqual(&
    message='ua_DivLat_xvb with Y_1^-1 cosφ and Y^2_1 cosφ',    &
    answer = xvb_ua(ua_DivLat_xvb(xvb_data1)),                    &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xv_data1 = -cos(xv_Lat)**2*sin(xv_Lon)              ! Y_1^-1 cosφ
  xv_ddata = -cos(xv_Lat)*cos(xv_Lon)                ! u_DivLon_xv
  call Assertequal(&
    message='u_DivLon_xv with Y_1^-1 cosφ and Y^2_1 cosφ',    &
    answer = xv_u(u_DivLon_xv(xv_data1)),                       &
    check = xv_ddata,                                           &
    significant_digits = check_digits, ignore_digits = ignore   &
    )
  xv_ddata = 3*sin(xv_Lat)*cos(xv_Lat)*sin(xv_Lon)        !u_DivLat_u
  call Assertequal(&
    message='u_DivLat_xv with Y_1^-1 cosφ and Y^2_1 cosφ',    &
    answer = xv_u(u_DivLat_xv(xv_data1)),                       &
    check = xv_ddata,                                           &
    significant_digits = check_digits, ignore_digits = ignore   &
    )
  
  xv_data1 = sin(xv_Lat)*cos(xv_Lat)**2 * cos(xv_Lon) ! Y_2^1 cosφ
  xv_ddata = -sin(xv_Lat)*cos(xv_Lat)*sin(xv_Lon)    ! u_DivLon_xv
  call Assertequal(&
    message='u_DivLon_xv with Y_1^-1 cosφ and Y^2_1 cosφ',    &
    answer = xv_u(u_DivLon_xv(xv_data1)),                       &
    check = xv_ddata,                                           &
    significant_digits = check_digits, ignore_digits = ignore   &
    )
  xv_ddata = cos(xv_Lat)*(1-4*sin(xv_Lat)**2)*cos(xv_Lon)     !u_DivLat_u
  call Assertequal(&
    message='u_DivLat_xv with Y_1^-1 cosφ and Y^2_1 cosφ',    &
    answer = xv_u(u_DivLat_xv(xv_data1)),                       &
    check = xv_ddata,                                           &
    significant_digits = check_digits, ignore_digits = ignore   &
    )

  
  !---- Jacobian のテスト ----
  if ( kp_kl(1) /= 0 ) &
       xvb_data1(:,:,kp_kl(1)) = -cos(xv_Lat)*sin(xv_Lon)          ! Y_1^{-1}

  if ( kp_kl(2) /= 0 ) &
       xvb_data1(:,:,kp_kl(2)) = -cos(xv_Lat)*sin(xv_Lon)          ! Y_1^{-1}

  if ( kp_kl(1) /= 0 ) &
       xvb_data2(:,:,kp_kl(1)) = -cos(xv_Lat)*sin(xv_Lon)         ! Y_1^{-1}
  if ( kp_kl(2) /= 0 ) &
       xvb_data2(:,:,kp_kl(2)) = sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon)   ! Y_2^1

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = 0.0
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = sin(xv_Lat)**2 - cos(xv_Lat)**2*cos(xv_Lon)**2

  call AssertEqual(&
    message='ua_JacobianMPI_ua_ua Y_1^-1 and Y^2_1',                 &
    answer = xvb_ua(ua_JacobianMPI_ua_ua(ua_xvb(xvb_data1),ua_xvb(xvb_data2))), &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xv_data1 = -cos(xv_Lat)*sin(xv_Lon)         ! Y_1^{-1}
  xv_data2 = -cos(xv_Lat)*sin(xv_Lon)         ! Y_1^{-1}
  xv_ddata = 0.0
  call AssertEqual(&
    message='u_JacobianMPI_u_u Y_1^-1 and Y^2_1',                    &
    answer = xv_u(u_JacobianMPI_u_u(u_xv(xv_data1),u_xv(xv_data2))), &
    check = xv_ddata,                                                &
    significant_digits = check_digits, ignore_digits = ignore        &
    )

  xv_data1 = -cos(xv_Lat)*sin(xv_Lon)                ! Y_1^{-1}
  xv_data2 = sin(xv_Lat)*cos(xv_Lat) * cos(xv_Lon)   ! Y_2^1
  xv_ddata = sin(xv_Lat)**2 - cos(xv_Lat)**2*cos(xv_Lon)**2
  call AssertEqual(&
    message='u_JacobianMPI_u_u Y_1^-1 and Y^2_1',                    &
    answer = xv_u(u_JacobianMPI_u_u(u_xv(xv_data1),u_xv(xv_data2))), &
    check = xv_ddata,                                                &
    significant_digits = check_digits, ignore_digits = ignore        &
    )
 

  !============== 微分計算 (λ,μ座標系用) のテスト ==============
  xv_mu = sin(xv_Lat)

  !----- Y_2^0, Y_1^1 のテスト -----
  if ( kp_kl(1) /= 0 ) &
       xvb_data1(:,:,kp_kl(1)) = 3*xv_mu**2-1                              ! Y_2^0
  if ( kp_kl(2) /= 0 ) &
       xvb_data1(:,:,kp_kl(2)) = sqrt(1-xv_mu**2)*cos(xv_Lon)              ! Y_1^1

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = 0.0
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = -sqrt(1-xv_mu**2)*sin(xv_Lon)
  call AssertEqual(&
    message='xvb_GradLambda_ua with Y_2^0 and Y^1_1',             &
    answer = xvb_GradLambda_ua(ua_xvb(xvb_data1)),                &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )
  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = 6*xv_mu*(1-xv_mu**2)
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = -xv_mu*sqrt(1-xv_mu**2)*cos(xv_Lon)
  call AssertEqual(&
    message='xvb_GradMu_ua with Y_2^0 and Y^1_1',                 &
    answer = xvb_GradMu_ua(ua_xvb(xvb_data1)),                    &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xv_data1 = 3*xv_mu**2-1                              ! Y_2^0
  xv_ddata = 0.0
  call AssertEqual(&
    message='xv_GradLambda_u with Y_2^0',                        &
    answer = xv_GradLambda_u(u_xv(xv_data1)),                    &
    check = xv_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore    &
    )
  xv_ddata = 6*xv_mu*(1-xv_mu**2)
  call AssertEqual(&
    message='xv_GradMu_u with Y_2^0',                             &
    answer = xv_GradMu_u(u_xv(xv_data1)),                         &
    check = xv_ddata,                                             &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xv_data1 = sqrt(1-xv_mu**2)*cos(xv_Lon)              ! Y_1^1
  xv_ddata = -sqrt(1-xv_mu**2)*sin(xv_Lon)
  call AssertEqual(&
    message='xv_GradLambda_u with Y^1_1',                        &
    answer = xv_GradLambda_u(u_xv(xv_data1)),                    &
    check = xv_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore    &
    )
  xv_ddata = -xv_mu*sqrt(1-xv_mu**2)*cos(xv_Lon)
  call AssertEqual(&
    message='xv_GradMu_u with Y_1^1',                             &
    answer = xv_GradMu_u(u_xv(xv_data1)),                         &
    check = xv_ddata,                                             &
    significant_digits = check_digits, ignore_digits = ignore     &
    )


  !----- Y_2^0(1-μ^2), Y_1^1 (1-μ^2) のテスト -----
  if ( kp_kl(1) /= 0 ) &
       xvb_data1(:,:,kp_kl(1)) = (3*xv_mu**2-1)*(1-xv_mu**2)       ! Y_2^0 (1-μ^2)
  if ( kp_kl(2) /= 0 ) &
       xvb_data1(:,:,kp_kl(2)) = (1-xv_mu**2)**(3.0/2)*cos(xv_Lon) ! Y_1^1 (1-μ^2)

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = 0.0
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = -(1-xv_mu**2)**(1.0d0/2)*sin(xv_Lon)
  call AssertEqual(&
    message='ua_DivLambda_xvb with Y_2^0 (1-μ^2) and Y^1_1 √(1-μ^2)',&
    answer = xvb_ua(ua_DivLambda_xvb(xvb_data1)),                 &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  if ( kp_kl(1) /= 0 ) &
       xvb_ddata(:,:,kp_kl(1)) = (2-3*xv_mu**2)*4*xv_mu
  if ( kp_kl(2) /= 0 ) &
       xvb_ddata(:,:,kp_kl(2)) = -3.0D0*xv_mu*(1-xv_mu**2)**(1.0D0/2)*cos(xv_Lon)
  call AssertEqual(&
    message='ua_DivMu_xvb with Y_2^0 (1-μ^2) and Y^1_1 √(1-μ^2)',&
    answer = xvb_ua(ua_DivMu_xvb(xvb_data1)),                     &
    check = xvb_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  
  xv_data1 = (3*xv_mu**2-1)*(1-xv_mu**2)       ! Y_2^0 (1-μ^2)
  xv_ddata = 0.0
  call AssertEqual(&
    message='u_DivLambda_xv with Y_2^0 (1-μ^2)',                &
    answer = xv_u(u_DivLambda_xv(xv_data1)),                     &
    check = xv_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore    &
    )
  xv_ddata = (2-3*xv_mu**2)*4*xv_mu
  call AssertEqual(&
    message='u_DivMu_xv with Y_2^0 (1-μ^2)',                    &
    answer = xv_u(u_DivMu_xv(xv_data1)),                         &
    check = xv_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore    &
    )

  xv_data1 = (1-xv_mu**2)**(3.0/2)*cos(xv_Lon) ! Y_1^1 (1-μ^2)
  xv_ddata = -(1-xv_mu**2)**(1.0d0/2)*sin(xv_Lon)
  call AssertEqual(&
    message='u_DivLambda_xv with Y^1_1 √(1-μ^2)',              &
    answer = xv_u(u_DivLambda_xv(xv_data1)),                     &
    check = xv_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore    &
    )
  xv_ddata = -3.0D0*xv_mu*(1-xv_mu**2)**(1.0D0/2)*cos(xv_Lon)
  call AssertEqual(&
    message='u_DivMu_xv with  Y^1_1 √(1-μ^2)',                 &
    answer = xv_u(u_DivMu_xv(xv_data1)),                         &
    check = xv_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore    &
    )

  
  call MessageNotify('M','ua_deriv_mpi_module_sypack_test', &
                         'ua_deriv_mpi_module_sypack function tests succeeded!')

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

  call MPI_FINALIZE(IERR)

end program ua_deriv_mpi_module_sypack_test
