!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
! サブルーチンのテスト
!    wt_VGradV
!  
program wttest12

  use wt_module
  implicit none

  integer,parameter  :: im=32, jm=16, km=8  ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: nm=10, lm=8         ! 切断波数の設定(水平, 動径)
  real(8),parameter  :: ri=7.0D0/13.0D0      ! 内外半径 \eta=0.35
  integer,parameter  :: np=8                ! OPENMP 最大スレッド数
  real(8),parameter  :: ro=20.0D0/13.0D0     ! 内外半径 \eta=0.35

  real(8), dimension(im,jm,0:km)     :: xyz_VLon       ! 速度(経度)
  real(8), dimension(im,jm,0:km)     :: xyz_VLat       ! 速度(緯度)
  real(8), dimension(im,jm,0:km)     :: xyz_Vrad       ! 速度(動径)
  real(8), dimension(im,jm,0:km)     :: xyz_VGradVLon  ! 慣性項(経度)
  real(8), dimension(im,jm,0:km)     :: xyz_VGradVLat  ! 慣性項(緯度)
  real(8), dimension(im,jm,0:km)     :: xyz_VGradVrad  ! 慣性項(動径)
  real(8), dimension(im,jm,0:km)     :: xyz_VGradV0Lon  ! 慣性項(正解, 経度)
  real(8), dimension(im,jm,0:km)     :: xyz_VGradV0Lat  ! 慣性項(正解, 緯度)
  real(8), dimension(im,jm,0:km)     :: xyz_VGradV0rad  ! 慣性項(正解,動径)

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

  integer :: i,j,k

  call wt_Initial(im,jm,km,nm,lm,ri,ro,np)

  write( 6,* ) 'Test for VGradV'
  write( 6,* ) 'Output is displayed if computational error is larger than',eps

! ----------------- 例 1 --------------------
  ! 剛体回転場
  xyz_Vlon = xyz_Rad * cos(xyz_Lat)
  xyz_Vlat = 0.0
  xyz_Vrad = 0.0

  xyz_VGradV0lon = 0.0
  xyz_VGradV0lat = xyz_Rad*sin(xyz_LAT)*cos(xyz_LAT)
  xyz_VGradV0rad = -xyz_Rad*cos(xyz_LAT)**2

  write(6,*)
  write(6,*)
  write(6,*)'Example 1 : rigid rotation'
  write(6,*)'    xyz_Vlon = xyz_Rad * cos(xyz_Lat)'
  write(6,*)'    xyz_Vlat = 0.0'
  write(6,*)'    xyz_VRad = 0.0'

  call checkresult

! ----------------- 例 2 --------------------
  ! 剛体回転場(南北流)
  xyz_VLon = -xyz_Rad*sin(xyz_Lat)*sin(xyz_Lon)
  xyz_VLat = -xyz_Rad*cos(xyz_Lon)
  xyz_VRad = 0.0

  xyz_VGradV0Lon = xyz_Rad*cos(xyz_Lat)*sin(xyz_Lon)*cos(xyz_Lon)
  xyz_VGradV0Lat = -xyz_Rad*sin(xyz_Lat)*cos(xyz_Lat)*sin(xyz_Lon)**2
  xyz_VGradV0Rad = -xyz_Rad*(sin(xyz_Lat)**2*sin(xyz_Lon)**2 + cos(xyz_lon)**2)

  write(6,*)
  write(6,*)
  write(6,*)'Example 2 : rigid rotation'
  write(6,*)'    xyz_Vlon=xyz_Rad*xyz_Rad*sin(xyz_Lat)*sin(xyz_Lon)'
  write(6,*)'    xyz_Vlat=-xyz_Rad*cos(xyz_Lon)'
  write(6,*)'    xyz_VRad = 0.0'

  call checkresult

! ----------------- 例 3 --------------------
 ! ベンチマーク case1 の初期磁場
  xyz_Vrad = 5.0D0/8.0D0 * ( 8*ro - 6*xyz_Rad - 2*ri**4/xyz_Rad**3 ) &
               * cos(pi/2-xyz_Lat)
  xyz_Vlat = - 5.0D0/8.0D0 * ( 9*xyz_Rad - 8*ro - ri**4/xyz_Rad**3 ) &
               * sin(pi/2-xyz_Lat)
  xyz_Vlon = 5 * sin(pi*(xyz_Rad-ri)) * sin(2*(pi/2-xyz_Lat))

  stop
contains

 !------- 結果比較 -------
  subroutine checkresult

    call wt_VGradV3(xyz_VGradVLon,xyz_VGradVLat,xyz_VGradVRad, &
                    xyz_Vlon,xyz_Vlat,xyz_Vrad )

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

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

    write(6,*)
    write(6,*)'Checking VGradV_Rad '
    do k=0,km
       do j=1,jm
          do i=1,im
             if ( abs(xyz_VGradVRad(i,j,k)-xyz_VGradV0Rad(i,j,k)) > eps ) then
                write(6,*) i,j,k, xyz_VGradVRad(i,j,k),xyz_VGradV0rad(i,j,k)
             endif
          end do
       end do
    end do


  end subroutine checkresult

 !------- VGradV 計算(その0) -------
  subroutine wt_VGradV0(xyz_VGRADV_LON,xyz_VGRADV_LAT,xyz_VGRADV_RAD, &
                       xyz_VLON,xyz_VLAT,xyz_VRAD )

    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_LON
    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_LAT
    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_RAD
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VLON
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VLAT
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VRAD

    xyz_VGRADV_LON = xyz_wt(    &
            wt_Div_xyz_xyz_xyz( &
                xyz_VLON * xyz_VLON, xyz_VLON*xyz_VLAT, xyz_VLON*xyz_VRAD )) &
          + xyz_VLON*xyz_VRAD/xyz_RAD              &
          - xyz_VLON*xyz_VLAT*tan(xyz_LAT)/xyz_RAD 

    xyz_VGRADV_LAT = xyz_wt(    &
            wt_Div_xyz_xyz_xyz( &
                xyz_VLAT*xyz_VLON, xyz_VLAT*xyz_VLAT, xyz_VLAT*xyz_VRAD )) &
          + xyz_VLAT*xyz_VRAD/xyz_RAD        &
          + xyz_VLON**2*tan(xyz_LAT)/xyz_RAD 

    xyz_VGRADV_RAD = xyz_wt(    &
            wt_Div_xyz_xyz_xyz( &
                xyz_VRAD*xyz_VLON, xyz_VRAD*xyz_VLAT, xyz_VRAD*xyz_VRAD )) &
          - (xyz_VLON**2 + xyz_VLAT**2)/xyz_RAD 

  end subroutine wt_VGradV0

 !------- VGradV 計算(その1) -------
  subroutine wt_VGradV1(xyz_VGRADV_LON,xyz_VGRADV_LAT,xyz_VGRADV_RAD, &
                       xyz_VLON,xyz_VLAT,xyz_VRAD )

    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_LON
    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_LAT
    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_RAD
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VLON
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VLAT
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VRAD

    real(8), dimension((nm+1)*(nm+1),0:lm)   :: wt_VLon
    real(8), dimension((nm+1)*(nm+1),0:lm)   :: wt_VLat
    real(8), dimension((nm+1)*(nm+1),0:lm)   :: wt_VRad

    wt_VLon=wt_xyz(xyz_Vlon)
    wt_VLat=wt_xyz(xyz_Vlat)
    wt_VRad=wt_xyz(xyz_VRad)

    xyz_VGRADV_LON = xyz_VLon*xyz_GradLon_wt(wt_VLon)     &
                   + xyz_VLat*xyz_GradLat_wt(wt_VLon)     &
                   + xyz_VRad*xyz_wt(wt_Drad_wt(wt_VLon)) &
                   + xyz_VLON*xyz_VRAD/xyz_RAD            &
                   - xyz_VLON*xyz_VLAT*tan(xyz_LAT)/xyz_Rad 

    xyz_VGRADV_LAT = xyz_VLon*xyz_GradLon_wt(wt_VLat)     &
                   + xyz_VLat*xyz_GradLat_wt(wt_VLat)     &
                   + xyz_VRad*xyz_wt(wt_Drad_wt(wt_VLat)) &
                   + xyz_VLAT*xyz_VRAD/xyz_RAD            &
                   + xyz_VLON**2*tan(xyz_LAT)/xyz_RAD 

    xyz_VGRADV_RAD = xyz_VLon*xyz_GradLon_wt(wt_VRad)     &
                   + xyz_VLat*xyz_GradLat_wt(wt_VRad)     &
                   + xyz_VRad*xyz_wt(wt_Drad_wt(wt_VLat)) &
                   - (xyz_VLON**2 + xyz_VLAT**2)/xyz_RAD 

  end subroutine wt_VGradV1

 !------- VGradV 計算(その2) -------
  subroutine wt_VGradV2(xyz_VGRADV_LON,xyz_VGRADV_LAT,xyz_VGRADV_RAD, &
                       xyz_VLON,xyz_VLAT,xyz_VRAD )
   ! V・▽V =  ▽(1/2V**2) - Vx▽xV

    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_LON
    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_LAT
    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_RAD
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VLON
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VLAT
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VRAD

    real(8), dimension(im,jm,0:km)           :: xyz_RotVlon
    real(8), dimension(im,jm,0:km)           :: xyz_RotVlat
    real(8), dimension(im,jm,0:km)           :: xyz_RotVrad
    real(8), dimension((nm+1)*(nm+1),0:lm)   :: wt_Ke

    wt_Ke = wt_xyz(1/2.0*(xyz_VLon**2 + xyz_VLat**2 + xyz_VRad**2))

    xyz_RotVLon=xyz_GradLat_wt(wt_xyz(xyz_VRad))-xyz_wt(wt_RotRad_wt(wt_xyz(xyz_VLat)))
    xyz_RotVLat=xyz_wt(wt_RotRad_wt(wt_xyz(xyz_VLon))-wt_DivLon_xyz(xyz_VRad))

    xyz_RotVRad=xyz_wt(wt_DivLon_xyz(xyz_VLat)-wt_DivLat_xyz(xyz_VLon))

    xyz_VGRADV_LON = xyz_GradLon_wt(wt_Ke) &
                     - (xyz_VLat*xyz_RotVRad- xyz_VRad*xyz_RotVLat)

    xyz_VGRADV_LAT = xyz_GradLat_wt(wt_Ke) &
                     - (xyz_VRad*xyz_RotVLon- xyz_VLon*xyz_RotVRad)

    xyz_VGRADV_RAD = xyz_wt(wt_DRad_wt(wt_Ke)) &
                     - (xyz_VLon*xyz_RotVLat- xyz_VLat*xyz_RotVLon)

  end subroutine wt_VGradV2

  subroutine wt_VGradV3(xyz_VGRADV_LON,xyz_VGRADV_LAT,xyz_VGRADV_RAD, &
                        xyz_VLON,xyz_VLAT,xyz_VRAD )

    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_LON
    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_LAT
    real(8), dimension(im,jm,0:km),intent(out)   :: xyz_VGRADV_RAD
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VLON
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VLAT
    real(8), dimension(im,jm,0:km),intent(in)    :: xyz_VRAD

    xyz_VGRADV_LON = &
            xyz_Div_xyz_xyz_xyz( &
                xyz_VLON * xyz_VLON, xyz_VLON*xyz_VLAT, xyz_VLON*xyz_VRAD ) &
          + xyz_VLON*xyz_VRAD/xyz_RAD              &
          - xyz_VLON*xyz_VLAT*tan(xyz_LAT)/xyz_RAD 

    xyz_VGRADV_LAT = &
            xyz_Div_xyz_xyz_xyz( &
                xyz_VLAT*xyz_VLON, xyz_VLAT*xyz_VLAT, xyz_VLAT*xyz_VRAD ) &
          + xyz_VLAT*xyz_VRAD/xyz_RAD        &
          + xyz_VLON**2*tan(xyz_LAT)/xyz_RAD 

    xyz_VGRADV_RAD = &
            xyz_Div_xyz_xyz_xyz( &
                xyz_VRAD*xyz_VLON, xyz_VRAD*xyz_VLAT, xyz_VRAD*xyz_VRAD ) &
          - (xyz_VLON**2 + xyz_VLAT**2)/xyz_RAD 

  end subroutine wt_VGradV3

end program wttest12

