!----------------------------------------------------------------------
!     Copyright (c) 2011 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wa_deriv_module_sjpack_cuda テストプログラム :: 微分関数のテスト
!
!履歴  2009/09/06  竹広真一 wa_deriv_module_sjpack_cuda_test.f90 より CUDA 用に改造
!
program wa_deriv_module_sjpack_cuda_test

  use dc_message, only : MessageNotify
  use dc_test, only : AssertEqual
  use wa_module_sjpack_cuda
  implicit none

  integer, parameter :: im=512, jm=256, nm=170, km=2

  real(8), dimension(0:im-1,1:jm,km)     ::  xya_data1              ! 元の関数
  real(8), dimension(0:im-1,1:jm,km)     ::  xya_data2              ! 元の関数
  real(8), dimension(0:im-1,1:jm,km)     ::  xya_ddata              ! 微分の正解
  real(8), dimension(0:im-1,1:jm)        ::  xy_mu                  ! μ=sinφ

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

  call MessageNotify('M','wa_deriv_module_sjpack_cuda_test', &
                         'wa_deriv_module_sjpack_cuda function tests') 

  call wa_Initial( nm, im, jm, km )

  !---- Y_1^-1 Y_2^1 のテスト ----
  xya_data1(:,:,1) = -cos(xy_Lat)*sin(xy_Lon)             ! Y_1^{-1}
  xya_data1(:,:,2)  = sin(xy_Lat)*cos(xy_Lat) * cos(xy_Lon) ! Y_2^1

  xya_ddata(:,:,1) = 2*cos(xy_Lat)*sin(xy_Lon)             ! wa_Lapla_wa
  xya_ddata(:,:,2) = -6*sin(xy_Lat)*cos(xy_Lat) * cos(xy_Lon) ! Y_2^1
  call AssertEqual(&
    message='wa_Lapla_wa with Y_1^-1 and Y^2_1',                  &
    answer = xya_wa(wa_Lapla_wa(wa_xya(xya_data1))),              &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = 1.0/2.0*cos(xy_Lat)*sin(xy_Lon)      ! wa_LaplaInv_wa
  xya_ddata(:,:,2) = -1.0D0/6.0*sin(xy_Lat)*cos(xy_Lat) * cos(xy_Lon)
  call AssertEqual(&
    message='wa_LaplaInv_wa with Y_1^-1 and Y^2_1',               &
    answer = xya_wa(wa_LaplaInv_wa(wa_xya(xya_data1))),           &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = -cos(xy_Lat)*cos(xy_Lon)         ! wa_DLon_wa
  xya_ddata(:,:,2) = -sin(xy_Lat)*cos(xy_Lat) * sin(xy_Lon)
  call AssertEqual(&
    message='wa_DLon_wa with Y_1^-1 and Y^2_1',                   &
    answer = xya_wa(wa_DLon_wa(wa_xya(xya_data1))),               &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = -cos(xy_Lon)                     ! xya_GradLon_wa
  xya_ddata(:,:,2) = -sin(xy_Lat) * sin(xy_Lon)
  call AssertEqual(&
    message='xya_GradLon_wa with Y_1^-1 and Y^2_1',               &
    answer = xya_GradLon_wa(wa_xya(xya_data1)),                   &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = sin(xy_Lat)*sin(xy_Lon)         ! xya_GradLat_wa
  xya_ddata(:,:,2) = cos(2*xy_Lat) * cos(xy_Lon)
  call AssertEqual(&
    message='xya_GradLat_wa with Y_1^-1 and Y^2_1',               &
    answer = xya_GradLat_wa(wa_xya(xya_data1)),                   &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  !---- Y_1^-1 cosφ, Y_2^1 cosφ のテスト ----
  xya_data1(:,:,1) = -cos(xy_Lat)**2*sin(xy_Lon)              ! Y_1^-1 cosφ
  xya_data1(:,:,2) = sin(xy_Lat)*cos(xy_Lat)**2 * cos(xy_Lon) ! Y_2^1 cosφ

  xya_ddata(:,:,1) = -cos(xy_Lat)*cos(xy_Lon)                ! wa_DivLon_xy
  xya_ddata(:,:,2) = -sin(xy_Lat)*cos(xy_Lat)*sin(xy_Lon)    ! wa_DivLon_xy
  call AssertEqual(&
    message='wa_DivLon_xya with Y_1^-1 cosφ and Y^2_1 cosφ',    &
    answer = xya_wa(wa_DivLon_xya(xya_data1)),                    &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = 3*sin(xy_Lat)*cos(xy_Lat)*sin(xy_Lon)        !wa_DivLat_wa
  xya_ddata(:,:,2) = cos(xy_Lat)*(1-4*sin(xy_Lat)**2)*cos(xy_Lon) !wa_DivLat_wa
  call AssertEqual(&
    message='wa_DivLat_xya with Y_1^-1 cosφ and Y^2_1 cosφ',    &
    answer = xya_wa(wa_DivLat_xya(xya_data1)),                    &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  !---- Jacobian のテスト ----
  xya_data1(:,:,1) = -cos(xy_Lat)*sin(xy_Lon)                ! Y_1^{-1}
  xya_data1(:,:,2) = -cos(xy_Lat)*sin(xy_Lon)                ! Y_1^{-1}
  xya_data2(:,:,1) = -cos(xy_Lat)*sin(xy_Lon)                ! Y_1^{-1}
  xya_data2(:,:,2) = sin(xy_Lat)*cos(xy_Lat) * cos(xy_Lon)   ! Y_2^1

  xya_ddata(:,:,1) = 0.0
  xya_ddata(:,:,2) = sin(xy_Lat)**2 - cos(xy_Lat)**2*cos(xy_Lon)**2
   call AssertEqual(&
    message='wa_Jacobian_wa_wa Y_1^-1 and Y^2_1',                 &
    answer = xya_wa(wa_Jacobian_wa_wa(wa_xya(xya_data1),wa_xya(xya_data2))), &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )


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

  !----- Y_2^0, Y_1^1 のテスト -----
  xya_data1(:,:,1) = 3*xy_mu**2-1                              ! Y_2^0
  xya_data1(:,:,2) = sqrt(1-xy_mu**2)*cos(xy_Lon)                  ! Y_1^1

  xya_ddata(:,:,1) = 0.0
  xya_ddata(:,:,2) = -sqrt(1-xy_mu**2)*sin(xy_Lon)
  call AssertEqual(&
    message='xya_GradLambda_wa with Y_2^0 and Y^1_1',             &
    answer = xya_GradLambda_wa(wa_xya(xya_data1)),                &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = 6*xy_mu*(1-xy_mu**2)
  xya_ddata(:,:,2) = -xy_mu*sqrt(1-xy_mu**2)*cos(xy_Lon)
  call AssertEqual(&
    message='xya_GradMu_wa with Y_2^0 and Y^1_1',                 &
    answer = xya_GradMu_wa(wa_xya(xya_data1)),                    &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  !----- Y_2^0(1-μ^2), Y_1^1 (1-μ^2) のテスト -----
  xya_data1(:,:,1) = (3*xy_mu**2-1)*(1-xy_mu**2)       ! Y_2^0 (1-μ^2)
  xya_data1(:,:,2) = (1-xy_mu**2)**(3.0/2)*cos(xy_Lon) ! Y_1^1 (1-μ^2)

  xya_ddata(:,:,1) = 0.0
  xya_ddata(:,:,2) = -(1-xy_mu**2)**(1.0d0/2)*sin(xy_Lon)
  call AssertEqual(&
    message='wa_DivLambda_xya with Y_2^0 (1-μ^2) and Y^1_1 √(1-μ^2)',&
    answer = xya_wa(wa_DivLambda_xya(xya_data1)),                 &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  xya_ddata(:,:,1) = (2-3*xy_mu**2)*4*xy_mu
  xya_ddata(:,:,2) = -3.0D0*xy_mu*(1-xy_mu**2)**(1.0D0/2)*cos(xy_Lon)
  call AssertEqual(&
    message='wa_DivMu_xya with Y_2^0 (1-μ^2) and Y^1_1 √(1-μ^2)',&
    answer = xya_wa(wa_DivMu_xya(xya_data1)),                     &
    check = xya_ddata,                                            &
    significant_digits = check_digits, ignore_digits = ignore     &
    )

  call MessageNotify('M','wa_deriv_module_sjpack_cuda_test', &
                         'wa_deriv_module_sjpack_cuda function tests succeeded!') 

end program wa_deriv_module_sjpack_cuda_test

