!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!  yt_radrot_ggg_ggg, yt_radrotrot_ggg_ggg_ggg  のテスト
!    r・(▽xv), r・(▽x▽xv)
!  
program spctest3

  use spcpack
  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      ! 内外半径

  real(8), dimension(im,jm,0:km)     :: ggg_vlon
  real(8), dimension(im,jm,0:km)     :: ggg_vlat
  real(8), dimension(im,jm,0:km)     :: ggg_vrad
  real(8), dimension(im,jm,0:km)     :: ggg_data
  real(8), dimension(im,jm,0:km)     :: ggg_radrot
  real(8), dimension(im,jm,0:km)     :: ggg_radrotrot
  real(8), dimension((nm+1)**2,0:lm) :: yt_psi
  real(8), dimension((nm+1)**2,0:km) :: yg_psi
  real(8), dimension(im,jm,0:km)     :: ggg_psi

  real(8), parameter :: eps = 1D-10

  integer :: i,j,k

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

  ggg_vlon = 0 ; ggg_vlat = 0 ; ggg_vrad = ggg_rad
  ggg_radrot = 0 ; ggg_radrotrot = 0
  write(6,*)
  write(6,*)'Simple example (v_r=r)'
  call checkresult

! 剛体回転解(東西流れ)
  ggg_vlon = 0
  ggg_vlat = ggg_rad*sin(ggg_lat)*cos(ggg_lat)
  ggg_vrad = -ggg_rad*cos(ggg_lat)**2
  ggg_radrot = 0 ; ggg_radrotrot = 0
  write(6,*)
  write(6,*)'Rigid rotation '
  call checkresult

! 剛体回転解(南北流れ)
  ggg_vlon = ggg_rad*cos(ggg_lat)*sin(ggg_lon)*cos(ggg_lon)
  ggg_vlat = -ggg_rad*sin(ggg_lat)*cos(ggg_lat)*sin(ggg_lon)**2
  ggg_vrad = -ggg_rad*(sin(ggg_lat)**2*sin(ggg_lon)**2 + cos(ggg_lon)**2)
  ggg_radrot = 0 ; ggg_radrotrot = 0

  write(6,*)
  write(6,*)'Rigid rotation'
  call checkresult

! 鉛直渦度を伴うベクトル場
  ggg_psi = ggg_rad**2 * cos(ggg_lat)*sin(ggg_lon)   ! r**2 P_1^1
  !ggg_psi = ggg_rad**2 * cos(ggg_lat)*sin(ggg_lat)*sin(ggg_lon)   ! r**2 P_2^1

  ggg_vlon =   ggg_gradlat_yt(yt_ggg(ggg_psi*ggg_rad))
  ggg_vlat = - ggg_gradlon_yt(yt_ggg(ggg_psi*ggg_rad))
  ggg_vrad = 0
  ggg_radrot = 2 * ggg_psi                           ! r・▽×▽×(ψr) = L_2ψ
  !ggg_radrot = 6 * ggg_psi                           ! r・▽×▽×(ψr) = L_2ψ
  ggg_radrotrot = 0

  write(6,*)
  write(6,*)'Vortical field'
  call checkresult

  ggg_vrad = ggg_yt(yt_divlon_ggg(ggg_vlat)-yt_divlat_ggg(ggg_vlon))
  ggg_data = ggg_vlat
  ggg_vlat = ggg_yt(yt_drad_yt(yt_ggg(ggg_vlon*ggg_rad)))/ggg_rad
  ggg_vlon = -ggg_yt(yt_drad_yt(yt_ggg(ggg_data*ggg_rad)))/ggg_rad
  ggg_radrot = 0
  ggg_radrotrot = ggg_yt(yt_l2_yt(yt_lapla_yt(yt_ggg(ggg_psi))))
                 ! r・▽×▽×▽×▽×(ψr) = L_2▽^2ψ
  write(6,*)
  write(6,*)'Non-Vortical field'
  call checkresult

  stop
contains

  subroutine checkresult
    write(6,*)'Checking r rot v...'

    ggg_data = ggg_yt(yt_radrot_ggg_ggg(ggg_vlon,ggg_vlat))

    do k=0,km
       do j=1,jm
          do i=1,im
             if ( abs(ggg_data(i,j,k)-ggg_radrot(i,j,k)) > eps ) then
                write(6,*) i,j,k, ggg_data(i,j,k), ggg_radrot(i,j,k)
             endif
          end do
       end do
    end do

    ggg_data = ggg_yt(yt_radrotrot_ggg_ggg_ggg(ggg_vlon,ggg_vlat,ggg_vrad))

    write(6,*)'Checking r rot rot v...'
    do k=0,km
       do j=1,jm
          do i=1,im
             if ( abs(ggg_data(i,j,k)-ggg_radrotrot(i,j,k)) > eps ) then
                write(6,*) i,j,k, ggg_data(i,j,k), ggg_radrotrot(i,j,k)
             endif
          end do
       end do
    end do
  end subroutine checkresult
end program spctest3
