!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
! 関数のテスト
!    ggg_gradlon_yt, ggg_gradlat_yt, yt_div_ggg_ggg_ggg
!  
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_gradlon
  real(8), dimension(im,jm,0:km)     :: ggg_gradlat
  real(8), dimension(im,jm,0:km)     :: ggg_div
  real(8), dimension(im,jm,0:km)     :: ggg_data
  real(8), dimension(im,jm,0:km)     :: ggg_psi

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

  integer :: i,j,k

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

! ----------------- 例 1 --------------------
  ggg_vrad = 0
  ggg_psi = ggg_rad**n * cos(ggg_lat)*sin(ggg_lon)   ! r**2 P_1^1

  ggg_gradlon =  ggg_rad**(n-1)*cos(ggg_lon)
  ggg_gradlat = -ggg_rad**(n-1)*sin(ggg_lat)*sin(ggg_lon)

  ggg_div = - 2* ggg_psi/ggg_rad**2

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

! ----------------- 例 2 --------------------
  ggg_vrad = 0
  ggg_psi = ggg_rad**n * cos(ggg_lat)*sin(ggg_lat) * sin(ggg_lon) ! P_2^1

  ggg_gradlon =  ggg_rad**(n-1)*sin(ggg_lat)*cos(ggg_lon)
  ggg_gradlat =  ggg_rad**(n-1)*cos(2*ggg_lat)*sin(ggg_lon)

  ggg_div = - 6* ggg_psi/ggg_rad**2

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

  stop
contains

  subroutine checkresult

    ggg_vlon =  ggg_gradlon_yt(yt_ggg(ggg_psi))
    write(6,*)'Checking gradlon (1/r cos(lat) d/dlon)'
    do k=0,km
       do j=1,jm
          do i=1,im
             if ( abs(ggg_vlon(i,j,k)-ggg_gradlon(i,j,k)) > eps ) then
                write(6,*) i,j,k, ggg_vlon(i,j,k), ggg_gradlon(i,j,k)
             endif
          end do
       end do
    end do

    ggg_vlat =  ggg_gradlat_yt(yt_ggg(ggg_psi))
    write(6,*)'Checking gradlat (1/r dlat)'
    do k=0,km
       do j=1,jm
          do i=1,im
             if ( abs(ggg_vlat(i,j,k)-ggg_gradlat(i,j,k)) > eps ) then
                write(6,*) i,j,k, ggg_vlat(i,j,k), ggg_gradlat(i,j,k)
             endif
          end do
       end do
    end do

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

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

  end subroutine checkresult

end program spctest3
