!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
! 関数のテスト
!    yt_divlon_ggg, yt_divlat_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_divlon
  real(8), dimension(im,jm,0:km)     :: ggg_divlat
  real(8), dimension(im,jm,0:km)     :: ggg_data

  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_vlon   = ggg_rad**n * cos(ggg_lat)**2*sin(ggg_lon)
  ggg_divlon = ggg_rad**(n-1)*cos(ggg_lat)*cos(ggg_lon)

  ggg_vlat   = ggg_rad**n * cos(ggg_lat)**2*sin(ggg_lon)
  ggg_divlat = -3*ggg_rad**(n-1)*cos(ggg_lat)*sin(ggg_lat)*sin(ggg_lon)

  write(6,*)
  write(6,*)'Example 1'
  call checkresult

  stop
contains

  subroutine checkresult

    write(6,*)'Checking divlon (1/r cos(lat) dlon)'
    ggg_data =  ggg_yt(yt_divlon_ggg(ggg_vlon))

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

    ggg_data =  ggg_yt(yt_divlat_ggg(ggg_vlat))
    write(6,*)'Checking divlat (1/r cos(lat) dlat cos(lat))'
    do k=0,km
       do j=1,jm
          do i=1,im
             if ( abs(ggg_data(i,j,k)-ggg_divlat(i,j,k)) > eps ) then
                write(6,*) i,j,k, ggg_data(i,j,k), ggg_divlat(i,j,k)
             endif
          end do
       end do
    end do
  end subroutine checkresult


end program spctest3
