function AvrLonLatRad_xyz(xyz) ! 緯度経度動径(全球)積分
      real(8), dimension(im,jm,0:km), intent(in) :: xyz    ! 3 次元格子点データ
      real(8)                     :: AvrLonLatRad_xyz      ! 全球平均

      AvrLonLatRad_xyz = IntLonLatRad_xyz(xyz)             /(sum(x_Lon_Weight)*sum(y_Lat_Weight) * sum(z_Rad_Weight))

    end function AvrLonLatRad_xyz

    function IntLonLatRad_xyz(xyz) ! 緯度経度動径(全球)積分
      real(8), dimension(im,jm,0:km), intent(in) :: xyz    ! 3 次元格子点データ
      real(8)                     :: IntLonLatRad_xyz      ! 全球積分値
      integer :: i, j, k

      IntLonLatRad_xyz = 0
      do k=0,km
         do j=1,jm
            do i=1,im
               IntLonLatRad_xyz = IntLonLatRad_xyz                     + xyz(i,j,k) * x_Lon_Weight(i)                          * y_Lat_Weight(j) * z_Rad_Weight(k)
            enddo
         enddo
      enddo
    end function IntLonLatRad_xyz

    function nmz_PoloidalEnergySpectrum_wt(wt_POLPOT)
      real(8), dimension((nm+1)*(nm+1),0:km), intent(in)            :: wt_POLPOT                      ! ポロイダルポテンシャル
      real(8), dimension(0:nm,-nm:nm,0:km)            :: nmz_PoloidalEnergySpectrum_wt   ! スペクトルトロイダル成分


      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA1   ! 作業領域
      real(8), dimension((nm+1)*(nm+1),0:km) ::wz_DATA2   ! 作業領域
      integer :: n, m

      nmz_PoloidalEnergySpectrum_wt = wt_VMiss

      wz_Data1 = wz_wt(wt_POLPOT)
      wz_Data2 = wz_Rad*wz_wt(wt_DRad_wt(wt_POLPOT))                + wz_wt(wt_POLPOT)                         ! = rdφ/dr+φ

      do n=0,nm
         do m=-n,n
            nmz_PoloidalEnergySpectrum_wt(n,m,:) =                  + 0.5* n*(n+1)* (4*pi)                  *( wz_Data2(l_nm(n,m),:)**2                     + n*(n+1)*wz_Data1(l_nm(n,m),:)**2 )
         enddo
      enddo

    end function nmz_PoloidalEnergySpectrum_wt

    function nmz_ToroidalEnergySpectrum_wt(wt_TORPOT)
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in)            :: wt_TORPOT                       ! トロイダルポテンシャル
      real(8), dimension(0:nm,-nm:nm,0:km)            :: nmz_ToroidalEnergySpectrum_wt   ! スペクトルトロイダル成分

      real(8), dimension((nm+1)*(nm+1),0:lm) ::wz_DATA   ! 作業領域
      integer :: n, m

      nmz_ToroidalEnergySpectrum_wt = wt_VMiss

      wz_DATA = wz_wt(wt_TORPOT)
      do n=0,nm
         do m=-n,n
            nmz_ToroidalEnergySpectrum_wt(n,m,:)               = 0.5 * n*(n+1)* (4*pi) * z_Rad**2                 * wz_DATA(l_nm(n,m),:)**2
         enddo
      enddo

    end function nmz_ToroidalEnergySpectrum_wt

    subroutine wt_BoundariesTau(wt,values,cond)   ! ディリクレ, ノイマン条件
      ! Chebyshev 空間での境界条件適用(タウ法)

      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)      :: wt
              ! 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension((nm+1)*(nm+1),2), intent(in), optional :: values
              ! 境界での 値/勾配 分布を水平スペクトル変換したものを与える. 
              ! 省略時は値/勾配 0 となる. 

      character(len=2), intent(in), optional             :: cond
              ! 境界条件. 省略時は 'DD'
              !   DD : 両端ディリクレ
              !   DN,ND : ディリクレ/ノイマン条件
              !   NN : 両端ノイマン

      if (.not. present(cond)) then
         if (present(values)) then
            call at_BoundariesTau_DD(wt,values)
         else
            call at_BoundariesTau_DD(wt)
         endif
         return
      endif

      select case(cond)
      case ('NN')
         if (present(values)) then
            call at_BoundariesTau_NN(wt,values)
         else
            call at_BoundariesTau_NN(wt)
         endif
      case ('DN')
         if (present(values)) then
            call at_BoundariesTau_DN(wt,values)
         else
            call at_BoundariesTau_DN(wt)
         endif
      case ('ND')
         if (present(values)) then
            call at_BoundariesTau_ND(wt,values)
         else
            call at_BoundariesTau_ND(wt)
         endif
      case ('DD')
         if (present(values)) then
            call at_BoundariesTau_DD(wt,values)
         else
            call at_BoundariesTau_DD(wt)
         endif
      case default
         call MessageNotify('E','wt_BoundariesTau','B.C. not supported')
      end select

    end subroutine wt_BoundariesTau

    subroutine wt_BoundariesGrid(wt,values,cond)   ! ディリクレ, ノイマン条件
      ! 実空間での境界条件適用

      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)      :: wt
              ! 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension((nm+1)*(nm+1),2), intent(in), optional :: values
              ! 境界での 値/勾配 分布を水平スペクトル変換したものを与える. 
              ! 省略時は値/勾配 0 となる. 

      character(len=2), intent(in), optional             :: cond
              ! 境界条件. 省略時は 'DD'
              !   DD : 両端ディリクレ
              !   DN,ND : ディリクレ/ノイマン条件
              !   NN : 両端ノイマン

      if (.not. present(cond)) then
         if (present(values)) then
            call at_boundariesGrid_DD(wt,values)
         else
            call at_boundariesGrid_DD(wt)
         endif
         return
      endif

      select case(cond)
      case ('NN')
         if (present(values)) then
            call at_BoundariesGrid_NN(wt,values)
         else
            call at_BoundariesGrid_NN(wt)
         endif
      case ('DN')
         if (present(values)) then
            call at_BoundariesGrid_DN(wt,values)
         else
            call at_BoundariesGrid_DN(wt)
         endif
      case ('ND')
         if (present(values)) then
            call at_BoundariesGrid_ND(wt,values)
         else
            call at_BoundariesGrid_ND(wt)
         endif
      case ('DD')
         if (present(values)) then
            call at_BoundariesGrid_DD(wt,values)
         else
            call at_BoundariesGrid_DD(wt)
         endif
      case default
         call MessageNotify('E','wt_BoundariesGrid','B.C. not supported')
      end select

    end subroutine wt_BoundariesGrid

    function wt_DRad_wt(wt)     ! 動径微分
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt
      real(8), dimension((nm+1)*(nm+1),0:lm)             :: wt_DRad_wt

      wt_DRad_wt = at_Dr_at(wt)

    end function wt_DRad_wt

    function wt_DivLon_xyz(xyz)   ! 格子に作用する発散型経度微分 
                                  ! 1/rcosφ・∂/∂λ

      real(8), dimension(im,jm,0:km), intent(in)   :: xyz
      real(8), dimension((nm+1)*(nm+1),0:lm)       :: wt_DivLon_xyz

      wt_DivLon_xyz = wt_wz(wa_DivLon_xya(xyz/xyz_Rad))
    end function wt_DivLon_xyz

    function wt_Div_xyz_xyz_xyz(xyz_Vlon,xyz_Vlat,xyz_Vrad) ! 発散

      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_Div_xyz_xyz_xyz

      wt_Div_xyz_xyz_xyz =   wt_DivLon_xyz(xyz_Vlon)                            + wt_DivLat_xyz(xyz_Vlat)                            + wt_DivRad_wt(wt_xyz(xyz_Vrad))

    end function wt_Div_xyz_xyz_xyz

    subroutine wt_Initial(i,j,k,n,l,r_in,r_out)

     integer,intent(in) :: i, j, k        ! 格子点の設定(経度, 緯度, 動径)
     integer,intent(in) :: n, l           ! 切断波数の設定(水平, 動径)

     real(8),intent(in) :: r_in, r_out    ! 球殻内外半径

     im = i  
 jm = j 
 km = k
     nm = n  
 lm = l
     ri = r_in 
 ro = r_out

     call wa_Initial(nm,im,jm,km+1)
     call at_Initial(km,lm,r_in,r_out)

     allocate(xyz_Lon(im,jm,0:km))
     allocate(xyz_Lat(im,jm,0:km))
     allocate(xyz_Rad(im,jm,0:km))

     allocate(wz_Rad((nm+1)*(nm+1),0:km))

     xyz_Lon = spread(xy_Lon,3,km+1)
     xyz_Lat = spread(xy_Lat,3,km+1)
     xyz_Rad = spread(spread(z_Rad,1,jm),1,im)

     wz_Rad = spread(z_Rad,1,(nm+1)*(nm+1))

     z_Rad_Weight = z_Rad_Weight * z_Rad**2       ! r^2 dr の積分重み
   end subroutine wt_initial

    function wt_KxRGrad_wt(wt)    ! k×r・▽ = ∂/∂λ
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt
      real(8), dimension((nm+1)*(nm+1),0:lm)             :: wt_KxRGrad_wt

      wt_KxRGrad_wt =  wa_Dlon_wa(wt)
    end function wt_KxRGrad_wt

    subroutine wt_PolmagBoundariesTau(wt_POL) ! 磁場ポロイダルポテンシャル境界
      ! Chebyshev 空間での境界条件適用
      ! 磁場ポロイダルポテンシャル \psi の境界条件(非電気伝導体)
      ! 外側
      !    \DP{wt_psi}{r} + (l+1) wt_psi/r = 0   at the outer boundary
      ! 内側
      !    \DP{wt_psi}{r} - l wt_psi/r = 0       at the inner boundary

      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)   :: wt_POL
              ! 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(8), dimension(:,:), allocatable    :: tt_I
      real(8), dimension(:,:), allocatable    :: tz_PSI
      real(8), dimension(:,:), allocatable    :: tz_DPSIDR

      logical :: first = .true.
      integer  :: l, n, nn(2)
      save     :: alu, kp, first

      if ( first ) then
         first = .false.

         allocate(alu((nm+1)*(nm+1),0:lm,0:lm),kp((nm+1)*(nm+1),0:lm))
         allocate(tt_I(0:lm,0:lm),tz_PSI(0:lm,0:km),tz_DPSIDR(0:lm,0:km))

         tt_I = 0.0
         do l=0,lm
            tt_I(l,l)=1
         enddo
         do n=1,(nm+1)*(nm+1)
            alu(n,:,:) = tt_I
         enddo

         ! 非電気伝導体
         tz_PSI = az_at(tt_I)
         tz_DPSIDR = az_at(at_dr_at(tt_I))

         do n=1,(nm+1)*(nm+1)
            nn=nm_l(n)
            alu(n,lm-1,:) = tz_DPSIDR(:,0) + (nn(1)+1) * tz_PSI(:,0)/z_RAD(0)
            alu(n,lm,:)   = tz_DPSIDR(:,km) - nn(1) * tz_PSI(:,km)/z_RAD(km)
         enddo
         call ludecomp(alu,kp)

         deallocate(tt_I,tz_PSI,tz_DPSIDR)
      endif

      wt_POL(:,lm-1) = 0.0
      wt_POL(:,lm)   = 0.0
      wt_POL = lusolve(alu,kp,wt_POL)

    end subroutine wt_PolmagBoundariesTau

    subroutine wt_Potential2Vector(         xyz_VLON,xyz_VLAT,xyz_VRAD,wt_TORPOT,wt_POLPOT)

      ! トロイダルポロイダルポテンシャルΨ,Φで表される非発散ベクトル場
      !     v = ▽x(Ψr) + ▽x▽x(Φr) 
      ! の各成分を計算する

      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((nm+1)*(nm+1),0:lm), intent(in)                                          :: wt_TORPOT ! トロイダルポテンシャル
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in)                                          :: wt_POLPOT ! ポロイダルポテンシャル

      xyz_VLON =   xyz_RAD * xyz_GradLat_wt(wt_TORPOT)                  + xya_GradLon_wa(wz_wt(wt_RotRad_wt(wt_POLPOT)))
      xyz_VLAT = - xyz_RAD * xyz_GradLon_wt(wt_TORPOT)                  + xya_GradLat_wa(wz_wt(wt_RotRad_wt(wt_POLPOT)))
      xyz_VRAD = xyz_wt(wt_L2_wt(wt_POLPOT))/xyz_RAD

    end subroutine wt_Potential2Vector

    function wt_RadRot_xyz_xyz(xyz_VLON,xyz_VLAT)  ! r・(▽×v)
      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((nm+1)*(nm+1),0:lm)     :: wt_RadRot_xyz_xyz

      wt_RadRot_xyz_xyz = wt_wz(wa_DivLon_xya(xyz_VLAT)                                 - wa_DivLat_xya(xyz_VLON))
      
    end function wt_RadRot_xyz_xyz

    subroutine wt_TorBoundariesTau(wt_TORPOT,values,cond) 
      ! トロイダルポテンシャル境界
      ! Chebyshev 空間での境界条件適用
      !  トロイダルポテンシャル \psi の境界条件を適用
      ! 粘着条件
      !    \psi = values  at boundaries (default 0)
      ! 応力なし条件
      !    \DP{\psi/r}{r} = 0  at boundaries

      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)      :: wt_TORPOT
              ! 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension((nm+1)*(nm+1),2), intent(in), optional :: values
              ! 両端境界でのトロイダルポテンシャル
              ! 粘着条件の時のみ有効

      character(len=2), intent(in), optional  :: cond
              ! 境界条件スイッチ. 省略時は 'RR'
              !   RR    : 両端粘着条件
              !   RF/FR : 粘着/応力なし条件
              !   FF    : 両端応力なし条件

      real(8), dimension(:,:), allocatable  :: alu
      integer, dimension(:), allocatable    :: kp
      real(8), dimension(0:lm,0:lm)         :: tt_data
      real(8), dimension(0:lm,0:km)         :: tg_data
      logical                               :: rigid1, rigid2   ! 境界条件

      logical :: first = .true.
      integer  :: l
      save     :: alu, kp, first

      if (.not. present(cond)) then
         rigid1=.TRUE. 
 rigid2=.TRUE.
      else
         select case (cond)
         case ('RR')
            rigid1 = .TRUE.  
 rigid2 = .TRUE.
         case ('RF')
            rigid1 = .TRUE.  
 rigid2 = .FALSE.
         case ('FR')
            rigid1 = .FALSE. 
 rigid2 = .TRUE.
         case ('FF')
            rigid1 = .FALSE. 
 rigid2 = .FALSE.
         case default
            call MessageNotify('E','wt_TorBoundariesTau','B.C. not supported')
         end select
      endif

      if ( first ) then
         first = .false.

         allocate(alu(0:lm,0:lm),kp(0:lm))

         tt_data = 0
         do l=0,lm
            tt_data(l,l)=1
         enddo
         alu = tt_data

         ! 力学的条件粘着条件 
         if ( rigid1 ) then
            tg_data = az_at(tt_data)
         else
            tg_data = az_at(at_dr_at(at_az(                  az_at(tt_data)/spread(z_rad,1,lm+1))))
         endif
         alu(lm-1,:) = tg_data(:,0)       ! 境界 k=0 での条件式代入

         if ( rigid2 ) then
            tg_data = az_at(tt_data)    
         else
            tg_data = az_at(at_dr_at(at_az(                  az_at(tt_data)/spread(z_rad,1,lm+1))))
         endif
         alu(lm,:)   = tg_data(:,km)      ! 境界 k=km での条件式代入

         call ludecomp(alu,kp)

         if ( rigid1 .AND. present(values) ) then 
            call MessageNotify('M','wt_TorBoundariesTau',                 'Toroidal potential at k=0 was given by the optional variable.')
         else if ( rigid1 .AND. (.NOT.present(values)) ) then
            call MessageNotify('M','wt_TorBoundariesTau',                 'Toroidal potential at k=0 was set to zero.')
         else if ( (.NOT. rigid1) .AND. present(values) ) then
            call MessageNotify('W','wt_TorBoundariesTau',                 'Boundary value k=0 cannot be set under stress-free condition.')
         endif

         if ( rigid2 .AND. present(values) ) then 
            call MessageNotify('M','wt_TorBoundariesTau',                 'Toroidal potential at k=0 was given by the optional variable.')
         else if ( rigid2 .AND. (.NOT.present(values)) ) then
            call MessageNotify('M','wt_TorBoundariesTau',                 'Toroidal potential at k=0 was set to zero.')
         else if ( (.NOT. rigid2) .AND. present(values) ) then
            call MessageNotify('W','wt_TorBoundariesTau',                 'Boundary value k=0 cannot be set under stress-free condition.')
         endif
      endif

      if ( rigid1 .AND. present(values) ) then
         wt_torpot(:,lm-1) = values(:,1)
      else
         wt_torpot(:,lm-1) = 0
      endif
      if ( rigid2 .AND. present(values) ) then
         wt_torpot(:,lm)   = values(:,2)
      else
         wt_torpot(:,lm) = 0
      endif

      wt_torpot = lusolve(alu,kp,wt_TORPOT)

    end subroutine wt_TorBoundariesTau

    subroutine wt_TormagBoundariesTau(wt_TOR) ! 磁場トロイダルポテンシャル境界
      ! Chebyshev 空間での境界条件適用
      ! 磁場トロイダルポテンシャル \psi の境界条件(非電気伝導体)
      ! 外側
      !    wt_psi = 0   at the outer boundary
      ! 内側
      !    wt_psi = 0       at the inner boundary
      ! 
      ! wt_Boundaries で対応可能だが, 将来のため別途作成しておく

      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)   :: wt_TOR
              ! 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(8), dimension(:,:), allocatable    :: tt_I
      real(8), dimension(:,:), allocatable    :: tz_PSI

      logical :: first = .true.
      integer  :: l, n
      save     :: alu, kp, first

      if ( first ) then
         first = .false.

         allocate(alu((nm+1)*(nm+1),0:lm,0:lm),kp((nm+1)*(nm+1),0:lm))
         allocate(tt_I(0:lm,0:lm),tz_PSI(0:lm,0:km))

         tt_I = 0.0
         do l=0,lm
            tt_I(l,l)=1
         enddo
         do n=1,(nm+1)*(nm+1)
            alu(n,:,:) = tt_I
         enddo

         ! 非電気伝導体
         tz_PSI = az_at(tt_I)

         do n=1,(nm+1)*(nm+1)
            alu(n,lm-1,:) = tz_PSI(:,0)
            alu(n,lm,:)   = tz_PSI(:,km)
         enddo
         call ludecomp(alu,kp)

         deallocate(tt_I,tz_PSI)
      endif

      wt_TOR(:,lm-1) = 0.0
      wt_TOR(:,lm)   = 0.0
      wt_TOR = lusolve(alu,kp,wt_TOR)

    end subroutine wt_TormagBoundariesTau

    subroutine wt_TormagBoundariesGrid(wt_TOR) ! 磁場トロイダルポテンシャル境界
      ! 鉛直実空間での境界条件適用
      ! 磁場トロイダルポテンシャル \psi の境界条件(非電気伝導体)
      ! 外側
      !    wt_psi = 0   at the outer boundary
      ! 内側
      !    wt_psi = 0       at the inner boundary
      ! 
      ! wt_Boundaries で対応可能だが, 将来のため別途作成しておく

      real(8), dimension((nm+1)*(nm+1),0:lm),intent(inout)   :: wt_TOR
              ! 境界条件を適用するデータ. 修正された値を返す. 

      real(8), dimension(:,:,:), allocatable  :: alu
      integer, dimension(:,:), allocatable    :: kp

      real(8), dimension(:,:), allocatable    :: tt_I
      real(8), dimension(:,:), allocatable    :: tz_PSI
      real(8), dimension((nm+1)*(nm+1),0:km)  :: wz_TOR

      logical :: first = .true.
      integer  :: l, n
      save     :: alu, kp, first

      if ( first ) then
         first = .false.

         if ( lm /= km ) then
            call MessageNotify('E','TorMagBoundariesGrid',              'Chebyshev truncation and number of grid points should be same.')
         endif

         allocate(alu((nm+1)*(nm+1),0:km,0:lm),kp((nm+1)*(nm+1),0:lm))
         allocate(tt_I(0:lm,0:lm),tz_PSI(0:lm,0:km))

         tt_I = 0.0
         do l=0,lm
            tt_I(l,l)=1
         enddo
         do n=1,(nm+1)*(nm+1)
            alu(n,:,:) = transpose(az_at(tt_I))   ! 内部領域は値そのまま.
         enddo

         ! 非電気伝導体
         tz_PSI = az_at(tt_I)

         do n=1,(nm+1)*(nm+1)
            alu(n,0,:) = tz_PSI(:,0)
            alu(n,km,:)   = tz_PSI(:,km)
         enddo
         call ludecomp(alu,kp)

         deallocate(tt_I,tz_PSI)
      endif
      
      wz_TOR       = wz_wt(wt_TOR)
      wz_TOR(:,0)  = 0.0
      wz_TOR(:,km) = 0.0
      wt_TOR = lusolve(alu,kp,wz_TOR)

    end subroutine wt_TormagBoundariesGrid

    subroutine wt_VGradV(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_VGradV

    function x_AvrLatRad_xyz(xyz)  ! 緯度動径(子午面)積分
      real(8), dimension(im,jm,0:km), intent(in) :: xyz    ! 3 次元格子点データ
      real(8), dimension(im)     :: x_AvrLatRad_xyz      ! 経度格子点データ

      x_AvrLatRad_xyz = x_IntLatRad_xyz(xyz)                    /( sum(y_Lat_Weight)*sum(z_Rad_Weight) )

    end function x_AvrLatRad_xyz

    function x_IntLatRad_xyz(xyz)  ! 緯度動径(子午面)積分
      real(8), dimension(im,jm,0:km), intent(in) :: xyz    ! 3 次元格子点データ
      real(8), dimension(im)     :: x_IntLatRad_xyz      ! 経度格子点データ
      integer :: j, k

      x_IntLatRad_xyz = 0
      do k=0,km
         do j=1,jm
            x_IntLatRad_xyz = x_IntLatRad_xyz                  + xyz(:,j,k) * y_Lat_Weight(j) * z_Rad_Weight(k)
         enddo
      enddo
    end function x_IntLatRad_xyz

    function xyz_GradLon_wt(wt) ! スペクトルに作用する勾配型経度微分
                                ! 1/rcosφ・∂/∂λ
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt
      real(8), dimension(im,jm,0:km)                     :: xyz_GradLon_wt

      xyz_GradLon_wt = xya_GradLon_wa(wz_wt(wt))/xyz_Rad

    end function xyz_GradLon_wt

    function xyz_RotLon_wt_wt(wt_Vrad,wt_Vlat) 
      ! ベクトル場の回転の経度成分
      ! 1/rcosφ・∂Vrad/∂λ-1/r ∂(r Vlat)/∂r

        real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt_Vrad ! 動径成分
        real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt_Vlat ! 緯度成分

        real(8), dimension(im,jm,0:km)                     :: xyz_RotLon_wt_wt

        xyz_RotLon_wt_wt =   xyz_GradLon_wt(wt_Vrad)                            - xyz_wt(wt_RotRad_wt(wt_Vlat))

    end function xyz_RotLon_wt_wt

    function xyz_wt(wt)  ! スペクトル -> 格子
      real(8), dimension(im,jm,0:km)                     :: xyz_wt
      real(8), dimension((nm+1)*(nm+1),0:lm), intent(in) :: wt

      xyz_wt = xya_wa(az_at(wt))

    end function xyz_wt

    function yz_AvrLon_xyz(xyz)  ! 経度(帯状)積分
      real(8), dimension(im,jm,0:km), intent(in) :: xyz    ! 3 次元格子点データ
      real(8), dimension(jm,0:km)  :: yz_AvrLon_xyz        ! 子午面格子点データ

      yz_AvrLon_xyz = yz_IntLon_xyz(xyz)/sum(x_Lon_Weight)

    end function yz_AvrLon_xyz

    function yz_IntLon_xyz(xyz)  ! 経度(帯状)積分
      real(8), dimension(im,jm,0:km), intent(in) :: xyz    ! 3 次元格子点データ
      real(8), dimension(jm,0:km)  :: yz_IntLon_xyz        ! 子午面格子点データ
      integer :: i

      yz_IntLon_xyz = 0.0d0
      do i=1,im
         yz_IntLon_xyz(:,:) = yz_IntLon_xyz(:,:)                        + xyz(i,:,:) * x_Lon_Weight(i)
      enddo
    end function yz_IntLon_xyz

    function z_AvrLat_yz(yz)  ! 緯度積分
      real(8), dimension(jm,0:km), intent(in) :: yz   ! 2 次元格子点データ
      real(8), dimension(0:km)  :: z_AvrLat_yz        ! 動径格子点データ

      z_AvrLat_yz = z_IntLat_yz(yz)/sum(y_Lat_Weight)
    end function z_AvrLat_yz

    function z_AvrLon_xz(xz)  ! 経度(帯状)積分
      real(8), dimension(im,0:km), intent(in) :: xz   ! 2 次元格子点データ
      real(8), dimension(0:km)  :: z_AvrLon_xz        ! 鉛直格子点データ

      z_AvrLon_xz = z_IntLon_xz(xz)/sum(x_Lon_Weight)

    end function z_AvrLon_xz

    function z_IntLat_yz(yz)  ! 緯度積分
      real(8), dimension(jm,0:km), intent(in) :: yz   ! 2 次元格子点データ
      real(8), dimension(0:km)  :: z_IntLat_yz        ! 動径格子点データ
      integer :: j

      z_IntLat_yz = 0.0d0
      do j=1,jm
         z_IntLat_yz(:) = z_IntLat_yz(:) + yz(j,:) * y_Lat_Weight(j)
      enddo
    end function z_IntLat_yz

    function z_IntLon_xz(xz)  ! 経度(帯状)積分
      real(8), dimension(im,0:km), intent(in) :: xz   ! 2 次元格子点データ
      real(8), dimension(0:km)  :: z_IntLon_xz        ! 鉛直格子点データ
      integer :: i

      z_IntLon_xz = 0.0d0
      do i=1,im
         z_IntLon_xz(:) = z_IntLon_xz(:) + xz(i,:) * x_Lon_Weight(i)
      enddo
    end function z_IntLon_xz