!----------------------------------------------------------------------
!     Copyright (c) 2009 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  wt_zonal_module テストプログラム
!
!      トロイダルポロイダル関係微分関数のテスト
!  
!履歴  2010/04/18  竹広真一   wt_module_sjpack_deriv_test2.f90 より改造
!      2012/08/26  竹広真一   wt_module_deriv_test2.f90 より改造
!      2013/08/20  竹広真一   スペクトル配列の大きさ修正
!
program wt_zonal_module_derivative2

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

  integer,parameter  :: im=1, jm=16, km=16   ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: nm=10, lm=10         ! 切断波数の設定(水平, 動径)
  real(8),parameter  :: ri=0.5, ro=1.5       ! 内外半径

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

  real(8), dimension(0:im-1,jm,0:km)     :: xyz_DData
  real(8), dimension(0:im-1,jm,0:km)     :: xyz_Psi
  real(8), dimension(nm+1,0:lm)          :: wt_Psi
  real(8), dimension(nm+1,0:lm)          :: wt_Phi

  real(8), dimension(0:im-1,jm,0:km)     :: xyz_VLon
  real(8), dimension(0:im-1,jm,0:km)     :: xyz_VLat
  real(8), dimension(0:im-1,jm,0:km)     :: xyz_VRad
  real(8), dimension(0:im-1,jm,0:km)     :: xyz_RadRot
  real(8), dimension(0:im-1,jm,0:km)     :: xyz_RadRotRot

  integer, parameter :: n=2

  call MessageNotify('M','wt_zonal_module_derivative_test2', &
       'wt_zonal_module derivativeative function test #2')

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

  !========== wt_KxRGrad_wt, xyz_KGrad_wt, wt_QOperator_wt ==========

  ! ----------------- 例 1 --------------------
  xyz_Psi = xyz_rad**n * sin(xyz_lat)   ! r**2 Y_1^0

  xyz_DData = 0.0
  call AssertEqual(&
       message='wt_KxRGrad_wt with r^2 Y_1^-1',                      &
       answer = xyz_DData,                                           &
       check = xyz_wt(wt_KxRGrad_wt(wt_xyz(xyz_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! k ・▽ r**n Y_1^0
  xyz_DData =  n*xyz_rad**(n-1)*sin(xyz_Lat)**2 &
            + xyz_rad**(n-1)*cos(xyz_Lat)**2
  call AssertEqual(&
       message='xyz_KGrad_wt with r^2 Y_1^0',                        &
       answer = xyz_DData,                                           &
       check = xyz_KGrad_wt(wt_xyz(xyz_Psi)),                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! Q r**n Y_1^0
  xyz_DData = -(n-1)*xyz_rad**(n-1) * (3*sin(xyz_Lat)**2 -1)
  call AssertEqual(&
       message='wt_QOperator_wt r^2 Y_1^0',                          &
       answer = xyz_DData,                                           &
       check = xyz_wt(wt_QOperator_wt(wt_xyz(xyz_Psi))),             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

! ----------------- 例 2 --------------------
  xyz_Psi = xyz_Rad*(3.0D0*sin(xyz_Lat)**2 - 1.0D0)/2.0D0 ! r Y_2^0

  xyz_DData = 0.0D0
  call AssertEqual(&
       message='wt_KxRGrad_wt with Y_2^1',                           &
       answer = xyz_DData,                                           &
       check = xyz_wt(wt_KxRGrad_wt(wt_xyz(xyz_Psi))),               &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! k・▽ Y_2^0
  xyz_DData = -3.0D0/2.0D0*sin(xyz_Lat)**3 + 5.0/2.0*sin(xyz_Lat)
  call AssertEqual(&
       message='xyz_KGrad_wt with Y_2^0',                            &
       answer = xyz_DData,                                           &
       check = xyz_KGrad_wt(wt_xyz(xyz_Psi)),                        &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  xyz_DData = -9.0/2.0*sin(xyz_Lat)*cos(xyz_Lat)**2 &
              +15.0/2.0*sin(xyz_Lat)**3 - 15.0/2.0*sin(xyz_Lat)
  call AssertEqual(&
       message='wt_QOperator_wt Y_2^0',                              &
       answer = xyz_DData,                                           &
       check = xyz_wt(wt_QOperator_wt(wt_xyz(xyz_Psi))),             &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !========== wt_RadRot_xyz_xyz, wt_RadRotRot_xyz_xyz_xyz ==========

  !---------- 単純流れ v_r=r ----------
  xyz_VLon = 0 ; xyz_VLat = 0 ; xyz_VRad = xyz_Rad
  xyz_RadRot = 0 ; xyz_RadRotRot = 0

  call AssertEqual(&
       message='wt_RadRot_xyz_xyz with v_r=r',                       &
       answer = xyz_RadRot,                                          &
       check = xyz_wt(wt_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wt_RadRotRot_xyz_xyz with v_r=r',                    &
       answer = xyz_RadRotRot,                                       &
       check = xyz_wt(wt_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !---------- 剛体回転解非線形項(東西流れ) ----------
  xyz_VLon = 0
  xyz_VLat = xyz_Rad*sin(xyz_Lat)*cos(xyz_Lat)
  xyz_VRad = -xyz_Rad*cos(xyz_Lat)**2

  xyz_RadRot = 0 ; xyz_RadRotRot = 0

  call AssertEqual(&
       message='wt_RadRot_xyz_xyz with rigid rotation(E-W)',         &
       answer = xyz_RadRot,                                          &
       check = xyz_wt(wt_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wt_RadRotRot_xyz_xyz with rigid rotation(E-W)',      &
       answer = xyz_RadRotRot,                                       &
       check = xyz_wt(wt_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !--------- 鉛直渦度を伴うベクトル場 ----------
  xyz_Psi = xyz_Rad**2 * sin(xyz_Lat)                              ! r**2 P_1^0
  !xyz_Psi = xyz_Rad**2 * cos(xyz_Lat)*sin(xyz_Lat)*sin(xyz_Lon)   ! r**2 P_2^1

  xyz_VLon =   xyz_GradLat_wt(wt_xyz(xyz_Psi*xyz_Rad))
  xyz_VLat = - xyz_GradLon_wt(wt_xyz(xyz_Psi*xyz_Rad))
  xyz_VRad = 0
  xyz_RadRot = 2 * xyz_Psi                           ! r・▽×▽×(ψr) = L_2ψ
  !xyz_RadRot = 6 * xyz_Psi                           ! r・▽×▽×(ψr) = L_2ψ
  xyz_RadRotRot = 0

  call AssertEqual(&
      message='wt_RadRot_xyz_xyz with r^2 Y_1^0 toroidal field',     &
       answer = xyz_RadRot,                                          &
       check = xyz_wt(wt_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wt_RadRotRot_xyz_xyz with r^2 Y_1^0 toroidal field', &
       answer = xyz_RadRotRot,                                       &
       check = xyz_wt(wt_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  !---------- 鉛直速度を伴うベクトル場 -----------
  xyz_VRad = xyz_wt(wt_l2_wt(wt_xyz(xyz_Psi/xyz_Rad)))
  xyz_VLat = xyz_gRadLat_wt(wt_dRad_wt(wt_xyz(xyz_Psi*xyz_Rad)))
  xyz_VLon = xyz_gRadLon_wt(wt_dRad_wt(wt_xyz(xyz_Psi*xyz_Rad)))

  xyz_RadRot = 0
  xyz_RadRotRot = -xyz_wt(wt_l2_wt(wt_Lapla_wt(wt_xyz(xyz_Psi))))
                 ! r・▽×▽×▽×▽×(ψr) = -L_2▽^2ψ

  call AssertEqual(&
      message='wt_RadRot_xyz_xyz with non-vortical field',           &
       answer = xyz_RadRot,                                          &
       check = xyz_wt(wt_RadRot_xyz_xyz(xyz_VLon,xyz_VLat)),         &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call AssertEqual(&
       message='wt_RadRotRot_xyz_xyz with r^2 Y_1^-1 toroidal field',&
       answer = xyz_RadRotRot,                                       &
       check = xyz_wt(wt_RadRotRot_xyz_xyz_xyz(xyz_VLon,xyz_VLat,xyz_VRad)), &
       significant_digits = check_digits, ignore_digits = ignore     &
       )


  !========== wt_Potential2Vector, wt_Potential2Rotation ==========

  !----------------- 剛体回転場 --------------------
  ! 
  wt_Psi = wt_xyz(xyz_Rad * sin(xyz_Lat))
  wt_Phi = 0.0D0

  call wt_Potential2Vector(&
       xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi )

  xyz_DData = xyz_Rad * cos(xyz_Lat)
  call AssertEqual(&
      message='wt_Potential2Vector (Lon) with rigid rotation field(E-W)', &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wt_Potential2Vector (Lat) with rigid rotation field(E-W)', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wt_Potential2Vector (Rad) with rigid rotation field(E-W)', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wt_Potential2Rotation(&
       xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi )

  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wt_Potential2Rotation (Lon) with rigid rotation field(E-W)', &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 2*cos(xyz_LAT)
  call AssertEqual(&
      message='wt_Potential2Rotation (Lat) with rigid rotation field(E-W)', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 2*sin(xyz_LAT)
  call AssertEqual(&
      message='wt_Potential2Rotation (Rad) with rigid rotation field(E-W)', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  ! ----------------- 渦無し場 --------------------
  wt_Psi = 0.0D0
  wt_Phi = wt_xyz(xyz_Rad * sin(xyz_Lat))

  call wt_Potential2Vector(&
       xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi )

  xyz_DData = 0
  call AssertEqual(&
       message='wt_Potential2Vector (Lon) with non-vortical field',  &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 2 * cos(xyz_Lat)
  call AssertEqual(&
       message='wt_Potential2Vector (Lat) with non-vortical field',  &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 2 * sin(xyz_Lat)
  call AssertEqual(&
      message='wt_Potential2Vector (Rad) with non-vortical field', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call wt_Potential2Rotation(&
       xyz_VLon,xyz_VLat,xyz_VRad, wt_Psi, wt_Phi )

  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wt_Potential2Rotation (Lon) with non-vortical field', &
       answer = xyz_VLon,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wt_Potential2Rotation (Lat) with non-vortical field', &
       answer = xyz_VLat,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )
  xyz_DData = 0.0D0
  call AssertEqual(&
      message='wt_Potential2Rotation (Rad) with non-vortical field', &
       answer = xyz_VRad,                                            &
       check = xyz_DData,                                            &
       significant_digits = check_digits, ignore_digits = ignore     &
       )

  call MessageNotify('M','wt_zonal_module_derivative_test2', &
       'wt_zonal_module derivative function test #2 succeeded!')

end program wt_zonal_module_derivative2
